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