使用 nls 将物理模型拟合到特定数据：过度参数化或无法识别的参数？

如何解决使用 nls 将物理模型拟合到特定数据：过度参数化或无法识别的参数？

我有一个有点复杂的物理模型，需要拟合五个未知参数，但目前没有成功。

我首先使用 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 点，因为它趋于无穷大）。