Finding \alpha and \beta of Beta-binomial model via method of moments

I am looking for a laymen step by step of how the process of finding the 1st and 2nd sample moments located:

http://en.wikipedia.org/wiki/Beta-binomial_distribution#Maximum_likelihood_estimation

Also it's my limited understanding that k-th sample moments are defined as $${\frac {\sum _{i=1}^{n}{x_{{i}}}^{k}}{n}}$$

For samples $x_1, x_2...x_n$ where $n$ = total number of samples. (source: http://en.wikipedia.org/wiki/Moment_(mathematics)#Sample_moments)

Given their example data:

Males       0   1   2   3   4   5     6     7     8    9    10  11  12
Families    3   24  104 286 670 1033  1343  1112  829  478  181 45  7


The first thing I don't understand is why they say $n=12$ when there are 13 data points. Wouldn't that imply $n=13$

I believe the sample moments are:

$m_1 = \frac{0+1+2+3+4+5+6+7+8+9+10+11+12}{13} = 6$

$m_2 = \frac{0^2+1^2+2^2+3^2+4^2+5^2+6^2+7^2+8^2+9^2+10^2+11^2+12^2}{13} = 50$

Yet they have

$$m_1 = 6.23$$ $$m_2=42.31$$

Even If I use $n=12$ and cut off either the first or last record I am left with different values.

Despite that, even using their values of $$m_1 = 6.23$$$$m_2=42.31$$$$n=12$$ going by the equation for the method of moments estimates:

$$\alpha= \frac{( nm_{{1}}-m_{{2}} ) }{n ( {\frac {m_{{2 }}}{m_{{1}}}}-m_{{1}}-1 ) +m_{{1}}} = 33.59257915$$

$$\beta= \frac{( n-m_{{1}} ) ( n-{\frac {m_{{2}}}{m_{{1}}}} )}{n ( {\frac {m_{{2}}}{m_{{1}}}}-m_{{1}}-1 ) +m_{{1}}} = 31.11222820$$

which do not match his values of:

$$\alpha= 34.1350$$ $$\beta = 31.6085$$

Edit: Given this question was spawned from Rating system incorporating experience; For purposes of record keeping for later googlers, I decided to reword this question to better suit the answers. A detail explanatin of Beta-binomial model and the MLE method of finding $\alpha$ and $\beta$ are located there.

-
To begin with, Binomial $(n,p)$ is a distribution on $\{0,1,\ldots,n\}$ hence $n=12$ is the correct parameter and $n+1=13$ is the number of possible results of each experiment. Second, the number of data points is not the number of possible results but the total number of families (in your case, this is 6115 the sum of 3, 24, 104, 286, and so on until 45, 7). – Did Aug 22 '12 at 8:02
OK, so how would I calculate the sample moments? – ParoX Aug 22 '12 at 8:52
For example the first moment is 0x3+1x24+2x104+3x286+...+11x45+12x7 divided by 6115. – Did Aug 22 '12 at 10:46
That's much closer to what he got, but it's still off by a little bit: i.minus.com/iblbDxkM2OlM9j.png (using maple you can see the values are ~6.22 and ~40.39 – ParoX Aug 22 '12 at 11:07
I added a more detailed explanation of the computations in the original question where the beta-binomial model came up: math.stackexchange.com/questions/183483/…. You should probably refer to this. – Einar Rødland Aug 23 '12 at 2:58

The moments should be $$m_k = \frac{ \sum_{i=0}^{12} f_i \times i^k}{\sum_{i=0}^{12} f_i}$$ where $f_i$ is the number of families with $i$ males.

The calculation of $\hat{\alpha}$ and $\hat{\beta}$ require the use of $n=12$, as @did says.

Do both and you will get the stated values.

-
Thank you for clarifying that to me. It should be noted that since I am using arrays that $f_0$ does not exist an returns and error and instead you should go i=1 to 13. I do have a question though, in this example you have $f_i \times i^k$ instead of my proposed $f_i * x_i^k$,where $x_1, x_2..x_n$ would be the number of males. Had the values of males not be a perfect count from 0..12 would you still use the counter (i) ? – ParoX Aug 22 '12 at 12:31
@BHare: Take the simple case of the mean number of males per family, which is $m_1$. You add up the number of males $3 \times 0 + 24 \times 1 + \cdots + 7 \times 12$ and divide by the number of families $3 + 24 + \cdots + 7$. If you can do that then the higher moments are similar. – Henry Aug 22 '12 at 17:07

The method of estimation that you are describing is called method of moments. It is not maximum likelihood estimation. To do maximum likelihood you have to write down the likelihood function for your observed data based on the parametric model. Then you search for a maximum value for that function (which is often unique). When the partial derivatives of the likelihood with respect to the parameters can be computed you take them and solve the equation for the local extrema which should turn out to be the global maximum.

-
Yes, this is what I had thought. It almost seems as if he derived the MLE from the method of moments. Can you go into detail on how to find this likelihood function with the chosen data, and finding the extrema? – ParoX Aug 22 '12 at 12:39
Suppose you have an unknown parameter or parameter vector θ and a set of observations X$_1$, X$_2$, ...,X$_n$ that are indeoendent and identically distributed with probability density f$_θ$(x) then the likelihood function is f$_θ$(X$_1$)f$_θ$(X$_2$)...f$_θ$(X$_n$). You can then usually maximize this likelihood by taking the partial derivates with respect to the components of θ and setting them equal to 0 and solving the system of equations that are obtained. – Michael Chernick Aug 22 '12 at 17:49

Just for historical purposes, here is how I got the $\alpha$ and $\beta$ using Maple 15:

Maple code:

males := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
fams := [3, 24, 104, 286, 670, 1033, 1343, 1112, 829, 478, 181, 45, 7];
n := 12;

k := 1;
m[1] := (sum(fams[i]*males[i]^k, i = 1 .. n+1))/(sum(fams[i], i = 1 .. n+1));

k := 2;
m[2] := (sum(fams[i]*males[i]^k, i = 1 .. n+1))/(sum(fams[i], i = 1 .. n+1));

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]);

printf("

n = %d
m_1 = %f
m_2 = %f
alpha = %f
beta = %f

", n, m[1], m[2], alpha, beta)


Output:

n = 12
m_1 = 6.230581
m_2 = 42.309403
alpha = 34.135021
beta = 31.608492

-