问题描述
我有一个有点复杂的物理模型,需要拟合五个未知参数,但目前没有成功。
我首先使用 nls2
来获得一些初始值的估计值,但随后 nls
、nlxb
和 nlsLM
都抛出了著名的“初始奇异梯度误差”参数估计”错误。
对于 nls2
的起始值,我是从文献中提取出来的,所以我认为至少对于 nls2
我有很好的起始值。从 nls2
中提取的参数估计在物理上也很有意义;但是,不要解决奇异梯度矩阵误差的问题。
既然是物理模型,每个系数都有物理意义,我宁愿不修正任何一个。
我还应该提到模型方程中的所有五个未知参数都是正的,形状参数m可以达到2。
阅读了许多帖子并尝试了不同的解决方案建议,我得出的结论是我存在参数过度化或无法识别的参数问题。
我的问题是我应该停止尝试在这个特定模型(有这么多未知参数)中使用 nls
还是有什么出路?
这是我的 MWE:
# Data
x <- c(0,1000,2000,2500,2750,3000,3250,3500,3750,4000,5000)
y <- c(1.0,0.99,0.98,0.95,0.795,0.59,0.35,0.295,0.175,0.14,0.095)
# Start values for nls2
bounds <- data.frame(a = c(0.8,1.5),b = c(1e+5,1e+7),c = c(0.4,1.4),n = c(0.1,2),m = c(0.1,2))
# Model equation function
mod <- function(x,a,b,c,n,m){
t <- b*85^n*exp(-c/0.0309)
(1 - exp(-(a/(t*x))^m))
}
# # Model equation
# mod <- y ~ (1 - exp(-(a/(b*85^n*exp(-c/0.0309)*x))^m))
# Model fit with nls2
fit2 <- nls2(y ~ mod(x,m),data = data.frame(x,y),start = bounds,algorithm = "brute-force")
# Model fit with nls
fit <- nls(y ~ mod(x,start = coef(fit2))
解决方法
我越看越困惑,但我会再试一次。
再看你的表达式,我们有指数里面的表达式
-(a/(b*85^n*exp(-c/0.0309)*x))^m
我们可以将其重写为
-( [a/(b*85^n*exp(-c/0.0309))] * 1/x )^m
(请检查我的代数!)
如果这是正确的,那么整个粗体 blob 不会影响 x
的函数形式——它全部折叠为等式中的单个常量。 (换句话说,{a,b,c,n}
都是共同无法识别。)将这些东西合并到单个参数 phi
中:
1 - exp(-(phi/x)^m)
-
phi
是一个形状参数(与x
具有相同的单位,应该与x
的典型值大致相同):让我们尝试使用 2500 的起始值(x
的平均值) -
m
是一个形状参数;从m==1
开始,我们不会出错
现在 nls
无需任何额外帮助即可正常工作:
n1 <- nls(y~1 - exp(-(phi/x)^m),start=list(phi=2500,m=1),data=data.frame(x,y))
并得到 phi
=2935,m
=6.49。
情节预测:
plot(x,y,ylim=c(0,1))
xvec <- seq(0,5000,length=101)
lines(xvec,predict(n1,newdata=data.frame(x=xvec)))
考虑这条曲线的作用的另一种方式:我们可以将方程转换为 -log(1-y) = phi^m*(1/x)^m
:也就是说,-log(1-y)
应该遵循关于 1/x
的幂律曲线。
这是它的样子:
plot(1/x,-log(1-y))
## curve() uses "x" as the current x-axis variable,i.e.
## read "x" as "1/x" below.
with(as.list(coef(n1)),curve(phi^m*x^m,add=TRUE))
在这种格式中,它似乎很好地拟合了中心数据,但对于较大的 1/x
值却失败了(此处缺少 x=0
点,因为它趋于无穷大)。