Observed data: Count vector $y : (y_1,y_2,y_3,y_4)^T$ postulated to be arising from a multinomial distribution with four cells and corresponding probabilities: $(p_1,p_2,p_3,p_4) = (\frac{1}{2}+\frac{1}{4}p, \frac{1}{4}(1-p), \frac{1}{4}(1-p), \frac{1}{4}p)$ such that $0 \leq p \leq 1$ and we need to estimate $p$

Probability density function $f(y;p) = \frac{n!}{y_1!y_2!y_3!y_4!}p_1^{y_1}p_2^{y_2}p_3^{y_3}p_4^{y_4}$

The log likelihood function is given by $\log(L(p|y)) = y_1\log(p_1) + y_2\log(p_2)+y_3\log(p_3)+y_4\log(p_4) + \log(n!) - \log(y_1!y_2!y_3!y_4!)$

MLE

\begin{eqnarray} \frac{\partial \log(L(p|y))}{\partial p} = 0\\ \frac{y_1}{p_1}(\frac{1}{4}) + \frac{y_2}{p_2}(-\frac{1}{4}) + \frac{y_3}{p_3}(-\frac{1}{4}) + \frac{y_4}{p_4}(\frac{1}{4}) = 0\\ \frac{y_1}{2+p} - \frac{y_2+y_3}{1-p} + \frac{y_4}{p} = 0\\ \end{eqnarray}

Analytical Solution

\begin{eqnarray} b = y_1-y_4-2(y_2+y_3)\\ a = -(y_1+y_2+y_3+y_4)\\ c = 2y_4\\ \hat{p} = \frac{-b\pm\sqrt{b^2-4ac}}{2a} \end{eqnarray}



In [15]:

    
import numpy as np
import scipy as sp
y1,y2,y3,y4 = 1997,906,904,32
n = 3839
b = y1-y4-2*(y2+y3)
a = -(y1+y2+y3+y4)
c = 2*y4
p = (-b-np.sqrt(b**2-4*a*c))/(2*a)



In [16]:

    
p









    Out[16]:





0.035712302240628144

MLE - Newton Raphson

\begin{eqnarray} S(y; p) = \frac{\partial \log(L(p|y))}{\partial p}\\ S(y;p) = 0\\ I(p;y) = -\frac{\partial^2 \log(L(p|y))}{\partial p^2}\\ I(p;y) = \frac{y_1}{(2+p)^2} + \frac{y_2+y_3}{(1-p)^2} + \frac{y_4}{p^2} \\ S(y;p) \approx S(y;p^k) + \frac{\partial S}{\partial p} I(p;y)(p-p^k) = S(y;p^k) - I(p;y)(p-p^k) \text{Taylor expansion}\\ p^{k+1}=p^k+I^{-1}(p^k;y)S(y;p^k) \end{eqnarray}

$$E[y_1+y_4-(y_2+y_3)] = np$$

Thus we choose starting point to be $p_0 = \frac{y_1+y_4-(y_2+y_3)}{n}$

NOTE: In order to reduce computational complexity we make this following inherent assumption: $$ E[\frac{\partial^2 S}{\partial p^2}] = E[\frac{\partial S}{\partial p} \times \frac{\partial S}{\partial p}] $$



In [101]:

    
pe = (y1+y4-(y2+y3))/n
pe









    Out[101]:





0.057046105756707474



In [134]:

    
pold = pe
pnew = pe
s = y1/(2+pe)-(y2+y3)/(1-pe)+y4/pe
print(s)









    



-387.74068038090877



In [135]:

    
def I(px):
    return y1/((2+px)**2)+(y2+y3)/((1-px)**2) + y4/(px**2)
def S(px):
    return y1/(2+px)-(y2+y3)/(1-px)+y4/px
iters = 0
while np.abs(s)>10e-9 and pnew>0 and pold>0:
    i = I(pnew)
    s = S(pnew) - i *(pnew-pold)
    pold = pnew
    pnew = pold + 1/(I(pold))*S(pold)
    print(pold, pnew, S(pold))
    iters+=1









    



0.057046105756707474 0.025626790824640956 -387.74068038090877
0.025626790824640956 0.03300085259227735 376.956468896675
0.03300085259227735 0.03552250304276715 80.19367816684894
0.03552250304276715 0.035711383960063195 5.248507070642177
0.035711383960063195 0.03571230221915427 0.025270958459486792
0.03571230221915427 0.03571230224062816 5.90944637224311e-07
0.03571230224062816 0.035712302240628144 -3.410605131648481e-13
0.035712302240628144 0.035712302240628144 0.0



In [ ]: