Beta-binomial model

In empirical Bayes methods, the Beta-binomial model is an analytic model where the likelihood function $L(x|\; heta)$ is specifed by a binomial distribution

:$L(x|\; heta)\; =\; operatorname\{Bin\}(x,\; heta),$

:::$=\; \{nchoose\; x\}\; heta^x(1-\; heta)^\{n-x\},$

and the conjugate prior is a Beta distribution

:

:::

The Beta-binomial is a two-dimensional multivariate Polya distribution, as the binomial and Beta distributions are special cases of the multinomial and Dirichlet distributions, respectively.

Derivation of the posterior and the marginal

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

:$pi(\; heta|mu,M)\; =\; Beta(mu,M),$

:::$=\; frac\{Gamma\; (M)\}\{Gamma(mu\; M)Gamma(M(1-mu))\}\; heta^\{Mmu-1\}(1-\; heta)^\{M(1-mu)-1\}$

where

::: and

so that

:::$E(\; heta|mu,M)\; =\; mu,$ :::$Var(\; heta|mu,M)\; =\; frac\{mu(1-mu)\}\{M+1\}$

The posterior distribution $ho(\; heta|x)$ is also a beta distribution

:$ho(\; heta|x)\; propto\; l(x|\; heta)pi(\; heta|mu,M)$

:::$=\; Beta(x+M\; mu\; ,\; n-x+M(1-\; mu)\; )\; ,$

:::$=\; frac\{Gamma\; (M)\; \}\{Gamma(Mmu)Gamma(M(1-mu))\}\{nchoose\; x\}\; heta^\{x+Mmu-1\}(1-\; heta)^\{n-x+M(1-mu)-1\},$

while the marginal distribution $m(x|mu,\; M)$ is given by

:$m(x|mu,M)\; =\; int\_\{0\}^\{1\}\; l(x|\; heta)pi(\; heta|mu,\; M)\; d\; heta$

::= $frac\{Gamma\; (M)\; \}\{Gamma(Mmu)Gamma(M(1-mu))\}\{nchoose\; x\}\; int\_\{0\}^\{1\}\; heta^\{x+Mmu-1\}(1-\; heta)^\{n-x+M(1-mu)-1\}\; d\; heta$

::= $frac\{Gamma\; (M)\; \}\{Gamma(Mmu)Gamma(M(1-mu))\}\{nchoose\; x\}\; frac\{Gamma\; (x+Mmu)Gamma(n-x+M(1-mu))\; \}\{Gamma(n+M)\}$

Moment estimates

Because the marginal is a complex, non-linear function of Gamma and Digamma functions, it is quite difficult to obtain a marginal maximum likelihood estimate (MMLE) for the mean and variance. Instead, we use the method of iterated expectations to find the expected value of the marginal moments.

Let us write our model as a two-stage compund sampling model (for each event "i" out of "n_i" possible events):

::$X\_\{i\}|\; heta\_\{i\}\; sim\; mathrm\{Bin\}(n\_\{i\},\; heta\_\{i\}),$

::$heta\_\{i\}\; sim\; mathrm\{Beta\}|(mu,M),\; i.i.d,$

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

::$Eleft(frac\{X\}\{n\}\; ight)\; =\; Eleft\; [Eleft(frac\{X\}\{n\}|\; heta\; ight)\; ight]\; =\; E(\; heta)\; =\; mu$

::$mathrm\{var\}left(frac\{X\}\{n\}\; ight)\; =\; Eleft\; [mathrm\{var\}left(frac\{X\}\{n\}|\; heta\; ight)\; ight]\; +\; mathrm\{var\}left\; [Eleft(frac\{X\}\{n\}|\; heta\; ight)\; ight]$

:::$=\; Eleft\; [left(frac\{1\}\{n\}\; ight)\; heta(1-\; heta)|mu,M\; ight]\; +\; mathrm\{var\}left(\; heta|mu,M\; ight)$

