geometric.js (2153B)
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 var ln = require( '@stdlib/math/base/special/ln' ); 25 26 27 // MAIN // 28 29 /** 30 * Returns a pseudorandom number drawn from a geometric distribution. 31 * 32 * ## Proof 33 * 34 * Consider 35 * 36 * ```tex 37 * N = \left \lfloor \ln (U) / \ln (1-p) \right \rfloor 38 * ``` 39 * 40 * where \\( U \\) is uniform on the interval \\((0,1)\\). Accordingly, \\(N\\) must be a nonnegative integer, and, for every \\( n \geq 0\\), the event \\(A_n = \left \{ N = n \right \}\\) is 41 * 42 * ```tex 43 * A_n = \left \{(n+1) \ln (1-p) < \ln (U) \leq n \ln (1-p) \right \} 44 * ``` 45 * 46 * where \\(\ln (1-p) < 0\\). Thus, 47 * 48 * ```tex 49 * A_n = \left \{(1-p)^{n+1} < U \leq (1-p)^n \right \} 50 * ``` 51 * 52 * For every \\(u < v\\) on the interval \\((0,1)\\), 53 * 54 * ```tex 55 * P\left \[u < U \leq v\right \] = v - u 56 * ``` 57 * 58 * Hence, 59 * 60 * ```tex 61 * P\left \[N = n \right \] = P\left \[A_n\right \] = (1-p)^n - (1-p)^{n+1} = (1-p)^n(1-(1-p)) = p(1-p)^n 62 * ``` 63 * 64 * which proves that \\(N\\) is a geometric random variable. 65 * 66 * 67 * @private 68 * @param {PRNG} rand - PRNG for uniformly distributed numbers 69 * @param {Probability} p - success probability 70 * @returns {NonNegativeInteger} pseudorandom number 71 */ 72 function geometric( rand, p ) { 73 var u = rand(); 74 if ( u === 0.0 ) { 75 // Drawing random variates from a PRNG (with period > 1) is effectively sampling without replacement. Thus, should not be possible to draw `0` twice in a row. 76 u = rand(); 77 } 78 return floor( ln( u ) / ln( 1.0-p ) ); 79 } 80 81 82 // EXPORTS // 83 84 module.exports = geometric;