MACS2/Prob.pyx

   1 # cython: language_level=3
   2 # cython: profile=True
   3 # Time-stamp: <2019-10-30 17:28:00 taoliu>
   4
   5 """Module Description: statistics functions to calculate p-values
   6
   7 This code is free software; you can redistribute it and/or modify it
   8 under the terms of the BSD License (see the file LICENSE included with
   9 the distribution).
  10 """
  11
  12 # ------------------------------------
  13 # python modules
  14 # ------------------------------------
  15 from libc.math cimport exp,log,log10, M_LN10 #,fabs,log1p
  16 from math import fabs
  17 from math import log1p #as py_log1p
  18 from math import sqrt
  19
  20 import numpy as np
  21 cimport numpy as np
  22
  23 from cpython cimport bool
  24 # ------------------------------------
  25 # constants
  26 # ------------------------------------
  27 cdef int LSTEP = 200
  28 cdef double EXPTHRES = exp(LSTEP)
  29 cdef double EXPSTEP  = exp(-1*LSTEP)
  30 cdef double bigx = 20
  31
  32 # ------------------------------------
  33 # Normal distribution functions
  34 # ------------------------------------
  35 # x is the value, u is the mean, v is the variance
  36 cpdef pnorm(int x, int u, int v):
  37     """The probability of X=x when X=Norm(u,v)
  38     """
  39     return 1.0/sqrt(2.0 * 3.141592653589793 * <float>v) * exp(-<float>(x-u)**2 / (2.0 * <float>v))
  40
  41 # ------------------------------------
  42 # Misc functions
  43 # ------------------------------------
  44
  45 cpdef factorial ( unsigned int n ):
  46     """Calculate N!.
  47
  48     """
  49     cdef double fact = 1
  50     cdef unsigned long i
  51     if n < 0:
  52         return 0
  53     for i in xrange( 2,n+1 ):
  54         fact = fact * i
  55     return fact
  56
  57 cdef double poz(double z):
  58     """ probability of normal z value
  59     """
  60     cdef:
  61         double y, x, w
  62         double Z_MAX = 6.0 # Maximum meaningful z value
  63     if z == 0.0:
  64         x = 0.0;
  65     else:
  66         y = 0.5 * fabs(z)
  67         if y >= (Z_MAX * 0.5):
  68             x = 1.0
  69         elif y < 1.0:
  70             w = y * y
  71             x = ((((((((0.000124818987 * w
  72                         - 0.001075204047) * w + 0.005198775019) * w
  73                         - 0.019198292004) * w + 0.059054035642) * w
  74                         - 0.151968751364) * w + 0.319152932694) * w
  75                         - 0.531923007300) * w + 0.797884560593) * y * 2.0
  76         else:
  77             y -= 2.0
  78             x = (((((((((((((-0.000045255659 * y
  79                              + 0.000152529290) * y - 0.000019538132) * y
  80                              - 0.000676904986) * y + 0.001390604284) * y
  81                              - 0.000794620820) * y - 0.002034254874) * y
  82                              + 0.006549791214) * y - 0.010557625006) * y
  83                              + 0.011630447319) * y - 0.009279453341) * y
  84                              + 0.005353579108) * y - 0.002141268741) * y
  85                              + 0.000535310849) * y + 0.999936657524
  86     if z > 0.0:
  87         return ((x + 1.0) * 0.5)
  88     else:
  89         return ((1.0 - x) * 0.5)
  90
  91 cdef double ex20 ( double x ):
  92     """Wrapper on exp function. It will return 0 if x is smaller than -20.
  93     """
  94     if x < -20.0:
  95         return 0.0
  96     else:
  97         return exp(x)
  98
  99 cpdef double chisq_pvalue_e ( double x, unsigned int df ):
 100     """Chisq CDF function for upper tail and even degree of freedom.
 101     Good for p-value calculation and designed for combining pvalues.
 102
 103     Note: we assume df to be an even number larger than 1! Do not
 104     violate this assumption and there is no checking.
 105
 106     df has to be even number! if df is odd, result will be wrong!
 107
 108     """
 109     cdef:
 110         double  a, y, s
 111         double  e, c, z
 112
 113     if x <= 0.0:
 114         return 1.0
 115
 116     a = 0.5 * x
 117     even = ((2*(df/2)) == df)
 118     y = ex20(-a)
 119     s = y
 120     if df > 2:
 121         x = 0.5 * (df - 1.0)
 122         z = 1.0
 123         if a > bigx:            # approximation for big number
 124             e = 0.0
 125             c = log (a)
 126             while z <= x:
 127                 e = log (z) + e
 128                 s += ex20(c*z-a-e)
 129                 z += 1.0
 130             return s
 131         else:
 132             e = 1.0
 133             c = 0.0
 134             while z <= x:
 135                 e = e * (a / z)
 136                 c = c + e
 137                 z += 1.0
 138             return c * y + s
 139     else:
 140         return s
 141
 142 cpdef double chisq_logp_e ( double x, unsigned int df, bool log10 = False ):
 143     """Chisq CDF function in log space for upper tail and even degree of freedom
 144     Good for p-value calculation and designed for combining pvalues.
 145
 146     Note: we assume df to be an even number larger than 1! Do not
 147     violate this assumption and there is no checking.
 148
 149     Return value is -logp. If log10 is set as True, return -log10p
 150     instead.
 151
 152     """
 153     cdef:
 154         double a, y
 155         double s                # s is the return value
 156         double e, c, z
 157
 158     if x <= 0.0:
 159         return 0.0
 160
 161     a = 0.5 * x
 162     y = exp(-a)             # y is for small number calculation
 163     # initialize s
 164     s = -a
 165     if df > 2:
 166         x = 0.5 * (df - 1.0)    # control number of iterations
 167         z = 1.0             # variable for iteration
 168         if a > bigx:            # approximation for big number
 169             e = 0.0         # constant
 170             c = log (a)         # constant
 171             while z <= x:       # iterations
 172                 e += log(z)
 173                 s = logspace_add(s, c*z-a-e)
 174                 z += 1.0
 175         else:                   # for small number
 176             e = 1.0             # not a constant
 177             c = 0.0             # not a constant
 178             while z <= x:
 179                 e = e * (a / z)
 180                 c = c + e
 181                 z += 1.0
 182             s = log(y+c*y) #logspace_add( s, log(c) ) - a
 183     # return
 184     if log10:
 185         return -s/log(10)
 186     else:
 187         return -s
 188
 189 cpdef double poisson_cdf ( unsigned int n, double lam, bool lower=False, bool log10=False ):
 190     """Poisson CDF evaluater.
 191
 192     This is a more stable CDF function. It can tolerate large lambda
 193     value. While the lambda is larger than 700, the function will be a
 194     little slower.
 195
 196     Parameters:
 197     n     : your observation
 198     lam   : lambda of poisson distribution
 199     lower : if lower is False, calculate the upper tail CDF, otherwise, to calculate lower tail; Default is False.
 200     log10 : if log10 is True, calculation will be in log space. Default is False.
 201     """
 202     assert lam > 0.0, "Lambda must > 0, however we got %d" % lam
 203
 204     if log10:
 205         if lower:
 206             # lower tail
 207             return log10_poisson_cdf_P_large_lambda(n, lam)
 208         else:
 209             # upper tail
 210             return log10_poisson_cdf_Q_large_lambda(n, lam)
 211
 212     if lower:
 213         if lam > 700: # may be problematic
 214             return __poisson_cdf_large_lambda (n, lam)
 215         else:
 216             return __poisson_cdf(n,lam)
 217     else:
 218         # upper tail
 219         if lam > 700: # may be problematic
 220             return __poisson_cdf_Q_large_lambda (n, lam)
 221         else:
 222             return __poisson_cdf_Q(n,lam)
 223
 224 cdef inline double __poisson_cdf ( unsigned int k, double a ):
 225     """Poisson CDF For small lambda. If a > 745, this will return
 226     incorrect result.
 227
 228     Parameters:
 229     k   : observation
 230     a   : lambda
 231     """
 232     cdef:
 233         double nextcdf
 234         double cdf
 235         unsigned int i
 236         double lastcdf
 237
 238     if k < 0:
 239         return 0.0                        # special cases
 240
 241     nextcdf = exp( -1 * a )
 242     cdf = nextcdf
 243
 244     for i in range( 1, k + 1 ):
 245         lastcdf = nextcdf
 246         nextcdf = lastcdf * a / i
 247         cdf = cdf + nextcdf
 248     if cdf > 1.0:
 249         return 1.0
 250     else:
 251         return cdf
 252
 253 cdef inline double __poisson_cdf_large_lambda ( unsigned int k, double a ):
 254     """Slower poisson cdf for large lambda ( > 700 )
 255
 256     Parameters:
 257     k   : observation
 258     a   : lambda
 259     """
 260     cdef:
 261         int num_parts
 262         double lastexp
 263         double nextcdf
 264         double cdf
 265         unsigned int i
 266         double lastcdf
 267
 268     assert a > 700
 269     if k < 0:
 270         return 0.0                        # special cases
 271
 272     num_parts = int( a / LSTEP )
 273     lastexp = exp( -1 * ( a % LSTEP ) )
 274     nextcdf = EXPSTEP
 275
 276     num_parts -= 1
 277
 278     for i in range( 1 , k + 1 ):
 279         lastcdf = nextcdf
 280         nextcdf = lastcdf * a / i
 281         cdf = cdf + nextcdf
 282         if nextcdf > EXPTHRES or cdf > EXPTHRES:
 283            if num_parts>=1:
 284                cdf *= EXPSTEP
 285                nextcdf *= EXPSTEP
 286                num_parts -= 1
 287            else:
 288                cdf *= lastexp
 289                lastexp = 1
 290
 291     for i in range( num_parts ):
 292         cdf *= EXPSTEP
 293     cdf *= lastexp
 294     return cdf
 295
 296 cdef inline double __poisson_cdf_Q ( unsigned int k, double a ):
 297     """internal Poisson CDF evaluater for upper tail with small
 298     lambda.
 299
 300     Parameters:
 301     k   : observation
 302     a   : lambda
 303     """
 304     cdef unsigned int i
 305
 306     if k < 0:
 307         return 1.0                        # special cases
 308     cdef double nextcdf
 309     nextcdf = exp( -1 * a )
 310     cdef double lastcdf
 311
 312     for i in xrange(1,k+1):
 313         lastcdf = nextcdf
 314         nextcdf = lastcdf * a / i
 315
 316     cdef double cdf = 0.0
 317     i = k+1
 318     while nextcdf >0.0:
 319         lastcdf = nextcdf
 320         nextcdf = lastcdf * a / i
 321         cdf += nextcdf
 322         i+=1
 323     return cdf
 324
 325 cdef inline double __poisson_cdf_Q_large_lambda ( unsigned int k, double a ):
 326     """Slower internal Poisson CDF evaluater for upper tail with large
 327     lambda.
 328
 329     Parameters:
 330     k   : observation
 331     a   : lambda
 332     """
 333     assert a > 700
 334     if k < 0:
 335         return 1.0                        # special cases
 336     cdef unsigned int num_parts = int(a/LSTEP)
 337     cdef double lastexp = exp(-1 * (a % LSTEP) )
 338     cdef double nextcdf = EXPSTEP
 339     cdef unsigned int i
 340     cdef double lastcdf
 341
 342     num_parts -= 1
 343
 344     for i in xrange(1,k+1):
 345         lastcdf = nextcdf
 346         nextcdf = lastcdf * a / i
 347         if nextcdf > EXPTHRES:
 348            if num_parts>=1:
 349                nextcdf *= EXPSTEP
 350                num_parts -= 1
 351            else:
 352                # simply raise an error
 353                raise Exception("Unexpected error")
 354                #cdf *= lastexp
 355                #lastexp = 1
 356     cdef double cdf = 0.0
 357     i = k+1
 358     while nextcdf >0.0:
 359         lastcdf = nextcdf
 360         nextcdf = lastcdf * a / i
 361         cdf += nextcdf
 362         i+=1
 363         if nextcdf > EXPTHRES or cdf > EXPTHRES:
 364            if num_parts>=1:
 365                cdf *= EXPSTEP
 366                nextcdf *= EXPSTEP
 367                num_parts -= 1
 368            else:
 369                cdf *= lastexp
 370                lastexp = 1
 371
 372     for i in xrange(num_parts):
 373         cdf *= EXPSTEP
 374     cdf *= lastexp
 375     return cdf
 376
 377 cdef inline double log10_poisson_cdf_P_large_lambda ( unsigned int k, double lbd ):
 378     """Slower Poisson CDF evaluater for lower tail which allow
 379     calculation in log space. Better for the pvalue < 10^-310.
 380
 381     Parameters:
 382     k   : observation
 383     lbd : lambda
 384
 385     ret = -lambda + \ln( \sum_{i=k+1}^{\inf} {lambda^i/i!} = -lambda + \ln( sum{ exp{ln(F)} } ), where F=lambda^m/m!
 386     \ln{F(m)} = m*ln{lambda} - \sum_{x=1}^{m}\ln(x)
 387     Calculate \ln( sum{exp{N} ) by logspace_add function
 388
 389     Return the log10(pvalue)
 390     """
 391     cdef double residue = 0
 392     cdef double logx = 0
 393     cdef double ln_lbd = log(lbd)
 394
 395     # first residue
 396     cdef int m = k
 397     cdef double sum_ln_m = 0
 398     cdef int i = 0
 399     for i in range(1,m+1):
 400         sum_ln_m += log(i)
 401     logx = m*ln_lbd - sum_ln_m
 402     residue = logx
 403
 404     while m > 1:
 405         m -= 1
 406         logy = logx-ln_lbd+log(m)
 407         pre_residue = residue
 408         residue = logspace_add(pre_residue,logy)
 409         if fabs(pre_residue-residue) < 1e-10:
 410             break
 411         logx = logy
 412
 413     return round((residue-lbd)/M_LN10,5)
 414
 415 cdef inline double log10_poisson_cdf_Q_large_lambda ( unsigned int k, double lbd ):
 416     """Slower Poisson CDF evaluater for upper tail which allow
 417     calculation in log space. Better for the pvalue < 10^-310.
 418
 419     Parameters:
 420     k   : observation
 421     lbd : lambda
 422
 423     ret = -lambda + \ln( \sum_{i=k+1}^{\inf} {lambda^i/i!} = -lambda + \ln( sum{ exp{ln(F)} } ), where F=lambda^m/m!
 424     \ln{F(m)} = m*ln{lambda} - \sum_{x=1}^{m}\ln(x)
 425     Calculate \ln( sum{exp{N} ) by logspace_add function
 426
 427     Return the log10(pvalue)
 428     """
 429     cdef double residue = 0
 430     cdef double logx = 0
 431     cdef double ln_lbd = log(lbd)
 432
 433     # first residue
 434     cdef int m = k+1
 435     cdef double sum_ln_m = 0
 436     cdef int i = 0
 437     for i in range(1,m+1):
 438         sum_ln_m += log(i)
 439     logx = m*ln_lbd - sum_ln_m
 440     residue = logx
 441
 442     while True:
 443         m += 1
 444         logy = logx+ln_lbd-log(m)
 445         pre_residue = residue
 446         residue = logspace_add(pre_residue,logy)
 447         if fabs(pre_residue-residue) < 1e-5:
 448             break
 449         logx = logy
 450
 451     return round((residue-lbd)/log(10),5)
 452
 453 cdef inline double logspace_add ( double logx, double logy ):
 454     """addition in log space.
 455
 456     Given two log values, such as logx and logy, return
 457     log(exp(logx)+exp(logy)).
 458
 459     """
 460     return max (logx, logy) + log1p (exp (-fabs (logx - logy)))
 461
 462 cpdef poisson_cdf_inv ( double cdf, double lam, int maximum=1000 ):
 463     """inverse poisson distribution.
 464
 465     cdf : the CDF
 466     lam : the lambda of poisson distribution
 467
 468     note: maxmimum return value is 1000
 469     and lambda must be smaller than 740.
 470     """
 471     assert lam < 740
 472     if cdf < 0 or cdf > 1:
 473         raise Exception ("CDF must >= 0 and <= 1")
 474     elif cdf == 0:
 475         return 0
 476     cdef double sum2 = 0
 477     cdef double newval = exp( -1*lam )
 478     sum2 = newval
 479
 480     cdef int i
 481     cdef double sumold
 482     cdef double lastval
 483
 484     for i in xrange(1,maximum+1):
 485         sumold = sum2
 486         lastval = newval
 487         newval = lastval * lam / i
 488         sum2 = sum2 + newval
 489         if sumold <= cdf and cdf <= sum2:
 490             return i
 491
 492     return maximum
 493
 494 cpdef poisson_cdf_Q_inv ( double cdf, double lam, int maximum=1000 ):
 495     """inverse poisson distribution.
 496
 497     cdf : the CDF
 498     lam : the lambda of poisson distribution
 499
 500     note: maxmimum return value is 1000
 501     and lambda must be smaller than 740.
 502     """
 503     assert lam < 740
 504     if cdf < 0 or cdf > 1:
 505         raise Exception ("CDF must >= 0 and <= 1")
 506     elif cdf == 0:
 507         return 0
 508     cdef double sum2 = 0
 509     cdef double newval = exp( -1 * lam )
 510     sum2 = newval
 511
 512     cdef int i
 513     cdef double lastval
 514     cdef double sumold
 515
 516     for i in xrange(1,maximum+1):
 517         sumold = sum2
 518         lastval = newval
 519         newval = lastval * lam / i
 520         sum2 = sum2 + newval
 521         if sumold <= cdf and cdf <= sum2:
 522             return i
 523
 524     return maximum
 525
 526 cpdef poisson_pdf ( unsigned int k, double a ):
 527     """Poisson PDF.
 528
 529     PDF(K,A) is the probability that the number of events observed in
 530     a unit time period will be K, given the expected number of events
 531     in a unit time as A.
 532     """
 533     if a <= 0:
 534         return 0
 535     return exp(-a) * pow (a, k, None) / factorial (k)
 536
 537
 538 cdef binomial_coef ( long n, long k ):
 539     """BINOMIAL_COEF computes the Binomial coefficient C(N,K)
 540
 541     n,k are integers.
 542     """
 543     cdef long mn = min (k, n-k)
 544     cdef long mx
 545     cdef double cnk
 546     cdef long i
 547     if mn < 0:
 548         return 0
 549     elif mn == 0:
 550         return 1
 551     else:
 552         mx = max(k,n-k)
 553         cnk = float(mx+1)
 554         for i in xrange(2,mn+1):
 555             cnk = cnk * (mx+i) / i
 556     return cnk
 557
 558 cpdef binomial_cdf ( long x, long a, double b, bool lower=True ):
 559     """ BINOMIAL_CDF compute the binomial CDF.
 560
 561     CDF(x)(A,B) is the probability of at most X successes in A trials,
 562     given that the probability of success on a single trial is B.
 563     """
 564     if lower:
 565         return _binomial_cdf_f (x,a,b)
 566     else:
 567         return _binomial_cdf_r (x,a,b)
 568
 569 cpdef binomial_sf ( long x, long a, double b, bool lower=True ):
 570     """ BINOMIAL_SF compute the binomial survival function (1-CDF)
 571
 572     SF(x)(A,B) is the probability of more than X successes in A trials,
 573     given that the probability of success on a single trial is B.
 574     """
 575     if lower:
 576         return 1.0 - _binomial_cdf_f (x,a,b)
 577     else:
 578         return 1.0 - _binomial_cdf_r (x,a,b)
 579
 580 cpdef pduplication (np.ndarray[np.float64_t] pmf, int N_obs):
 581     """return the probability of a duplicate fragment given a pmf
 582     and a number of observed fragments N_obs
 583     """
 584     cdef:
 585         n = pmf.shape[0]
 586         float p, sf = 0.0
 587     for p in pmf:
 588         sf += binomial_sf(2, N_obs, p)
 589     return sf / <float>n
 590
 591 cdef _binomial_cdf_r ( long x, long a, double b ):
 592     """ Binomial CDF for upper tail.
 593
 594     """
 595     if x < 0:
 596         return 1
 597     elif a < x:
 598         return 0
 599     elif b == 0:
 600         return 0
 601     elif b == 1:
 602         return 1
 603
 604     cdef long argmax=int(a*b)
 605     cdef double seedpdf
 606     cdef double cdf
 607     cdef double pdf
 608     cdef long i
 609
 610     if x<argmax:
 611         seedpdf=binomial_pdf(argmax,a,b)
 612         pdf=seedpdf
 613         cdf = pdf
 614         for i in xrange(argmax-1,x,-1):
 615             pdf/=(a-i)*b/(1-b)/(i+1)
 616             if pdf==0.0: break
 617             cdf += pdf
 618
 619         pdf = seedpdf
 620         i = argmax
 621         while True:
 622             pdf*=(a-i)*b/(1-b)/(i+1)
 623             if pdf==0.0: break
 624             cdf+=pdf
 625             i+=1
 626         cdf=min(1,cdf)
 627         cdf = float("%.10e" %cdf)
 628         return cdf
 629     else:
 630         pdf=binomial_pdf(x+1,a,b)
 631         cdf = pdf
 632         i = x+1
 633         while True:
 634             pdf*=(a-i)*b/(1-b)/(i+1)
 635             if pdf==0.0: break
 636             cdf += pdf
 637             i+=1
 638         cdf=min(1,cdf)
 639         cdf = float("%.10e" %cdf)
 640         return cdf
 641
 642
 643 cdef _binomial_cdf_f ( long x, long a, double b ):
 644     """ Binomial CDF for lower tail.
 645
 646     """
 647     if x < 0:
 648         return 0
 649     elif a < x:
 650         return 1
 651     elif b == 0:
 652         return 1
 653     elif b == 1:
 654         return 0
 655
 656     cdef long argmax=int(a*b)
 657     cdef double seedpdf
 658     cdef double cdf
 659     cdef double pdf
 660     cdef long i
 661
 662     if x>argmax:
 663         seedpdf=binomial_pdf(argmax,a,b)
 664         pdf=seedpdf
 665         cdf = pdf
 666         for i in xrange(argmax-1,-1,-1):
 667             pdf/=(a-i)*b/(1-b)/(i+1)
 668             if pdf==0.0: break
 669             cdf += pdf
 670
 671         pdf = seedpdf
 672         for i in xrange(argmax,x):
 673             pdf*=(a-i)*b/(1-b)/(i+1)
 674             if pdf==0.0: break
 675             cdf+=pdf
 676         cdf=min(1,cdf)
 677         cdf = float("%.10e" %cdf)
 678         return cdf
 679     else:
 680         pdf=binomial_pdf(x,a,b)
 681         cdf = pdf
 682         for i in xrange(x-1,-1,-1):
 683             pdf/=(a-i)*b/(1-b)/(i+1)
 684             if pdf==0.0: break
 685             cdf += pdf
 686         cdf=min(1,cdf)
 687         cdf = float("%.10e" %cdf)
 688         return cdf
 689
 690 cpdef binomial_cdf_inv ( double cdf, long a, double b ):
 691     """BINOMIAL_CDF_INV inverts the binomial CDF. For lower tail only!
 692
 693     """
 694     if cdf < 0 or cdf >1:
 695         raise Exception("CDF must >= 0 or <= 1")
 696
 697     cdef double cdf2 = 0
 698     cdef long x
 699
 700     for x in xrange(0,a+1):
 701         pdf = binomial_pdf (x,a,b)
 702         cdf2 = cdf2 + pdf
 703         if cdf < cdf2:
 704             return x
 705     return a
 706
 707 cpdef binomial_pdf( long x, long a, double b ):
 708     """binomial PDF by H. Gene Shin
 709
 710     """
 711
 712     if a<1:
 713         return 0
 714     elif x<0 or a<x:
 715         return 0
 716     elif b==0:
 717         if x==0:
 718             return 1
 719         else:
 720             return 0
 721     elif b==1:
 722         if x==a:
 723             return 1
 724         else:
 725             return 0
 726
 727     cdef double p
 728     cdef long mn
 729     cdef long mx
 730     cdef double pdf
 731     cdef long t
 732     cdef long q
 733
 734     if x>a-x:
 735         p=1-b
 736         mn=a-x
 737         mx=x
 738     else:
 739         p=b
 740         mn=x
 741         mx=a-x
 742     pdf=1
 743     t = 0
 744     for q in xrange(1,mn+1):
 745         pdf*=(a-q+1)*p/(mn-q+1)
 746         if pdf < 1e-100:
 747             while pdf < 1e-3:
 748                 pdf /= 1-p
 749                 t-=1
 750         if pdf > 1e+100:
 751             while pdf > 1e+3 and t<mx:
 752                 pdf *= 1-p
 753                 t+=1
 754
 755     for i in xrange(mx-t):
 756         pdf *= 1-p
 757
 758     pdf=float("%.10e" % pdf)
 759     return pdf
 760
 761 # cdef normal_01_cdf ( double x ):
 762 #     """NORMAL_01_CDF evaluates the Normal 01 CDF.
 763 #     """
 764 #     cdef double a1 = 0.398942280444
 765 #     cdef double a2 = 0.399903438504
 766 #     cdef double a3 = 5.75885480458
 767 #     cdef double a4 = 29.8213557808
 768 #     cdef double a5 = 2.62433121679
 769 #     cdef double a6 = 48.6959930692
 770 #     cdef double a7 = 5.92885724438
 771 #     cdef double b0 = 0.398942280385
 772 #     cdef double b1 = 3.8052E-08
 773 #     cdef double b2 = 1.00000615302
 774 #     cdef double b3 = 3.98064794E-04
 775 #     cdef double b4 = 1.98615381364
 776 #     cdef double b5 = 0.151679116635
 777 #     cdef double b6 = 5.29330324926
 778 #     cdef double b7 = 4.8385912808
 779 #     cdef double b8 = 15.1508972451
 780 #     cdef double b9 = 0.742380924027
 781 #     cdef double b10 = 30.789933034
 782 #     cdef double b11 = 3.99019417011
 783 #     cdef double cdf
 784
 785 #     if abs ( x ) <= 1.28:
 786 #         y = 0.5 * x * x
 787 #         q = 0.5 - abs ( x ) * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5 + a6 / ( y + a7 ) ) ) )
 788 #     elif abs ( x ) <= 12.7:
 789 #         y = 0.5 * x * x
 790
 791 #         q = exp ( - y ) * b0 / ( abs ( x ) - b1 \
 792 #                                 + b2 / ( abs ( x ) + b3 \
 793 #                                 + b4 / ( abs ( x ) - b5 \
 794 #                                 + b6 / ( abs ( x ) + b7 \
 795 #                                 - b8 / ( abs ( x ) + b9 \
 796 #                                 + b10 / ( abs ( x ) + b11 ) ) ) ) ) )
 797 #     else :
 798 #         q = 0.0
 799 #     #
 800 #     #  Take account of negative X.
 801 #     #
 802 #     if x < 0.0:
 803 #         cdf = q
 804 #     else:
 805 #         cdf = 1.0 - q
 806
 807 #     return cdf
 808
 809 # def normal_cdf_inv(p, mu = None, sigma2 = None, lower=True):
 810
 811 #     upper = not lower
 812 #     if p < 0 or p > 1:
 813 #         raise Exception("Illegal argument %f for qnorm(p)." % p)
 814
 815 #     split = 0.42
 816 #     a0 = 2.50662823884
 817 #     a1 = -18.61500062529
 818 #     a2 = 41.39119773534
 819 #     a3 = -25.44106049637
 820 #     b1 = -8.47351093090
 821 #     b2 = 23.08336743743
 822 #     b3 = -21.06224101826
 823 #     b4 = 3.13082909833
 824 #     c0 = -2.78718931138
 825 #     c1 = -2.29796479134
 826 #     c2 = 4.85014127135
 827 #     c3 = 2.32121276858
 828 #     d1 = 3.54388924762
 829 #     d2 = 1.63706781897
 830 #     q = p - 0.5
 831
 832 #     r = 0.0
 833 #     ppnd = 0.0
 834
 835 #     if abs(q) <= split:
 836 #         r = q * q
 837 #         ppnd = q * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + 1)
 838 #     else:
 839 #         r = p
 840 #         if q > 0:
 841 #             r = 1 - p
 842
 843 #         if r > 0:
 844 #             r = math.sqrt(- math.log(r))
 845 #             ppnd = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1)
 846
 847 #             if q < 0:
 848 #                 ppnd = - ppnd
 849 #         else:
 850 #             ppnd = 0
 851
 852 #     if upper:
 853 #         ppnd = - ppnd
 854
 855 #     if mu != None and sigma2 != None:
 856 #         return ppnd * math.sqrt(sigma2) + mu
 857 #     else:
 858 #         return ppnd
 859
 860 # def normal_cdf (z, mu = 0.0, sigma2 = 1.0, lower=True):
 861
 862 #     upper = not lower
 863
 864 #     z = (z - mu) / math.sqrt(sigma2)
 865
 866 #     ltone = 7.0
 867 #     utzero = 18.66
 868 #     con = 1.28
 869 #     a1 = 0.398942280444
 870 #     a2 = 0.399903438504
 871 #     a3 = 5.75885480458
 872 #     a4 = 29.8213557808
 873 #     a5 = 2.62433121679
 874 #     a6 = 48.6959930692
 875 #     a7 = 5.92885724438
 876 #     b1 = 0.398942280385
 877 #     b2 = 3.8052e-8
 878 #     b3 = 1.00000615302
 879 #     b4 = 3.98064794e-4
 880 #     b5 = 1.986153813664
 881 #     b6 = 0.151679116635
 882 #     b7 = 5.29330324926
 883 #     b8 = 4.8385912808
 884 #     b9 = 15.1508972451
 885 #     b10 = 0.742380924027
 886 #     b11 = 30.789933034
 887 #     b12 = 3.99019417011
 888
 889 #     y = 0.0
 890 #     alnorm = 0.0
 891
 892 #     if z < 0:
 893 #         upper = not upper
 894 #         z = - z
 895
 896 #     if z <= ltone or upper and z <= utzero:
 897 #         y = 0.5 * z * z
 898 #         if z > con:
 899 #             alnorm = b1 * math.exp(- y) / (z - b2 + b3 / (z + b4 + b5 / (z - b6 + b7 / (z + b8 - b9 / (z + b10 + b11 / (z + b12))))))
 900 #         else:
 901 #             alnorm = 0.5 - z * (a1 - a2 * y / (y + a3 - a4 / (y + a5 + a6 / (y + a7))))
 902 #     else:
 903 #         alnorm = 0.0
 904 #     if not upper:
 905 #         alnorm = 1.0 - alnorm
 906 #     return alnorm