beta 二项式和 beta 分布的 alpha 和 beta 估计答案

【问题标题】：alpha and beta estimates for beta binomial and beta distributionsbeta 二项式和 beta 分布的 alpha 和 beta 估计
【发布时间】：2015-04-09 23:06:25
【问题描述】：

我正在尝试将我的数据拟合到 beta 二项分布并估计 alpha 和 beta 形状参数。对于此分布，先验取自 beta 分布。 Python 没有适用于 beta-binomial 的 fit 函数，但它适用于 beta。 python beta 拟合和 R beta 二项式拟合很接近，但系统性地偏离了。

R：

library("VGAM")
x = c(222,909,918,814,970,346,746,419,610,737,201,865,573,188,450,229,629,708,250,508)
y = c(2,18,45,11,41,38,22,7,40,24,34,21,49,35,31,44,20,28,39,17)
fit=vglm(cbind(y, x) ~ 1, betabinomialff, trace = TRUE)
Coef(fit)
   shape1    shape2 
  1.736093 26.870768

蟒蛇：

import scipy.stats
import numpy as np
x = np.array([222,909,918,814,970,346,746,419,610,737,201,865,573,188,450,229,629,708,250,508], dtype=float)
y = np.array([2,18,45,11,41,38,22,7,40,24,34,21,49,35,31,44,20,28,39,17])
scipy.stats.beta.fit((y)/(x+y), floc=0, fscale=1)
    (1.5806623978910086, 24.031893492546242, 0, 1)

我已经这样做了很多次，似乎 python 系统地比 R 的结果要低一点。我想知道这是我的输入错误还是只是计算方式的不同？

【问题讨论】：

标签： python r scipy beta parameterization

【解决方案1】：

您的问题是拟合 Beta 二项式模型与拟合 Beta 模型的值等于比率不同。我将在此处使用bbmle 包进行说明，该包将适合与VGAM 类似的模型（但我更熟悉）。

预赛：

library("VGAM")  ## for dbetabinom.ab
x <- c(222,909,918,814,970,346,746,419,610,737,
       201,865,573,188,450,229,629,708,250,508)
y <- c(2,18,45,11,41,38,22,7,40,24,34,21,49,35,31,44,20,28,39,17)

library("bbmle")

拟合 beta 二项式模型：

mle2(y~dbetabinom.ab(size=x+y,shape1,shape2),
     data=data.frame(x,y),
     start=list(shape1=2,shape2=30))
## Coefficients:
##    shape1    shape2 
##  1.736046 26.871526

这与您引用的VGAM 结果或多或少完全一致。

现在使用相同的框架来拟合 Beta 模型：

mle2(y/(x+y) ~ dbeta(shape1,shape2),
     data=data.frame(x,y),
     start=list(shape1=2,shape2=30))
## Coefficients:
##    shape1    shape2 
## 1.582021 24.060570

这符合您的 Python 测试结果。（我敢肯定，如果您使用 VGAM 来适应 Beta，您也会得到相同的答案。）

【讨论】：

如果有人不确定 bbmle 包是谁编写的（以及为什么 Ben 更熟悉它），他们只需在 R 控制台会话中输入 maintainer('bbmle')。
那么，我想我知道bbmle中的'bb'代表什么了！
如果我知道从长远来看它会很有用/流行的话，我可能会用不同的名字来命名它......

【解决方案2】：

您可以将conjugate_prior 包用于python

请参阅掷硬币示例的代码：

from conjugate_prior import BetaBinomial
heads = 95
tails = 105
prior_model = BetaBinomial() #Uninformative prior
updated_model = prior_model.update(heads, tails)
credible_interval = updated_model.posterior(0.45, 0.55)
print ("There's {p:.2f}% chance that the coin is fair".format(p=credible_interval*100))
predictive = updated_model.predict(50, 50)
print ("The chance of flipping 50 Heads and 50 Tails in 100 trials is {p:.2f}%".format(p=predictive*100))

代码取自here

【讨论】：