Beta-binomial distribution

The first three raw moments are

${\displaystyle {\begin{aligned}\mu _{1}&={\frac {n\alpha }{\alpha +\beta }}\\[8pt]\mu _{2}&={\frac {n\alpha [n(1+\alpha )+\beta ]}{(\alpha +\beta )(1+\alpha +\beta )}}\\[8pt]\mu _{3}&={\frac {n\alpha [n^{2}(1+\alpha )(2+\alpha )+3n(1+\alpha )\beta +\beta (\beta -\alpha )]}{(\alpha +\beta )(1+\alpha +\beta )(2+\alpha +\beta )}}\end{aligned}}}$

and the kurtosis is

${\displaystyle \beta _{2}={\frac {(\alpha +\beta )^{2}(1+\alpha +\beta )}{n\alpha \beta (\alpha +\beta +2)(\alpha +\beta +3)(\alpha +\beta +n)}}\left[(\alpha +\beta )(\alpha +\beta -1+6n)+3\alpha \beta (n-2)+6n^{2}-{\frac {3\alpha \beta n(6-n)}{\alpha +\beta }}-{\frac {18\alpha \beta n^{2}}{(\alpha +\beta )^{2}}}\right].}$

Letting ${\displaystyle \pi ={\frac {\alpha }{\alpha +\beta }}\!}$ we note, suggestively, that the mean can be written as

${\displaystyle \mu ={\frac {n\alpha }{\alpha +\beta }}=n\pi \!}$

and the variance as

${\displaystyle \sigma ^{2}={\frac {n\alpha \beta (\alpha +\beta +n)}{(\alpha +\beta )^{2}(\alpha +\beta +1)}}=n\pi (1-\pi ){\frac {\alpha +\beta +n}{\alpha +\beta +1}}=n\pi (1-\pi )[1+(n-1)\rho ]\!}$

where ${\displaystyle \rho ={\tfrac {1}{\alpha +\beta +1}}\!}$. The parameter ${\displaystyle \rho \!}$ is known as the "intra class" or "intra cluster" correlation. It is this positive correlation which gives rise to overdispersion.

It is convenient to reparameterize the distributions so that the expected mean of the prior is a single parameter: Let

${\displaystyle {\begin{aligned}\pi (\theta \mid \mu ,M)&=\operatorname {Beta} (M\mu ,M(1-\mu ))\\[6pt]&={\frac {\Gamma (M)}{\Gamma (M\mu )\Gamma (M(1-\mu ))}}\theta ^{M\mu -1}(1-\theta )^{M(1-\mu )-1}\end{aligned}}}$

where

${\displaystyle {\begin{aligned}\mu &={\frac {\alpha }{\alpha +\beta }}\\[6pt]M&=\alpha +\beta \end{aligned}}}$

so that

${\displaystyle {\begin{aligned}\operatorname {E} (\theta \mid \mu ,M)&=\mu \\[6pt]\operatorname {Var} (\theta \mid \mu ,M)&={\frac {\mu (1-\mu )}{M+1}}.\end{aligned}}}$

The posterior distribution ρ(θ | k) is also a beta distribution:

${\displaystyle {\begin{aligned}\rho (\theta \mid k)&\propto \ell (k\mid \theta )\pi (\theta \mid \mu ,M)\\[6pt]&=\operatorname {Beta} (k+M\mu ,n-k+M(1-\mu ))\\[6pt]&={\frac {\Gamma (M)}{\Gamma (M\mu )\Gamma (M(1-\mu ))}}{n \choose k}\theta ^{k+M\mu -1}(1-\theta )^{n-k+M(1-\mu )-1}\end{aligned}}}$

And

${\displaystyle \operatorname {E} (\theta \mid k)={\frac {k+M\mu }{n+M}}.}$

while the marginal distribution m(k|μ, M) is given by

${\displaystyle {\begin{aligned}m(k\mid \mu ,M)&=\int _{0}^{1}\ell (k\mid \theta )\pi (\theta \mid \mu ,M)\,d\theta \\[6pt]&={\frac {\Gamma (M)}{\Gamma (M\mu )\Gamma (M(1-\mu ))}}{n \choose k}\int _{0}^{1}\theta ^{k+M\mu -1}(1-\theta )^{n-k+M(1-\mu )-1}\,d\theta \\[6pt]&={\frac {\Gamma (M)}{\Gamma (M\mu )\Gamma (M(1-\mu ))}}{n \choose k}{\frac {\Gamma (k+M\mu )\Gamma (n-k+M(1-\mu ))}{\Gamma (n+M)}}.\end{aligned}}}$

Substituting back M and μ, in terms of ${\displaystyle \alpha }$ and ${\displaystyle \beta }$, this becomes:

${\displaystyle m(k\mid \alpha ,\beta )={\frac {\Gamma (n+1)}{\Gamma (k+1)\Gamma (n-k+1)}}{\frac {\Gamma (k+\alpha )\Gamma (n-k+\beta )}{\Gamma (n+\alpha +\beta )}}{\frac {\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}}.}$

which is the expected beta-binomial distribution with parameters ${\displaystyle n,\alpha }$ and ${\displaystyle \beta }$.

We can also use the method of iterated expectations to find the expected value of the marginal moments. Let us write our model as a two-stage compound sampling model. Let k_i be the number of success out of n_i trials for event i:

