From 1e3e7ffe1e233b7a58cbd3b7245f39446868b3b2 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 17 Aug 2020 11:52:51 -0500 Subject: simplify likelihood calculation in chi2 --- utils/chi2 | 13 +++++-------- utils/sddm/plot_energy.py | 2 ++ 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/utils/chi2 b/utils/chi2 index 4515b91..3b519e2 100755 --- a/utils/chi2 +++ b/utils/chi2 @@ -110,14 +110,11 @@ def make_nll(data, mc_hists): nll = 0 for id in data_hists: - N = data_hists[id].sum() - nll -= poisson.logpmf(N,mc_hists[id].sum()*x[0]) - if N > 0: - p = mc_hists[id]/mc_hists[id].sum() - # Fix a bug in scipy(). See https://github.com/scipy/scipy/issues/8235 (I think). - p += 1e-10 - p /= p.sum() - nll -= multinomial.logpmf(data_hists[id],N,p) + oi = data_hists[id].sum() + ei = mc_hists[id] + EPSILON + N = oi.sum() + n = ei.sum() + nll -= -N - np.sum(gammaln(oi+1)) + np.sum(oi*np.log(ei)) return nll - norm.logpdf(x[1],ENERGY_SCALE_MEAN,ENERGY_SCALE_UNCERTAINTY) return nll diff --git a/utils/sddm/plot_energy.py b/utils/sddm/plot_energy.py index 84a7094..a4b7005 100755 --- a/utils/sddm/plot_energy.py +++ b/utils/sddm/plot_energy.py @@ -23,6 +23,8 @@ from . import grouper, print_warning, AV_RADIUS, PSUP_RADIUS, read_hdf import pandas as pd from scipy.misc import logsumexp +EPSILON = 1e-10 + def lnfact(n): return gammaln(n+1) -- cgit