raofalsi <- function(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 (as a list):
#
#   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)
#   
# *************************************************************************

# Setup vectors

low <- lower
upp <- upper
mid <-rep(0,n+1)
l <- lower
u <- upper

# Iterate n times

for(i in 1:n) {

  gl <- rao(x,l)
  gu <- rao(x,u)
  xint <- l - (u-l)*gl/(gu-gl)
  mid[i] <- xint
  gxint <- rao(x,xint)

# Update interval

  if((gl > 0) && (gu <0) && (gxint > 0)) {
    low <- c(low, xint)
    upp <- c(upp, u)
    l <- xint
  }
  else if((gl > 0) && (gu <0) && (gxint < 0)) {
    low <- c(low, l)
    upp <- c(upp, xint)
    u <- xint
  }
  else if((gl < 0) && (gu >0) && (gxint > 0)) {
    low <- c(low, l)
    upp <- c(upp, xint)
    u <- xint
  }
  else {
    low <- c(low, xint)
    upp <- c(upp, u)
    l <- xint
  }
}

# Get final estimate

gl <- rao(x,l)
gu <- rao(x,u)
xint <- l - (u-l)*gl/(gu-gl)
mid[n+1] <- xint

list(low=low, upp=upp, mid=mid)

}