还给出 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)
只有这样不行,我该如何解决这个问题?
我认为第一个问题是您在内循环中
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
一些注意事项:
Y
是从 1 开始的数字 + 顺序,并且仅处理单个数字 X
预测器,正确的拟合例程应该允许更大的灵活性(例如,在内部将 Y
的级别转换为此有序序列) .