# ul_lima_bayes_1d.py - 2013-05-25 SJF
# Bayesian upper limit in Li & Ma problem
from math import *
import scipy.stats, scipy.optimize, scipy.integrate, sys
# non, noff, alpha, T = (2808, 4959, 1.0/3, 27.2)
# non, noff, alpha, T = (15, 24, 1.0/3, 10.0)
non, noff, alpha, T = (4, 36, 1.0/3, 10.0)
C = 0.95; # Use 95% confidence region (C must be >0.5)
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
sig_S    = sqrt(non+noff*alpha**2)/T
logL_max = profileLogL(S_hat)
def logPrior(S):
    return log(1);
def logPosterior(S):
    return logPrior(S)+profileLogL(S)-logL_max
def integralPosterior(Smax):
    integrand = lambda S: exp(logPosterior(S))
    y, err = scipy.integrate.quad(integrand,0,Smax)
    return y
total_integral = integralPosterior(S_hat+100*sig_S);
root_fn  = lambda S: integralPosterior(S) - total_integral*C
S_ul = scipy.optimize.brentq(root_fn, 0, S_hat+100*sig_S)
print S_ul, integralPosterior(S_ul)/total_integral, total_integral