1 | /* ================================================================
|
2 | * probability-distributions by Matt Asher (me[at]mattasher.com)
|
3 | * Originally created for StatisticsBlog.com
|
4 | *
|
5 | * first created at : Sat Oct 10 2015
|
6 | *
|
7 | * ================================================================
|
8 | * Copyright 2015 Matt Asher
|
9 | *
|
10 | * Licensed under the MIT License
|
11 | * You may not use this file except in compliance with the License.
|
12 | *
|
13 | * ================================================================ */
|
14 |
|
15 | var crypto = require('crypto');
|
16 |
|
17 | // Shortcuts
|
18 | var exp = Math.exp;
|
19 | var ln = Math.log;
|
20 | var PI = Math.PI;
|
21 | var pow = Math.pow;
|
22 |
|
23 | module.exports = {
|
24 |
|
25 | /**
|
26 | * This is the core function for generating entropy
|
27 | *
|
28 | * @param len number of bytes of entropy to create
|
29 | * @returns {number} A pseduo random number between 0 and 1
|
30 | *
|
31 | */
|
32 | prng: function(len) {
|
33 | if(len === undefined) len=16;
|
34 |
|
35 | var entropy = crypto.randomBytes(len);
|
36 | var result = 0;
|
37 |
|
38 | for(var i=0; i<len; i++) {
|
39 | result = result + Number(entropy[i])/Math.pow(256,(i+1))
|
40 | }
|
41 | return result
|
42 | },
|
43 |
|
44 |
|
45 |
|
46 |
|
47 | /**
|
48 | *
|
49 | * @param n The number of random variates to create. Must be a positive integer.
|
50 | * @param alpha First shape parameter
|
51 | * @param beta Second shape parameter
|
52 | * @param loc Location or Non-centrality parameter
|
53 | */
|
54 | rbeta: function(n, alpha, beta, loc) {
|
55 | // Uses relationship with gamma to calculate
|
56 |
|
57 | // Validations
|
58 | n = this._v(n, "n");
|
59 | alpha = this._v(alpha, "nn", 1);
|
60 | beta = this._v(beta, "nn", 1);
|
61 | loc = this._v(loc, "r", 0);
|
62 |
|
63 | var toReturn = [];
|
64 |
|
65 | for(var i=0; i<n; i++) {
|
66 | var g1 = this.rgamma(1, alpha, 1)[0];
|
67 | var g2 = this.rgamma(1, beta, 1)[0];
|
68 |
|
69 |
|
70 | toReturn[i] = loc + g1/(g1+g2);
|
71 | }
|
72 | return toReturn
|
73 |
|
74 | },
|
75 |
|
76 |
|
77 | /**
|
78 | *
|
79 | * @param n Number of variates to return.
|
80 | * @param size Number of Bernoulli trials to be summed up. Defaults to 1
|
81 | * @param p Probability of a "success". Defaults to 0.5
|
82 | * @returns {Array} Random variates array
|
83 | */
|
84 | rbinom: function(n, size, p) {
|
85 | n = this._v(n, "n");
|
86 | size = this._v(size, "nni", 1);
|
87 | p = this._v(p, "p", 0.5);
|
88 |
|
89 | var toReturn = [];
|
90 |
|
91 | for(var i=0; i<n; i++) {
|
92 | var result = 0;
|
93 | for(var j=0; j<size; j++) {
|
94 | if(this.prng() < p) {
|
95 | result++
|
96 | }
|
97 | }
|
98 | toReturn[i] = result;
|
99 | }
|
100 | return toReturn
|
101 | },
|
102 |
|
103 |
|
104 | /**
|
105 | *
|
106 | * @param n The number of variates to create
|
107 | * @param loc Location parameter
|
108 | * @param scale Scale parameter
|
109 | * @returns {Array} Random variates array
|
110 | */
|
111 | rcauchy: function(n, loc, scale) {
|
112 | n = this._v(n, "n");
|
113 | loc = this._v(loc, "r", 0);
|
114 | scale = this._v(scale, "nn", 1);
|
115 |
|
116 | var toReturn = [];
|
117 | for(var i=0; i<n; i++) {
|
118 | var x = scale * Math.tan(PI * (this.prng()-0.5))+loc;
|
119 |
|
120 | toReturn[i] = x;
|
121 | }
|
122 |
|
123 | return toReturn
|
124 | },
|
125 |
|
126 | /**
|
127 | *
|
128 | * @param n The number of variates to create
|
129 | * @param df Degrees of freedom for the distribution
|
130 | * @param ncp Non-centrality parameter
|
131 | * @returns {Array} Random variates array
|
132 | */
|
133 | rchisq: function(n, df, ncp) {
|
134 | n = this._v(n, "n");
|
135 | df = this._v(df, "nn");
|
136 | ncp = this._v(ncp, "r", 0);
|
137 |
|
138 | var toReturn = [];
|
139 | for(var i=0; i<n; i++) {
|
140 | // Start at ncp
|
141 | var x = ncp;
|
142 | for(var j=0; j<df; j++) {
|
143 | x = x + Math.pow(this.rnorm(1)[0],2);
|
144 | }
|
145 | toReturn[i] = x;
|
146 | }
|
147 | return toReturn
|
148 | },
|
149 |
|
150 | /**
|
151 | *
|
152 | * @param x Where to sample the density
|
153 | * @param rate The rate parameter. Must be a positive number
|
154 | * @returns {Number} The density given the parameter values
|
155 | */
|
156 | dexp: function(x, rate) {
|
157 | x = this._v(x, "r");
|
158 | rate = this._v(rate, "pos", 1);
|
159 | if(x < 0) return 0
|
160 |
|
161 | return rate * exp(-rate * x)
|
162 | },
|
163 |
|
164 | /**
|
165 | *
|
166 | * @param n The number of random variates to create. Must be a positive integer.
|
167 | * @param rate The rate parameter. Must be a positive number
|
168 | */
|
169 | rexp: function(n, rate) {
|
170 | n = this._v(n, "n");
|
171 | rate = this._v(rate, "pos", 1);
|
172 |
|
173 | var toReturn = [];
|
174 |
|
175 | for(var i=0; i<n; i++) {
|
176 |
|
177 | toReturn[i] = -ln(this.prng())/rate;
|
178 | }
|
179 |
|
180 | return toReturn
|
181 | },
|
182 |
|
183 | /**
|
184 | *
|
185 | * @param n The number of variates to create
|
186 | * @param df1 Degrees of freedom for the first parameter
|
187 | * @param df2 Degrees of freedom for the first parameter
|
188 | * @returns {Array} Random variates array
|
189 | */
|
190 | rf: function(n, df1, df2) {
|
191 | n = this._v(n, "n");
|
192 | df1 = this._v(df1, "nn");
|
193 | df2 = this._v(df2, "nn");
|
194 |
|
195 | var toReturn = [];
|
196 | for(var i=0; i<n; i++) {
|
197 | var num = this.rchisq(1, df1)[0]/df1;
|
198 | var denom = this.rchisq(1, df2)[0]/df2;
|
199 |
|
200 | toReturn[i] = num/denom;
|
201 | }
|
202 |
|
203 | return toReturn
|
204 |
|
205 | },
|
206 |
|
207 |
|
208 | /**
|
209 | *
|
210 | * @param n The number of random variates to create. Must be a positive integer
|
211 | * @param alpha
|
212 | * @param rate
|
213 | * @returns {Array} Random variates array
|
214 | */
|
215 | rgamma: function(n, alpha, rate) {
|
216 | // Adapted from https://github.com/mvarshney/simjs-source/ & scipy
|
217 | n = this._v(n, "n");
|
218 | alpha = this._v(alpha, "nn");
|
219 | rate = this._v(rate, "pos", 1);
|
220 |
|
221 | var LOG4 = ln(4.0);
|
222 | var SG_MAGICCONST = 1.0 + ln(4.5);
|
223 | var beta = 1/rate;
|
224 |
|
225 | var toReturn = [];
|
226 | for(var i = 0; i<n; i++) {
|
227 |
|
228 | /* Based on Python 2.6 source code of random.py.
|
229 | */
|
230 |
|
231 | if (alpha > 1.0) {
|
232 | var ainv = Math.sqrt(2.0 * alpha - 1.0);
|
233 | var bbb = alpha - LOG4;
|
234 | var ccc = alpha + ainv;
|
235 |
|
236 | while (true) {
|
237 | var u1 = this.prng();
|
238 | if ((u1 < 1e-7) || (u > 0.9999999)) {
|
239 | continue;
|
240 | }
|
241 | var u2 = 1.0 - this.prng();
|
242 | var v = ln(u1 / (1.0 - u1)) / ainv;
|
243 | var x = alpha * exp(v);
|
244 | var z = u1 * u1 * u2;
|
245 | var r = bbb + ccc * v - x;
|
246 | if ((r + SG_MAGICCONST - 4.5 * z >= 0.0) || (r >= ln(z))) {
|
247 | var result = x * beta;
|
248 | break;
|
249 | }
|
250 | }
|
251 | } else if (alpha == 1.0) {
|
252 | var u = this.prng();
|
253 | while (u <= 1e-7) {
|
254 | u = this.prng();
|
255 | }
|
256 | var result = - ln(u) * beta;
|
257 | } else {
|
258 | while (true) {
|
259 | var u = this.prng();
|
260 | var b = (Math.E + alpha) / Math.E;
|
261 | var p = b * u;
|
262 | if (p <= 1.0) {
|
263 | var x = Math.pow(p, 1.0 / alpha);
|
264 | } else {
|
265 | var x = - ln((b - p) / alpha);
|
266 | }
|
267 | var u1 = this.prng();
|
268 | if (p > 1.0) {
|
269 | if (u1 <= Math.pow(x, (alpha - 1.0))) {
|
270 | break;
|
271 | }
|
272 | } else if (u1 <= exp(-x)) {
|
273 | break;
|
274 | }
|
275 | }
|
276 | var result = x * beta;
|
277 | }
|
278 |
|
279 | toReturn[i] = result;
|
280 | }
|
281 |
|
282 | return toReturn;
|
283 |
|
284 | },
|
285 |
|
286 |
|
287 | /**
|
288 | *
|
289 | * @param n The number of random variates to create. Must be a positive integer
|
290 | * @param min Minimum value
|
291 | * @param max Maximum value
|
292 | * @param inclusive By default the minimum and maximum are inclusive. To make exclusive, set to false
|
293 | * @returns {Array}
|
294 | */
|
295 | rint: function(n, min, max, inclusive) {
|
296 | n = this._v(n, "n");
|
297 | min = this._v(min, "int");
|
298 | max = this._v(max, "int");
|
299 | if(inclusive === false) {
|
300 | min++;
|
301 | max--;
|
302 | }
|
303 |
|
304 | if(min > max) throw "Minimum value cannot be greater than maximum value. For non-inclusive, minimum and maximum must be separated by at least 2.";
|
305 |
|
306 | var toReturn = [];
|
307 |
|
308 | var raw = this.runif(n, min, max);
|
309 |
|
310 | for(var i=0; i<n; i++) {
|
311 | toReturn[i] = Math.round(raw[i]);
|
312 | }
|
313 |
|
314 | return toReturn
|
315 | },
|
316 |
|
317 | // Syntax as in R library VGAM
|
318 | /**
|
319 | *
|
320 | * @param n The number of random variates to create. Must be a positive integer
|
321 | * @param loc Mean
|
322 | * @param scale Scale parameter
|
323 | * @returns {Array} Random variates array
|
324 | */
|
325 | rlaplace: function(n, loc, scale) {
|
326 | n = this._v(n, "n");
|
327 | loc = this._v(loc, "r", 0);
|
328 | scale = this._v(scale, "nn", 1);
|
329 |
|
330 | var toReturn = [];
|
331 |
|
332 | for(var i=0; i<n; i++) {
|
333 | var core = this.sample([-1,1])[0] * ln(this.prng());
|
334 |
|
335 | var x = loc - scale * core;
|
336 |
|
337 | toReturn[i] = x;
|
338 | }
|
339 |
|
340 | return toReturn
|
341 | },
|
342 |
|
343 |
|
344 | /**
|
345 | *
|
346 | * @param n The number of random variates to create. Must be a positive integer.
|
347 | * @param meanlog The mean log.
|
348 | * @param sdlog Log SD. Must be greater than 0.
|
349 | * @returns {Array} Random variates array
|
350 | */
|
351 | rlnorm: function(n, meanlog, sdlog) {
|
352 | n = this._v(n, "n");
|
353 | meanlog = this._v(meanlog, "r", 0);
|
354 | sdlog = this._v(sdlog, "nn", 1);
|
355 |
|
356 | var toReturn = [];
|
357 |
|
358 | for(var i=0; i<n; i++) {
|
359 | var x = this.rnorm(1, meanlog, sdlog)[0];
|
360 |
|
361 | toReturn[i] = exp(x);
|
362 | }
|
363 |
|
364 | return toReturn
|
365 | },
|
366 |
|
367 | /**
|
368 | *
|
369 | * @param n The number of random variates to create. Must be a positive integer.
|
370 | * @param size Number of hits required
|
371 | * @param p Hit probability
|
372 | * @param mu Optional way to specify hit probability
|
373 | * @returns {Array} Random variates array
|
374 | */
|
375 | rnbinom: function(n, size, p, mu) {
|
376 | n = this._v(n, "n");
|
377 | if(size === undefined) size=1;
|
378 | if(Math.round(size) != size) throw "Size must be a whole number";
|
379 | if(size < 1) throw "Size must be one or greater";
|
380 | if(p !== undefined && mu !== undefined) throw "You must specify probability or mean, not both";
|
381 | if(mu !== undefined) p = size/(size+mu);
|
382 | p = this._v(p, "p");
|
383 |
|
384 |
|
385 | var toReturn = [];
|
386 |
|
387 | for(var i=0; i<n; i++) {
|
388 |
|
389 | // Core distribution
|
390 | var result = 0;
|
391 | var leftToFind = size;
|
392 | while(leftToFind > 0) {
|
393 | result++;
|
394 | if(this.prng() < p) leftToFind--;
|
395 | }
|
396 |
|
397 | toReturn[i] = result - 1;
|
398 | }
|
399 |
|
400 | return toReturn
|
401 |
|
402 | },
|
403 |
|
404 | /**
|
405 | *
|
406 | * @param x Where to sample the density
|
407 | * @param mean Mean of the distribution
|
408 | * @param sd Standard deviation for the distribution
|
409 | * @returns {Number} The density given the parameter values
|
410 | */
|
411 | dnorm: function(x, mean, sd) {
|
412 | x = this._v(x, "r");
|
413 | mean = this._v(mean, "r", 0);
|
414 | sd = this._v(sd, "nn", 1);
|
415 |
|
416 | // Check for degeneracy
|
417 | if(sd === 0) {
|
418 | if(x === mean) return Infinity;
|
419 | return 0
|
420 | }
|
421 |
|
422 | var a = sd*(Math.sqrt(2*PI));
|
423 | var b = -(x-mean)*(x-mean);
|
424 | var c = 2*sd*sd;
|
425 |
|
426 | return (1/a)*exp(b/c)
|
427 | },
|
428 |
|
429 |
|
430 | /**
|
431 | *
|
432 | * @param n The number of random variates to create. Must be a positive integer.
|
433 | * @param mean Mean of the distribution
|
434 | * @param sd Standard Deviation of the distribution
|
435 | * @returns {Array} Random variates array
|
436 | */
|
437 | rnorm: function(n, mean, sd) {
|
438 | // Adapted from http://blog.yjl.im/2010/09/simulating-normal-random-variable-using.html
|
439 |
|
440 | n = this._v(n, "n");
|
441 | mean = this._v(mean, "r", 0);
|
442 | sd = this._v(sd, "nn", 1);
|
443 |
|
444 | var toReturn = [];
|
445 |
|
446 | for(var i=0; i<n; i++) {
|
447 | var V1, V2, S, X;
|
448 |
|
449 | do {
|
450 | var U1 = this.prng();
|
451 | var U2 = this.prng();
|
452 | V1 = (2 * U1) - 1;
|
453 | V2 = (2 * U2) - 1;
|
454 | S = (V1 * V1) + (V2 * V2);
|
455 | } while (S > 1);
|
456 |
|
457 | X = Math.sqrt(-2 * ln(S) / S) * V1;
|
458 | X = mean + sd * X;
|
459 | toReturn.push(X);
|
460 | }
|
461 |
|
462 | return toReturn
|
463 | },
|
464 |
|
465 |
|
466 | /**
|
467 | *
|
468 | * @param x Where to sample the density
|
469 | * @param lambda Mean/variance
|
470 | * @returns {Number} The density given the parameter values
|
471 | */
|
472 | dpois: function(x, lambda) {
|
473 | x = this._v(x, "nni");
|
474 | lambda = this._v(lambda, "nn");
|
475 |
|
476 | // Check for degeneracy
|
477 | if(lambda === 0) {
|
478 | if(x === 0) return 1;
|
479 | return 0
|
480 | }
|
481 |
|
482 | var a = pow(lambda, x);
|
483 | var b = exp(-lambda);
|
484 | var c = this._factorial(x);
|
485 |
|
486 | return a*b/c
|
487 | },
|
488 |
|
489 |
|
490 | /**
|
491 | *
|
492 | * @param n The number of random variates to create. Must be a positive integer.
|
493 | * @param lambda Mean/Variance of the distribution
|
494 | * @returns {Array} Random variates array
|
495 | */
|
496 | rpois: function(n, lambda) {
|
497 | n = this._v(n, "n");
|
498 | lambda = this._v(lambda, "pos");
|
499 |
|
500 | var toReturn = [];
|
501 |
|
502 | for(var i=0; i<n; i++) {
|
503 |
|
504 | // Adapted from http://wiki.q-researchsoftware.com/wiki/How_to_Generate_Random_Numbers:_Poisson_Distribution
|
505 | if (lambda < 30) {
|
506 |
|
507 | var L = exp(-lambda);
|
508 | var p = 1;
|
509 | var k = 0;
|
510 | do {
|
511 | k++;
|
512 | p *= this.prng();
|
513 | } while (p > L);
|
514 | toReturn.push(k - 1);
|
515 |
|
516 | } else {
|
517 |
|
518 | // Roll our own
|
519 | // Fix total number of samples
|
520 | var samples = 10000;
|
521 | var p = lambda/samples;
|
522 | var k = 0;
|
523 | for(var j=0; j<samples; j++) {
|
524 | if(this.prng() < p) {
|
525 | k++
|
526 | }
|
527 | }
|
528 | toReturn[i] = k;
|
529 | }
|
530 | }
|
531 |
|
532 | return toReturn
|
533 | },
|
534 |
|
535 | dunif: function(x, min, max) {
|
536 | x = this._v(x, "r");
|
537 | min = this._v(min, "r", 0);
|
538 | max = this._v(max, "r", 1);
|
539 | if(min > max) throw "Minimum value cannot be greater than maximum value";
|
540 |
|
541 | if(x < min || x > max) return 0;
|
542 | if(min === max) return Infinity;
|
543 |
|
544 |
|
545 | return 1/(max-min);
|
546 | },
|
547 |
|
548 |
|
549 | /**
|
550 | *
|
551 | * @param n Number of variates to return
|
552 | * @param min Lower bound
|
553 | * @param max Upper bound
|
554 | * @returns {Array} Random variates array
|
555 | */
|
556 | runif: function(n, min, max) {
|
557 | n = this._v(n, "n");
|
558 | min = this._v(min, "r", 0);
|
559 | max = this._v(max, "r", 1);
|
560 | if(min > max) throw "Minimum value cannot be greater than maximum value";
|
561 |
|
562 | var toReturn = [];
|
563 |
|
564 | for(var i=0; i<n; i++) {
|
565 | var raw = this.prng();
|
566 | var scaled = min + raw*(max-min);
|
567 | toReturn.push(scaled)
|
568 | }
|
569 | return toReturn
|
570 | },
|
571 |
|
572 |
|
573 | /**
|
574 | *
|
575 | * @param collection Array of items to sample from
|
576 | * @param n Number of items to sample. If missing, n will be set to the length of the collection and it will shuffle
|
577 | * @param replace Sample with replacement? False by default
|
578 | * @param ratios Ratios to weight items. Can be any non-negative number. By default all items are given equal weight
|
579 | * @returns {Array} Array of sampled items
|
580 | */
|
581 | sample: function(collection, n, replace, ratios) {
|
582 |
|
583 | // Validations
|
584 | collection = this._v(collection, "a");
|
585 | n = this._v(n, "n", collection.length); // If n is undefined, sample the full array
|
586 | if(replace === undefined) replace = false;
|
587 | if(!replace && collection.length < n)
|
588 | throw "You cannot select " + n + " items from an array of length " + collection.length + " without replacement";
|
589 |
|
590 | if(ratios === undefined) {
|
591 | ratios = [];
|
592 | for(var m=0; m<collection.length; m++) { ratios[m] = 1 }
|
593 | }
|
594 |
|
595 | var cumulativeProbs = this._getCumulativeProbs(ratios, collection.length);
|
596 |
|
597 | // Main loop
|
598 | var toReturn = [];
|
599 |
|
600 | for(var i=0; i<n; i++) {
|
601 |
|
602 | var chosen = this._sampleOneIndex(cumulativeProbs);
|
603 |
|
604 | if(replace) {
|
605 | toReturn[i] = collection[chosen];
|
606 | } else {
|
607 |
|
608 | // Remove from collection and ratios
|
609 | toReturn[i] = collection.splice(chosen, 1)[0];
|
610 | ratios.splice(chosen, 1);
|
611 |
|
612 | // Make sure we aren't at the end
|
613 | if(ratios.length) {
|
614 | cumulativeProbs = this._getCumulativeProbs(ratios, collection.length);
|
615 | }
|
616 | }
|
617 | }
|
618 |
|
619 | return toReturn;
|
620 |
|
621 | },
|
622 |
|
623 |
|
624 | // HELPERS
|
625 |
|
626 | /**
|
627 | *
|
628 | * @param ratios Array of non-negative numbers to be turned into CDF
|
629 | * @param len length of the collection
|
630 | * @returns {Array}
|
631 | * @private
|
632 | */
|
633 | _getCumulativeProbs: function(ratios, len) {
|
634 | if(len === undefined) throw "An error occurred: len was not sent to _getCumulativeProbs";
|
635 | if(ratios.length !== len) throw "Probabilities for sample must be same length as the array to sample from";
|
636 |
|
637 | var toReturn = [];
|
638 |
|
639 | if(ratios !== undefined) {
|
640 | ratios = this._v(ratios, "a");
|
641 | if(ratios.length !== len) throw "Probabilities array must be the same length as the array you are sampling from";
|
642 |
|
643 | var sum = 0;
|
644 | ratios.map(function(ratio) {
|
645 | ratio = this._v(ratio, "nn"); // Note validating as ANY non-negative number
|
646 | sum+= ratio;
|
647 | toReturn.push(sum);
|
648 | }.bind(this));
|
649 |
|
650 | // Divide by total to normalize
|
651 | for(var k=0; k<toReturn.length; k++) { toReturn[k] = toReturn[k]/sum }
|
652 | return toReturn
|
653 | }
|
654 | },
|
655 |
|
656 | _sampleOneIndex: function(cumulativeProbs) {
|
657 |
|
658 | var toTake = this.prng();
|
659 |
|
660 | // Find out where this lands in weights
|
661 | var cur = 0;
|
662 | while(toTake > cumulativeProbs[cur]) cur++;
|
663 |
|
664 | return cur;
|
665 | },
|
666 |
|
667 | _factorial: function(n) {
|
668 | var toReturn=1;
|
669 | for (var i = 2; i <= n; i++)
|
670 | toReturn = toReturn * i;
|
671 |
|
672 | return toReturn;
|
673 | },
|
674 |
|
675 | // Return default if undefined, otherwise validate
|
676 | _v: function(param, type, defaultParam) {
|
677 | if(param == null && defaultParam != null)
|
678 | return defaultParam;
|
679 |
|
680 | switch(type) {
|
681 |
|
682 | // Array of 1 item or more
|
683 | case "a":
|
684 | if(!Array.isArray(param) || !param.length) throw "Expected an array of length 1 or greater";
|
685 | return param;
|
686 |
|
687 | // Integer
|
688 | case "int":
|
689 | if(param !== Number(param)) throw "A required parameter is missing or not a number";
|
690 | if(param !== Math.round(param)) throw "Parameter must be a whole number";
|
691 | if(param === Infinity) throw 'Sent "infinity" as a parameter';
|
692 | return param;
|
693 |
|
694 | // Natural number
|
695 | case "n":
|
696 | if(param === undefined) throw "You must specify how many values you want";
|
697 | if(param !== Number(param)) throw "The number of values must be numeric";
|
698 | if(param !== Math.round(param)) throw "The number of values must be a whole number";
|
699 | if(param < 1) throw "The number of values must be a whole number of 1 or greater";
|
700 | if(param === Infinity) throw "The number of values cannot be infinite ;-)";
|
701 | return param;
|
702 |
|
703 | // Valid probability
|
704 | case "p":
|
705 | if(Number(param) !== param) throw "Probability value is missing or not a number";
|
706 | if(param > 1) throw "Probability values cannot be greater than 1";
|
707 | if(param < 0) throw "Probability values cannot be less than 0";
|
708 | return param;
|
709 |
|
710 | // Positive numbers
|
711 | case "pos":
|
712 | if(Number(param) !== param) throw "A required parameter is missing or not a number";
|
713 | if(param <= 0) throw "Parameter must be greater than 0";
|
714 | if(param === Infinity) throw 'Sent "infinity" as a parameter';
|
715 | return param;
|
716 |
|
717 | // Look for numbers (reals)
|
718 | case "r":
|
719 | if(Number(param) !== param) throw "A required parameter is missing or not a number";
|
720 | if(param === Infinity) throw 'Sent "infinity" as a parameter';
|
721 | return param;
|
722 |
|
723 | // Non negative real number
|
724 | case "nn":
|
725 | if(param !== Number(param)) throw "A required parameter is missing or not a number";
|
726 | if(param < 0) throw "Parameter cannot be less than 0";
|
727 | if(param === Infinity) throw 'Sent "infinity" as a parameter';
|
728 | return param;
|
729 |
|
730 | // Non negative whole number (integer)
|
731 | case "nni":
|
732 | if(param !== Number(param)) throw "A required parameter is missing or not a number";
|
733 | if(param !== Math.round(param)) throw "Parameter must be a whole number";
|
734 | if(param < 0) throw "Parameter cannot be less than zero";
|
735 | if(param === Infinity) throw 'Sent "infinity" as a parameter';
|
736 | return param;
|
737 |
|
738 | }
|
739 | },
|
740 |
|
741 |
|
742 | // ________ _______ ______ _____ _____ __ __ ______ _ _ _______ _
|
743 | // | ____\ \ / / __ \| ____| __ \|_ _| \/ | ____| \ | |__ __|/\ | |
|
744 | // | |__ \ V /| |__) | |__ | |__) | | | | \ / | |__ | \| | | | / \ | |
|
745 | // | __| > < | ___/| __| | _ / | | | |\/| | __| | . ` | | | / /\ \ | |
|
746 | // | |____ / . \| | | |____| | \ \ _| |_| | | | |____| |\ | | |/ ____ \| |____
|
747 | // |______/_/ \_\_| |______|_| \_\_____|_| |_|______|_| \_| |_/_/ \_\______|
|
748 |
|
749 | /**
|
750 | *
|
751 | * @param n Number of variates to return
|
752 | * @param loc Starting point. Must be a non-negative integer. 0 for degenerate distribution of 0.
|
753 | * @param p Probability of moving towards finish
|
754 | * @param cap Maximum steps before giving up
|
755 | * @param trace Variable to track progress
|
756 | * @returns {Array} Random variates array
|
757 | *
|
758 | * The FML distribution is a is based on the number of steps taken to return to the origin
|
759 | * from a given position, with transition probabilities set at the beginning by picking a
|
760 | * random variate from U(0,1).
|
761 | */
|
762 | rfml: function (n, loc, p, cap, trace) {
|
763 | n = this._v(n, "n");
|
764 | loc = this._v(loc, "nni", 1);
|
765 | if(p === undefined) p=this.prng;
|
766 | cap = this._v(cap, "n", 10000);
|
767 | if(trace === undefined) trace={};
|
768 |
|
769 | var toReturn = [];
|
770 |
|
771 | for(var i=0; i<n; i++) {
|
772 | var x = 0;
|
773 | var s = loc;
|
774 | var currP = p();
|
775 | if(loc === 0) {
|
776 |
|
777 | toReturn[i] = 0;
|
778 | } else {
|
779 |
|
780 | do {
|
781 |
|
782 | var trial = this.prng();
|
783 | if(trial < currP) {
|
784 | s++;
|
785 | trace[String(i) + "_" + String(x)] = { problems: s, p: currP, result: 1 }
|
786 | } else {
|
787 | s--;
|
788 | trace[String(i) + "_" + String(x)] = { problems: s, p: currP, result: -1 }
|
789 | }
|
790 | x++
|
791 | } while(s > 0 && x < cap);
|
792 |
|
793 | if(x === cap) x = -1; // Indicate we failed to do it in time.
|
794 | toReturn[i] = x;
|
795 | }
|
796 | }
|
797 | return toReturn
|
798 | },
|
799 |
|
800 | // http://www.statisticsblog.com/2013/05/uncovering-the-unreliable-friend-distribution-a-case-study-in-the-limits-of-mc-methods/
|
801 | /**
|
802 | *
|
803 | * The Unrelaible Friend distribution
|
804 | * @param n
|
805 | * @returns {Array} Random variates array
|
806 | */
|
807 | ruf: function(n) {
|
808 | n = this._v(n, "n");
|
809 |
|
810 | var toReturn = [];
|
811 |
|
812 | for(var i=0; i<n; i++) {
|
813 | toReturn[i] = this.rexp(1, this.prng())[0];
|
814 | }
|
815 |
|
816 | return toReturn
|
817 | }
|
818 | };
|
819 |
|
820 | // TODO: Add "perfect fake" functions: http://www.statisticsblog.com/2010/06/the-perfect-fake/
|
821 | // NOTES
|
822 | // Potential config options:
|
823 | // default entropy amount
|
824 | // Need pathway to make ready for secure applications (NIST/diehard?)
|
825 | // Always return a vector unless number is 1? This could be config option or put "1" at end of fcn to get 1 only
|
826 | // Separate out core random variate creation from number to create loop
|
827 | // TODO: To test out quality of randomness, stub in specific values for this.prng and make sure correct stuff is returned. |
\ | No newline at end of file |