1 files changed, 42 insertions, 0 deletions
diff --git a/misc.c b/misc.c
new file mode 100644
index 0000000..74c729e
--- /dev/null
+++ b/misc.c
@@ -0,0 +1,42 @@
+#include "misc.h"
+#include <math.h>
+#include <stdlib.h> /* for size_t */
+
+double logsumexp(double *a, size_t n)
+{
+    /* Returns the log of the sum of the exponentials of the array `a`.
+     *
+     * This function is designed to reduce underflow when the exponentials of
+     * `a` are very small, for example when computing probabilities. */
+    size_t i;
+    double amax, sum;
+
+    amax = a[0];
+    for (i = 0; i < n; i++) {
+        if (a[i] > amax) amax = a[i];
+    }
+
+    sum = 0.0;
+
+    for (i = 0; i < n; i++) {
+        sum += exp(a[i]-amax);
+    }
+
+    sum = log(sum);
+
+    return amax + sum;
+}
+
+double norm(double x, double mu, double sigma)
+{
+    /* Returns the PDF for a gaussian random variable with mean `mu` and
+     * standard deviation `sigma`. */
+    return exp(-pow(x-mu,2)/(2*pow(sigma,2)))/(sqrt(2*M_PI)*sigma);
+}
+
+double norm_cdf(double x, double mu, double sigma)
+{
+    /* Returns the CDF for a gaussian random variable with mean `mu` and
+     * standard deviation `sigma`. */
+    return erfc(-(x-mu)/(sqrt(2)*sigma))/2.0;
+}