0

I am attempting to work through the method of moments estimation of alpha and beta for a beta binomial distribution. Taking the steps found at: http://en.wikipedia.org/wiki/Beta-binomial_distribution#Maximum_likelihood_estimation

I have think I have been able to code this accurately in python:

import numpy as np

males = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], dtype=float)
fams = np.array([3, 24, 104, 286, 670, 1033, 1343, 1112, 829, 478, 181, 45, 7], dtype=float)
n = 12

k = 1
m_1 = (sum(fams*males**k))/(sum(fams)) 

k = 2
m_2 = (sum(fams*males**k))/(sum(fams)) 

alpha = (n*m_1-m_2)/(n*(m_2/m_1-m_1-1)+m_1)
beta = (n-m_1)*(n-m_2/m_1)/(n*(m_2/m_1-m_1-1)+m_1)

print "n =", n
print "m_1 =", m_1
print "m_2 =", m_2
print "alpha =", alpha
print "beta =", beta

which outputs:

n = 12
m_1 = 6.23058053966
m_2 = 42.3094031071
alpha = 34.135021177
beta = 31.6084920506

This are the same results as the example. However, if I use the R package VGAM which uses maximum likelihood the estimates for alpha and beta are completely different

x = c(0,1,2,3,4,5,6,7,8,9,10,11,12)
y = c(3,24,104,286,670,1033,1343,1112,829,478,181,45,7)
library("VGAM")
fit=vglm(cbind(x, y) ~ 1, betabinomialff, trace = TRUE)
Coef(fit)
  shape1    shape2 
0.4241806 4.9069560

Am I doing something wrong?

user3266890
  • 465
  • 5
  • 15
  • 2
    Method of moments and maximum likelihood are two completely different methods of parameter estimation, you are comparing apples to oranges. – mlegge Feb 10 '15 at 18:51
  • OP Is comparing to the Wikipedia estimates via mle. clearly these do not match. likely to be some r convention OP is not getting right. See if you can correct that. – pythOnometrist Aug 27 '15 at 14:55

0 Answers0