R_using 有约束的 nlsLM()
R_using nlsLM() with constraints
我正在使用 nlsLM {minpack.lm}
来查找函数 myfun
的参数 a
和 b
的值,它们最适合数据集,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
设置约束。但帮助不大。有人可以帮我解决这个问题吗?
据我所知,在 minpack.lm 包的 nlsLM 中,下限和上限参数边界是唯一可用的约束。
您的问题的一种可能解决方案是基于包 nloptr
的使用,它是免费开源库 NLopt 的 R 接口,用于非线性优化。
在下面的代码中,我使用了 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
希望对您有所帮助。
似乎nlsLM
不接受约束,一种可能的方法是使用constrOptim
函数代替,但缺点是你必须将问题转化为[=]的形式12=] 接受。顺便说一句,我目前正在开发一个 R 包,它为 optim
和 constrOptim
提供了一个很好的接口,并自动进行必要的转换,但它还没有准备好使用。下面是如何使用 constrOptim
来解决这个问题。
第一步是constrOptim
只接受线性约束,
a * b < 1000
=> log(a) + log(b) < log(1000)
,因此您需要使用 log(a)
和 log(b)
重新参数化问题,约束变为:-log(a) - log(b) > -log(1000)
、log(a) >= log(1000)
、[ =24=、-log(b) >= 0
、ui
和 ci
是这些约束的矩阵形式。 myfun1
是您的函数的重新参数化版本,然后您可以使用 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)
我正在使用 nlsLM {minpack.lm}
来查找函数 myfun
的参数 a
和 b
的值,它们最适合数据集,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
设置约束。但帮助不大。有人可以帮我解决这个问题吗?
据我所知,在 minpack.lm 包的 nlsLM 中,下限和上限参数边界是唯一可用的约束。
您的问题的一种可能解决方案是基于包 nloptr
的使用,它是免费开源库 NLopt 的 R 接口,用于非线性优化。
在下面的代码中,我使用了 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
希望对您有所帮助。
似乎nlsLM
不接受约束,一种可能的方法是使用constrOptim
函数代替,但缺点是你必须将问题转化为[=]的形式12=] 接受。顺便说一句,我目前正在开发一个 R 包,它为 optim
和 constrOptim
提供了一个很好的接口,并自动进行必要的转换,但它还没有准备好使用。下面是如何使用 constrOptim
来解决这个问题。
第一步是constrOptim
只接受线性约束,
a * b < 1000
=> log(a) + log(b) < log(1000)
,因此您需要使用 log(a)
和 log(b)
重新参数化问题,约束变为:-log(a) - log(b) > -log(1000)
、log(a) >= log(1000)
、[ =24=、-log(b) >= 0
、ui
和 ci
是这些约束的矩阵形式。 myfun1
是您的函数的重新参数化版本,然后您可以使用 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)