我试图通过在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
我尝试了几个初始值，但仍然得到相同的错误。你能解决这个问题吗？或者可能还有其他代码可以解决这个问题？
谢谢。
问题来自于计算一个大数字的阶乘，然后取其对数。阶乘数太高，R无法识别为有限数，但它的对数不是。在这种情况下，使用lgamma函数可以得到与log(factorial(y))相同的结果。
这不是黑客攻击；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
无法进行优化，因为您有非常大的值，这会导致无穷大或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可能为负，这是一个问题。您可能必须使用约束优化来解决此问题。
注意，最大化的lambda值
lnL <- sum(y*log(lambda) - log(factorial(y)) - lambda)
是最大化的相同值
lnL_2 <- sum(y*log(lambda) - lambda)
因此您可以优化lnL_ 2而不是lnL。例如，请参阅数学堆栈交换中的此答案以获取推导。