diff options
| -rw-r--r-- | src/likelihood.c | 4 | ||||
| -rw-r--r-- | src/misc.c | 7 | ||||
| -rw-r--r-- | src/misc.h | 1 | ||||
| -rw-r--r-- | src/test.c | 31 |
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; } @@ -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 @@ -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); @@ -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) |