Merge pull request #830 from null-a/poisson

Fix Poisson sampler

src/dists.ad.js

@@ -1264,31 +1264,38 @@ function lnfactExact(x) {
   return t;
 }
 
+// This method comes from Ahrens and Dieters' 1974 paper "Computer
+// Methods for Sampling from Gamma, Beta, Poisson and Binomial
+// Distributions". The method is called Method PG, see page 240.
+
+function poissonSample(mu) {
+  var k = 0;
+  while (mu > 10) {
+    var m = Math.floor(7 / 8 * mu);
+    var x = gammaSample(m, 1);
+    if (x >= mu) {
+      return k + binomialSample(mu / x, m - 1);
+    } else {
+      mu -= x;
+      k += m;
+    }
+  }
+  var emu = Math.exp(-mu);
+  var p = 1;
+  do {
+    p *= util.random();
+    k++;
+  } while (p > emu);
+  return k - 1;
+}
+
 var Poisson = makeDistributionType({
   name: 'Poisson',
   desc: 'Distribution over integers.',
   params: [{name: 'mu', desc: 'mean', type: types.positiveReal}],
   wikipedia: true,
   sample: function() {
-    var k = 0;
-    var mu = ad.value(this.params.mu);
-    while (mu > 10) {
-      var m = 7 / 8 * mu;
-      var x = gammaSample(m, 1);
-      if (x > mu) {
-        return (k + binomialSample(mu / x, m - 1)) || 0;
-      } else {
-        mu -= x;
-        k += m;
-      }
-    }
-    var emu = Math.exp(-mu);
-    var p = 1;
-    do {
-      p *= util.random();
-      k++;
-    } while (p > emu);
-    return (k - 1) || 0;
+    return poissonSample(ad.value(this.params.mu));
   },
   score: function(val) {
     'use ad';
@@ -1614,6 +1621,7 @@ module.exports = _.assign({
   tensorLaplaceSample: tensorLaplaceSample,
   gammaSample: gammaSample,
   dirichletSample: dirichletSample,
+  poissonSample: poissonSample,
   // helpers
   serialize: serialize,
   deserialize: deserialize,

tests/test-data/sampler/poisson.js

@@ -0,0 +1,62 @@
+var assert = require('assert');
+var dists = require('../../../src/dists');
+var util = require('../../../src/util');
+var statistics = require('../../../src/math/statistics');
+
+var ln = Math.log,
+    pow = Math.pow,
+    sqrt = Math.sqrt,
+    abs = Math.abs;
+
+module.exports = {
+  name: 'poisson',
+  sampler: dists.poissonSample,
+  type: 'integer',
+  inSupport: function(params, x) {
+    return Number.isInteger(x) && x >= 0;
+  },
+  settings: [
+    {params: [0.5], n: 1e05, reltol: {mode: 0.05}},
+    {params: [1], n: 1e05, skip: ['mode']},
+    {params: [4], n: 1e05, skip: ['mode']},
+    {params: [4.5], n: 1e05, reltol: {mode: 0.05}},
+    {params: [10], n: 1e05, skip: ['mode']},
+    {params: [11], n: 1e05, skip: ['mode']},
+    {params: [11.5], n: 1e05, reltol: {mode: 0.05}},
+    {params: [20], n: 1e05, skip: ['mode']},
+    {params: [200], n: 1e05, skip: ['mode']},
+    {params: [200.5], n: 1e05, reltol: {mode: 0.05}},
+  ],
+  moment: function(params, N) {
+    assert.ok(N === 4, "Don't know how to compute moment N=" + N);
+    var mu = params[0];
+    // http://mathworld.wolfram.com/PoissonDistribution.html
+    return mu * (1 + 3 * mu);
+  },
+  populationStatisticFunctions: {
+    // https://en.wikipedia.org/wiki/Poisson_distribution
+    mean: function(params) {
+      var mu = params[0];
+      return mu;
+    },
+    mode: function(params) {
+      var mu = params[0];
+      var mode1 = Math.ceil(mu) - 1;
+      var mode2 = Math.floor(mu);
+      assert.ok(mode1 === mode2, "Don't know how to test multimodal distributions.");
+      return mode1;
+    },
+    variance: function(params) {
+      var mu = params[0];
+      return mu;
+    },
+    skew: function(params) {
+      var mu = params[0];
+      return 1 / Math.sqrt(mu);
+    },
+    kurtosis: function(params) {
+      var mu = params[0];
+      return 3 + (1 / mu);
+    }
+  }
+};

tests/test-samplers.js

@@ -46,7 +46,8 @@ var distMetadataList = [
   require('./test-data/sampler/gamma'),
   require('./test-data/sampler/binomial'),
   require('./test-data/sampler/beta'),
-  require('./test-data/sampler/gaussian')
+  require('./test-data/sampler/gaussian'),
+  require('./test-data/sampler/poisson')
 ];
 
 var generateSettingTest = function(seed, distMetadata, settings) {