1 Week 6: EM Solutions

Write down the likelihood $L(p;𝐱)$ of $𝐱=(x_1,x_2,x_3)$ given $p$, where $x_i$ denotes the total number of households with $i$ people infected.

The data is multinomial, in that, there are three possible outcomes for each household and the outcome in each household are independent and identically distributed. Hence,

	$L(p;𝐱)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}}{({(1-p)}^2)}^{x_1}{(2p{(1-p)}^2)}^{x_2}{\{2p^2(1-p)+p^2\}}^{x_3}$
		$=$	${\displaystyle \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}}{({(1-p)}^2)}^{x_1}{(2p{(1-p)}^2)}^{x_2}{\{p^2(3-2p)\}}^{x_3}$
		$=$	${\displaystyle \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}}{2}^{x_2}{(1-p)}^{2x_1+2x_2}{p}^{x_2+2x_3}{(3-2p)}^{x_3}$

The log-likelihood is

$$l(p;𝐱)=K+2(x_1+x_2)\log (1-p)+(x_2+2x_3)\log p+x_3\log (3-2p).$$

(Leave this question to the end, if you are short on time.) Find $\frac{d}{dp}l(p;𝐱)=\frac{d}{dp}\log L(p;𝐱)$. Solve $\frac{d}{dp}l(p;𝐱)=0$ to find the MLE, $\widehat{p}$.

This forms a useful check for the EM algorithm and should hopefully convince you of the advantages of using the EM algorithm for this problem.

	${\displaystyle \frac{d}{dp}}l(p;𝐱)$	$=$	${\displaystyle \frac{-2(x_1+x_2)}{1-p}}+{\displaystyle \frac{x_2+2x_3}{p}}-{\displaystyle \frac{2x_3}{3-2p}}$
		$=$	${\displaystyle \frac{-118}{1-p}}+{\displaystyle \frac{575}{p}}-{\displaystyle \frac{550}{3-2p}}$

Setting $\frac{d}{d\mathit{}p}\mathit{}l\mathit{}\mathrm{(}p\mathrm{;}\mathrm{x}\mathrm{)}\mathrm{=}\mathrm{0}$, we get

	$-118p(3-2p)+575(1-p)(3-2p)-550p(1-p)$	$=$	$0$
	$-354p+236p^2+1725-2875p+1150p^2-550p+550p^2$	$=$	$0$
	$1936p^2-3779p+1725$	$=$	$0.$

Solving the quadratic gives

$$\widehat{p}=\frac{3779\pm \sqrt{{(-3779)}^2-4\times 1936\times 1725}}{2*1936}=1.2240\mathit{\text{ or }}0.7279.$$

Therefore $\widehat{p}\mathrm{=}\mathrm{0.7279}$ since .

Let $y$ denote the total number of households where the initial infective infects both the other individuals in the household. i.e. Outcome $\{1,2\}$ occurs.

Write down the likelihood $L(p;𝐱,y)$ of $𝐱$ and $y$ given $p$.

The data is again multinomial, but there are now four outcomes correspond as final size 3 is split into two. Hence,

	$L(p;𝐱,y)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3)!}{x_1!x_2!(x_3-y)!y!}}{({(1-p)}^2)}^{x_1}{(2p{(1-p)}^2)}^{x_2}{(2p^2(1-p))}^{x_3-y}{(p^2)}^y$
		$=$	${\displaystyle \frac{(x_1+x_2+x_3)!}{x_1!x_2!(x_3-y)!y!}}{2}^{x_2+x_3-y})(1-p){}^{2x_1+2x_2+x_3-y}p{}^{x_2+2x_3}$

Simplify $l(p;𝐱,y)=\log L(p;𝐱,y)$ and compute the MLE, $\widehat{p}$ for $(𝐱,y)$.
(This will help form the M-step of the EM algorithm.)

From the previous answer,

$l(p;𝐱,y)$

$=$

$K+(2x_1+2x_2+x_3-y)\log (1-p)+(x_2+2x_3)\log p.$

Therefore

$$\frac{d}{dp}l(p;𝐱,y)=-\frac{2x_1+2x_2+x_3-y}{1-p}+\frac{x_2+2x_3}{p}.$$

Setting $\frac{d}{d\mathit{}p}\mathit{}l\mathit{}\mathrm{(}p\mathrm{;}\mathrm{x}\mathrm{,}y\mathrm{)}\mathrm{=}\mathrm{0}$ gives

	$-(2x_1+2x_2+x_3-y)p+(1-p)(x_2+2x_3)$	$=$	$0$
	$x_2+2x_3$	$=$	$p\{2x_1+2x_2+x_3-y+x_2+2x_3\}.$

Hence

$$\widehat{p}=\frac{x_2+2x_3}{2x_1+3x_2+3x_3-y}.$$

Notice the similarities with the genetic linkage example with

