# lima.py - 2013-05-15 SJF
# Example of Li & Ma significance calculation
import math, scipy.stats

def ts_lima(non,noff,alpha):
    opa  = 1.0+alpha
    ntot = non+noff
    return 2.0*(non*math.log(opa*non/alpha/ntot) \
                + noff*math.log(opa*noff/ntot))

non    = 2808
noff   = 4959
alpha  = 1.0/3
T      = 27.2

S_hat  = (non - noff*alpha)/T
sig2_S = (non + noff*alpha**2)/T**2
ts     = ts_lima(non,noff,alpha)
signif = math.sqrt(ts)
Pval   = scipy.stats.chi2.sf(ts,1)

print S_hat, math.sqrt(sig2_S), ts, signif, Pval