拟合非线性方程的参数以匹配实验数据

Question

我尝试修改我在这里找到的代码，使用 nls 函数来处理这个更简单的方程，但到目前为止没有成功。

new_f <- function(a, t, t1, dt, m){
  1-exp(-a*((x-t1)/dt)^(m+1))
}

fit_d <- nls(y~new_f(a, t, -30, dt, m), start=list(a=1, dt=46, m=1))

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

有人可以通过傻瓜分步指南帮助解决这个问题吗？任何建议将不胜感激！

我的数据：t=θ 和 y=Y(θ)

t <- c(-30, -29, -28, -27, -26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89)

y <- c(0.0000, -0.0013, -0.0027, -0.0041, -0.0055, -0.0069, -0.0085, -0.0101, -0.0119, -0.0137, -0.0155, -0.0173, -0.0192, -0.0212, -0.0234, -0.0258, -0.0285, -0.0291, -0.0250, -0.0224, -0.0268, -0.0298, -0.0297, -0.0324, -0.0347, -0.0353, -0.0377, -0.0392, -0.0392, -0.0383, -0.0309, -0.0166, 0.0006, 0.0233, 0.0521, 0.0843, 0.1200, 0.1597, 0.1997, 0.2381, 0.2762, 0.3137, 0.3501, 0.3864, 0.4229, 0.4590, 0.4951, 0.5312, 0.5655, 0.5967, 0.6245, 0.6494, 0.6726, 0.6939, 0.7132, 0.7316, 0.7487, 0.7641, 0.7784, 0.7917, 0.8037, 0.8146, 0.8248, 0.8343, 0.8432, 0.8516, 0.8596, 0.8671, 0.8740, 0.8805, 0.8868, 0.8931, 0.8992, 0.9050, 0.9105, 0.9155, 0.9202, 0.9247, 0.9290, 0.9330, 0.9369, 0.9405, 0.9439, 0.9470, 0.9501, 0.9529, 0.9557, 0.9583, 0.9608, 0.9633, 0.9658, 0.9682, 0.9704, 0.9726, 0.9747, 0.9768, 0.9788, 0.9808, 0.9827, 0.9844, 0.9858, 0.9869, 0.9880, 0.9890, 0.9899, 0.9909, 0.9918, 0.9927, 0.9938, 0.9947, 0.9956, 0.9964, 0.9970, 0.9975, 0.9980, 0.9985, 0.9989, 0.9993, 0.9997, 1.0000)

Answer 1

存在几个问题：

new_f 中的 x 应该是 t
问题是过度参数化，导致参数无法识别。有无数种参数组合可以满足该问题。特别是，a / dt^(m+1) 可以用单个参数代替。要解决这个问题，例如从起始值中删除 dt 并将其固定为 46。

通过这些更改，代码会产生一个结果。

new_f <- function(a, t, t1, dt, m){
  1-exp(-a*((t-t1)/dt)^(m+1))
}

dt <- 46
fit_d <- nls(y~new_f(a, t, -30, dt, m), start = list(a = 1, m = 1))
fit_d
## Nonlinear regression model
##   model: y ~ new_f(a, t, -30, dt, m)
##    data: parent.frame()
##     a     m 
## 0.560 3.315 
## residual sum-of-squares: 0.3099
##
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 4.165e-06

plot(y ~ t)
lines(fitted(fit_d) ~ t, col = "red")

Answer 2

从观察曲线的形状开始可以观察到：

对于 x 的高值，形状显示为指数类型，因为 y 的对数标度非常线性：

对于较低的 x 值，形状显示为线性类型：

如果我们尝试将曲线的两部分与 Heaviside 阶跃函数 H(x) 结合起来结果不是很令人满意，因为交界处不光滑，如下图所示：

为了获得平滑的曲线，可以用逻辑函数代替 Heaviside 函数：

以上参数数值未优化。为了更好地拟合，必须进行非线性回归（这还没有完成）。非线性回归的一个主要困难是设置好的参数“猜测”值以启动迭代演算。这就是为什么上述近似值可能非常有用。

拟合非线性方程的参数以匹配实验数据

Fitting parameters of a nonlinear equation to match experimental data

r

curve-fitting

nonlinear-optimization

model-fitting

non-linear-regression