Find the error in the mean of the ratio of two means from normal distributions

mean_error_ratio.py

from scipy.stats import t
import math

# Method from "Intuitive Biostatistics" pp. 285-286
#
#@book{harvey1995intuitive,
#   title={Intuitive Biostatistics},
#   author={Harvey Motulsky},
#   isbn={9780195086072},
#   lccn={95008166},
#   url={https://books.google.co.uk/books?id=PzaimmdNRogC},
#   year={1995},
#   publisher={Oxford University Press}
# }

def standard_error_quotient(mean_a, mean_b, stdev_a, stdev_b, samples, confint=0.95):
  #Student T distribution, n=samples-1, p<1-confint, 2-tail
  # t_intervals = stats.t.ppf(confint, samples-1)
  # alpha=confint, df = samples-1
  lower, upper = t.interval(confint, samples-1)
  sem_a = stdev_a/math.sqrt(samples)  
  sem_b = stdev_b/math.sqrt(samples)
  quotient = mean_a/mean_b
  std_err_q = quotient * math.sqrt(((sem_a**2)/(mean_a**2))+((sem_b**2)/(mean_b**2)))
  # print "std err q = "+str(std_err_q)
  return (quotient - (upper*std_err_q), quotient + (upper*std_err_q))

# Method from Fieller's Theorem.
# Method assumes that a and b are jointly normally distributed, and b is not too near 0
# Fieller, EC. (1954). "Some problems in interval estimation."
#
# @article{1954,
#      jstor_articletype = {research-article},
#      title = {Some Problems in Interval Estimation},
#      author = {Fieller, E. C.},
#      journal = {Journal of the Royal Statistical Society. Series B (Methodological)},
#      jstor_issuetitle = {},
#      volume = {16},
#      number = {2},
#      jstor_formatteddate = {1954},
#      pages = {pp. 175-185},
#      url = {http://www.jstor.org/stable/2984043},
#      ISSN = {00359246},
#      year = {1954},
#      publisher = {Wiley for the Royal Statistical Society},
#      copyright = {Copyright 1954 Royal Statistical Society},
#     }

def standard_error_quotient_fieller(mean_a, mean_b, stdev_a, stdev_b, samples, confint=0.95):
  #Student T distribution, n=samples-1, p<1-confint, 2-tail
  # t_intervals = stats.t.ppf(confint, samples-1)
  # alpha=confint, df = samples-1
  lower, upper = t.interval(confint, samples-1)
  quotient = mean_a/mean_b
  var_mean_quot = (quotient**2)*(((stdev_a**2)/(mean_a**2))+((stdev_b**2)/(mean_b**2)))
  # print "var mean quot = "+str(var_mean_quot)
  return (quotient - (upper*var_mean_quot), quotient + (upper*var_mean_quot))

ma = 88.1034
mb = 175.403
sa = 3.95574
sb = 9.34458
s = 100

mean = ma/mb
print "Biostatistics values:"
low, high = standard_error_quotient(ma, mb, sa, sb, s)
print "\tMean: "+str(mean)
print "\tLow: "+str(low)
print "\tHigh: "+str(high)
print "Fieller values:"
fieller_low, fieller_high = standard_error_quotient_fieller(ma, mb, sa, sb, s)
print "\tMean: "+str(mean)
print "\tLow: "+str(fieller_low)
print "\tHigh: "+str(fieller_high)
print "Method deltas:"
print "\tLow: " +str(abs(low-fieller_low))+"  "+str(100*abs(low-fieller_low)/low)+"%  "+str(100*abs(low-fieller_low)/fieller_low)+"%"
print "\tHigh: " +str(abs(high-fieller_high))+"  "+str(100*abs(high-fieller_high)/high)+"%  "+str(100*abs(high-fieller_high)/fieller_high)+"%"