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