vose.js (2502B)
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 floor = require( '@stdlib/math/base/special/floor' ); 24 25 26 // MAIN // 27 28 /** 29 * Samples with replacement and non-uniform probabilities using Vose's [alias method][alias-method]. 30 * 31 * ## References 32 * 33 * - Vose, Michael D. 1991. "A linear algorithm for generating random numbers with a given distribution." _IEEE Transactions on Software Engineering_ 17 (9): 972–75. doi:[10.1109/32.92917][@vose:1991]. 34 * 35 * [alias-method]: http://keithschwarz.com/darts-dice-coins/ 36 * [@vose:1991]: https://doi.org/10.1109/32.92917 37 * 38 * 39 * @private 40 * @param {ArrayLike} x - array-like object from which to sample 41 * @param {NonNegativeInteger} size - sample size 42 * @param {Function} rand - PRNG for generating uniformly distributed numbers 43 * @param {ProbabilityArray} probabilities - element probabilities 44 * @returns {Array} sample 45 */ 46 function vose( x, size, rand, probabilities ) { 47 var small; 48 var large; 49 var probs; 50 var alias; 51 var out; 52 var N; 53 var p; 54 var g; 55 var i; 56 var l; 57 58 probs = probabilities.slice(); 59 N = x.length; 60 61 small = []; 62 large = []; 63 for ( i = 0; i < N; i++ ) { 64 probs[ i ] *= N; 65 if ( probs[ i ] < 1.0 ) { 66 small.push( i ); 67 } else { 68 large.push( i ); 69 } 70 } 71 alias = new Array( N ); 72 p = new Array( N ); 73 while ( small.length !== 0 && large.length !== 0 ) { 74 l = small.shift(); 75 g = large.shift(); 76 p[ l ] = probs[ l ]; 77 alias[ l ] = g; 78 probs[ g ] = probs[ g ] + probs[ l ] - 1.0; 79 if ( probs[ g ] < 1.0 ) { 80 small.push( g ); 81 } else { 82 large.push( g ); 83 } 84 } 85 for ( i = 0; i < large.length; i++ ) { 86 p[ large[ i ] ] = 1.0; 87 } 88 for ( i = 0; i < small.length; i++ ) { 89 p[ small[ i ] ] = 1.0; 90 } 91 out = new Array( size ); 92 for ( i = 0; i < size; i++ ) { 93 l = floor( N*rand() ); 94 if ( rand() < p[ l ] ) { 95 out[ i ] = x[ l ]; 96 } else { 97 out[ i ] = x[ alias[ l ] ]; 98 } 99 } 100 return out; 101 } 102 103 104 // EXPORTS // 105 106 module.exports = vose;