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

# *************************************************************************
#
# FILE: rwmcemlinkage.r
#
# Date Created: April 22, 2004
# Author: Mark Irwin
#
# Contents: Finds the MLE  by reweighted 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
#           ([34 20 18 125]).
#
# 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

# Generate sample
p0 <- lambdan / (2 + lambdan)
x4samp <- rbinom(niter, y[4], p0)
prob0 <- dbinom(x4samp, y[4], p0)


# Run algorithm n times

while (n < Nmax && diff > TOL) {
  n <- n + 1
  lambdaold <- lambdan
  
  # update weights
  pn <- lambdan / (2 + lambdan)
  weight <- dbinom(x4samp, y[4], pn) / prob0
  weight <- weight / sum(weight)
  x4 <- sum(x4samp * weight)
  lambdan <- (y[1] + x4) / (sum(y[1:3]) + x4)
  lambda <- c(lambda, lambdan)
  diff <-abs(lambdan - lambdaold)
}

lambda
}