fork download
copy 
  1. import math
  2.  
  3. def binconf(p, n, c=0.95):
  4.   '''
  5.   Calculate binomial confidence interval based on the number of positive and
  6.   negative events observed.  Uses Wilson score and approximations to inverse
  7.   of normal cumulative density function.
  8.  
  9.   Parameters
  10.   ----------
  11.   p: int
  12.       number of positive events observed
  13.   n: int
  14.       number of negative events observed
  15.   c : optional, [0,1]
  16.       confidence percentage. e.g. 0.95 means 95% confident the probability of
  17.       success lies between the 2 returned values
  18.  
  19.   Returns
  20.   -------
  21.   theta_low  : float
  22.       lower bound on confidence interval
  23.   theta_high : float
  24.       upper bound on confidence interval
  25.   '''
  26.   p, n = float(p), float(n)
  27.   N    = p + n
  28.  
  29.   if N == 0.0: return (0.0, 1.0)
  30.  
  31.   p = p / N
  32.   z = normcdfi(1 - 0.5 * (1-c))
  33.  
  34.   a1 = 1.0 / (1.0 + z * z / N)
  35.   a2 = p + z * z / (2 * N)
  36.   a3 = z * math.sqrt(p * (1-p) / N + z * z / (4 * N * N))
  37.  
  38.   return (a1 * (a2 - a3), a1 * (a2 + a3))
  39.  
  40.  
  41. def erfi(x):
  42.   """Approximation to inverse error function"""
  43.   a  = 0.147  # MAGIC!!!
  44.   a1 = math.log(1 - x * x)
  45.   a2 = (
  46.     2.0 / (math.pi * a)
  47.     + a1 / 2.0
  48.   )
  49.  
  50.   return (
  51.     sign(x) *
  52.     math.sqrt( math.sqrt(a2 * a2 - a1 / a) - a2 )
  53.   )
  54.  
  55.  
  56. def sign(x):
  57.   if x  < 0: return -1
  58.   if x == 0: return  0
  59.   if x  > 0: return  1
  60.  
  61.  
  62. def normcdfi(p, mu=0.0, sigma2=1.0):
  63.   """Inverse CDF of normal distribution"""
  64.   if mu == 0.0 and sigma2 == 1.0:
  65.     return math.sqrt(2) * erfi(2 * p - 1)
  66.   else:
  67.     return mu + math.sqrt(sigma2) * normcdfi(p)
  68.  
  69. print binconf(14224,8935)
  70. print "\n"
  71. print binconf(6553,1181)
Success #stdin #stdout 0.01s 8968KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
(0.6079039577456685, 0.6204359386460506)


(0.839112731337121, 0.8551380514201897)
https://ideone.com/bQDO1N
language:
Python (cpython 2.7.16)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!