Animals that are crosses of AB/ab x AB/ab are classified in 4 categories: \[\begin{array}{*{20}{c}} Y&{Genotype}&{\Pr (category|\theta )}&{\Pr (category|\pi )} \\ {125}&{AB}&{(3 - 2\theta + {\theta ^2})/4}&{\frac{1}{2} + \frac{\pi }{4}} \\ {18}&{Ab}&{(2\theta - {\theta ^2})/4}&{\frac{1}{4} - \frac{\pi }{4}} \\ {20}&{aB}&{(2\theta - {\theta ^2})/4}&{\frac{1}{4} - \frac{\pi }{4}} \\ {34}&{ab}&{(1 - 2\theta + {\theta ^2})/4}&{\frac{\pi }{4}} \end{array}\] The data are:

y <- y.rao <- c(125,18,20,34)  #Original data.

The likelihood function is: \[\begin{gathered} ODL \propto [Y|\pi ] \propto {\left( {\frac{1}{2} + \frac{\pi }{4}} \right)^{{Y_1}}}{\left( {\frac{1}{4} - \frac{\pi }{4}} \right)^{{Y_2} + {Y_3}}}{\left( {\frac{\pi }{4}} \right)^{{Y_4}}} \hfill \\ \quad \quad \quad \quad \quad \quad \propto \sum\limits_{{X_1} = 0}^{{Y_1}} {\frac{{n!}}{{{X_1}!{X_2}!{Y_2}!{Y_3}!{Y_4}!}}{{\left( {\frac{1}{2}} \right)}^{{X_1}}}{{\left( {\frac{1}{4} - \frac{\pi }{4}} \right)}^{{Y_2} + {Y_3}}}{{\left( {\frac{\pi }{4}} \right)}^{{X_2} + {Y_4}}}} \hfill \\ \quad \quad \quad \quad = \sum\limits_{{X_1} = 0}^{{Y_1}} {\frac{{n!}}{{{X_1}!{X_2}!{Y_2}!{Y_3}!{Y_4}!}}CDL(\pi ;{X_1},{X_2},{Y_2},{Y_3},{Y_4})} \hfill \\ \end{gathered} \]

The code for the likelihood function and its log:

odl <- odl.rao <- function(y,p)
{
  (1/2+p/4)^y[1] * ((1-p)/4)^(y[2]+y[3]) * (p/4)^y[4] 
}

log.odl <- function(y,p) {
  log(odl(y,p))
}

Plots of the ODL and its log:

p.delta <- 0.01
p.range <- seq(.01,.99,by=p.delta)
plot (p.range, odl(y,p.range), type="l",
      main="likelihood function")

plot (p.range, log.odl(y,p.range), type="l",
      main="log(likelihood)")

Newton-Raphson method uses the first and second derivatives:

logodlDot = function(y,p)
  y[1]/(2+p) -(y[2]+y[3])/(1-p) + y[4]/p
logodlDotDot = function(y,p)
  -(y[1]/(2+p)^2 + (y[2]+y[3])/(1-p)^2 + y[4]/p^2)

plot (p.range, log.odl(y,p.range), type="l",
      main="log(likelihood);  Newton-Raphson")
p0 = 0.01
for(iter in 1:9) {
  points(p0, log(odl.rao(y=y.rao, p0)), col="red", cex=4, pch=as.character(iter))
  p0 = p0 - logodlDot(y.rao, p0)/logodlDotDot(y.rao, p0)
}
p0 = 0.99
for(iter in 1:9) {
  points(p0, log(odl.rao(y=y.rao, p0)), col="blue", cex=4, pch=as.character(iter))
  p0 = p0 - logodlDot(y.rao, p0)/logodlDotDot(y.rao, p0)
}

For the EM algorithm, we need the complete data likelihood. \[\begin{gathered} CDL = [X|\pi ] = \frac{{n!}}{{{X_1}!{X_2}!{Y_2}!{Y_3}!{Y_4}!}}{\left( {\frac{1}{2}} \right)^{{Y_1} - {X_2}}}{\left( {\frac{\pi }{4}} \right)^{{X_2}}}{\left( {\frac{1}{4} - \frac{\pi }{4}} \right)^{{Y_2} + {Y_3}}}{\left( {\frac{\pi }{4}} \right)^{{Y_4}}} \hfill \\ \quad \quad \propto {\pi ^{{X_2} + {Y_4}}}{(1 - \pi )^{{Y_2} + {Y_3}}} \\ \quad \quad \quad \quad \\ \end{gathered} \]

and the the expected sufficient statistic,

\[E({X_2}|Y,\pi ) = {Y_1}\frac{\pi }{{2 + \pi }}\]

nsteps <-4
plot (p.range, log.odl(y,p.range), type="l",
      main="log(likelihood);  EM algorithm")
p0 = 0.01
for(iter in 1:4) {
  points(p0, log(odl.rao(y=y.rao, p0)), col="red", cex=4, pch=as.character(iter))
  ExpX2 =  y[1]*p0/(2+p0) + y[4] 
  p0 <- ExpX2/(ExpX2 + y[2]+y[3])
}
p0 = 0.99
for(iter in 1:4) {
  points(p0, log(odl.rao(y=y.rao, p0)), col="blue", cex=4, pch=as.character(iter))
  ExpX2 =  y[1]*p0/(2+p0) + y[4] 
  p0 <- ExpX2/(ExpX2 + y[2]+y[3])
}