function [lambda, niter] = raoem(y, lambda0, TOL, Nmax)

% *************************************************************************
%
% FILE: raoemt.m
%
% Date Created: Feb 4, 2004
% Author: Mark E. Irwin
%
% Contents:  Finds the MLE 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
%   TOL: convergence tolerance
%   Nmax: number of functional iterations
%
% Output:
%
%   lambda: MLE
%   niter: number of EM iterations
%
% *************************************************************************

% Initial vectors

n = 0;
diff = 2 * TOL;
lambdan = lambda0;

% Run algorithm n times

while (n < Nmax && diff > TOL)
  n = n + 1;
  lambdaold = lambdan;
  x4 = y(4) * lambdan / (2 + lambdan);
  lambdan = (y(1) + x4) / (sum(y(1:3)) + x4);
  
  diff = abs(lambdan - lambdaold);
  
end

if (diff > TOL)
  warning('Method did not converage after Nmax steps');
end

niter = n;
lambda = lambdan;