binomial.py 2.63 KB
Newer Older
Mikael Boden's avatar
Mikael Boden committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
from math import log, exp
import sys


MAXIT = 100
EPS = 3.0e-7
FPMIN = 1.0e-300
gamma_c = [76.18009172947146,
           -86.50532032941677,
           24.01409824083091,
           -1.23173957245,
           0.1208650973866179e-2,
           -0.5395239384953e-5]


def log_binomial_ncdf(N, k, p):
    """
    Log of one minus the cumulative distribution function of the binomial dist.

    The binomial density gives the probability of k successes in N independent
    trials each with probability p of success.
    """
    if (k==0):
        return 0
    else:
        return log_betai(k, N-k+1, p)


def betai (a, b, x):
    """
    Incomplete beta function
    """
    if (x<0 or x>1): die("Bad x=`" + str(x) + "' in routine betai")
    if (x==0 or x==1):
        bt = 0
    else:
        bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b)+a*log(x)+b*log(1-x))

    thresh = (a+1)/(a+b+2.0)
    if (x<thresh):
        return(bt*betacf(a,b,x)/a)
    else:
        return(1.0-bt*betacf(b,a,1.0-x)/b)


def log_betai(a, b, x):
    """
    log incomplete beta function
    """
    if (x<0 or x>1): die("Bad x=`" + str(x) + "' in routine betai")
    if (x==0 or x==1):
        log_bt = -1e300           # log(0)
    else:
        log_bt = gammaln(a+b)-gammaln(a)-gammaln(b)+a*log(x)+b*log(1.0-x)

    thresh = (a+1.0)/(a+b+2.0)
    if (x<thresh):
        return(log_bt + log(betacf(a,b,x)/a))
    else:
        return(log(1.0 - exp(log_bt)*betacf(b,a,1.0-x)/b))


def betacf(a, b, x):
    """
    used by betai
    """
    qab = a+b
    qap = a+1.0
    qam = a-1.0
    c = 1.0
    d = 1.0-qab*x/qap

    if (abs(d) < FPMIN): d = FPMIN
    d = 1.0/d
    h = d

    for m in range(1, MAXIT+1):
        m2 = 2.0*m
        aa = m*(b-m)*x/((qam+m2)*(a+m2))
        d=1.0+aa*d
        if (abs(d) < FPMIN): d=FPMIN
        c=1.0+aa/c
        if (abs(c) < FPMIN): c=FPMIN
        d = 1.0/d
        h *= d*c
        aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2))

        d=1.0+aa*d
        if (abs(d) < FPMIN): d=FPMIN
        c=1.0+aa/c
        if (abs(c) < FPMIN): c=FPMIN
        d = 1.0/d

        delta = d*c
        h *= delta
        if (abs(delta-1.0) < EPS): break

Mikael Boden's avatar
Mikael Boden committed
98 99
    if (m > MAXIT):  print(("a or b too big or MAXIT too small "
                                           "in betacf"), file=sys.stderr)
Mikael Boden's avatar
Mikael Boden committed
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
    return h


def gammaln(x):
    """
    Compute log gamma function
    """
    xx = x
    s = 1.000000000190015
    for i in range(0, 6):
        xx += 1
        s += gamma_c[i]/xx

    res = ((x+0.5) * log(x+5.5)) - (x+5.5) + log(2.5066282746310005*s/x)
    if (res >= 0):
        return res
    else:
        return 0					# avoid roundoff error


def die(string):
Mikael Boden's avatar
Mikael Boden committed
121
    print(string, file=sys.stderr)
Mikael Boden's avatar
Mikael Boden committed
122