我正在嘗試使用 R 中的 optim 包手動優化負二項式回歸模型，嘗試使用以下代碼使用y因子矩陣預測計數變數X：
# generating some fake data
n <- 1000
X <- matrix(NA, ncol = 5, nrow = n)
X[,1] <- 1
X[,2] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,3] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,4] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,5] <- sample(size = n, x = c(0,1), replace = TRUE)
beta0 <- 3
beta1 <- -2
beta2 <- -2
beta3 <- -4
beta4 <- -0.9
k <- 0.9
## draws from negative binomial distribution
mu <- exp(beta0 beta1 * X[,2] beta2 * X[,3] beta3 * X[,4] beta4 * X[,5])
theta <- mu mu ^2 / k
# dependent variable
y <- rnegbin(n, mu = mu, theta = theta)
# function to be optimised
negbin_ll <- function(y, X, theta){
beta <- theta[1:ncol(X)]
alpha <- theta[ncol(X) 1]
logll <- y * log(alpha) y *( beta %*% t(X) ) - (y (1 / alpha ) ) * log( 1 alpha * exp(beta %*% t(X))) lgamma(y (1 / alpha)) - lgamma ( y 1) - lgamma ( 1 / alpha)
logll <- sum( logll )
return(logll)
}
stval <- rep(0, ncol(X) 1)
res <-
optim(
stval,
negbin_ll,
y = y,
X = X,
control = list(fnscale = -1),
hessian = TRUE,
method = "BFGS"
)
代碼應該從優化程序中產生點估計，但在執行優化函式時失敗 error in optim(stval, negbin_ll, y = y, X = X, control = list(fnscale = -1), : initial value in 'vmmin' is not finite.
我已經嘗試在似然函式中更改log(gamma(...))為lgamma(...)并嘗試了許多其他方法，但無法獲得估計值。
更改 optim 的起始值也無濟于事。
您是否知道似然函式是否有任何特殊性導致以任何奇怪的方式處理值？
幫助將不勝感激。
uj5u.com熱心網友回復：
optim嘗試幾個點以達到最小值，在您的情況下，它在日志中的引數中遇到了一些非正值。一種方法是通過回傳負（在您的情況下）大數，例如-lenght(series)*10^6. 重新制作對數似然函式，像這樣它有點作業：
negbin_ll <- function(y, X, theta){
beta <- theta[1:ncol(X)]
alpha <- theta[ncol(X) 1]
if(any(alpha<=0)) return(-length(y)*10^6)
if(any(1 alpha * exp(beta %*% t(X))<=0)) return(-length(y)*10^6)
logll <- y * log(alpha) y *( beta %*% t(X) ) - (y (1 / alpha ) ) * log( 1 alpha * exp(beta %*% t(X))) lgamma(y (1 / alpha)) - lgamma ( y 1) - lgamma ( 1 / alpha)
logll <- sum( logll )
return(logll)
}
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/388783.html
上一篇：如何在R代碼中為這種情況設定條件
下一篇：洗掉R中的單引號