\documentstyle[11pt,fullpage]{article}

\begin{document}

\parindent=0pt
\parskip=0.6em

\begin{center}\bf
Example of E-M: Genetic Linkage Problem
\end{center}

Rao (1973) considers 197 animals randomly divided into four categories
as follows:
\[
                Y = (y_1, y_2, y_3, y_4)^T = (125, 18, 20, 34)^T
\]
with cell probabilities
\[
                ( 1/2 + \theta/4, (1-\theta)/4, (1-\theta)/4, \theta/4 )^T
\]
The observed-data likelihood is 
\[
        g(Y|\theta) 
        \propto \cdot (2 + \theta)^{y_1}(1-\theta)^{y_2+y_3} \theta^{y_4}
\]
If we ``augment'' the data by splitting the first category into two, 
\begin{eqnarray*}
         X & = & (x_1, x_2, x_3, x_4, x_5)^T \\
\mbox{with probabilities} & & (1/2, \theta/4, (1-\theta)/4,
                              (1-\theta)/4, \theta/4 )^T 
\end{eqnarray*}
the complete-data likelihood is now simpler:
\[
        f(X|\theta) \propto \theta^{x_2+x_5} (1-\theta)^{x_3+x_4}
\]
The ``missing data'' $Z$ is like a coin toss for which category each
observation of $y_1$ gets split into:
\[
        (Z|Y,\theta) \sim Bin(y_1, 2/(2+\theta))
\]
Now we compute the $Q$-function,
\begin{eqnarray*}
        Q(\theta,\theta^k) & = & E[(x_2+x_5)\log(\theta) + (x_3 +
        x_4)\log(1-\theta) ~|~ Y,\theta^k] \\
        & = & (y_1\theta^k/(2+\theta^k) + x_5)\log(\theta) +
        (x_3+x_4)\log(1-\theta) 
\end{eqnarray*}
This is of the form $a\log(\theta) + b\log(1-\theta)$, a binomial
whose MLE is $\hat\theta=a/(a+b)$; translating back to the particulars
of $Q(\cdot,\cdot)$ we get
\[
        \theta^{k+1} = \frac{y_1\theta^k/(2+\theta^k) +
        x_5}{y_1\theta^k/(2+\theta^k) + x_5 + x_3+x_4}
\]

Now we attempt to find the MLE of $\theta$, two ways:

\begin{itemize}

\item By Newton-Raphson
\begin{verbatim}
> Y _ c(125, 18, 20, 34)
> g _ function(th) {
+ (2 + th)^(125) * (1-th)^(18+20) * th^(34)
+ }
> ng _ function(th) {-log(g(th))}
> summary(ms( ~ ng(th),start=c(th=0.5)))
Final value: -67.3841
Solution:
        Par. 
th 0.6268215

Convergence: RELATIVE FUNCTION CONVERGENCE.

Computations done:
 Iterations Function Gradient 
          6        8        9
> source("infomat.s")
> info.mtx.fun(ng,0.6268215)
          [,1] 
[1,] -377.5169
> sqrt(1/377.5169)
[1] 0.05146735
\end{verbatim}

\item By E-M
\begin{verbatim}
> emfun <- 
function(th0, maxit = 1000, releps = 1e-09)
{
        i <- 0
        th1 <- th0
        th0 <- th1 + 1
        while((i != maxit) && (abs(th1 - th0) > releps * abs(th0))) {
                i <- i + 1
                th0 <- th1
                th1 <- ((125 * th1)/(2 + th1) + 34)/((125 * th1)/(2 + th1) +
                        34 + 18 + 20)
                print(c(th0,th1,i))
        }
        return(i, th1)
}
> emfun (0.5)
[1] 0.5000000 0.6082474 1.0000000
[1] 0.6082474 0.6243211 2.0000000
[1] 0.6243211 0.6264889 3.0000000
[1] 0.6264889 0.6267773 4.0000000
[1] 0.6267773 0.6268156 5.0000000
[1] 0.6268156 0.6268207 6.0000000
[1] 0.6268207 0.6268214 7.0000000
[1] 0.6268214 0.6268215 8.0000000
[1] 0.6268215 0.6268215 9.0000000
[1]  0.6268215  0.6268215 10.0000000
[1]  0.6268215  0.6268215 11.0000000
$i:
[1] 11

$th1:
[1] 0.6268215
\end{verbatim}

\end{itemize}

\end{document}

