使用 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 中的每一行对应一个数据集。通过运行代码，第 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 使缩放比例变得更糟，因此问题中给出的该问题的参数化似乎具有固有的不良数值缩放比例。

请注意，其中 none 需要运行 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

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

Find initial conditions for nonlinear models using the nlsLM function

r

non-linear-regression