# conf_lima_1d.py - 2013-05-25 SJF
# 1-D 2-sided confidence interval in Li & Ma problem
from math import *
import scipy.stats, scipy.optimize, sys
# non, noff, alpha, T = (2808, 4959, 1.0/3, 27.2)
non, noff, alpha, T = (15, 24, 1.0/3, 10.0)
C = 0.68; # Use 1-sigma confidence region
d2logL = scipy.stats.chi2.ppf(C,1)
def logL(S,B):
    return non*log(max((S+alpha*B)*T,sys.float_info.min)) + \
    noff*log(max(B*T,sys.float_info.min))-(S+(1+alpha)*B)*T
def profileLogL(S):
    opt_fn = lambda B: -logL(S,B)
    opt_res = scipy.optimize.minimize(opt_fn, 1)
    return -opt_res.fun
S_hat    = (non-noff*alpha)/T
B_hat    = noff/T
logL_max = logL(S_hat,B_hat)
sig_S    = sqrt(non+noff*alpha**2)/T
TS       = -2.0*(profileLogL(0)-logL_max)
root_fn  = lambda S: 2.0*(profileLogL(S)-logL_max)+d2logL
S_lo     = scipy.optimize.brentq(root_fn, 1e-8, S_hat)
S_hi     = scipy.optimize.brentq(root_fn, S_hat, 1e8)
print S_hat, S_lo-S_hat, S_hi-S_hat, sig_S, TS, sqrt(TS)