$$l(p;𝐱,y)=K+B\log (1-p)+A\log p$$

and $\widehat{p}\mathrm{=}A\mathrm{/}\mathrm{(}A\mathrm{+}B\mathrm{)}$.

What is the distribution of $y$ given $𝐱$ and $p$?
Hint: Think about similarities between this problem and the genetics example.

There are $x_{\mathrm{3}}\mathrm{=}\mathrm{275}$ households where all three individuals are infected. For each of these households there are two possibilities $\mathrm{\{}\mathrm{1}\mathrm{,}\mathrm{1}\mathrm{,}\mathrm{1}\mathrm{\}}$ and $\mathrm{\{}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{\}}$. The conditional probability of $\mathrm{\{}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{\}}$ given that three people are infected in the household is

$$P(\{1,2\}|\{1,1,1\}\cup \{1,2\})=\frac{P(\{1,2\})}{P(\{1,1,1\})+P(\{1,2\})}=\frac{p^2}{2p^2(1-p)+p^2}=\frac{1}{3-2p}.$$

Hence, $y\mathrm{|}\mathrm{x}\mathrm{,}p\mathrm{\sim }\mathrm{Bin}\mathit{}\mathrm{(}{x}_{\mathrm{3}}\mathrm{(}\mathrm{=}\mathrm{275}\mathrm{)}\mathrm{,}\mathrm{1}\mathrm{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{}p\mathrm{)}\mathrm{)}$.

Write down $Q(p,p^{*})=E_{p^{*}}[l(p;𝐱,Y)|𝐱]$.
(This will help form the E-step of the EM algorithm.)

	$Q(p,p^{*})$	$=$	$E_{p^{*}}[l(p;𝐱,Y)|𝐱]$
		$=$	$E_{p^{*}}[K+(2x_1+2x_2+x_3-Y)\log (1-p)+(x_2+2x_3)\log p|𝐱]$
		$=$	$E_{p^{*}}[K|𝐱]+(2x_1+2x_2+x_3-E_{p^{*}}[Y|𝐱])\log (1-p)+(x_2+2x_3)\log p.$

Note that $E_{p^{\mathrm{*}}}\mathit{}\mathrm{[}K\mathrm{|}\mathrm{x}\mathrm{]}$ is not important since it is not a function of p. Therefore we only need to find $E_{p^{\mathrm{*}}}\mathit{}\mathrm{[}Y\mathrm{|}\mathrm{x}\mathrm{]}$ which using $y\mathrm{|}\mathrm{x}\mathrm{,}p\mathrm{\sim }\mathrm{Bin}\mathit{}\mathrm{(}{x}_{\mathrm{3}}\mathrm{(}\mathrm{=}\mathrm{275}\mathrm{)}\mathrm{,}\mathrm{1}\mathrm{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{}p\mathrm{)}\mathrm{)}$, gives $E_{p^{\mathrm{*}}}\mathit{}\mathrm{[}Y\mathrm{|}\mathrm{x}\mathrm{]}\mathrm{=}\mathrm{275}\mathrm{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{}p\mathrm{)}$.

Write an EM algorithm in R to find the MLE $\widehat{p}$.

See R code solutions.

Calculate the standard error of the MLE $\widehat{p}$.

Firstly note that

	$Q(p,p^{*})$	$=$	$E_{p^{*}}[l(p;𝐱,Y)|𝐱]$
		$=$	$E_{p^{*}}[K+(2x_1+2x_2+x_3-Y)\log (1-p)+(x_2+2x_3)\log p|𝐱]$
		$=$	$E_{p^{*}}[K|𝐱]+(2x_1+2x_2+x_3-E_{p^{*}}[Y|𝐱])\log (1-p)+(x_2+2x_3)\log p$
		$=$	$E_{p^{*}}[K|𝐱]+(393-E_{p^{*}}[Y|𝐱])\log (1-p)+575\log p.$

Therefore $E_{p^{\mathrm{*}}}\mathit{}\mathrm{[}Y\mathrm{|}\mathrm{x}\mathrm{]}\mathrm{=}\mathrm{275}\mathrm{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{}{p}^{\mathrm{*}}\mathrm{)}$, we have that

$-{{\displaystyle \frac{{\partial }^{2}Q(p,p^{*})}{\partial p^2}}|}_{p=p^{\ast }=\widehat{p}}$

$=$

${\displaystyle \frac{393-275/(3-2\widehat{p})}{{(1-\widehat{p})}^2}}+{\displaystyle \frac{575}{{\widehat{p}}^2}}=3987.977.$

