python - beta 二项分布和 beta 分布的 alpha 和 beta 估计

Question

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

回复:

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

python :

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 的结果低一点。我想知道这是我的输入错误还是它们的计算方式不同？

Answer 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，beta-fit 结果。 (我敢肯定，如果您使用 VGAM 来拟合 Beta，您也会得到相同的答案。)