# missing counts
mm <- c(NA,NA,NA,NA,3,17)
# data
nn <- c(4,9,3,17)

# Initial guess
#theta0 <- c(.2,.3,.5) # guess of p,q,r
theta0 <- rep(1/3,3)

err <- function(theta0,theta1)
	{
	tmp <- .5*abs(theta0-theta1)/(theta0+theta1)
	return( mean(tmp) )
	}

# report the loglikelihood (up to a constant)

ll <- function(theta,nn)
	{
	 pp <- rep(NA,4)
	 pp[1] <- theta[1]^2 + 2*theta[1]*theta[3]
	 pp[2] <- theta[2]^2 + 2*theta[2]*theta[3]
	 pp[3] <- 2*theta[1]*theta[2]
	 pp[4] <- theta[3]^2
	 loglik <- sum( nn*log(pp) )
	 return(loglik)
	}
# iterate

notok <- T
theta <- theta0
tol <- 1e-5
maxiter <- 100
niter <- 1
while( notok )
 {
  print(theta)
  print( ll(theta,nn) )
  # E step
  mm[1]  <-  nn[1] * theta[1]^2/( theta[1]^2+2*theta[1]*theta[3] )
  mm[2] <- nn[1] - mm[1]
  mm[3] <- nn[2] * theta[2]^2/( theta[2]^2 + 2*theta[2]*theta[3] )
  mm[4] <- nn[2] - mm[3]

  print( round(mm,3) )
  # M step
  theta[1] <-  (2*mm[1]+mm[2]+mm[5])
  theta[2] <- (2*mm[3]+mm[4]+mm[5])
  theta[3] <- (2*mm[6]+mm[2]+mm[4])
  theta <- theta/sum(theta)
 
  # convergence?
  if( err(theta,theta0) <= tol | niter == maxiter )
     notok <- F

  niter <- niter+1
  theta0 <- theta  # update old theta
 }