fix the likelihood function to return the *negative* log likelihood of the path coefficients

Previously I was adding the log likelihood of the path coefficients instead of
the *negative* log likelihood! When fitting electrons this would sometimes
cause the fit to become unstable and continue increasing the path coefficients
without bound since the gain in the likelihood caused by increasing the
coefficients was more than the loss caused by a worse fit to the PMT data.

Doh!

4 files changed, 40 insertions, 3 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 7c3d10f..f426271 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -704,7 +704,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
 
     logp_path = 0.0;
     for (i = 0; i < n; i++)
-        logp_path += log(norm(z1[i],0,1)) + log(norm(z2[i],0,1));
+        logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1);
 
-    return kahan_sum(nll,nhit) + logp_path;
+    return kahan_sum(nll,nhit) - logp_path;
 }
diff --git a/src/misc.c b/src/misc.c
index 6f50033..ba0a090 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -384,6 +384,13 @@ double logsumexp(double *a, size_t n)
     return amax + sum;
 }
 
+double log_norm(double x, double mu, double sigma)
+{
+    /* Returns the log of the PDF for a gaussian random variable with mean `mu`
+     * and standard deviation `sigma`. */
+    return -pow(x-mu,2)/(2*pow(sigma,2)) - log(sqrt(2*M_PI)*sigma);
+}
+
 double norm(double x, double mu, double sigma)
 {
     /* Returns the PDF for a gaussian random variable with mean `mu` and
diff --git a/src/misc.h b/src/misc.h
index 46931cb..dadbbce 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -14,6 +14,7 @@ double interp1d(double x, double *xp, double *yp, size_t n);
 int isclose(double a, double b, double rel_tol, double abs_tol);
 int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol);
 double logsumexp(double *a, size_t n);
+double log_norm(double x, double mu, double sigma);
 double norm(double x, double mu, double sigma);
 double norm_cdf(double x, double mu, double sigma);
 double mean(const double *x, size_t n);
diff --git a/src/test.c b/src/test.c
index 2645982..20a66ec 100644
--- a/src/test.c
+++ b/src/test.c
@@ -823,7 +823,7 @@ err:
 
 int test_std(char *err)
 {
-    /* Tests the mean() function. */
+    /* Tests the std() function. */
     size_t i, j;
     double x[100];
     double sigma, expected, mu;
@@ -851,6 +851,34 @@ err:
     return 1;
 }
 
+int test_log_norm(char *err)
+{
+    /* Tests the log_norm() function. */
+    size_t i;
+    double x, mu, sigma, logp, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        x = randn();
+        mu = randn();
+        sigma = genrand_real2() + 1.0;
+
+        logp = log_norm(x,mu,sigma);
+        expected = log(norm(x,mu,sigma));
+
+        if (!isclose(logp, expected, 0, 1e-9)) {
+            sprintf(err, "log_norm(%.2g,%.2g,%.2g) returned %.5g, but expected %.5g", x, mu, sigma, logp, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -877,6 +905,7 @@ struct tests {
     {test_proton_get_dEdx, "test_proton_get_dEdx"},
     {test_proton_get_range, "test_proton_get_range"},
     {test_proton_get_energy, "test_proton_get_energy"},
+    {test_log_norm, "test_log_norm"},
 };
 
 int main(int argc, char **argv)