:::$=\; frac\{1\}\{n\}left(mu(1-mu)\; ight)+\; frac\{n\_\{i\}-1\}\{n\_\{ifrac\{(mu(1-mu))\}\{M+1\}$

:::$=\; frac\{mu(1-mu)\}\{n\}left(1+frac\{n-1\}\{M+1\}\; ight).$

We now seek a point estimate $ilde\{\; heta\}$ as a weighted average of the sample estimate $hat\{\; heta\_\{i$ and an estimate for $hat\{mu\}$. The sample estimate $hat\{\; heta\_\{i$ is given by

::$hat\{\; heta\_\{i\; =\; frac\{x\_\{i\{n\_\{i,$

Therefore we need point estimates for $mu$ and $M$. The estimated mean $hat\{mu\}$ is given as a weighted average

::$hat\{mu\}\; =\; frac\{sum\_\{i=1\}^\{N\}\; n\_\{i\}\; hat\{\; heta\_\{i\; \}\{sum\_\{i=1\}^\{N\}\; n\_\{i\}\; \}\; =\; frac\{sum\_\{i=1\}^\{N\}\; x\_\{i\}\; \}\{sum\_\{i=1\}^\{N\}\; n\_\{i\}\; \}$

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

::$s^\{2\}\; =\; frac\{1\}\{N\}\; sum\_\{i=1\}^\{N\}\; mathrm\{var\}left(frac\{x\_\{i\{n\_\{i\; ight)\; =\; frac\{1\}\{N\}\; sum\_\{i=1\}^\{N\}\; frac\{hat\{mu\}(1-hat\{mu\})\}\{n\_\{i\; left\; [1+frac\{n\_\{i\}-1\}\{hat\{M\}+1\}\; ight]\; =$

::$s^\{2\}\; =\; frac\{N\; sum\_\{i=1\}^\{N\}\; n\_\{i\}\; (hat\{\; heta\_\{i\; -\; hat\{mu\})^\{2\}\; \}\{(N-1)sum\_\{i=1\}^\{N\}\; n\_\{i\}\; \}$

We can now solve for $hat\{M\}$:

::$hat\{M\}\; =\; frac\{hat\{mu\}(1-hat\{mu\})\; -\; s^\{2\{s^\{2\}\; -\; frac\{hat\{mu\}(1\; -\; hat\{mu\})\}\{N\}sum\_\{i=1\}^\{N\}\; 1/n\_\{i\}\; \}.$

Given our point estimates for the prior, we may now plug in these values to find a point estimate for the posterior

:$ilde\{\; heta\_\{i\; =\; E(Theta|hat\{mu\},hat\{M\})\; =\; frac\{x\_\{i\}\; +\; hat\{M\}hat\{mu\{n\_\{i\}+hat\{M\; =\; frac\{hat\{M\{n\_\{i\}+hat\{M\_\{i\}hat\{mu\}\; +\; frac\{n\_\{i\{n\_\{i\}+hat\{M\_\{i\}frac\{x\_\{i\{n\_\{i$

hrinkage factors

We may write the posterior estimate as a weighted average:

::$ilde\{\; heta\_\{i\; =\; hat\{B\_\{ihat\{mu\}\; +\; (1-hat\{B\_\{i)hat\{\; heta\_\{i$

where $hat\{B\_\{i$ is called the "Shrinkage Factor".

::$hat\{B\_\{i\; =\; frac\{hat\{M\_\{i\}\{hat\{M\_\{i+n\_\{i$

Example

Maximum likelihood estimation

Maximum likelihood estimates from empirical data can be computed using general methods for fitting multinomial Polya distributions, methods for which are described in (Minka 2003).

Improved estimates: James-Stein estimator

ee also

* multivariate Polya distribution

External links

* [http://www.cs.ubc.ca/~murphyk/Teaching/Stat406-Spring07/reading/ebHandout.pdf Empirical Bayes for Beta-Binomial model]
* [http://it.stlawu.edu/~msch/biometrics/papers.htm Using the Beta-binomial distribution to assess performance of a biometric identification device]
* [http://www.emse.fr/g2i/publications/rapports/RR_2005-500-012.pdf Extended Beta-Binomial Model for Demand Forecasting of Multiple Slow-Moving Items with Low Consumption and Short Request History]
* [http://research.microsoft.com/~minka/software/fastfit/ Fastfit] contains Matlab code for fitting Beta-Binomial distributions (in the form of two-dimensional Polya distributions) to data.

References

* Minka, Thomas P. (2003). [http://research.microsoft.com/~minka/papers/dirichlet/ Estimating a Dirichlet distribution] . Microsoft Technical Report.