2012-07-08 51 views
4

我一直在使用Excel的求解器處理以下問題錯誤試圖約束優化使用Optim()

解決AB和C的方程中時:

y = a*b*c*x/((1 - c*x)(1 - c*x + b*c*x)) 

受到限制

0 < a < 100 
0 < b < 100 
0 < c < 100 

f(x[1]) < 10 
f(x[2]) > 20 
f(x[3]) < 40 

其中我有大約10(x,y)值對。我最小化abs(y - f(x))的和。我可以在每個x上限制我的函數結果的係數和值範圍。

我試過nls(沒有試圖強制約束),而Excel提供了幾乎任何我關心的初始值的估計值,nls幾乎從未返回答案。

我切換到使用優化,但我無法應用約束。

這是我已經變得如此之遠

best = function(p,x,y){sum(abs(y - p[1]*p[2]*p[3]*x/((1 - p[3]*x)*(1 - p[3]*x + p[2]*p[3]*x))))} 
p = c(1,1,1) 
x = c(.1,.5,.9) 
y = c(5,26,35) 
optim(p,best,x=x,y=y) 

我這樣做是爲了增加第一套constraints-

optim(p,best,x=x,y=y,method="L-BFGS-B",lower=c(0,0,0),upper=c(100,100,100)) 

我得到的錯誤 「」 錯誤:ABNORMAL_TERMINATION_IN_LNSRCH」

並最終得到一個更高的錯誤值($ value)。所以看起來我做錯了什麼,我無法弄清楚如何應用我的其他約束。

有人能給我提供一個基本的想法如何解決這個非統計學家能理解的問題嗎?我看了很多帖子,看了幾本R書。 R書停止了最簡單的使用優化。

+0

一般洞察關於優化與約束[這裏](http://stackoverflow.com/questions/9817001/optimization-with-constraints) – Chase 2012-07-08 23:31:21

回答

10

絕對值引入了一個奇點: 您可能希望使用方形代替 ,尤其是對於基於梯度的方法(例如L-BFGS)。

函數的分母可以是零。

參數出現在產品 中並且您允許它們(任意接近)零 也會導致問題。

您可以嘗試使用其他優化程序 (optimization task view上的完整列表), ,直到找到優化收斂的一個爲止。

x0 <- c(.1,.5,.9) 
y0 <- c(5,26,35) 
p <- c(1,1,1) 
lower <- 0*p 
upper <- 100 + lower 
f <- function(p,x=x0,y=y0) sum( 
    (
    y - p[1]*p[2]*p[3]*x/((1 - p[3]*x)*(1 - p[3]*x + p[2]*p[3]*x)) 
)^2 
) 

library(dfoptim) 
nmkb(p, f, lower=lower, upper=upper) # Converges 

library(Rvmmin) 
Rvmmin(p, f, lower=lower, upper=upper) # Does not converge 

library(DEoptim) 
DEoptim(f, lower, upper) # Does not converge 

library(NMOF) 
PSopt(f, list(min=lower, max=upper))[c("xbest", "OFvalue")] # Does not really converge 
DEopt(f, list(min=lower, max=upper))[c("xbest", "OFvalue")] # Does not really converge 

library(minqa) 
bobyqa(p, f, lower, upper) # Does not really converge 

作爲最後的手段,您始終可以使用網格搜索。

library(NMOF) 
r <- gridSearch(f, 
    lapply(seq_along(p), function(i) seq(lower[i],upper[i],length=200)) 
) 
+0

當我使用Excel,我試圖儘量減少兩者之和的問題絕對差異和絕對差異的平方和。我所得到的數據有時具有高度不確定性的數據點,廣場正在努力適應它。我在Excel中也能做的是使用個人估計值作爲約束條件。例如,y0 [3]應該是35。我可以應用約束y0 [3] <40.我看到如何設置參數的界限。我想我需要修改我最小化的函數來處理這種情況,但不知道該怎麼做。 – 2012-07-10 15:54:13

+0

您可以將約束違規作爲懲罰 添加到目標函數(「軟約束」), ,例如 'pmax(y0 [3] -40,0)^ 2'。 您還可以嘗試其他成本函數: 絕對值, 平滑絕對值, 'cost < - function(x)ifelse(abs(x)<。5,x^2,abs(x) - 。25 )' 或函數增長速度更慢,例如, 'cost < - function(x)ifelse(abs(x)<。5,x^2,sqrt(abs(x)) - sqrt(.5)+ 0.25)'。 如果您知道哪些觀測值是噪聲, 您可以在目標函數 中使用權重 來減少它們的影響。 – 2012-07-11 00:05:03