无法使用 optim 优化负二项式模型

Question

我正在尝试使用 R 中的 optim 包手动优化负二项式回归模型，尝试使用以下代码使用因子矩阵 X 预测计数变量 y：

# 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 的起始值也无济于事。

您是否知道似然函数是否有任何特殊性导致值以任何奇怪的方式处理？

我们将不胜感激。

Answer 1

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)
  
}