如何在 DEoptim (R) 中使用 "fnmap" 设置整数约束



我在R中使用DEoptim,以便最大化二变量学生t分布的似然函数。

这种分布要求自由度是整数(我已经可以使用DEoptim中的下界和上界来设置正逆)。

该pdf详细解释了DEoptim,并指出

[fnmap]是一个可选函数,它将在创建每个填充之后运行,但在将总体传递给目标函数。这允许用户强制整数/基数约束。

它没有详细说明如何使用参数"fnmap",也没有给出任何细节。在网上,我发现了一个由包的开发者给出的施加肉欲约束的示例函数,但他也没有解释其中的合理性。

那么,如果我有一个参数(例如参数(x,y,z)中的"z"),我应该把什么作为"fnmap"来将z限制为整数值?

我想这就是您想要的。

fnmap_f <- fuction(x) c(x[1], x[2], round(x[3]))
DEoptim(..., fnMap = fnmap_f)

框约束可以设置为目标函数本身的惩罚约束,如下所示。

lambda乘数的值决定了目标函数和约束之间的相对重要性。lambda的值越高,求解器将优先考虑约束而不是目标函数。

在这里,我们求解一个简单的方程,其整数解为x1^2+x2^2=25,其中x2<x1作为约束,具有c(4,3)的预期输出

library("DEoptim")
#run DEoptim and set a seed first for replicability
#note the results are sensitive to seed values and parameters (lambda,NP,itermax,F,CR)
set.seed(1234)
#create a vector/grid of lambda values 
lambdaVec = sapply(0:12,function(x) 10^x)
lambdaVec
#[1] 1e+00 1e+01 1e+02 1e+03 1e+04 1e+05 1e+06 1e+07 1e+08 1e+09 1e+10

#For each value of lambda compute the output of optimization function and combine the results
optimSummary = do.call(rbind,lapply(lambdaVec, function(lambdaParam)  {

fnCustom = function(x,lambda=lambdaParam) {
x1 <- x[1]
x2 <- x[2]
#integer param constraints  
firstIntPenalty  <- (x1-round(x1))^2
secondIntPenalty  <-  (x2-round(x2))^2
# x2 < x1, note the sign is opposite of desired constraint
inEqualityConstraint <- sum(x2>x1)

100 * (x1^2 + x2^2 - 25)^2 + lambda * ( firstIntPenalty + secondIntPenalty + inEqualityConstraint ) 
}

lower <- c(0,0)
upper <- c(5,5)

#you will have to tinker with values of NP,F and CR and monitor convergence, see ?DEoptim last paragraph
outDEoptim <- DEoptim(fnCustom, lower, upper, DEoptim.control(NP = 80, itermax = 1000, F = 1.2, CR = 0.7,trace=TRUE))
#output data.frame of optimization result
optimRes <- data.frame(lambda = lambdaParam ,param1 = outDEoptim$optim$bestmem[1],param2 = outDEoptim$optim$bestmem[2])

rownames(optimRes) <- NULL
return(optimRes)
} ))

浮点表示法:

由于浮点表示,在大多数情况下,结果不会完全等于您期望的整数输出。根据您的域,必须定义一个可接受的阈值,低于该阈值,我们将数字视为整数。

有关更多详细信息,请参阅?.Machine和这个R浮点精度

输出收敛和验证:

threshold = 1e-6
expectedOut = c(4,3)
#optimization summary
optimSummary
#   lambda   param1       param2
#1   1e+00 4.999996 0.0002930537
#2   1e+01 4.000000 3.0000000000
#3   1e+02 4.000000 3.0000000000
#4   1e+03 4.000000 3.0000000000
#5   1e+04 4.000000 3.0000000000
#6   1e+05 4.000000 3.0000000000
#7   1e+06 4.000000 3.0000000000
#8   1e+07 4.000000 3.0000000000
#9   1e+08 4.000000 2.9999999962
#10  1e+09 3.999999 2.9999998843
#11  1e+10 4.000000 2.9999999569
#12  1e+11 4.000000 3.0000000140
#13  1e+12 4.000000 3.0000000194


#Verify output
#1)With constraintBreach1 and constraintBreach2 we check if difference between output and expected result
#has breached threshold 
#2)With constraintBreach3 we check if x1 > x2 condition is violated
#3)Columns with TRUE observations indicate breach of respective constraints at particular lambda values

verifyDF = data.frame(lambdaVec,constraintBreach1 = abs(optimSummary$param1-expectedOut[1]) > threshold
, constraintBreach2 = abs(optimSummary$param2-expectedOut[2]) > threshold
,  constraintBreach3 =  optimSummary$param1 < optimSummary$param1)
verifyDF
#   lambdaVec constraintBreach1 constraintBreach2 constraintBreach3
#1      1e+00              TRUE              TRUE             FALSE
#2      1e+01             FALSE             FALSE             FALSE
#3      1e+02             FALSE             FALSE             FALSE
#4      1e+03             FALSE             FALSE             FALSE
#5      1e+04             FALSE             FALSE             FALSE
#6      1e+05             FALSE             FALSE             FALSE
#7      1e+06             FALSE             FALSE             FALSE
#8      1e+07             FALSE             FALSE             FALSE
#9      1e+08             FALSE             FALSE             FALSE
#10     1e+09             FALSE             FALSE             FALSE
#11     1e+10             FALSE             FALSE             FALSE
#12     1e+11             FALSE             FALSE             FALSE
#13     1e+12             FALSE             FALSE             FALSE

对于较低级别的lambda,解算器会忽略约束,并且随着lambda的增加,解算程序会为其指定更多权重满足相对于目标函数的约束,从而满足约束。

对于您的特定问题,您将不得不修改lambda,NP,itermax,F,CR的值。

最新更新