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 from 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
  31 from cpython cimport bool
  32 # ------------------------------------
  33 # constants
  34 # ------------------------------------
  35 cdef int32_t LSTEP = 200
  36 cdef float64_t EXPTHRES = exp(LSTEP)
  37 cdef float64_t EXPSTEP  = exp(-1*LSTEP)
  38 cdef float64_t bigx = 20
  39
  40 # ------------------------------------
  41 # Normal distribution functions
  42 # ------------------------------------
  43 # x is the value, u is the mean, v is the variance
  44 cpdef 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
  53 cpdef 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
  65 cdef 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
  99 cdef 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
 107 cpdef 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
 150 cpdef 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
 197 cpdef 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
 232 cdef 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
 261 cdef 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
 304 cdef 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
 333 cdef 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
 384 cdef 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
 422 cdef 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
 460 cdef 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
 469 cpdef 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
 501 cpdef 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
 533 cpdef 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
 545 cdef 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
 565 cpdef 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
 576 cpdef 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
 587 cpdef 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
 598 cdef _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
 650 cdef _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
 697 cpdef 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
 714 cpdef 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