Polya's Distribution (Beta-Binomial) Parameter Estimation

This project deals with estimation of the parameters of the Polya’s distribution. The PMF of the observations $x_1,\dots ,x_m$ is given by: $$p(x_1,x_2,\dots ,x_m)=\prod _{i=1}^m(\frac{n(x_i)!}{\prod _k{n}_k(x_i)!}\frac{\mathrm{\Gamma }(\sum _k{\alpha }_k)}{\mathrm{\Gamma }(n(x_i)+\sum _k{\alpha }_k)}\prod _k\frac{\mathrm{\Gamma }(n_k(x_i)+{\alpha }_k)}{\mathrm{\Gamma }({\alpha }_k)})$$

For our problem, we focus on the beta binomial model, with $k\in 1,2$.

Defining parameter vector $\theta$.

In general form: $\theta =[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k{]}^T$

Beta binomial form: $\theta =[{\alpha }_1,{\alpha }_2{]}^T$

Compute the FIM and CRLB.

For our problem, we are given the following assumptions:

• $m=20$
• $n(x_i)=20$ for $1\le i\le m$
• $k\in 1,2$

Note: since $k\in 1,2$, then $\sum _k{\alpha }_k={\alpha }_1+{\alpha }_2$

We now find $\frac{\partial \log p}{\partial \theta }=\left[\begin{array}{c}\frac{\partial \log p}{\partial {\alpha }_1}\frac{\partial \log p}{\partial {\alpha }_2}\end{array}\right]$ from the above log likelihood function.

We now find $\frac{{\partial }^2\log p}{\partial \theta \partial {\theta }^T}=\left[\begin{array}{ccc}\frac{{\partial }^2\log p}{\partial {\alpha }_1\partial {\alpha }_1}& \frac{{\partial }^2\log p}{\partial {\alpha }_1\partial {\alpha }_2}\frac{{\partial }^2\log p}{\partial {\alpha }_2\partial {\alpha }_1}& \frac{{\partial }^2\log p}{\partial {\alpha }_2\partial {\alpha }_2}\end{array}\right]$ from the above first-order partial derivatives.

We now find $FIM=-E[\frac{{\partial }^2\log p}{\partial \theta \partial {\theta }^T}]$ and $CRLB=FI{M}^{-1}$ from the above second-order partial derivatives.

For simplicity and spacing, let $v=m({\psi }^{\prime }({\alpha }_1+{\alpha }_2)-{\psi }^{\prime }(20+{\alpha }_1+{\alpha }_2))$ since it occurs in each partial derivative. $$\begin{array}{rl}FIM& =-E\left[\begin{array}{ccc}\frac{{\partial }^2\log p}{\partial {\alpha }_1\partial {\alpha }_1}& \frac{{\partial }^2\log p}{\partial {\alpha }_1\partial {\alpha }_2}\frac{{\partial }^2\log p}{\partial {\alpha }_2\partial {\alpha }_1}& \frac{{\partial }^2\log p}{\partial {\alpha }_2\partial {\alpha }_2}\end{array}\right]\\ & =-E\left[\begin{array}{cc}v-m{\psi }^{\prime }({\alpha }_1)+\sum _{i=1}^m{\psi }^{\prime }(n_1(x_i)+{\alpha }_1)& v\\ v& v-m{\psi }^{\prime }({\alpha }_2)+\sum _{i=1}^m{\psi }^{\prime }(n_1(x_i)+{\alpha }_2)\end{array}\right]\\ & =\left[\begin{array}{cc}-E[v-m{\psi }^{\prime }({\alpha }_1)+\sum _{i=1}^m{\psi }^{\prime }(n_1(x_i)+{\alpha }_1)]& -E[v]\\ -E[v]& -E[v-m{\psi }^{\prime }({\alpha }_2)+\sum _{i=1}^m{\psi }^{\prime }(n_2(x_i)+{\alpha }_2)]\end{array}\right]\\ CRLB& =FI{M}^{-1}\\ & ={\left[\begin{array}{cc}-E[v-m{\psi }^{\prime }({\alpha }_1)+\sum _{i=1}^m{\psi }^{\prime }(n_1(x_i)+{\alpha }_1)]& -E[v]\\ -E[v]& -E[v-m{\psi }^{\prime }({\alpha }_2)+\sum _{i=1}^m{\psi }^{\prime }(n_2(x_i)+{\alpha }_2)]\end{array}\right]}^{-1}\end{array}$$

Maximum Likelihood Estimation

