Skip to content
Permalink
Browse files

Fix wrong poisson sample with large lambda

Add a separate function to accurately generate a poisoon distribution
sample when the lambda is large.

PR-URL: #246
Fixes: #243
Ref: https://pdfs.semanticscholar.org/00f1/8468dd53e4e3f0c28a5d504f8cd10376a673.pdf
Reviewed-by: Trevor Norris <trev.norris@gmail.com>
  • Loading branch information...
Divyanshu Kalra authored and trevnorris committed Aug 12, 2019
1 parent b8fda65 commit a22a7f259e21dcdc3cdec675dc95b4d4e754827b
Showing with 171 additions and 3 deletions.
  1. +78 −1 dist/jstat.js
  2. +1 −1 dist/jstat.min.js
  3. +39 −1 src/distribution.js
  4. +39 −0 src/special.js
  5. +14 −0 test/distribution/poisson-test.js
@@ -1261,6 +1261,45 @@ jStat.gammaln = function gammaln(x) {
return Math.log(2.5066282746310005 * ser / xx) - tmp;
};

/*
* log-gamma function to support poisson distribution sampling. The
* algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their
* book "Computation of Special Functions", 1996, John Wiley & Sons, Inc.
*/
jStat.loggam = function loggam(x) {
var x0, x2, xp, gl, gl0;
var k, n;

var a = [8.333333333333333e-02, -2.777777777777778e-03,
7.936507936507937e-04, -5.952380952380952e-04,
8.417508417508418e-04, -1.917526917526918e-03,
6.410256410256410e-03, -2.955065359477124e-02,
1.796443723688307e-01, -1.39243221690590e+00];
x0 = x;
n = 0;
if ((x == 1.0) || (x == 2.0)) {
return 0.0;
}
if (x <= 7.0) {
n = Math.floor(7 - x);
x0 = x + n;
}
x2 = 1.0 / (x0 * x0);
xp = 2 * Math.PI;
gl0 = a[9];
for (k = 8; k >= 0; k--) {
gl0 *= x2;
gl0 += a[k];
}
gl = gl0 / x0 + 0.5 * Math.log(xp) + (x0 - 0.5) * Math.log(x0) - x0;
if (x <= 7.0) {
for (k = 1; k <= n; k++) {
gl -= Math.log(x0 - 1.0);
x0 -= 1.0;
}
}
return gl;
}

