如何在 R 中编码,这个对数似然函数

问题描述 投票:0回答:1

给出了有序logit模型;

考虑 Logit 模型的以下对数似然函数:

还给出 k = {0,1,2,3,4},其中对于 k = 0 和 k = 4,我们的 alpha 分别为 -Inf 和 +Inf。现在我的问题是如何编写一个优化对数似然函数的函数,即找到 B = B1 的 MLE 和 k = {1,2,3} 的 alpha。

这是我的尝试;

J <- length(Poss)

log_likelihood <- function(params, Poss = Poss){
  beta <- params[1]
  alpha <- c(-Inf, params[2:4], +Inf)
  
  for (k in 2:4) {
    for (j in J) {
      delta[k] <- alpha[k] - Poss[j] * beta
      ll <- delta[k] - log(1 + exp(delta[k]) - delta[k] + log(1 + exp(delta[k])))
      return(ll)
    }
  }
  return(sum(ll))
}

guess_optim <- c(0,0,0,0) #for beta1 and alpha1,2,3

optim(guess_optim, log_likelihood, Poss = Poss)

只有这样不行,我该如何解决这个问题?

r numerical-methods mle
1个回答
0
投票

我认为第一个问题是您在内循环中

return(ll)
,这将导致第一个记录的第一个值——您的过程不会迭代,直到它到达
return(sum(ll))
。第二个问题是,观察结果应该只对它们实际所属类别的可能性做出贡献,你不做这种区分。最后,
optim
最小化,所以你不应该返回你想要最大化的值。

为了测试实现,我们首先获取一些示例数据(来自here):

## Example data
dat <- foreign::read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")

Y <- as.integer(cut(dat$gpa, 4))  ## Create 4 quartiles
X <- dat$public                   ## Dummy predictor, could be slope too

## This is the result we'll be reproducing
MASS::polr(as.factor(Y) ~ X)

#> Coefficients:
#>        X 
#> 1.178314 
#> 
#> Intercepts:
#>         1|2         2|3         3|4 
#> -2.43862931 -0.06816316  2.05272249 

现在,这是一个可以为您提供大致相同的参数估计值的实现:

## Logit link inversion
expit <- function(x) 1/(1+exp(-x))

## -log-likelihood function
ll <- function(params, X, Y) {
  alpha <- params[-1]
  beta <- params[1]
  
  ll <- unlist(mapply(\(x, y) {
    eta <- alpha - x*beta
    mu <- expit(eta)
    ## This line calculates expit(delta) - expit(delta_-1) in probability scale
    probs <- pmin(pmax(diff(c(0, mu, 1)), .Machine$double.eps), 1)
    ## Only use the probability of the observed category!
    log(probs[y])
  }, x=X, y=Y))
  ## Return inverse so value is maximized
  -sum(ll)
}

guess <- c(1, -2, 0, 2)
optim(guess, ll, X=X, Y=Y)
#> $par
#> [1]  1.17991554 -2.43853589 -0.06775819  2.05298068

一些注意事项:

  • 我计算了概率尺度的可能性,然后记录下来,而不是直接使用您给出的公式(但简单的代数会表明它们是相同的)。我没有将对数赔率设置为(负)无穷大,而是将概率设置为 0/1。
  • 您需要提供合适的起始值才能使其发挥作用。特别是,截距应该已经是不同的并且顺序正确(总是根据您想要对类别进行排序的方式增加或减少),从所有它们相等开始会导致问题。
  • 这依赖于
    Y
    是从 1 开始的数字 + 顺序,并且仅处理单个数字
    X
    预测器,正确的拟合例程应该允许更大的灵活性(例如,在内部将
    Y
    的级别转换为此有序序列) .
© www.soinside.com 2019 - 2024. All rights reserved.