From 1e3e7ffe1e233b7a58cbd3b7245f39446868b3b2 Mon Sep 17 00:00:00 2001
From: tlatorre <tlatorre@uchicago.edu>
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