2017-02-13 41 views
2

我在R中使用DEoptim以最大化雙變量學生t分佈的似然函數。如何在DEoptim(R)中使用「fnmap」設置整數約束

此分佈要求自由度爲整數(我可以使用DEoptim中的下限和上限設置正面禁忌)。

This pdf詳細解釋DEoptim並指出

[fnmap]是在創建每個種羣之後將要運行的可選功能,但是 之前的人口被傳遞給目標函數。這允許用戶強加 整數/基數consulatnts。

它沒有詳細說明如何使用參數「fnmap」,也沒有提供任何細節。在網上查找時,我發現了一個示例函數,用來強制開發者對軟件包進行強制約束,但他並沒有解釋那裏的理性。因此,如果我有一個參數(例如參數(x,y,z)中的'z'),我將「fnmap」作爲參數來限制z爲整數值?

回答

1

我想這是你想要的。

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

這工作完美。 – Patty

1

方框約束可以設置爲目標函數本身的懲罰約束,如下所示。

lambda乘數的值表示目標函數和約束之間的相對重要性。 對於lambda的更高值,求解器將優先考慮目標函數的約束條件。

在這裏,我們解決對於x1^2 + X 2^2 = 25×2 < 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 floating point precision

輸出收斂和驗證:

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,NP,itermax,F,CR的值。