使用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个数据集(或第3行Y)有一条错误消息: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个参数,如斜率和截距但在quesiton尝试式定义三个:a,c和e.

我们需要修复其中一个变量,或者通常添加一个额外的约束.现在if 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