使用nlsLM函数查找非线性模型的初始条件

Question

我正在使用nlsLM包中的minpack.lm函数来查找参数a, e和c的值，它们最适合数据out。这是我的代码：

n <- seq(0, 70000, by = 1)
TR <- 0.946
b <- 2000
k <- 50000
nr <- 25
na <- 4000
nd <- 3200
d <- 0.05775

y <- d + ((TR*b)/k)*(nr/(na + nd + nr))*n
## summary(y)
out <- data.frame(n = n, y = y)
plot(out$n, out$y)

## Estimate the parameters of a nonlinear model
library(minpack.lm)
k1 <- 50000
k2 <- 5000

fit_r <- nlsLM(y ~ a*(e*n + k1*k2 + c), data=out,
               start=list(a = 2e-10,
                          e = 6e+05, 
                          c = 1e+07), lower = c(0, 0, 0), algorithm="port")
print(fit_r)
## summary(fit_r)

df_fit <- data.frame(n = seq(0, 70000, by = 1))
df_fit$y <- predict(fit_r, newdata = df_fit)
plot(out$n, out$y, type = "l", col = "red", ylim = c(0,10))
lines(df_fit$n, df_fit$y, col="green")
legend(0,ceiling(max(out$y)),legend=c("observed","predicted"), col=c("red","green"), lty=c(1,1), ncol=1)

对数据的拟合似乎对初始条件非常敏感。例如：

• 与list(a = 2e-10, e = 6e+05, c = 1e+07)，这是一个很好的契合：

Nonlinear regression model
  model: y ~ a * (e * n + k1 * k2 + c)
   data: out
        a         e         c 
2.221e-10 5.895e+05 9.996e+06 
 residual sum-of-squares: 3.225e-26

Algorithm "port", convergence message: Relative error between `par' and the solution is at most `ptol'.

• 与list(a = 2e-01, e = 100, c = 2)，这是一个不合适的：

Nonlinear regression model
  model: y ~ a * (e * n + k1 * k2 + c)
   data: out
        a         e         c 
1.839e-08 1.000e+02 0.000e+00 
 residual sum-of-squares: 476410

