123 lines
2.8 KiB
JavaScript
123 lines
2.8 KiB
JavaScript
|
/**
|
|||
|
|
* @license Apache-2.0
|
|||
|
|
*
|
|||
|
|
* Copyright (c) 2018 The Stdlib Authors.
|
|||
|
|
*
|
|||
|
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
|||
|
|
* you may not use this file except in compliance with the License.
|
|||
|
|
* You may obtain a copy of the License at
|
|||
|
|
*
|
|||
|
|
* http://www.apache.org/licenses/LICENSE-2.0
|
|||
|
|
*
|
|||
|
|
* Unless required by applicable law or agreed to in writing, software
|
|||
|
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
|||
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|||
|
|
* See the License for the specific language governing permissions and
|
|||
|
|
* limitations under the License.
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
'use strict';
|
|||
|
|
|
|||
|
|
// MODULES //
|
|||
|
|
|
|||
|
|
var factorialln = require( '@stdlib/math/base/special/factorialln' );
|
|||
|
|
var floor = require( '@stdlib/math/base/special/floor' );
|
|||
|
|
var sign = require( '@stdlib/math/base/special/signum' );
|
|||
|
|
var sqrt = require( '@stdlib/math/base/special/sqrt' );
|
|||
|
|
var abs = require( '@stdlib/math/base/special/abs' );
|
|||
|
|
var ln = require( '@stdlib/math/base/special/ln' );
|
|||
|
|
var LN_SQRT_TWO_PI = require( '@stdlib/constants/float64/ln-sqrt-two-pi' );
|
|||
|
|
|
|||
|
|
|
|||
|
|
// VARIABLES //
|
|||
|
|
|
|||
|
|
var ONE_12 = 1.0 / 12.0;
|
|||
|
|
var ONE_360 = 1.0 / 360.0;
|
|||
|
|
|
|||
|
|
|
|||
|
|
// MAIN //
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* Returns a pseudorandom number drawn from a Poisson distribution with parameter `lambda`.
|
|||
|
|
*
|
|||
|
|
* ## References
|
|||
|
|
*
|
|||
|
|
* - Hörmann, W. 1993. "The transformed rejection method for generating Poisson random variables." _Insurance: Mathematics and Economics_ 12 (1): 39–45. doi:[10.1016/0167-6687(93)90997-4][@hormann:1993b].
|
|||
|
|
*
|
|||
|
|
* [@hormann:1993b]: http://dx.doi.org/10.1016/0167-6687(93)90997-4
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* @private
|
|||
|
|
* @param {PRNG} rand - PRNG for generating uniformly distributed numbers
|
|||
|
|
* @param {PositiveNumber} lambda - mean
|
|||
|
|
* @returns {NonNegativeInteger} pseudorandom number
|
|||
|
|
*/
|
|||
|
|
function poisson( rand, lambda ) {
|
|||
|
|
var slambda;
|
|||
|
|
var ainv;
|
|||
|
|
var urvr;
|
|||
|
|
var us;
|
|||
|
|
var vr;
|
|||
|
|
var a;
|
|||
|
|
var b;
|
|||
|
|
var k;
|
|||
|
|
var u;
|
|||
|
|
var v;
|
|||
|
|
|
|||
|
|
slambda = sqrt( lambda );
|
|||
|
|
|
|||
|
|
b = (2.53*slambda) + 0.931;
|
|||
|
|
a = (0.02483*b) - 0.059;
|
|||
|
|
|
|||
|
|
ainv = (1.1328/(b-3.4)) + 1.1239;
|
|||
|
|
vr = (-3.6224/(b-2.0)) + 0.9277;
|
|||
|
|
urvr = 0.86 * vr;
|
|||
|
|
|
|||
|
|
while ( true ) {
|
|||
|
|
v = rand();
|
|||
|
|
if ( v <= urvr ) {
|
|||
|
|
u = (v / vr) - 0.43;
|
|||
|
|
u *= (2.0*a / (0.5-abs(u))) + b;
|
|||
|
|
u += lambda + 0.445;
|
|||
|
|
return floor( u );
|
|||
|
|
}
|
|||
|
|
if ( v >= vr ) {
|
|||
|
|
u = rand() - 0.5;
|
|||
|
|
} else {
|
|||
|
|
u = (v / vr) - 0.93;
|
|||
|
|
u = (sign( u )*0.5) - u;
|
|||
|
|
v = vr * rand();
|
|||
|
|
}
|
|||
|
|
us = 0.5 - abs( u );
|
|||
|
|
if (
|
|||
|
|
us >= 0.013 ||
|
|||
|
|
us >= v
|
|||
|
|
) {
|
|||
|
|
k = floor( (((2.0*a/us) + b)*u) + lambda + 0.445 );
|
|||
|
|
v *= ainv / ( (a/(us*us)) + b );
|
|||
|
|
u = (k+0.5) * ln( lambda/k );
|
|||
|
|
u += -lambda - LN_SQRT_TWO_PI + k;
|
|||
|
|
u -= ( ONE_12 - (ONE_360/(k*k)) ) / k;
|
|||
|
|
if (
|
|||
|
|
k >= 10 &&
|
|||
|
|
u >= ln( v*slambda )
|
|||
|
|
) {
|
|||
|
|
return k;
|
|||
|
|
}
|
|||
|
|
u = (k*ln( lambda )) - lambda - factorialln( k );
|
|||
|
|
if (
|
|||
|
|
k >= 0 &&
|
|||
|
|
k <= 9 &&
|
|||
|
|
u >= ln( v )
|
|||
|
|
) {
|
|||
|
|
return k;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
// EXPORTS //
|
|||
|
|
|
|||
|
|
module.exports = poisson;
|