a few statistical helpers i'm using to figure some algorithms out

2 files changed, 61 insertions, 0 deletions

diff --git a/utils/m.h b/utils/m.h
index 5e45efee..dc881b36 100644
--- a/utils/m.h
+++ b/utils/m.h
@@ -4,6 +4,10 @@
 #include <cassert>
 #include <cmath>
 #include <boost/math/special_functions/digamma.hpp>
+#include <boost/math/constants/constants.hpp>
+
+// TODO right now I sometimes assert that x is in the support of the distributions
+// should be configurable to return -inf instead
 
 template <typename F>
 struct M {
@@ -59,6 +63,47 @@ struct M {
     return (alpha-1)*std::log(x)+(beta-1)*std::log(1-x) - log_beta_fn(alpha, beta);
   }
 
+  // support x \in R
+  static inline F log_laplace_density(const F& x, const F& mu, const F& b) {
+    assert(b > 0.0);
+    return -std::log(2*b) - std::fabs(x - mu) / b;
+  }
+
+  // support x \in R
+  // this is NOT the "log normal" density, it is the log of the "normal density at x"
+  static inline F log_gaussian_density(const F& x, const F& mu, const F& var) {
+    assert(var > 0.0);
+    return -0.5 * std::log(var * 2 * boost::math::constants::pi<F>()) - (x - mu)*(x - mu) / (2 * var);
+  }
+
+  // (x1,x2) \in R^2
+  // parameterized in terms of two means, a two "variances", a correlation < 1
+  static inline F log_bivariate_gaussian_density(const F& x1, const F& x2,
+                                                 const F& mu1, const F& mu2,
+                                                 const F& var1, const F& var2,
+                                                 const F& cor) {
+    assert(var1 > 0);
+    assert(var2 > 0);
+    assert(std::fabs(cor) < 1.0);
+    const F cor2 = cor*cor;
+    const F var1var22 = var1 * var2;
+    const F Z = 0.5 * std::log(var1var22 * (1 - cor2)) + std::log(2 * boost::math::constants::pi<F>());
+    return -Z -1.0 / (2 * (1 - cor2)) * ((x1 - mu1)*(x1-mu1) / var1 + (x2-mu2)*(x2-mu2) / var2 - 2*cor*(x1 - mu1)*(x2-mu2) / std::sqrt(var1var22));
+  }
+
+  // support x \in [a,b]
+  static inline F log_triangle_density(const F& x, const F& a, const F& b, const F& c) {
+    assert(a < b);
+    assert(a <= c);
+    assert(c <= b);
+    assert(x >= a);
+    assert(x <= b);
+    if (x <= c)
+      return std::log(2) + std::log(x - a) - std::log(b - a) - std::log(c - a);
+    else
+      return std::log(2) + std::log(b - x) - std::log(b - a) - std::log(b - c);
+  }
+
   // note: this has been adapted so that 0 is in the support of the distribution
   // support [0, 1, 2 ...)
   static inline F log_yule_simon(unsigned x, const F& rho) {
diff --git a/utils/m_test.cc b/utils/m_test.cc
index fca8f895..c4d6a166 100644
--- a/utils/m_test.cc
+++ b/utils/m_test.cc
@@ -14,6 +14,22 @@ class MTest : public testing::Test {
   virtual void TearDown() { }
 };
 
+TEST_F(MTest, Densities) {
+  double px1 = Md::log_gaussian_density(1.0, 0.0, 1.0);
+  double px2 = Md::log_gaussian_density(-1.0, 0.0, 1.0);
+  double py1 = Md::log_laplace_density(1.0, 0.0, 1.0);
+  double py2 = Md::log_laplace_density(1.0, 0.0, 1.0);
+  double pz1 = Md::log_triangle_density(1.0, -2.0, 2.0, 0.0);
+  double pz2 = Md::log_triangle_density(1.0, -2.0, 2.0, 0.0);
+  cerr << px1 << " " << py1 << " " << pz2 << endl;
+  EXPECT_FLOAT_EQ(px1, px2);
+  EXPECT_FLOAT_EQ(py1, py2);
+  EXPECT_FLOAT_EQ(pz1, pz2);
+  double b1 = Md::log_bivariate_gaussian_density(1.0, -1.0, 0.0, 0.0, 1.0, 1.0, -0.8);
+  double b2 = Md::log_bivariate_gaussian_density(-1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -0.8);
+  cerr << b1 << " " << b2 << endl;
+}
+
 TEST_F(MTest, Poisson) {
   double prev = 1.0;
   double tot = 0;