time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

rejection.js (2853B)

      1 /**
      2 * @license Apache-2.0
      3 *
      4 * Copyright (c) 2018 The Stdlib Authors.
      5 *
      6 * Licensed under the Apache License, Version 2.0 (the "License");
      7 * you may not use this file except in compliance with the License.
      8 * You may obtain a copy of the License at
      9 *
     10 *    http://www.apache.org/licenses/LICENSE-2.0
     11 *
     12 * Unless required by applicable law or agreed to in writing, software
     13 * distributed under the License is distributed on an "AS IS" BASIS,
     14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     15 * See the License for the specific language governing permissions and
     16 * limitations under the License.
     17 */
     18 
     19 'use strict';
     20 
     21 // MODULES //
     22 
     23 var factorialln = require( '@stdlib/math/base/special/factorialln' );
     24 var floor = require( '@stdlib/math/base/special/floor' );
     25 var sign = require( '@stdlib/math/base/special/signum' );
     26 var sqrt = require( '@stdlib/math/base/special/sqrt' );
     27 var abs = require( '@stdlib/math/base/special/abs' );
     28 var ln = require( '@stdlib/math/base/special/ln' );
     29 var LN_SQRT_TWO_PI = require( '@stdlib/constants/float64/ln-sqrt-two-pi' );
     30 
     31 
     32 // VARIABLES //
     33 
     34 var ONE_12 = 1.0 / 12.0;
     35 var ONE_360 = 1.0 / 360.0;
     36 
     37 
     38 // MAIN //
     39 
     40 /**
     41 * Returns a pseudorandom number drawn from a Poisson distribution with parameter `lambda`.
     42 *
     43 * ## References
     44 *
     45 * -   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].
     46 *
     47 * [@hormann:1993b]: http://dx.doi.org/10.1016/0167-6687(93)90997-4
     48 *
     49 *
     50 * @private
     51 * @param {PRNG} rand - PRNG for generating uniformly distributed numbers
     52 * @param {PositiveNumber} lambda - mean
     53 * @returns {NonNegativeInteger} pseudorandom number
     54 */
     55 function poisson( rand, lambda ) {
     56 	var slambda;
     57 	var ainv;
     58 	var urvr;
     59 	var us;
     60 	var vr;
     61 	var a;
     62 	var b;
     63 	var k;
     64 	var u;
     65 	var v;
     66 
     67 	slambda = sqrt( lambda );
     68 
     69 	b = (2.53*slambda) + 0.931;
     70 	a = (0.02483*b) - 0.059;
     71 
     72 	ainv = (1.1328/(b-3.4)) + 1.1239;
     73 	vr = (-3.6224/(b-2.0)) + 0.9277;
     74 	urvr = 0.86 * vr;
     75 
     76 	while ( true ) {
     77 		v = rand();
     78 		if ( v <= urvr ) {
     79 			u = (v / vr) - 0.43;
     80 			u *= (2.0*a / (0.5-abs(u))) + b;
     81 			u += lambda + 0.445;
     82 			return floor( u );
     83 		}
     84 		if ( v >= vr ) {
     85 			u = rand() - 0.5;
     86 		} else {
     87 			u = (v / vr) - 0.93;
     88 			u = (sign( u )*0.5) - u;
     89 			v = vr * rand();
     90 		}
     91 		us = 0.5 - abs( u );
     92 		if (
     93 			us >= 0.013 ||
     94 			us >= v
     95 		) {
     96 			k = floor( (((2.0*a/us) + b)*u) + lambda + 0.445 );
     97 			v *= ainv / ( (a/(us*us)) + b );
     98 			u = (k+0.5) * ln( lambda/k );
     99 			u += -lambda - LN_SQRT_TWO_PI + k;
    100 			u -= ( ONE_12 - (ONE_360/(k*k)) ) / k;
    101 			if (
    102 				k >= 10 &&
    103 				u >= ln( v*slambda )
    104 			) {
    105 				return k;
    106 			}
    107 			u = (k*ln( lambda )) - lambda - factorialln( k );
    108 			if (
    109 				k >= 0 &&
    110 				k <= 9 &&
    111 				u >= ln( v )
    112 			) {
    113 				return k;
    114 			}
    115 		}
    116 	}
    117 }
    118 
    119 
    120 // EXPORTS //
    121 
    122 module.exports = poisson;