R估计二项分布的参数

Question

我正在尝试根据 R 中的最大似然估计二项分布的参数 n 和 p。

我正在使用 stats 包中的函数 optim，但出现错误。

这是我的代码:

xi = rbinom(100, 20, 0.5) # Sample
n = length(xi) # Sample size

# Log-Likelihood
lnlike <- function(theta){
log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
(n*theta[1] - sum(xi))*log(1-theta[2])
}

# Optimizing 
optim(theta <- c(10,.3), lnlike, hessian=TRUE)

Error in optim(theta <- c(10, 0.3), lnlike, hessian = TRUE) : function cannot be evaluated at initial parameters

有人做过吗？使用了哪个函数？

Answer 1

tl;dr 如果响应变量大于二项式 N(这是理论上的 < em>响应的最大值)。在大多数实际问题中，N 被认为是已知的，只是估计概率。如果确实要估计 N，则需要 (1) 将其约束为 >= 样本中的最大值； (2) 做一些特殊的事情来优化必须离散的参数(这是一个高级/棘手的问题)。

此答案的第一部分展示了用于识别问题的调试策略，第二部分展示了同时优化 N 和 p 的策略(通过在合理的 N 范围内的蛮力)。

设置:

set.seed(101)
n <- 100
xi <- rbinom(n, size=20, prob=0.5) # Sample

对数似然函数:

lnlike <- function(theta){
    log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
       (n*theta[1] - sum(xi))*log(1-theta[2])
}

让我们分解一下。

theta <- c(10,0.3)  ## starting values
lnlike(c(10,0.3))  ## -Inf

好的，初始值的对数似然为 -Inf。 optim() 无法使用它也就不足为奇了。

让我们研究一下这些条款。

log(prod(choose(theta[1],xi))) ## -Inf

好吧，我们在第一学期就遇到了麻烦。

prod(choose(theta[1],xi)) ## 0

乘积为零……为什么？

choose(theta[1],xi)
##  [1] 120 210  10   0   0  10 120 210   0   0  45 210   1   0

很多零。为什么？有问题的 xi 值是什么？

## [1]  7  6  9 12 11  9  7  6

啊哈!我们对 7、6、9 没问题……但对 12 有麻烦。

badvals <- (choose(theta[1],xi)==0)
all(badvals==(xi>10))  ## TRUE

如果你真的想这样做，你可以通过对 n 的合理值进行暴力枚举来实现......

## likelihood function
llik2 <- function(p,n) {
    -sum(dbinom(xi,prob=p,size=n,log=TRUE))
}
## possible N values (15 to 50)
nvec <- max(xi):50
Lvec <- numeric(length(nvec))
for (i in 1:length(nvec)) {
    ## optim() wants method="Brent"/lower/upper for 1-D optimization
    Lvec[i] <- optim(par=0.5,fn=llik2,n=nvec[i],method="Brent",
                     lower=0.001,upper=0.999)$val
}
nvec[which.min(Lvec)]  ## 20
par(las=1,bty="l")
plot(nvec,Lvec,type="b")