${\displaystyle {\begin{aligned}k_{i}&\sim \operatorname {Bin} (n_{i},\theta _{i})\\[6pt]\theta _{i}&\sim \operatorname {Beta} (\mu ,M),\ \mathrm {i.i.d.} \end{aligned}}}$

We can find iterated moment estimates for the mean and variance using the moments for the distributions in the two-stage model:

${\displaystyle \operatorname {E} \left({\frac {k}{n}}\right)=\operatorname {E} \left[\operatorname {E} \left(\left.{\frac {k}{n}}\right|\theta \right)\right]=\operatorname {E} (\theta )=\mu }$

${\displaystyle {\begin{aligned}\operatorname {var} \left({\frac {k}{n}}\right)&=\operatorname {E} \left[\operatorname {var} \left(\left.{\frac {k}{n}}\right|\theta \right)\right]+\operatorname {var} \left[\operatorname {E} \left(\left.{\frac {k}{n}}\right|\theta \right)\right]\\[6pt]&=\operatorname {E} \left[\left(\left.{\frac {1}{n}}\right)\theta (1-\theta )\right|\mu ,M\right]+\operatorname {var} \left(\theta \mid \mu ,M\right)\\[6pt]&={\frac {1}{n}}\left(\mu (1-\mu )\right)+{\frac {n-1}{n}}{\frac {(\mu (1-\mu ))}{M+1}}\\[6pt]&={\frac {\mu (1-\mu )}{n}}\left(1+{\frac {n-1}{M+1}}\right).\end{aligned}}}$

(Here we have used the law of total expectation and the law of total variance.)

We want point estimates for ${\displaystyle \mu }$ and ${\displaystyle M}$. The estimated mean ${\displaystyle {\widehat {\mu }}}$ is calculated from the sample

${\displaystyle {\widehat {\mu }}={\frac {\sum _{i=1}^{N}k_{i}}{\sum _{i=1}^{N}n_{i}}}.}$

The estimate of the hyperparameter M is obtained using the moment estimates for the variance of the two-stage model:

${\displaystyle s^{2}={\frac {1}{N}}\sum _{i=1}^{N}\operatorname {var} \left({\frac {k_{i}}{n_{i}}}\right)={\frac {1}{N}}\sum _{i=1}^{N}{\frac {{\widehat {\mu }}(1-{\widehat {\mu }})}{n_{i}}}\left[1+{\frac {n_{i}-1}{{\widehat {M}}+1}}\right]}$

Solving:

${\displaystyle {\widehat {M}}={\frac {{\widehat {\mu }}(1-{\widehat {\mu }})-s^{2}}{s^{2}-{\frac {{\widehat {\mu }}(1-{\widehat {\mu }})}{N}}\sum _{i=1}^{N}1/n_{i}}},}$

where

${\displaystyle s^{2}={\frac {N\sum _{i=1}^{N}n_{i}({\widehat {\theta _{i}}}-{\widehat {\mu }})^{2}}{(N-1)\sum _{i=1}^{N}n_{i}}}.}$

Since we now have parameter point estimates, ${\displaystyle {\widehat {\mu }}}$ and ${\displaystyle {\widehat {M}}}$, for the underlying distribution, we would like to find a point estimate ${\displaystyle {\tilde {\theta }}_{i}}$ for the probability of success for event i. This is the weighted average of the event estimate ${\displaystyle {\widehat {\theta _{i}}}=k_{i}/n_{i}}$ and ${\displaystyle {\widehat {\mu }}}$. Given our point estimates for the prior, we may now plug in these values to find a point estimate for the posterior

${\displaystyle {\tilde {\theta _{i}}}=\operatorname {E} (\theta \mid k_{i})={\frac {k_{i}+{\widehat {M}}{\widehat {\mu }}}{n_{i}+{\widehat {M}}}}={\frac {\widehat {M}}{n_{i}+{\widehat {M}}}}{\widehat {\mu }}+{\frac {n_{i}}{n_{i}+{\widehat {M}}}}{\frac {k_{i}}{n_{i}}}.}$

We may write the posterior estimate as a weighted average:

${\displaystyle {\tilde {\theta }}_{i}={\widehat {B}}_{i}\,{\widehat {\mu }}+(1-{\widehat {B}}_{i}){\widehat {\theta }}_{i}}$

where ${\displaystyle {\widehat {B}}_{i}}$ is called the shrinkage factor.

${\displaystyle {\widehat {B_{i}}}={\frac {\widehat {M}}{{\widehat {M}}+n_{i}}}}$

• ${\displaystyle BB(1,1,n)\sim U(0,n)\,}$ where ${\displaystyle U(a,b)\,}$ is the discrete uniform distribution.

• Dirichlet-multinomial distribution

• Minka, Thomas P. (2003). Estimating a Dirichlet distribution. Microsoft Technical Report.

• Using the Beta-binomial distribution to assess performance of a biometric identification device
• Fastfit contains Matlab code for fitting Beta-Binomial distributions (in the form of two-dimensional Pólya distributions) to data.
• Interactive graphic: Univariate Distribution Relationships
• Beta-binomial functions in VGAM R package
• Beta-binomial distribution in Sandia National Labs Cognitive Foundry Java library