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)

优化错误（theta

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

Answer 1

tl;dr 如果响应变量大于二项式 N（这是理论 响应的最大值）。在大多数实际问题中，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")

Answer 2

你为什么会惹上麻烦？

如果你执行lnlike(c(10, 0.3))，你会得到-Inf。这就是为什么您的错误消息是在抱怨lnlike，而不是optim。

通常，n 是已知的，并且只需要估计 p。在这种情况下，矩估计器或最大似然估计器都是封闭形式，不需要数值优化。所以，估计n真的很奇怪。

如果您确实想估计，您必须意识到它是受约束的。检查

range(xi)  ## 5  15

您的观测值范围为[5, 15]，因此，需要n >= 15。你怎么能传递一个初始值 10？ n的搜索方向，应该从一个大的起始值开始，然后逐渐向下搜索，直到到达max(xi)。所以，你可以尝试30作为@987654334的初始值@。

另外，您不需要以当前方式定义lnlike。这样做：

lnlike <- function(theta, x) -sum(dbinom(x, size = theta[1], prob = theta[2], log = TRUE))

• optim 通常用于最小化（尽管它可以最大化）。我在函数中放了一个减号以获得负对数似然。通过这种方式，您可以最小化 lnlike w.r.t。 theta。
• 您还应该将您的观察xi 作为附加参数传递给lnlike，而不是从全局环境中获取。

天真地尝试optim：

在我的评论中，我已经说过我不相信使用optim 来估计n 会起作用，因为n 必须是整数，而optim 用于连续变量。这些错误和警告会让您信服。

optim(c(30,.3), fn = lnlike, x = xi, hessian = TRUE)

Error in optim(c(30, 0.3), fn = lnlike, x = xi, hessian = TRUE) : 
  non-finite finite-difference value [1]
In addition: There were 15 or more warnings (use warnings() to see the 
first 15

> warnings()
Warning messages:
1: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced
2: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced
3: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced
4: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced
5: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced

解决方案？

Ben 为您提供了一种方法。我们没有让optim 估计n，而是手动对n 进行网格搜索。对于每个候选 n，我们执行 w.r.t 的单变量优化。 p。 （哎呀，其实这里不需要做数值优化。）这样，你得到的profile可能性为n。然后，我们在网格上找到n 以最小化这个轮廓的可能性。

Ben 已为您提供了完整的详细信息，我将不再重复。干得好（而且很快），本！