0
$\begingroup$

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.

  • 0
    For better help on how to run the analyses, e.g. R routines that do the analyses, the statistics forum http://stats.stackexchange.com/ might be better suited.2012-08-23

3 Answers 3

2

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.

  • 0
    @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.2012-08-22
0

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.

  • 0
    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.2012-08-22
0

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

enter image description here

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