3

I am trying to fit my data to a beta-binomial distribution and estimate the alpha and beta shape parameters. For this distribution, the prior is taken from a beta distribution. Python does not have a fit function for beta-binomial but it does for beta. The python beta fitting and R beta binomial fitting is close but systematically off.

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

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)

I have done this many times and it seems like python is systematically a little bit lower than the R results. I was wondering if this is an input error on my part or just a difference in the way they are calculated?

4

2 に答える 2

5

あなたの問題は、ベータ二項モデルをフィッティングすることは、比率に等しい値でベータモデルをフィッティングすることと同じではないということです。ここでは、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")

ベータ二項モデルを当てはめる：

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これは、引用した結果とほぼ完全に一致します。

代わりに、同じフレームワークを使用してベータ モデルに適合させます。

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ベータ版に慣れていれば、同じ答えが得られるはずです。)

于 2015-02-10T02:14:27.750 に答える