# ul_lima_1d.py - 2013-05-25 SJF # Frequentest upper limit 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) non, noff, alpha, T = (4, 36, 1.0/3, 10.0) C = 0.95; # Use 95% confidence region (C must be >0.5) d2logL = scipy.stats.chi2.ppf(2*C-1,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) TS = -2.0*(profileLogL(0)-logL_max) root_fn = lambda S: 2.0*(profileLogL(S)-logL_max)+d2logL S_ul = scipy.optimize.brentq(root_fn, S_hat, 1e8) print S_hat, S_ul, TS, sqrt(TS)