Recall from the previous section that we derived the partial derivatives of the log likelihood function as follows: $$\begin{array}{rl}\frac{\partial \log p}{\partial {\alpha }_1}& =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)-\psi ({\alpha }_1)\\ & =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2)-\psi ({\alpha }_1))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)\\ \frac{\partial \log p}{\partial {\alpha }_2}& =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)-\psi ({\alpha }_2)\\ & =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2)-\psi ({\alpha }_2))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)\end{array}$$ We can then set each of the derivatives equal to 0 and solve for our ${\alpha }_k$ parameter in order to find the Maximum Likelihood Estimator for each. $$\begin{array}{rl}\frac{\partial \log p}{\partial {\alpha }_1}& =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2)-\psi ({\alpha }_1))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)=0\\ \frac{\partial \log p}{\partial {\alpha }_2}& =m(\psi ({\alpha }_1+{\alpha }_2)-\psi (20+{\alpha }_1+{\alpha }_2)-\psi ({\alpha }_2))+\sum _{i=1}^m\psi (n_1(x_i)+{\alpha }_1)=0\end{array}$$

However, no closed form solution exists for this equation. As proposed in the project guidelines, we used the method of fixed point iteration proposed by Minka et al. in order to update each of the ${\alpha }_k$ values: $${\alpha }_k^{(t)}={\alpha }_k^{(t-1)}\frac{\sum _i\psi (n_{ik}+{\alpha }_k^{(t-1)})-\psi ({\alpha }_k^{(t-1)})}{\sum _i\psi (n_i+\sum _k{\alpha }_k^{(t-1)})-\psi (\sum _k{\alpha }_k^{(t-1)})}$$

Specifically, for our problem of $k\in 1,2$, we implement: $$\begin{array}{rl}{\alpha }_1^{(t)}& ={\alpha }_1^{(t-1)}\frac{\sum _i\psi (n_{i1}+{\alpha }_1^{(t-1)})-\psi ({\alpha }_1^{(t-1)})}{\sum _i\psi (n_i+{\alpha }_1^{(t-1)}+{\alpha }_2^{(t-1)})-\psi ({\alpha }_1^{(t-1)}+{\alpha }_2^{(t-1)})}\\ {\alpha }_2^{(t)}& ={\alpha }_2^{(t-1)}\frac{\sum _i\psi (n_{i2}+{\alpha }_2^{(t-1)})-\psi ({\alpha }_2^{(t-1)})}{\sum _i\psi (n_i+{\alpha }_1^{(t-1)}+{\alpha }_2^{(t-1)})-\psi ({\alpha }_1^{(t-1)}+{\alpha }_2^{(t-1)})}\end{array}$$

Method of Moments

For the beta-binomial distribution, via Wikipedia, the first moment (${\mu }_1$) and second moment (${\mu }_2$) based on $n_1$ are given as follows: $${\mu }_1=\frac{n_1{\alpha }_1}{{\alpha }_1+{\alpha }_2}$$ $${\mu }_2=\frac{n_1{\alpha }_1(n_1(1+{\alpha }_1)+{\alpha }_2)}{({\alpha }_1+{\alpha }_2)({\alpha }_1+{\alpha }_2+1)}$$

For solving the first and second order moments, we set the following: $$\bar{x}=\frac{1}{N}\sum _{i=1}^N{x}_i$$ $$\overline{x^2}=\frac{1}{N}\sum _{i=1}^N{x}_i^2$$

Let $v=\frac{n_1}{\bar{x}}-1$, and $v+1=\frac{n_1}{\bar{x}}$, such that ${\alpha }_2={\alpha }_{1}v$.

We can now plug this back into our first order moment to solve for $\beta$. $$\begin{array}{rl}{\alpha }_2& ={\alpha }_1(\frac{n_1}{\bar{x}}-1)\\ & =(\frac{n_1\bar{x}-\overline{x^2}}{n_1(\frac{\overline{x^2}}{\bar{x}}-\bar{x}-1)+\bar{x}})(\frac{n_1}{\bar{x}}-1)\\ & =\frac{n_1^2-n_1\bar{x}-\frac{n_1\overline{x^2}}{\bar{x}}+\overline{x^2}}{n_1(\frac{\overline{x^2}}{\bar{x}}-\bar{x}-1)+\bar{x}}\\ {\widehat{{\alpha }_2}}_{mom}& =\frac{(n_1-\bar{x})(n_1-\frac{\overline{x^2}}{\bar{x}})}{n_1(\frac{\overline{x^2}}{\bar{x}}-\bar{x}-1)+\bar{x}}\end{array}$$