gamma-poisson word length model

1 files changed, 33 insertions, 0 deletions

diff --git a/utils/gamma_poisson.h b/utils/gamma_poisson.h
new file mode 100644
index 00000000..fec763f6
--- /dev/null
+++ b/utils/gamma_poisson.h
@@ -0,0 +1,33 @@
+#ifndef _GAMMA_POISSON_H_
+#define _GAMMA_POISSON_H_
+
+#include <m.h>
+
+// http://en.wikipedia.org/wiki/Conjugate_prior
+struct GammaPoisson {
+  GammaPoisson(double shape, double rate) :
+    a(shape), b(rate), n(), marginal() {}
+
+  double prob(unsigned x) const {
+    return exp(Md::log_negative_binom(x, a + marginal, 1.0 - (b + n) / (1 + b + n)));
+  }
+
+  void increment(unsigned x) {
+    ++n;
+    marginal += x;
+  }
+
+  void decrement(unsigned x) {
+    --n;
+    marginal -= x;
+  }
+
+  double log_likelihood() const {
+    return 0;
+  }
+
+  double a, b;
+  int n, marginal;
+};
+
+#endif