R_using nlsLM（）與約束

Question

我使用nlsLM {minpack.lm}找到的參數a和功能myfunb其數值給出的數據集，mydata最合適的。

mydata=data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45)) 

myfun=function(a,b,r,t){ 
    prd=a*b*(1-exp(-b*r*t)) 
    return(prd) 
} 

，並使用nlsLM

myfit=nlsLM(y~myfun(a,b,r=2,t=x),data=mydata,start=list(a=2000,b=0.05), 
        lower = c(1000,0), upper = c(3000,1)) 

它的工作原理。但是現在我想介紹一個約束條件，即a*b<1000。我查看了nlsLM中的可用選項，通過nls.lm.control設置約束條件。但這沒什麼幫助。有人可以幫助我，或者提出一個不同的方法來解決這個問題嗎？

Answer 1

據我所知，在minpack.lm包的nlsLM中，下限和上限參數是唯一可用的約束條件。
針對您的問題的一種可能的解決方案是基於使用包nloptr，R接口免費開源庫NLopt進行非線性優化。
在下面我使用cobyla，對於具有非線性不等式和等式約束自由導數的優化算法的代碼：

mydata <- data.frame(x=c(0,5,9,13,17,20), y=c(0,11,20,29,38,45)) 

myfun=function(a,b,r,t){ 
    prd=a*b*(1-exp(-b*r*t)) 
    return(prd) 
} 

library(nloptr) 

sse <- function(ab,r,x,y){ 
    sum((y - myfun(ab[1],ab[2],r,x))^2) 
} 

nonlincons <- function(ab,r,x,y) { 
    h <- numeric(1) 
    h[1] <- 1000 - ab[1]*ab[2] 
    return(h) 
} 

x0 <- c(2000,0.05) 
optpar <- cobyla(x0=x0, fn=sse, lower = c(1000,0), upper = c(3000,1), 
       hin=nonlincons, r=2, x=mydata$x, y=mydata$y, 
       control = list(xtol_rel = 1e-12, maxeval = 20000)) 
(optpar <- optpar$par) 

sum((mydata$y-myfun(optpar[1],optpar[2],r=2,mydata$x))^2) 

的最佳參數的值是：

[1] 3.000000e+03 2.288303e-02 

和的總和平方誤差爲：

[1] 38.02078 

希望它能幫助你。

Answer 2

看來，nlsLM不接受約束，一種可能的方法是使用constrOptim函數來代替，但缺點是，你必須把這個問題成constrOptim接受的形式。順便說一下，我目前正在研究一個R包，它爲optim和constrOptim提供了一個很好的界面，並自動進行必要的轉換，但尚未準備好使用。以下是如何使用constrOptim解決此問題。

第一步是constrOptim只接受線性約束， a * b < 1000 =>log(a) + log(b) < log(1000)，所以需要reparametrize使用log(a)和log(b)問題，約束變成：-log(a) - log(b) > -log(1000)，log(a) >= log(1000)，-log(a) >= -log(3000)，-log(b) >= 0，ui和ci是矩陣形式爲這些約束。 myfun1是您的功能的reparametrize版本，然後您可以使用constrOptim獲得您的解決方案log(a)和log(b)。

mydata=data.frame(x=c(0,5,9,13,17,20),y = c(0,11,20,29,38,45)) 

myfun=function(a,b,r,t){ 
    prd=a*b*(1-exp(-b*r*t)) 
    return(prd) 
} 

myfun1 <- function(logab, data){ 
    sum((data$y - myfun(exp(logab[1]), exp(logab[2]), r = 2, t = data$x))^2) 
} 

ui <- rbind(c(-1,-1), c(1,0), c(-1,0), c(0,-1)) 
ci <- c(-log(1000), log(1000), -log(3000), 0) 
init <- log(c(2000, 0.05)) 
r <- constrOptim(theta = init, f = myfun1, grad = NULL, ui = ui, ci = ci, data = mydata)