Algorithm "port", convergence message: Relative error in the sum of squares is at most `ftol'.

那么，有没有一种有效的方法来找到能够很好地拟合数据的初始条件？

编辑：

我添加了以下代码来更好地解释问题。该代码用于查找a，e和c的值，这些值最适合来自多个数据集的数据。 Y中的每一行对应一个数据集。通过运行代码，第3个数据集（或Y中的第3行）有一条错误消息：singular gradient matrix at initial parameter estimates.这是代码：

TR <- 0.946
b <- 2000
k <- 50000
nr <- 25
na <- 4000
nd <- 3200
d <- 0.05775

Y <- data.frame(k1 = c(114000, 72000, 2000, 100000), k2 = c(47356, 30697, 214, 3568), n = c(114000, 72000, 2000, 100000), 
           na = c(3936, 9245, 6834, 2967), nd = c(191, 2409, 2668, 2776), nr = c(57, 36, 1, 50), a = NA, e = NA, c = NA)

## Create a function to round values
roundDown <- function(x) {
  k <- floor(log10(x))
  out <- floor(x*10^(-k))*10^k
  return(out)
}

ID_line_NA <- which(is.na(Y[,c("a")]), arr.ind=TRUE)
## print(ID_line_NA)

for(i in ID_line_NA){

  print(i)

  ## Define the variable y
  seq_n <- seq(0, Y[i, c("n")], by = 1)
  y <- d + (((TR*b)/(Y[i, c("k1")]))*(Y[i, c("nr")]/(Y[i, c("na")] + Y[i, c("nd")] + Y[i, c("nr")])))*seq_n
  ## summary(y)
  out <- data.frame(n = seq_n, y = y)
  ## plot(out$n, out$y)

  ## Build the linear model to find the values of parameters that give the best fit 
  mod <- lm(y ~ n, data = out)
  ## print(mod)

  ## Define initial conditions
  test_a <- roundDown(as.vector(coefficients(mod)[1])/(Y[i, c("k1")]*Y[i, c("k2")]))
  test_e <- as.vector(coefficients(mod)[2])/test_a
  test_c <- (as.vector(coefficients(mod)[1])/test_a) - (Y[i, c("k1")]*Y[i, c("k2")])

  ## Build the nonlinear model
  fit <- tryCatch( nlsLM(y ~ a*(e*n + Y[i, c("k1")]*Y[i, c("k2")] + c), data=out,
                                   start=list(a = test_a,
                                              e = test_e,
                                              c = test_c), lower = c(0, 0, 0)),
                             warning = function(w) return(1), error = function(e) return(2))
  ## print(fit)

  if(is(fit,"nls")){

    ## Plot
    tiff(paste("F:/Sources/Test_", i, ".tiff", sep=""), width = 10, height = 8, units = 'in', res = 300)
    par(mfrow=c(1,2),oma = c(0, 0, 2, 0))
    df_fit <- data.frame(n = seq_n)
    df_fit$y <- predict(fit, newdata = df_fit)
    plot(out$n, out$y, type = "l", col = "red", ylim = c(0, ceiling(max(out$y))))
    lines(df_fit$n, df_fit$y, col="green")
    dev.off()

    ## Add the parameters a, e and c in the data frame
    Y[i, c("a")] <- as.vector(coef(fit)[c("a")])
    Y[i, c("e")] <- as.vector(coef(fit)[c("e")])
    Y[i, c("c")] <- as.vector(coef(fit)[c("c")])

  } else{

    print("Error in the NLM")

  }
}

因此，使用约束a > 0, e > 0, and c > 0，是否有一种有效的方法来查找nlsLM函数的初始条件，该函数可以很好地拟合数据并避免错误消息？

我添加了一些条件来定义参数a，e和c的初始条件：

使用线性模型lm(y ~ n)的结果：

c = intercept/a - k1*k2 > 0 and 
e = slope/a > 0
0 < a < intercept/(k1*k2)

其中intercept和slope分别是lm(y ~ n)的截距和斜率。

Answer 1

问题不在于如何找到参数的初始值。问题是这是一个带参数的重新参数化线性模型。线性模型的参数是斜率a*e和截距a*(k1 * k2 + c)所以只能有2个参数，如斜率和截距，但问题中的公式试图定义三个：a，c和e。

我们需要修复其中一个变量，或者通常添加一个额外的约束。现在，如果co是一个向量，其第一个元素是截距，第二个元素是线性模型的斜率

fm <- lm(y ~ n)
co <- coef(fm)

然后我们有方程式：

co[[1]] = a*e
co[[2]] = a*(k1*k2+c)

co，k1和k2是众所周知的，如果我们认为c是固定的，那么我们可以解决a和e给出：

a = co[[2]] / (k1*k2 + c)
e = (k1 * k2 + c) * co[[1]] / co[[2]]

既然co[[1]]和co[[2]]都是正面的，c必须太a和e必然是正面的，一旦我们任意修复c给我们解决方案。这给出了无数个a，e对，它们最小化模型，每个非负值c。请注意，我们不需要为此调用nlsLM。

例如，对于c = 1e-10，我们有：

fm <- lm(y ~ n)
co <- coef(fm)

c <- 1e-10
a <- co[[2]] / (k1*k2 + c)
e <- (k1 * k2 + c) * co[[1]] / co[[2]]
a; e
## [1] 5.23737e-13
## [1] 110265261628

注意，由于系数之间的大小差异，可能存在数值问题;然而，如果我们增加c将增加e并减少a使得缩放更加糟糕，因此问题中给出的这个问题的参数化似乎具有继承性的不良数值缩放。

请注意，这一切都不需要运行nlsLM来获得最佳系数;但是，由于缩放不良，可能仍有可能在某种程度上改善答案。

co <- coef(lm(y ~ n))
c <- 1e-10
a <- co[[2]] / (k1*k2 + c)
e <- (k1 * k2 + c) * co[[1]] / co[[2]]
nlsLM(y ~ a * (e * n + k1 * k2 + c), start = list(a = a, e = e), lower = c(0, 0))

这使：

Nonlinear regression model
  model: y ~ a * (e * n + k1 * k2 + c)
   data: parent.frame()
        a         e 
2.310e-10 5.668e+05 
 residual sum-of-squares: 1.673e-26

Number of iterations to convergence: 12 
Achieved convergence tolerance: 1.49e-08