# errors_poisson.py - 2013-05-07 SJF
# Evaluate the errors on the Poisson mean
import math, scipy.optimize
n_meas   = 2
logL     = lambda lam: n_meas*math.log(lam)-lam
opt_fn   = lambda lam: -logL(lam)
opt_res  = scipy.optimize.minimize(opt_fn, 1e-8)
lam_est  = opt_res.x[0]
logL_max = logL(lam_est)
root_fn  = lambda lam: 2.0*(logL(lam)-logL_max)+1.0
lam_lo   = scipy.optimize.brentq(root_fn, 1e-8, lam_est)
lam_hi   = scipy.optimize.brentq(root_fn, lam_est, 1e8)
print lam_est, lam_lo-lam_est, lam_hi-lam_est