时间戳

首页 > 解决方案 > 优化 R 中的错​​误 - 无法在初始参数处进行评估

r - 优化 R 中的错​​误 - 无法在初始参数处进行评估

问题描述

我正在尝试通过在 R 中使用 optim 来估计参数 a、b、c 和 s。这是我的代码。

age <- c(0,30,60,90)
Dx <- c(49294.57, 2975.1, 11456.38, 2977.08)
Ex <- c(1572608.38, 1531956.05, 650404.58, 9728.47)

log_lik <- function(par,x,y,z){
  a <- par[1]
  b <- par[2]
  c <- par[3]
  s <- par[4]
  mu <- (a*exp(b*x))/(1+s * (a)/(b) * (exp(b*x)-1)) + c
  lambda <- mu * z
  
  lnL <- sum(y*log(lambda) - log(factorial(y)) - lambda)
  -lnL
}

optim(c(1,1,1,1),log_lik, x = age, y = Dx, z = Ex)

但是，我得到一个错误

Error in optim(c(1, 1, 1, 1), log_lik, x = age, y = Dx, z = Ex) : 
  function cannot be evaluated at initial parameters

我尝试了几个初始值，但仍然得到相同的错误。你能解决这个问题吗？或者也许有另一个代码来解决这个问题？

谢谢你。

标签： roptimizationnonlinear-optimization

解决方案

问题来自计算一个大数的阶乘然后取其对数。阶乘数太高，R 无法将其识别为有限数，但其对数不是。log(factorial(y))在这种情况下，我们可以得到与使用lgamma函数相同的结果。

这不是黑客攻击；R 中的factorial函数只是该函数的一个薄包装器gamma：

factorial
#> function (x) 
#> gamma(x + 1)

所以我们可以得到一个函数，它log(factorial(y))不需要实际执行计算极高数字然后获取它们的日志的步骤，如下所示：

log_factorial <- function(x) lgamma(x + 1)

我们可以看到给了我们正确的结果：

log(factorial(21))
#> [1] 45.38014

log_factorial(21)
#> [1] 45.38014

但是允许我们输入更高的数字而不会产生无穷大。

log(factorial(200))
#> [1] Inf

log_factorial(200)
#> [1] 863.232

因此，我们可以将您的代码稍微更改为：

log_lik <- function(par,x,y,z){
  a <- par[1]
  b <- par[2]
  c <- par[3]
  s <- par[4]
  mu <- (a*exp(b*x))/(1+s * (a)/(b) * (exp(b*x)-1)) + c
  lambda <- mu * z
  
  lnL <- sum(y*log(lambda) - lgamma(y + 1) - lambda)
  -lnL
}

现在我们得到：

optim(c(1,1,1,1), log_lik, x = age, y = Dx, z = Ex)
#> $par
#> [1]  0.6114036  1.1267546 -0.5800334  1.9163744
#> 
#> $value
#> [1] 15828.8
#> 
#> $counts
#> function gradient 
#>      161       NA 
#> 
#> $convergence
#> [1] 0

$message
NULL

推荐阅读