Secondly, since $y\mathrm{|}\mathrm{x}\mathrm{,}{p}^{\mathrm{\ast }}\mathrm{\sim }\mathrm{Bin}\mathit{}\mathrm{(}\mathrm{275}\mathrm{,}\mathrm{1}\mathrm{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{}{p}^{\mathrm{\ast }}\mathrm{)}\mathrm{)}$, we have that

	$\mathrm{Var}\left({\displaystyle \frac{\partial \log f(𝐱,Y|p)}{dp}}|𝐱,\widehat{p}\right)$	$=$	$\mathrm{Var}\left(-{\displaystyle \frac{2x_1+2x_2+x_3-Y}{1-p}}+{\displaystyle \frac{x_2+2x_3}{p}}|𝐱,\widehat{p}\right)$
		$=$	${\displaystyle \frac{1}{{(1-\widehat{p})}^2}}\mathrm{Var}(Y|𝐱,\widehat{p})$
		$=$	${\displaystyle \frac{1}{{(1-\widehat{p})}^2}}\times 275\times {\displaystyle \frac{1}{3-2\widehat{p}}}\times {\displaystyle \frac{2(1-\widehat{p})}{3-2\widehat{p}}}=847.6705$

Therefore the standard error of $\widehat{p}$ is $\sqrt{\mathrm{1}\mathrm{/}\mathrm{(}\mathrm{3987.977}\mathrm{-}\mathrm{847.6705}\mathrm{)}}\mathrm{=}\mathrm{0.01784}$.

1. (a)

   The likelihood of $p$ and $\lambda$ is

   $$L(p,\lambda |𝐱)=\prod _{i=1}^n\left\{p{\delta }_{x_{i}0}+(1-p)\frac{{\lambda }^{x_i}}{x_i!}{e}^{-\lambda }\right\},$$

   where ${\delta }_{ij}=1$ if $i=j$ and ${\delta }_{ij}=0$ if $i\ne j$

2. (b)

   The log-likelihood of $p$ and $\lambda$ given $𝐮$ is

   	$l(p,\lambda |𝐱,𝐮)$	$=$	${\displaystyle \sum _{i;u_i=1}}\{\log p+\log f(x_i)\}+{\displaystyle \sum _{i;u_i=0}}\{\log (1-p)+\log g(x_i)\}$
   		$=$	${\displaystyle \sum _{i;u_i=1}}\{\log p+\log f(x_i)\}+{\displaystyle \sum _{i;u_i=0}}\{\log (1-p)+x_i\log \lambda -\log x_i!-\lambda \}$

   Now $m={\sum }_{i;u_i=1}1$ and $n-m={\sum }_{i;u_i=0}1$. Note that if $u_i=1$ then $x_i=0$, so ${\sum }_{i;u_i=0}{x}_i={\sum }_{i=1}^n{x}_i=n\mu$. (Since $U_i=1$ only if $x_i=0$.)

   Therefore

   	$l(p,\lambda |𝐱,𝐮)$	$=$	$m\log p+(n-m)\log (1-p)+n\mu \log \lambda +{\displaystyle \sum _{i;u_i=0}}\log x_i!-(n-m)\lambda$		(1.1)
   		$=$	$m\log p+(n-m)\log (1-p)+n\mu \lambda -(n-m)\lambda +A$

   as required.

3. (c)

   To find the MLEs simply differentiate (1.1) with respect to $p$ and $\lambda$, set the derivative equal to 0 and solve.

   Firstly,

   ${\displaystyle \frac{\partial l}{\partial p}}$

   $=$

   ${\displaystyle \frac{m}{p}}-{\displaystyle \frac{n-m}{1-p}}$

   Therefore $\widehat{p}$ satisfies

   	${\displaystyle \frac{m}{p}}-{\displaystyle \frac{n-m}{1-p}}$	$=$	$0$
   	$m(1-p)-(n-m)p$	$=$	$0$
   	$m-np$	$=$	$0.$

   Therefore $\widehat{p}=\frac{m}{n}$.

   Secondly,

   ${\displaystyle \frac{\partial l}{\partial \lambda }}$

   $=$

   ${\displaystyle \frac{n\mu }{\lambda }}-(n-m)$

   Therefore $\widehat{\lambda }$ satisfies

   	${\displaystyle \frac{n\mu }{\lambda }}-(n-m)$	$=$	$0$
   	$n\mu -(n-m)\lambda$	$=$	$0.$

   Therefore $\widehat{\lambda }=\frac{n}{n-m}\mu$.

4. (d)

   This involves find $E[U_i|p,x_i]$ for $i=1,2,\mathrm{\dots },n$ which is

   	$E[U_i|p,x_i]$	$=$	${\displaystyle \frac{P(U_i=1,X_i=0|p,\lambda )}{P(X_i=0|p,\lambda )}}$
   		$=$	${\displaystyle \frac{pf(x_i)}{pf(x_i)+(1-p)g(x_i)}}.$

   For $x_i>0$, $f(x_i)=0$ and so $E[U_i|p,x_i]=0$.

   For $x_i=0$, $f(x_i)=1$ and $g(x_i)=\exp (-\lambda )$. Therefore

   $E[U_i|p,x_i]={\displaystyle \frac{p\times 1}{p\times 1+(1-p)\times \exp (-\lambda )}}={\displaystyle \frac{p}{p+(1-p)\exp (-\lambda )}}$

   as required.