// gamma of x
jStat.gammafn = function gammafn(x) {
@@ -2739,13 +2778,51 @@ jStat.extend(jStat.poisson, {
return l;
},

sample: function sample(l) {
sampleSmall: function sampleSmall(l) {
var p = 1, k = 0, L = Math.exp(-l);
do {
k++;
p *= jStat._random_fn();
} while (p > L);
return k - 1;
},

sampleLarge: function sampleLarge(l) {
var lam = l;
var k;
var U, V, slam, loglam, a, b, invalpha, vr, us;

slam = Math.sqrt(lam);
loglam = Math.log(lam);
b = 0.931 + 2.53 * slam;
a = -0.059 + 0.02483 * b;
invalpha = 1.1239 + 1.1328 / (b - 3.4);
vr = 0.9277 - 3.6224 / (b - 2);

while (1) {
U = Math.random() - 0.5;
V = Math.random();
us = 0.5 - Math.abs(U);
k = Math.floor((2 * a / us + b) * U + lam + 0.43);
if ((us >= 0.07) && (V <= vr)) {
return k;
}
if ((k < 0) || ((us < 0.013) && (V > us))) {
continue;
}
/* log(V) == log(0.0) ok here */
/* if U==0.0 so that us==0.0, log is ok since always returns */
if ((Math.log(V) + Math.log(invalpha) - Math.log(a / (us * us) + b)) <= (-lam + k * loglam - jStat.loggam(k + 1))) {
return k;
}
}
},

sample: function sample(l) {
if (l < 10)
return this.sampleSmall(l);
else
return this.sampleLarge(l);
}
});

Large diffs are not rendered by default.

@@ -1027,13 +1027,51 @@ jStat.extend(jStat.poisson, {
return l;
},

sample: function sample(l) {
sampleSmall: function sampleSmall(l) {
var p = 1, k = 0, L = Math.exp(-l);
do {
k++;
p *= jStat._random_fn();
} while (p > L);
return k - 1;
},

sampleLarge: function sampleLarge(l) {
var lam = l;
var k;
var U, V, slam, loglam, a, b, invalpha, vr, us;

slam = Math.sqrt(lam);
loglam = Math.log(lam);
b = 0.931 + 2.53 * slam;
a = -0.059 + 0.02483 * b;
invalpha = 1.1239 + 1.1328 / (b - 3.4);
vr = 0.9277 - 3.6224 / (b - 2);

while (1) {
U = Math.random() - 0.5;
V = Math.random();
us = 0.5 - Math.abs(U);
k = Math.floor((2 * a / us + b) * U + lam + 0.43);
if ((us >= 0.07) && (V <= vr)) {
return k;
}
if ((k < 0) || ((us < 0.013) && (V > us))) {
continue;
}
/* log(V) == log(0.0) ok here */
/* if U==0.0 so that us==0.0, log is ok since always returns */
if ((Math.log(V) + Math.log(invalpha) - Math.log(a / (us * us) + b)) <= (-lam + k * loglam - jStat.loggam(k + 1))) {
return k;
}
}
},

sample: function sample(l) {
if (l < 10)
return this.sampleSmall(l);
else
return this.sampleLarge(l);
}
});

@@ -17,6 +17,45 @@ jStat.gammaln = function gammaln(x) {
return Math.log(2.5066282746310005 * ser / xx) - tmp;
};

/*
* log-gamma function to support poisson distribution sampling. The
* algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their
* book "Computation of Special Functions", 1996, John Wiley & Sons, Inc.
*/
jStat.loggam = function loggam(x) {
var x0, x2, xp, gl, gl0;
var k, n;

var a = [8.333333333333333e-02, -2.777777777777778e-03,
7.936507936507937e-04, -5.952380952380952e-04,
8.417508417508418e-04, -1.917526917526918e-03,
6.410256410256410e-03, -2.955065359477124e-02,
1.796443723688307e-01, -1.39243221690590e+00];
x0 = x;
n = 0;
if ((x == 1.0) || (x == 2.0)) {
return 0.0;
}
if (x <= 7.0) {
n = Math.floor(7 - x);
x0 = x + n;
}
x2 = 1.0 / (x0 * x0);
xp = 2 * Math.PI;
gl0 = a[9];
for (k = 8; k >= 0; k--) {
gl0 *= x2;
gl0 += a[k];
}
gl = gl0 / x0 + 0.5 * Math.log(xp) + (x0 - 0.5) * Math.log(x0) - x0;
if (x <= 7.0) {
for (k = 1; k <= n; k++) {
gl -= Math.log(x0 - 1.0);
x0 -= 1.0;
}
}
return gl;
}

// gamma of x
jStat.gammafn = function gammafn(x) {
@@ -40,6 +40,20 @@ suite.addBatch({

assert.epsilon(tol, mean(1), 1);
assert.epsilon(tol, mean(3.5), 3.5);
},
'check sample': function(jStat) {
var samplingTol = 10;
var lambdaToTest = [10, 100, 1000, 10000, 100000];
var sampleSize = 10000;
for(var i in lambdaToTest) {
var lambda = lambdaToTest[i];
var sum = 0;
for(let i = 0; i < sampleSize; i++) {
sum += jStat.poisson.sample(lambda);
}
var mean = sum/sampleSize;
assert.epsilon(samplingTol, mean, lambda);
}
}
}
});

0 comments on commit a22a7f2

Please sign in to comment.
You can’t perform that action at this time.