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

% *************************************************************************
%
% FILE: raobisect.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 bisection 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: midpoint of intervals (mid(i) = (low(i) + upp(i))/2 (length n+1)
%        The last entry should be treated as the zero of the function.
%
% *************************************************************************

% Initial vectors

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

% Run algorithm n times

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

% Determine whether midpoint should me new right or left endpoint and
% 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

mid = [mid (u+l)/2];