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

Question

我正在使用 minpack.lm 包中的 nlsLM 函数来查找参数 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 中的每一行对应一个数据集。通过运行代码，第三个数据集(或 Y 中的第 3 行)出现错误消息:初始参数估计处的奇异梯度矩阵。以下是代码:

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 和 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)

，其中截距和斜率分别是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