nls fit 中起始参数的线性模型

Question

我读到了关于使用日志作为创建线性方程的替代方法，以便我可以从 R 中的线性方程中提取 nls 拟合的起始值； 因此，

对于等式： Y=q/(1+bDX)^(1/b) 其中 Y 和 X 是我的数据； q,b,D 是我要估计的参数。我创建了线性模型：

  X<-1:45
  Y <- c(35326L, 30339L, 23379L, 21877L, 18629L, 17627L, 15691L, 15435L, 
   14205L, 11732L, 10560L, 10592L, 9939L, 7491L, 4928L, 3427L, 8123L, 
   9027L, 8733L, 9599L, 8737L, 9135L, 8548L, 7279L, 8940L, 8459L, 
   8460L, 7700L, 6817L, 7167L, 7089L, 7091L, 7538L, 9206L, 9680L, 
   5876L, 7799L, 8384L, 10586L, 8623L, 7848L, 5534L, 6610L, 6539L, 
   6650L)
lmodel <- coef(lm(log(Y)~X+I(X^2)))   
q0 <- exp(lmodel[1])
D0 <- -lmodel[2]
b0 <- lmodel[3]*2/D0^2 
Start1=list(q=q0,b=b0,D=D0)
theta_hat2 <- nls(y~q*(1+b*D*x)^(-1/b),start=Start1) 

这每次都能为我提供所需的结果。但是，我想将 nls 拟合到下面提到的其他一些对我来说相当困难的方程。如果有人可以帮助我制作一个适合这些的线性模型，我将不胜感激。

等式 2：

Y=q*(-D+(b/n)*X^n here q,b,D,n are to be estimated.

等式 3：

Y=q*exp⁡(-(X/D)^b here q,b,D are to be estimated.

等式 4：

Y=q*X^((-b) )*exp⁡(D/((1-b) )*(X^(1-b)-1) here q,b,D are to be estimated.

我正在附加此问题的示例数据集here：

Answer 1

主要问题是方程 (2) 和 (4) 参数化过大，因此无论使用哪种优化算法都会出现问题。在 (3) 的情况下，存在语法错误，我们需要更好的起始值。测试输入应包含在问题中，我们已在最后的注释中提供。

等式（1）

x 和 y 应该是 X 和 Y。我们可以使用plinear 算法来避免必须为线性参数q 指定起始值，在这种情况下，其他参数的起始值1 似乎就足够了。

fo1 <- Y ~ (1+b*D*X)^(-1/b)
fm1 <- nls(fo1, start = list(D = 1, b = 1), alg = "plinear")

等式（2）

等式 (2) 被过度参数化，因为如果你将 q 乘以任意数字 a 并同时将 D 和 b 除以 a 右手边不变。去掉 q 我们注意到 D 和 b 是线性输入的，只有 n 实际上是非线性的，所以我们可以使用 plinear 算法来避免线性参数的初始值：

fo2 <- Y ~ cbind(-1,(1/n)*X^n)
fm2 <- nls(fo2, start = list(n = 1), alg = "plinear")

等式（3）

问题中给出的公式中缺少右括号。如果我们解决了这个问题，那么问题是我们需要更好的起始值。首先尝试将 b 固定为 1，然后将结果用作整个方程的起始值。再次使用plinear 以避免必须为线性参数提供起始值。

fo3 <- Y ~ exp(-(X/D)^b)
b <- 1
fm <- nls(fo3, start = list(D = 1), alg = "plinear")
fm3 <- nls(fo3, start = list(D = coef(fm)[["D"]], b = 1), alg = "plinear")

等式 (4)

等式 (4) 也被过度参数化，因此将 D 设置为 1-b 在这种情况下 exp(D/(1-b)) 是常数 exp(1) 所以：

fo4 <- Y ~ X^((-b))* exp(1) * (X^(1-b)-1)
fm4 <- nls(fo4, start = list(b = .5), alg = "plinear")

比较

我们可以绘制各种解决方案：

plot(Y ~ X, pch = 20)
lines(fitted(fm1) ~ X)
lines(fitted(fm2) ~ X, col = "red")
lines(fitted(fm3) ~ X, col = "blue")
lines(fitted(fm4) ~ X, col = "green")
legend("topright", legend = 1:4, lty = 1,  col = c("black", "red", "blue", "green"))

注意

测试数据应该包含在问题中，由于它没有以那种形式提供，我们这次将在此处提供：

X <- 1:45
Y <- c(35326L, 30339L, 23379L, 21877L, 18629L, 17627L, 15691L, 15435L, 
14205L, 11732L, 10560L, 10592L, 9939L, 7491L, 4928L, 3427L, 8123L, 
9027L, 8733L, 9599L, 8737L, 9135L, 8548L, 7279L, 8940L, 8459L, 
8460L, 7700L, 6817L, 7167L, 7089L, 7091L, 7538L, 9206L, 9680L, 
5876L, 7799L, 8384L, 10586L, 8623L, 7848L, 5534L, 6610L, 6539L, 
6650L)