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

【问题标题】：Error in optim R - cannot be evaluated at initial parameters优化 R 中的错误 - 无法在初始参数处进行评估
【发布时间】：2021-11-05 17:14:11
【问题描述】：

我正在尝试通过在 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

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

谢谢。

【问题讨论】：

问题已经出在factorial(y) 部分，这从一开始就是无限的，你的值太大了。
那么，如何解决这个问题呢？因为方程由factorial(y)
不要使用log(factorial(x))；而是使用lfactorial(x)。

标签： r optimization nonlinear-optimization

【解决方案1】：

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

这不是 hack； 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

【讨论】：

感谢您的回答。我忘了说所有参数都必须是正数。从这个结果来看，参数 c 是负数。如何修改代码？
@mathz 函数中的公式使用年龄、Dx 和 Ex 的值给出这些结果。我不知道您要计算什么，也不知道您是从哪张纸上得出的，所以除此之外我无法告诉您任何其他信息。如果结果没有意义，那么要么是您的数据、公式存在问题，要么是您将此模型应用于这些数据。在不知道您要做什么的情况下，我无能为力。
既然有开箱即用的lfactorial，为什么还要定义log_factorial 函数？
@nicola 我不知道它存在，所以感谢您指出这一点。有趣的是，它的代码与我定义的函数完全相同。嗯，一直在学习！

【解决方案2】：

无法进行优化，因为您的值非常大，这会导致无限或 NA 值。一种选择可能是重新调整您的变量，例如，如果您的变量自然在 100 万左右的范围内，则将所有值除以 100 万。例如。

age=age/1e2
Dx=Dx/1e5
Ex=Ex/1e6

现在优化工作并返回

$par
[1]  1.418161 37.235806 -1.104942 31.443860

$value
[1] 1.421373

$counts
function gradient 
     479       NA 

$convergence
[1] 0

$message
NULL

Warning messages:
1: In log(lambda) : NaNs produced
2: In log(lambda) : NaNs produced
3: In log(lambda) : NaNs produced
4: In log(lambda) : NaNs produced
5: In log(lambda) : NaNs produced
6: In log(lambda) : NaNs produced
7: In log(lambda) : NaNs produced
8: In log(lambda) : NaNs produced
9: In log(lambda) : NaNs produced

log(lambda) 部分仍然存在问题，因为 lambda 可能是负数，这是一个问题。您可能必须使用约束优化来解决此问题。

【讨论】：

【解决方案3】：

注意，最大化的 lambda 值

  lnL <- sum(y*log(lambda) - log(factorial(y)) - lambda)

是最大化的相同值

  lnL_2 <- sum(y*log(lambda) - lambda)

因此您可以优化 lnL_2 而不是 lnL。参见，例如，this answer 在 Math Stackexchange 上进行推导。

【讨论】：