mcemlinkage <- function(y, lambda0, niter = 1000, TOL = 1e-6, Nmax=20) {

# *************************************************************************
#
# FILE:mcemlinkage.r
#
# Date Created: April 22, 2004
# Author: Mark Irwin
#
# Contents: Finds the MLE  by MCEM for the AB|ab x AB|ab linkage analysis
#           from Rao, 1973 by EM.  The functions allows for different data 
#           sets to be used, not just the standard one used in most papers
#           ([125 18 20 34]).
#
# Revision History
# Date          Name            Changes/Reasons
#
# *************************************************************************
#
# Input:
#
#   y: the vector of counts for ab, Ab, aB, and AB phenotypes
#   lambda0: initial estimate of lambda
#   niter: number of MC imputations for each EM step
#   TOL: convergence tolerance
#   Nmax: number of EM steps
#
# Output:
#
#   lambda: MCMLE estimate sequence
#
# *************************************************************************

# Initial vectors

n <- 0
diff <- 2 * TOL
lambdan <- lambda0
lambda <- lambda0

# Run algorithm n times

while (n < Nmax && diff > TOL) {
  n <- n + 1
  lambdaold <- lambdan
  prob <- lambdan / (2 + lambdan)
  x4samp <- rbinom(niter, y[4], prob)
  x4 <- mean(x4samp)
  lambdan <- (y[1] + x4) / (sum(y[1:3]) + x4)
  lambda <- c(lambda, lambdan)
  diff <-abs(lambdan - lambdaold)
}

lambda
}