5. (e)

   The EM algorithm alternates between:-

   • •

     E-step: $E[l(p|𝐱,𝐔)]$ given the current estimate of p.

     Note that

     $$l(p|𝐱,𝐔)=\left(\sum _{i=1}^n{U}_i\right)\log p+\left(n-\sum _{i=1}^n{U}_i\right)\log (1-p)+n\mu \log \lambda +\left(n-\sum _{i=1}^n{U}_i\right)\lambda +A.$$

     Therefore it suffices to compute $E[U_i|p,x_i]$ $(i=1,2,\mathrm{\dots },n)$ which is derived in $(d)$.

   • •

     M-step: Maximizing $l(p,\lambda |𝐱,𝐮)$ with respect to $p$ and $\lambda$ where $𝐮$ is replaced by the corresponding expected values from the E-step (part $(d)$) and the MLEs are given in the part $(c)$.

1. (a)

   The likelihood is

   	$L(\beta ;𝐲,𝐳)$	$=$	$f(𝐲,𝐳|\beta )$
   		$=$	${\displaystyle \prod _{i=1}^n}f(y_i,z_i|\beta )$
   		$=$	${\displaystyle \prod _{i=1}^n}f(y_i|z_i,\beta )f(z_i|\beta )$
   		$=$	${\displaystyle \prod _{i=1}^n}\left\{1_{\{z_i>0\}}^{y_i}{1}_{\{z_i\le 0\}}^{1-y_i}\right\}\times {\displaystyle \frac{1}{\sqrt{2\pi }}}\exp (-1/2{(z_i-x_i\beta )}^2)$
   		$=$	${\displaystyle \prod _{i=1}^n}\left\{1_{\{z_i>0\}}^{y_i}{1}_{\{z_i\le 0\}}^{1-y_i}\right\}{(2\pi )}^{-n/2}{\displaystyle \prod _{i=1}^n}\exp (-1/2{(z_i-x_i\beta )}^2)$
   		$=$	${\displaystyle \prod _{i=1}^n}{1}_{\{z_i>0\}}^{y_i}{1}_{\{z_i\le 0\}}^{1-y_i}{(2\pi )}^{-n/2}\exp \left(-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{(z_i-x_i\beta )}^2\right).$

2. (b)

   Note that the log-likelihood satisfies

   $$l(\beta ;𝐲,𝐳)=K-\frac{1}{2}\sum _{i=1}^n{(z_i-x_i\beta )}^2,$$

   where $K$ is a constant not depending upon $\beta$.
   Hence

   	${\displaystyle \frac{dl(\beta ;𝐲,𝐳)}{d\beta }}$	$=$	$-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}2(z_i-x_i\beta )(-x_i)$
   		$=$	$-{\displaystyle \sum _{i=1}^n}{z}_i{x}_i+\beta {\displaystyle \sum _{i=1}^n}{x}_i^2.$

   Setting $\frac{dl(\beta ;𝐲,𝐳)}{d\beta }=0$ yields

   $$\tilde{\beta }=\frac{{\sum }_{i=1}^n{z}_i{x}_i}{{\sum }_{i=1}^n{x}_i^2}.$$

3. (c)

   Note that $y_i=1$ implies that $z_i>0$. Therefore $z_i|y_i=1,x_i,\beta$ follows $N(x_i\beta ,1)$ conditioned to be positive.
   Therefore

   	$E[z_i|y_i=1,x_i,\beta ]$	$=$	$E[z_i|z_i>0]$
   		$=$	$x_i\beta +{\displaystyle \frac{\varphi (-x_i\beta )}{1-\mathrm{\Phi }(-x_i\beta )}}$
   		$=$	$x_i\beta +{\displaystyle \frac{\varphi (x_i\beta )}{\mathrm{\Phi }(x_i\beta )}}.$

4. (d)
   • •

     Choose an initial value for $\beta$, say $\beta =0$.

   • •

     E-step: For each $i$, compute $E[z_i|y_i,x_i,\beta ]$. For $y_i=1$, $E[z_i|y_i,x_i,\beta ]$ is given by (c) and $E[z_i|y_i=0,x_i,\beta ]$ can be computed similarly.

   • •

     M-Step: Compute $\tilde{\beta }=\frac{{\sum }_{i=1}^n{z}_i{x}_i}{{\sum }_{i=1}^n{x}_i^2}$.

   • •

     Stop when two consecutive estimates of $\tilde{\beta }$ agree to a predefined precision, say ${10}^{-4}$.