function [low, upp , mid] = raofalsi(x,lower,upper,n)

% *************************************************************************
%
% FILE: raofalsi.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 the regula falsi method.  The first derivative of
%            log likelihood is expected to be in a function rao (in rao.m).
%            The functions allows for different data sets to be used, not
%            just the standard one used in most papers ([125 18 20 34]).
%            Note that this code does not to error checking to verify that 
%            the initial interval contains a root.
%
% Revision History
% Date          Name            Changes/Reasons
%
% *************************************************************************
%
% Input:
%
%   x: the vector of counts for AB, Ab, aB, and ab phenotypes
%   lower: lower bound for root
%   upper: upper bound for root
%   n: number in iterations of bisection algorithm
%
% Output:
%
%   low: sequence of lower bounds of root (length n+1)
%   upp: sequence of upper bounds of root (length n+1)
%   mid: x intercept of line joining (low,g(low)) and (upp,g(upp)) (length n+1)
%        The last entry should be treated as the zero of the function.
%   
% *************************************************************************

% Set up vectors

low = lower;
upp = upper;
mid = [];
l = lower;
u = upper;

% Iterate n times

for i = 1:n
  gl = rao(x,l);
  gu = rao(x,u);
  c = l - (u-l)*gl/(gu-gl);
  mid = [mid c];
  gc = rao(x,c);

% Update interval

  if((gl > 0) && (gu <0) && (gc > 0))
    low = [low  c];
    upp = [upp  u];
    l = c;
  elseif((gl > 0) && (gu <0) && (gc < 0))
    low = [low  l];
    upp = [upp  c];
    u = c;
  elseif((gl < 0) && (gu >0) && (gc > 0))
    low = [low  l];
    upp = [upp  c];
    u = c;
  else
    low = [low  c];
    upp = [upp  u];
    l = c;
  end
end

gl = rao(x,l);
gu = rao(x,u);
c = l - (u-l)*gl/(gu-gl);
mid = [mid c];
