-
Notifications
You must be signed in to change notification settings - Fork 7
Expand file tree
/
Copy pathrandomnotes.html
More file actions
815 lines (682 loc) · 61.4 KB
/
randomnotes.html
File metadata and controls
815 lines (682 loc) · 61.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
<!DOCTYPE html><html xmlns:dc="http://purl.org/dc/terms/" xmlns:og="http://ogp.me/ns#" ><head><meta http-equiv=Content-Type content="text/html; charset=utf-8"><title>More Random Sampling Methods</title><meta name="citation_pdf_url" content="https://peteroupc.github.io/randomnotes.pdf"><meta name="citation_url" content="https://peteroupc.github.io/randomnotes.html"><meta name="citation_title" content="More Random Sampling Methods"><meta name="dc.date" content="2025-09-30"><meta name="citation_date" content="2025/09/30"><meta name="citation_publication_date" content="2025/09/30"><meta name="citation_online_date" content="2025/09/30"><meta name="og:title" content="More Random Sampling Methods"><meta name="og:type" content="article"><meta name="og:url" content="https://peteroupc.github.io/randomnotes.html"><meta name="og:site_name" content="peteroupc.github.io"><meta name="dc.format" content="text/html"><meta name="dc.language" content="en"><meta name="title" content="More Random Sampling Methods"><meta name="dc.title" content="More Random Sampling Methods"><meta name="twitter:title" content="More Random Sampling Methods"><meta name="dc.creator" content="Peter Occil"/><meta name="author" content="Peter Occil"/><meta name="citation_author" content="Peter Occil"/><meta name="viewport" content="width=device-width"><link rel=stylesheet type="text/css" href="/style.css">
<script type="text/x-mathjax-config"> MathJax.Hub.Config({"HTML-CSS": { availableFonts: ["STIX","TeX"], linebreaks: { automatic:true }, preferredFont: "TeX" },
tex2jax: { displayMath: [ ["$$","$$"], ["\\[", "\\]"] ], inlineMath: [ ["$", "$"], ["\\\\(","\\\\)"] ], processEscapes: true } });
</script><script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-AMS_HTML-full"></script></head><body> <div class="header">
<nav><p><a href="#navigation">Menu</a> - <a href="#top">Top</a> - <a href="/">Home</a></nav></div>
<div class="mainarea" id="top">
<h1 id="more-random-sampling-methods">More Random Sampling Methods</h1><p>This version of the document is dated 2025-09-30. <a href='https://github.com/peteroupc/peteroupc.github.io/issues'>Post an issue or comment on this document</a>.</p>
<p><a href="mailto:poccil14@gmail.com"><strong>Peter Occil</strong></a></p>
<p><a id="Contents"></a></p>
<h2 id="contents">Contents</h2>
<ul>
<li><a href="#Contents"><strong>Contents</strong></a>
<ul>
<li><a href="#About_This_Document"><strong>About This Document</strong></a></li>
<li><a href="#Specific_Distributions"><strong>Specific Distributions</strong></a>
<ul>
<li><a href="#Normal_Gaussian_Distribution"><strong>Normal (Gaussian) Distribution</strong></a></li>
<li><a href="#Gamma_Distribution"><strong>Gamma Distribution</strong></a></li>
<li><a href="#Beta_Distribution"><strong>Beta Distribution</strong></a></li>
<li><a href="#Uniform_Partition_with_a_Positive_Sum"><strong>Uniform Partition with a Positive Sum</strong></a></li>
<li><a href="#Noncentral_Hypergeometric_Distributions"><strong>Noncentral Hypergeometric Distributions</strong></a></li>
<li><a href="#von_Mises_Distribution"><strong>von Mises Distribution</strong></a></li>
<li><a href="#Stable_Distribution"><strong>Stable Distribution</strong></a></li>
<li><a href="#Phase_Type_Distributions"><strong>Phase-Type Distributions</strong></a></li>
<li><a href="#Multivariate_Normal_Multinormal_Distribution"><strong>Multivariate Normal (Multinormal) Distribution</strong></a></li>
<li><a href="#Gaussian_and_Other_Copulas"><strong>Gaussian and Other Copulas</strong></a></li>
<li><a href="#Multivariate_Phase_Type_Distributions"><strong>Multivariate Phase-Type Distributions</strong></a></li>
</ul>
</li>
</ul>
</li>
<li><a href="#Notes"><strong>Notes</strong></a></li>
<li><a href="#Appendix"><strong>Appendix</strong></a>
<ul>
<li><a href="#Exact_Error_Bounded_and_Approximate_Algorithms"><strong>Exact, Error-Bounded, and Approximate Algorithms</strong></a></li>
</ul>
</li>
<li><a href="#License"><strong>License</strong></a></li>
</ul>
<p><a id="About_This_Document"></a></p>
<h3 id="about-this-document">About This Document</h3>
<p><strong>This is an open-source document; for an updated version, see the</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/raw/master/randomnotes.md"><strong>source code</strong></a> <strong>or its</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/blob/master/randomnotes.md"><strong>rendering on GitHub</strong></a><strong>. You can send comments on this document on the</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/issues"><strong>GitHub issues page</strong></a><strong>.</strong></p>
<p>My audience for this article is <strong>computer programmers with mathematics knowledge, but little or no familiarity with calculus</strong>.</p>
<p>I encourage readers to implement any of the algorithms given in this page, and report their implementation experiences. In particular, <a href="https://github.com/peteroupc/peteroupc.github.io/issues/18"><strong>I seek comments on the following aspects</strong></a>:</p>
<ul>
<li>Are the algorithms in this article (in conjunction with “<a href="https://peteroupc.github.io/randomfunc.html"><strong>Randomization and Sampling Methods</strong></a>”) easy to implement? Is each algorithm written so that someone could write code for that algorithm after reading the article?</li>
<li>Does this article have errors that should be corrected?</li>
<li>Are there ways to make this article more useful to the target audience?</li>
</ul>
<p>Comments on other aspects of this document are welcome.</p>
<p><a id="Specific_Distributions"></a></p>
<h3 id="specific-distributions">Specific Distributions</h3>
<p><strong>Requires random real numbers.</strong> This section shows algorithms to sample several popular nonuniform distributions. The algorithms are exact unless otherwise noted, and applications should choose algorithms with either no error (including rounding error) or a user-settable error bound. See the <a href="#Exact_Error_Bounded_and_Approximate_Algorithms"><strong>appendix</strong></a> for more information.</p>
<p><a id="Normal_Gaussian_Distribution"></a></p>
<h4 id="normal-gaussian-distribution">Normal (Gaussian) Distribution</h4>
<p>The <a href="https://en.wikipedia.org/wiki/Normal_distribution"><strong><em>normal distribution</em></strong></a> (also called the Gaussian distribution) takes the following two parameters:
- <code>mu</code> (μ) is the mean (average), or where the peak of the distribution’s “bell curve” is.
- <code>sigma</code> (σ), the standard deviation, affects how wide the “bell curve” appears. The
probability that a number sampled from the normal distribution will be within one standard deviation from the mean is about 68.3%; within two standard deviations (2 times <code>sigma</code>), about 95.4%; and within three standard deviations, about 99.7%. (Some publications give σ<sup>2</sup>, or variance, rather than standard deviation, as the second parameter. In this case, the standard deviation is the variance’s square root.)</p>
<p>There are a number of methods for sampling the normal distribution. An application can combine some or all of these.</p>
<ol>
<li>The ratio-of-uniforms method (given as <code>NormalRatioOfUniforms</code> below).</li>
<li>In the <em>Box–Muller transformation</em>, <code>mu + radius * cos(angle)</code> and <code>mu + radius * sin(angle)</code>, where <code>angle = RNDRANGEMinMaxExc(0, 2 * pi)</code> and <code>radius = sqrt(Expo(0.5)) * sigma</code>, are two independent values sampled from the normal distribution. The polar method (given as <code>NormalPolar</code> below) likewise produces two independent values sampled from that distribution at a time.</li>
<li>Karney’s algorithm to sample from the normal distribution, in a manner that minimizes approximation error and without using floating-point numbers (Karney 2016)<sup id="fnref:1"><a href="#fn:1" class="footnote" rel="footnote" role="doc-noteref">1</a></sup>.</li>
</ol>
<p>For surveys of Gaussian samplers, see (Thomas et al. 2007)<sup id="fnref:2"><a href="#fn:2" class="footnote" rel="footnote" role="doc-noteref">2</a></sup>, and (Malik and Hemani 2016)<sup id="fnref:3"><a href="#fn:3" class="footnote" rel="footnote" role="doc-noteref">3</a></sup>.</p>
<pre>METHOD NormalRatioOfUniforms(mu, sigma)
while true
a=RNDRANGEMinMaxExc(0,1)
bv = sqrt(2.0/exp(1.0))
// Or bv = 858/1000.0, which is also correct
b=RNDRANGEMinMaxExc(0,bv)
if b*b <= -a * a * 4 * ln(a)
return (RNDINT(1) * 2 - 1) *
(b * sigma / a) + mu
end
end
END METHOD
METHOD NormalPolar(mu, sigma)
while true
a = RNDRANGEMinMaxExc(0,1)
b = RNDRANGEMinMaxExc(0,1)
if RNDINT(1) == 0: a = 0 - a
if RNDINT(1) == 0: b = 0 - b
c = a * a + b * b
if c != 0 and c <= 1
c = sqrt(-ln(c) * 2 / c)
return [a * sigma * c + mu, b * sigma * c + mu]
end
end
END METHOD
</pre>
<blockquote>
<p><strong>Notes:</strong></p>
<ol>
<li>The <em>standard normal distribution</em> is implemented as <code>Normal(0, 1)</code>.</li>
<li>Methods implementing a variant of the normal distribution, the <em>discrete Gaussian distribution</em>, generate <em>integers</em> that closely follow the normal distribution. Examples include the one in (Karney 2016)<sup id="fnref:1:1"><a href="#fn:1" class="footnote" rel="footnote" role="doc-noteref">1</a></sup>, an improved version in (Du et al. 2021)<sup id="fnref:4"><a href="#fn:4" class="footnote" rel="footnote" role="doc-noteref">4</a></sup>, as well as so-called “constant-time” methods such as (Micciancio and Walter 2017)<sup id="fnref:5"><a href="#fn:5" class="footnote" rel="footnote" role="doc-noteref">5</a></sup> that are used above all in <em>lattice-based cryptography</em>.</li>
<li>The following are some approximations to the normal distribution that papers have suggested:
<ul>
<li>The sum of twelve <code>RNDRANGEMinMaxExc(0, sigma)</code> numbers, subtracted by 6 * <code>sigma</code>, to generate an approximate normal variate with mean 0 and standard deviation <code>sigma</code>. (Kabal 2000/2019)<sup id="fnref:6"><a href="#fn:6" class="footnote" rel="footnote" role="doc-noteref">6</a></sup> “warps” this sum in the following way (before adding the mean <code>mu</code>) to approximate the normal distribution better: <code>ssq = sum * sum; sum = ((((0.0000001141*ssq - 0.0000005102) * ssq + 0.00007474) * ssq + 0.0039439) * ssq + 0.98746) * sum</code>. See also <a href="https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution"><strong>“Irwin–Hall distribution”</strong></a>, namely the sum of <code>n</code> many <code>RNDRANGEMinMaxExc(0, 1)</code> numbers, on Wikipedia. D. Thomas (2014)<sup id="fnref:7"><a href="#fn:7" class="footnote" rel="footnote" role="doc-noteref">7</a></sup>, describes a more general approximation called CLT<sub>k</sub>, which combines <code>k</code> numbers in [0, 1] sampled from the uniform distribution as follows: <code>RNDRANGEMinMaxExc(0, 1) - RNDRANGEMinMaxExc(0, 1) + RNDRANGEMinMaxExc(0, 1) - ...</code>.</li>
<li>Numerical <a href="#Inverse_Transform_Sampling"><strong>inversions</strong></a> of the normal distribution’s cumulative distribution function (CDF, or the probability of getting X or less at random), including those by Wichura, by Acklam, and by Luu (Luu 2016)<sup id="fnref:8"><a href="#fn:8" class="footnote" rel="footnote" role="doc-noteref">8</a></sup>. See also <a href="https://www.johndcook.com/blog/normal_cdf_inverse/"><strong>“A literate program to compute the inverse of the normal CDF”</strong></a>.</li>
</ul>
</li>
<li>A pair of <em>q-Gaussian</em> random variates with parameter <code>q</code> less than 3 can be generated using the Box–Muller transformation, except <code>radius</code> is <code>radius=sqrt(-2*(pow(u,1-qp)-1)/(1-qp))</code> (where <code>qp=(1+q)/(3-q)</code> and <code>u=RNDRANGEMinMaxExc(0, 1)</code>), and the two variates are not statistically independent (Thistleton et al. 2007)<sup id="fnref:9"><a href="#fn:9" class="footnote" rel="footnote" role="doc-noteref">9</a></sup>.</li>
<li>A well-known result says that adding <code>n</code> many <code>Normal(0, 1)</code> variates, and dividing by <code>sqrt(n)</code>, results in a new <code>Normal(0, 1)</code> variate.</li>
</ol>
</blockquote>
<p><a id="Gamma_Distribution"></a></p>
<h4 id="gamma-distribution">Gamma Distribution</h4>
<p>The following method samples a number from a <em>gamma distribution</em> and is based on Marsaglia and Tsang’s method from 2000<sup id="fnref:10"><a href="#fn:10" class="footnote" rel="footnote" role="doc-noteref">10</a></sup> and (Liu et al. 2015)<sup id="fnref:11"><a href="#fn:11" class="footnote" rel="footnote" role="doc-noteref">11</a></sup>. Usually, the number expresses either—</p>
<ul>
<li>the lifetime (in days, hours, or other fixed units) of a random component with an average lifetime of <code>meanLifetime</code>, or</li>
<li>a random amount of time (in days, hours, or other fixed units) that passes until as many events as <code>meanLifetime</code> happen.</li>
</ul>
<p>Here, <code>meanLifetime</code> must be an integer or noninteger greater than 0.</p>
<pre>METHOD GammaDist(meanLifetime)
// Needs to be greater than 0
if meanLifetime <= 0: return error
// Exponential distribution special case if
// `meanLifetime` is 1 (see also (Devroye 1986), p. 405)
if meanLifetime == 1: return Expo(1)
if meanLifetime < 0.3 // Liu, Martin, Syring 2015
lamda = (1.0/meanLifetime) - 1
w = meanLifetime / (1-meanLifetime) * exp(1)
r = 1.0/(1+w)
while true
z = 0
x = RNDRANGEMinMaxExc(0, 1)
if x <= r: z = -ln(x/r)
else: z = -Expo(lamda)
ret = exp(-z/meanLifetime)
eta = 0
if z>=0: eta=exp(-z)
else: eta=w*lamda*exp(lamda*z)
if RNDRANGEMinMaxExc(0, eta) < exp(-ret-z): return ret
end
end
d = meanLifetime
v = 0
if meanLifetime < 1: d = d + 1
d = d - (1.0 / 3) // NOTE: 1.0 / 3 must be a fractional number
c = 1.0 / sqrt(9 * d)
while true
x = 0
while true
x = Normal(0, 1)
v = c * x + 1;
v = v * v * v
if v > 0: break
end
u = RNDRANGEMinMaxExc(0,1)
x2 = x * x
if u < 1 - (0.0331 * x2 * x2): break
if ln(u) < (0.5 * x2) + (d * (1 - v + ln(v))): break
end
ret = d * v
if meanLifetime < 1
ret = ret * pow(RNDRANGEMinMaxExc(0, 1), 1.0 / meanLifetime)
end
return ret
END METHOD
</pre>
<blockquote>
<p><strong>Notes:</strong></p>
<ol>
<li>The following is a useful identity for the gamma distribution: <code>GammaDist(a) = BetaDist(a, b - a) * GammaDist(b)</code> (Stuart 1962)<sup id="fnref:12"><a href="#fn:12" class="footnote" rel="footnote" role="doc-noteref">12</a></sup>.</li>
<li>The gamma distribution is usually defined to have a second parameter (called <code>theta</code> here), which is unfortunately defined differently in different works. For example, the gamma variate can be either multiplied or divided by <code>theta</code> depending on the work.</li>
<li>For other algorithms to sample from the gamma distribution, see Luengo (2022)<sup id="fnref:13"><a href="#fn:13" class="footnote" rel="footnote" role="doc-noteref">13</a></sup>.</li>
</ol>
<p><strong>Example:</strong> <strong>Moment exponential</strong> distribution (Dara and Ahmad 2012): <code>GammaDist(2)*beta</code> (or <code>(Expo(1)+Expo(1))*beta</code>), where <code>beta > 0</code>.</p>
</blockquote>
<p><a id="Beta_Distribution"></a></p>
<h4 id="beta-distribution">Beta Distribution</h4>
<p>The beta distribution takes on values on the interval (0, 1). Its two parameters, <code>a</code> and <code>b</code>, are both greater than 0 and describe the distribution’s shape. Depending on <code>a</code> and <code>b</code>, the shape can be a smooth peak or a smooth valley.</p>
<p>The following method samples a number from a <em>beta distribution</em>, in the interval [0, 1).</p>
<pre>METHOD BetaDist(a, b)
if b==1 and a==1: return RNDRANGEMinMaxExc(0, 1)
// Min-of-uniform
if a==1: return 1.0-pow(RNDRANGEMinMaxExc(0, 1),1.0/b)
// Max-of-uniform. Use only if a is small to
// avoid accuracy problems, as pointed out
// by Devroye 1986, p. 675.
if b==1 and a < 10: return pow(RNDRANGEMinMaxExc(0, 1),1.0/a)
x=GammaDist(a)
return x/(x+GammaDist(b))
END METHOD
</pre>
<p>I give an <a href="https://peteroupc.github.io/exporand.html"><strong>error-bounded sampler</strong></a> for the beta distribution (when <code>a</code> and <code>b</code> are both 1 or greater) in a separate page.</p>
<p><a id="Uniform_Partition_with_a_Positive_Sum"></a></p>
<h4 id="uniform-partition-with-a-positive-sum">Uniform Partition with a Positive Sum</h4>
<p>The following algorithm chooses at random a uniform partition of the number <code>sum</code> into <code>n</code> parts, and returns an <code>n</code>-item list of the chosen numbers, which sum to <code>sum</code> assuming no rounding error. In this algorithm, <code>n</code> must be an integer greater than 0, and <code>sum</code> must be greater than 0. The method was described in Bini and Buttazzo (2005)<sup id="fnref:14"><a href="#fn:14" class="footnote" rel="footnote" role="doc-noteref">14</a></sup> and Mai et al. (2022)<sup id="fnref:15"><a href="#fn:15" class="footnote" rel="footnote" role="doc-noteref">15</a></sup>.</p>
<p><code>
METHOD UniformSum(n, sum):
if n<=0 or sum<=0: return error
w=1; nn=n-1;ret=NewList()
while nn>0
v=w*(1-pow(RNDU01MinMaxExc(),1.0/nn))
ret.append(v*sum)
w=w-v; nn=nn-1
end
AddItem(ret, w*sum); return ret
END METHOD
</code></p>
<p><a id="Noncentral_Hypergeometric_Distributions"></a></p>
<h4 id="noncentral-hypergeometric-distributions">Noncentral Hypergeometric Distributions</h4>
<p>The following variants of the hypergeometric distribution are described in detail by Agner Fog in “<a href="https://cran.r-project.org/web/packages/BiasedUrn/vignettes/UrnTheory.pdf"><strong>Biased Urn Theory</strong></a>”.</p>
<p>Let there be <em>m</em> balls that each have one of two or more colors. For each color, assign each ball of that color the same weight (a real number 0 or greater). Then:</p>
<ol>
<li><strong>Wallenius’s hypergeometric distribution:</strong> Choose one ball not yet chosen, with probability equal to its weight divided by the sum of weights of balls not yet chosen. Repeat until exactly <em>n</em> items are chosen this way. Then for each color, count the number of items of that color chosen this way.</li>
<li><strong>Fisher’s hypergeometric distribution:</strong> For each ball, choose it with probability equal to its weight divided by the sum of weights of all balls. (Thus, each ball is independently chosen or not chosen depending on its weight.) If exactly <em>n</em> items were chosen this way, stop. Otherwise, start over. Then among the last <em>n</em> items chosen this way, count the number of items of each color.</li>
</ol>
<p>For both distributions, if there are two colors, there are four parameters: <em>m</em>, <em>ones</em>, <em>n</em>, <em>weight</em>, such that—</p>
<ul>
<li>for the first color, there are <em>ones</em> many balls each with weight <em>weight</em>;</li>
<li>for the second color, there are (<em>m</em>−<em>ones</em>) many balls each with weight 1; and</li>
<li>the random variate is the number of chosen balls of the first color.</li>
</ul>
<p><a id="von_Mises_Distribution"></a></p>
<h4 id="von-mises-distribution">von Mises Distribution</h4>
<p>The <em>von Mises distribution</em> describes a distribution of circular angles and uses two parameters: <code>mean</code> is the mean angle and <code>kappa</code> is a shape parameter. The distribution is uniform at <code>kappa = 0</code> and approaches a normal distribution with increasing <code>kappa</code>.</p>
<p>The algorithm below samples a number from the von Mises distribution, and is based on the Best–Fisher algorithm from 1979 (as described in (Devroye 1986)<sup id="fnref:16"><a href="#fn:16" class="footnote" rel="footnote" role="doc-noteref">16</a></sup> with errata incorporated).</p>
<pre>METHOD VonMises(mean, kappa)
if kappa < 0: return error
if kappa == 0
return RNDRANGEMinMaxExc(mean-pi, mean+pi)
end
r = 1.0 + sqrt(4 * kappa * kappa + 1)
rho = (r - sqrt(2 * r)) / (kappa * 2)
s = (1 + rho * rho) / (2 * rho)
while true
u = RNDRANGEMinMaxExc(-pi, pi)
v = RNDRANGEMinMaxExc(0, 1)
z = cos(u)
w = (1 + s*z) / (s + z)
y = kappa * (s - w)
if y*(2 - y) - v >=0 or ln(y / v) + 1 - y >= 0
if angle<-1: angle=-1
if angle>1: angle=1
// NOTE: Inverse cosine replaced here
// with `atan2` equivalent
angle = atan2(sqrt(1-w*w),w)
if u < 0: angle = -angle
return mean + angle
end
end
END METHOD
</pre>
<p><a id="Stable_Distribution"></a></p>
<h4 id="stable-distribution">Stable Distribution</h4>
<p>As more and more numbers, sampled independently at random in the same way, are added together, their distribution tends to a <a href="https://en.wikipedia.org/wiki/Stable_distribution"><strong><em>stable distribution</em></strong></a>, which resembles a curve with a single peak, but with generally “fatter” tails than the normal distribution. (Here, the stable distribution means the “alpha-stable distribution”.) The pseudocode below uses the Chambers–Mallows–Stuck algorithm. The <code>Stable</code> method, implemented below, takes two parameters:</p>
<ul>
<li><code>alpha</code> is a stability index in the interval (0, 2].</li>
<li><code>beta</code> is an asymmetry parameter in the interval [-1, 1]; if <code>beta</code> is 0, the curve is symmetric.</li>
</ul>
<p> </p>
<pre>METHOD Stable(alpha, beta)
if alpha <=0 or alpha > 2: return error
if beta < -1 or beta > 1: return error
halfpi = pi * 0.5
unif=RNDRANGEMinMaxExc(-halfpi, halfpi)
c=cos(unif)
expo=Expo(1)
if alpha == 1
s=sin(unif)
if beta == 0: return s/c
return 2.0*((unif*beta+halfpi)*s/c -
beta * ln(halfpi*expo*c/(unif*beta+halfpi)))/pi
else
z=-tan(alpha*halfpi)*beta
ug=unif+atan2(-z, 1)/alpha
cpow=pow(c, -1.0 / alpha)
return pow(1.0+z*z, 1.0 / (2*alpha))*
(sin(alpha*ug)*cpow)*
pow(cos(unif-alpha*ug)/expo, (1.0 - alpha) / alpha)
end
END METHOD
</pre>
<p>Methods implementing the strictly geometric stable and general geometric stable distributions are shown next (Kozubowski 2000)<sup id="fnref:17"><a href="#fn:17" class="footnote" rel="footnote" role="doc-noteref">17</a></sup>. Here, <code>alpha</code> is in (0, 2], <code>lamda</code> is greater than 0, and <code>tau</code>’s absolute value is not more than min(1, 2/<code>alpha</code> − 1). The result of <code>GeometricStable</code> is a symmetric Linnik distribution if <code>tau = 0</code>, or a Mittag-Leffler distribution if <code>tau = 1</code> and <code>alpha < 1</code>.</p>
<pre>METHOD GeometricStable(alpha, lamda, tau)
rho = alpha*(1-tau)/2
sign = -1
if tau==1 or RNDINT(1)==0 or RNDRANGEMinMaxExc(0, 1) < tau
rho = alpha*(1+tau)/2
sign = 1
end
w = 1
if rho != 1
rho = rho * pi
cotparam = RNDRANGEMinMaxExc(0, rho)
w = sin(rho)*cos(cotparam)/sin(cotparam)-cos(rho)
end
return Expo(1) * sign * pow(lamda*w, 1.0/alpha)
END METHOD
METHOD GeneralGeoStable(alpha, beta, mu, sigma)
z = Expo(1)
if alpha == 1: return mu*z+Stable(alpha, beta)*sigma*z+
sigma*z*beta*2*pi*ln(sigma*z)
else: return mu*z+
Stable(alpha, beta)*sigma*pow(z, 1.0/alpha)
END METHOD
</pre>
<p><a id="Phase_Type_Distributions"></a></p>
<h4 id="phase-type-distributions">Phase-Type Distributions</h4>
<p>A <em>phase-type distribution</em> models a sum of exponential random variates driven by a <a href="https://peteroupc.github.io/randomnotes.html"><strong>Markov chain</strong></a>. The Markov chain has <code>n</code> normal states and one “absorbing” or terminating state. This distribution has two parameters:</p>
<ul>
<li><code>alpha</code>, an <code>n</code>-item array showing the probability of starting the chain at each normal state.</li>
<li><code>s</code>, an <code>n</code>×<code>n</code> <em>subgenerator matrix</em>, a list of <code>n</code> lists of <code>n</code> values each. The values in each list (each normal state of the Markov chain) must sum to 0 or less, and for each state <code>i</code>, <code>s[i][i]</code> is 0 minus the rate of that state’s exponential random variate, and each entry <code>s[i][j]</code> with <code>i!=j</code> is the relative probability for moving to state <code>j</code>.</li>
</ul>
<p>The method <code>PhaseType</code>, given later, samples from a phase-type distribution given the two parameters above. (The pseudocode assumes each number in <code>alpha</code> and <code>s</code> is a rational number, because it uses <code>NormalizeRatios</code>.)</p>
<p>```
METHOD GenToTrans(s)
// Converts a subgenerator matrix to a
// more intuitive transition matrix.
m=[];
for j in 0…size(s)
m[j]=[]; for i in 0…size(s)+1: AddItem(m[j],0)
end
for i in 0…size(s)
isum=Sum(s[i])
if isum<0: m[i][size(s)]=isum/s[i][i]
for j in 0…size(s)
if j!=i: m[i][j]=-s[i][j]/s[i][i]
end
end
return m
END METHOD</p>
<p>METHOD PhaseType(alpha, s)
// Setup
trans=GenToTrans(s)
// Sampling
state=WeightedChoice(NormalizeRatios(alpha))
ret=0
while state<size(s)
ret=ret+Expo(-s[state][state])
state=WeightedChoice(NormalizeRatios(trans[state]))
end
return ret
END METHOD
```</p>
<blockquote>
<p><strong>Note:</strong> An <strong>inhomogeneous phase-type</strong> random variate has the form <code>G(PhaseType(alpha, s))</code>, where <code>G(x)</code> is a function designed to control the heaviness of the distribution’s tail (Bladt 2021)<sup id="fnref:18"><a href="#fn:18" class="footnote" rel="footnote" role="doc-noteref">18</a></sup>. For example, <code>G(x) = pow(x, 1.0/beta)</code>, where <code>beta>0</code>, leads to a tail as heavy as a Weibull distribution.</p>
</blockquote>
<p><a id="Multivariate_Normal_Multinormal_Distribution"></a></p>
<h4 id="multivariate-normal-multinormal-distribution">Multivariate Normal (Multinormal) Distribution</h4>
<p>The following pseudocode generates a random vector (list of numbers) that follows a <a href="https://en.wikipedia.org/wiki/Multivariate_normal_distribution"><strong><em>multivariate normal (multinormal) distribution</em></strong></a>. The method <code>MultivariateNormal</code> takes the following parameters:</p>
<ul>
<li>A list, <code>mu</code> (μ), which indicates the means to add to the random vector’s components. <code>mu</code> can be <code>nothing</code>, in which case each component will have a mean of zero.</li>
<li>A list of lists <code>cov</code>, that specifies a <em>covariance matrix</em> (Σ), a symmetric positive definite N×N matrix, where N is the number of components of the random vector. (An N×N matrix is <em>positive definite</em> if its determinant [overall scale] is greater than 0 and if either the matrix is 1 × 1 or a smaller matrix formed by removing the last row and column is positive definite.)</li>
</ul>
<p> </p>
<pre>METHOD Decompose(matrix)
numrows = size(matrix)
if size(matrix[0])!=numrows: return error
// Does a Cholesky decomposition of a matrix
// assuming it's positive definite and invertible
ret=NewList()
for i in 0...numrows
submat = NewList()
for j in 0...numrows: AddItem(submat, 0)
AddItem(ret, submat)
end
s1 = sqrt(matrix[0][0])
if s1==0: return ret // For robustness
for i in 0...numrows
ret[0][i]=matrix[0][i]*1.0/s1
end
for i in 0...numrows
msum=0.0
for j in 0...i: msum = msum + ret[j][i]*ret[j][i]
sq=matrix[i][i]-msum
if sq<0: sq=0 // For robustness
ret[i][i]=math.sqrt(sq)
end
for j in 0...numrows
for i in (j + 1)...numrows
// For robustness
if ret[j][j]==0: ret[j][i]=0
if ret[j][j]!=0
msum=0
for k in 0...j: msum = msum + ret[k][i]*ret[k][j]
ret[j][i]=(matrix[j][i]-msum)*1.0/ret[j][j]
end
end
end
return ret
END METHOD
METHOD VecAdd(a, b)
c=[]; for j in 0...size(a): c[j]=a[j]+b[j]
return c
END METHOD
METHOD VecScale(a, scalar)
c=[]; for j in 0...size(a): c[j]=a[j]*scalar
return c
END METHOD
METHOD MultivariateNormal(mu, cov)
vars=NewList()
for j in 0...mulen: AddItem(vars, Normal(0, 1))
return MultivariateCov(mu,cov,vars)
END METHOD
METHOD MultivariateCov(mu, cov, vars)
// Returns mu + cov^(1/2)*vars
mulen=size(cov)
if mu != nothing
mulen = size(mu)
if mulen!=size(cov): return error
if mulen!=size(cov[0]): return error
end
// NOTE: If multiple random points will
// be generated using the same covariance
// matrix, an implementation can consider
// precalculating the decomposed matrix
// in advance rather than calculating it here.
cho=Decompose(cov)
i=0
ret=NewList()
while i<mulen
msum = 0
for j in 0...mulen: msum=cho[j][i]*vars[j]
AddItem(ret, msum)
i=i+1
end
if mu!=nothing: ret=VecAdd(ret, mu)
return ret
end
</pre>
<blockquote>
<p><strong>Note:</strong> The <a href="https://peteroupc.github.io/randomgen.zip"><strong>Python sample code</strong></a> contains a variant of this
method for generating multiple random vectors in one call.</p>
<p><strong>Examples:</strong></p>
<ol>
<li>A vector that follows a <strong>binormal distribution</strong> (two-variable multinormal distribution) is a vector of two numbers from the normal distribution, and can be sampled using the following idiom: <code>MultivariateNormal([mu1, mu2], [[s1*s1, s1*s2*rho], [rho*s1*s2, s2*s2]])</code>, where <code>mu1</code> and <code>mu2</code> are the means of the vector’s two components, <code>s1</code> and <code>s2</code> are their standard deviations, and <code>rho</code> is a <em>correlation coefficient</em> greater than -1 and less than 1 (0 means no correlation).</li>
<li><strong>Log-multinormal distribution</strong>: Generate a multinormal random vector, then apply <code>exp(n)</code> to each component <code>n</code>.</li>
<li>A <strong>Beckmann distribution</strong>: Generate a random binormal vector <code>vec</code>, then apply <code>PNorm(vec, 2)</code> to that vector. (<code>PNorm</code> is given in the main page’s section “<a href="https://peteroupc.github.io/randomnotes.md#Random_Points_on_a_Sphere"><strong>Random Points on a Sphere</strong></a>.”)</li>
<li>A <strong>Rice (Rician) distribution</strong> is a Beckmann distribution in which the binormal random pair is generated with <code>m1 = m2 = a / sqrt(2)</code>, <code>rho = 0</code>, and <code>s1 = s2 = b</code>, where <code>a</code> and <code>b</code> are the parameters to the Rice distribution.</li>
<li><strong>Rice–Norton distribution</strong>: Generate <code>vec = MultivariateNormal(</code> <code>[v,v,v], [[w,0,0], [0,w,0],[0,0,w]])</code> (where <code>v = a/sqrt(m*2)</code>, <code>w = b*b/m</code>, and <code>a</code>, <code>b</code>, and <code>m</code> are the parameters to the Rice–Norton distribution), then apply <code>PNorm(vec, 2)</code> to that vector.</li>
<li>A <strong>standard</strong> <a href="https://en.wikipedia.org/wiki/Complex_normal_distribution"><strong>complex normal distribution</strong></a> is a binormal distribution in which the binormal random pair is generated with <code>s1 = s2 = sqrt(0.5)</code> and <code>mu1 = mu2 = 0</code> and treated as the real and imaginary parts of a complex number.</li>
<li><strong>Multivariate Linnik distribution</strong>: Generate a multinormal random vector, then multiply each component by <code>x</code>, where <code>x = GeometricStable(alpha/2.0, 1, 1)</code>, where <code>alpha</code> is a parameter in (0, 2] (Kozubowski 2000)<sup id="fnref:17:1"><a href="#fn:17" class="footnote" rel="footnote" role="doc-noteref">17</a></sup>.</li>
<li><strong>Multivariate exponential power distribution</strong> (Solaro 2004)<sup id="fnref:27"><a href="#fn:27" class="footnote" rel="footnote" role="doc-noteref">19</a></sup>: <code>MultivariateCov(mu, cov, vec)</code>, where <code>vec = RandomPointOnSphere(m, pow(Gamma(m/s,1)*2,1.0/s), 2)</code>, <code>m</code> is the dimension, <code>s > 0</code> is a shape parameter, <code>mu</code> is the mean as an <code>m</code>-dimensional vector (<code>m</code>-item list), and <code>cov</code> is a covariance matrix.</li>
<li><strong>Elliptical distribution</strong>: <code>MultivariateCov(mu, cov, RandomPointOnSphere(dims, z, 2))</code>, where <code>z</code> is an arbitrary random variate,<code>dims</code> is the number of dimensions, <code>mu</code> is a <code>dims</code>-dimensional location vector, and <code>cov</code> is a <code>dims</code>×<code>dims</code> covariance matrix. See, for example, Fang et al. (1990)<sup id="fnref:19"><a href="#fn:19" class="footnote" rel="footnote" role="doc-noteref">20</a></sup></li>
<li><strong>Mean-variance mixture of normal distributions</strong> (Barndorff-Nielsen et al. 1982)<sup id="fnref:3:1"><a href="#fn:3" class="footnote" rel="footnote" role="doc-noteref">3</a></sup>: <code>VecAdd(mu, VecAdd(VecScale(delta, v), VecScale(MultivariateNormal(nothing, cov), sqrt(z))))</code>, where <code>mu</code> and <code>delta</code> are<code>n</code>-dimensional vectors, <code>cov</code> is a covariance matrix, and <code>v</code> is an arbitrary random variate 0 or greater.</li>
<li><strong>Mean mixture of normal distributions</strong> (Bhagwat and Marchand 2022)<sup id="fnref:4:1"><a href="#fn:4" class="footnote" rel="footnote" role="doc-noteref">4</a></sup>: <code>MultivariateNormal(VecAdd(theta,VecScale(a,v)), cov)</code> where <code>theta</code> is an <code>n</code>-dimensional location vector, <code>a</code> is an <code>n</code>-dimensional “perturbation vector”, <code>cov</code> is a covariance matrix, and <code>v</code> is an arbitrary random variate.</li>
</ol>
</blockquote>
<p><a id="Gaussian_and_Other_Copulas"></a></p>
<h4 id="gaussian-and-other-copulas">Gaussian and Other Copulas</h4>
<p>A <em>copula</em> is a way to describe the dependence between randomly sampled numbers.</p>
<p>One example is a <em>Gaussian copula</em>; this copula is sampled by sampling from a <a href="#Multivariate_Normal_Multinormal_Distribution"><strong>multinormal distribution</strong></a>, then converting the resulting numbers to <em>dependent</em> uniform random values. In the following pseudocode, which implements a Gaussian copula:</p>
<ul>
<li>The parameter <code>covar</code> is the covariance matrix for the multinormal distribution.</li>
<li><code>erf(v)</code> is the <a href="https://en.wikipedia.org/wiki/Error_function"><strong>error function</strong></a> of the number <code>v</code>.</li>
</ul>
<p> </p>
<pre>METHOD GaussianCopula(covar)
mvn=MultivariateNormal(nothing, covar)
for i in 0...size(covar)
// Apply the normal distribution's CDF
// to get uniform numbers
mvn[i] = (erf(mvn[i]/(sqrt(2)*sqrt(covar[i][i])))+1)*0.5
end
return mvn
END METHOD
</pre>
<p>Each of the resulting uniform random values will be in the interval [0, 1], and each one can be further transformed to any other probability distribution (which is called a <em>marginal distribution</em> or <em>marginal</em> here) by taking the quantile of that uniform number for that distribution (see “<a href="https://peteroupc.github.io/randomfunc.html#Inverse_Transform_Sampling"><strong>Inverse Transform Sampling</strong></a>”, and see also (Cario and Nelson 1997)<sup id="fnref:20"><a href="#fn:20" class="footnote" rel="footnote" role="doc-noteref">21</a></sup>.)</p>
<blockquote>
<p><strong>Note:</strong> The Gaussian copula is also known as the <em>normal-to-anything</em> method.</p>
<p><strong>Examples:</strong></p>
<ol>
<li>To generate two correlated uniform random values with a Gaussian copula, generate <code>GaussianCopula([[1, rho], [rho, 1]])</code>, where <code>rho</code> is the Pearson correlation coefficient, in the interval [-1, 1]. (Other correlation coefficients besides <code>rho</code> exist. For example, for a two-variable Gaussian copula, the <a href="https://en.wikipedia.org/wiki/Rank_correlation"><strong>Spearman correlation coefficient</strong></a> <code>srho</code> can be converted to <code>rho</code> by <code>rho = sin(srho * pi / 6) * 2</code>. Other correlation coefficients, and other measures of dependence between randomly sampled numbers, are not further discussed in this document.)</li>
<li>
<p>The following example generates a two-dimensional random vector that follows a Gaussian copula with exponential marginals (<code>rho</code> is the Pearson correlation coefficient, and <code>rate1</code> and <code>rate2</code> are the rates of the two exponential marginals).</p>
<pre> METHOD CorrelatedExpo(rho, rate1, rate2)
copula = GaussianCopula([[1, rho], [rho, 1]])
// Transform to exponentials using that
// distribution's quantile function
return [-log1p(-copula[0]) / rate1,
-log1p(-copula[1]) / rate2]
END METHOD
</pre>
</li>
<li>
<p>The <em><strong>T</strong>–Poisson hierarchy</em> (Knudson et al. 2021)<sup id="fnref:21"><a href="#fn:21" class="footnote" rel="footnote" role="doc-noteref">22</a></sup> is a way to generate N-dimensional Poisson-distributed random vectors via copulas. Each of the N dimensions is associated with—</p>
<ul>
<li>a parameter <code>lamda</code>, and</li>
<li>a marginal distribution that may not be discrete and takes on only nonnegative values.</li>
</ul>
<p>To sample from the <strong>T</strong>–Poisson hierarchy—</p>
<ol>
<li>sample an N-dimensional random vector via a copula (such as <code>GaussianCopula</code>), producing an N-dimensional vector of correlated uniform numbers; then</li>
<li>for each component in the vector, replace it with that component’s quantile for the corresponding marginal; then</li>
<li>for each component in the vector, replace it with <code>Poisson(lamda * c)</code>, where <code>c</code> is that component and <code>lamda</code> is the <code>lamda</code> parameter for the corresponding dimension.</li>
</ol>
<p>The following example implements the T-Poisson hierarchy using a Gaussian copula and exponential marginals.</p>
<pre> METHOD PoissonH(rho, rate1, rate2, lambda1, lambda2)
vec = CorrelatedExpo(rho, rate1, rate2)
return [Poisson(lambda1*vec[0]),Poisson(lambda2*vec[1])]
END METHOD
</pre>
</li>
</ol>
</blockquote>
<p>Other kinds of copulas describe different kinds of dependence between randomly sampled numbers. Examples of other copulas are—</p>
<ul>
<li>the <strong>Fréchet–Hoeffding upper bound copula</strong> <em>[x, x, …, x]</em> (for example, <code>[x, x]</code>), where <code>x = RNDRANGEMinMaxExc(0, 1)</code>,</li>
<li>the <strong>Fréchet–Hoeffding lower bound copula</strong> <code>[x, 1.0 - x]</code> where <code>x = RNDRANGEMinMaxExc(0, 1)</code>,</li>
<li>the <strong>product copula</strong>, where each number is a separately generated <code>RNDRANGEMinMaxExc(0, 1)</code> (indicating no dependence between the numbers), and</li>
<li>the <strong>Archimedean copulas</strong>, described by M. Hofert and M. Mächler (2011)<sup id="fnref:22"><a href="#fn:22" class="footnote" rel="footnote" role="doc-noteref">23</a></sup>.</li>
</ul>
<p><a id="Multivariate_Phase_Type_Distributions"></a></p>
<h4 id="multivariate-phase-type-distributions">Multivariate Phase-Type Distributions</h4>
<p>The following pseudocode generates a random vector (of <code>d</code> coordinates) following a <em>multivariate phase-type distribution</em> called MPH*. In addition to parameters <code>alpha</code> and <code>s</code>, there is also a <em>reward matrix</em> <code>r</code>, such that <code>r[i][j]</code> is the probability of adding to coordinate <code>j</code> when state <code>i</code> is visited. (The pseudocode assumes each number in <code>alpha</code>, <code>s</code>, and <code>r</code> is a rational number, because it uses <code>NormalizeRatios</code>.)</p>
<p><code>
METHOD MPH(alpha, s, r)
if len(r[0])<1 or len(r)!=len(s): return error
// Setup
trans=GenToTrans(s)
ret=[]; for i in 0...size(r[0]): AddItem(ret,0)
// Sampling
state=WeightedChoice(NormalizeRatios(alpha))
ret=0
while state<size(s)
rs=WeightedChoice(NormalizeRatios(r[state]))
ret[rs]=ret[rs]+Expo(-s[state][state])
state=WeightedChoice(NormalizeRatios(trans[state]))
end
return ret
END METHOD
</code></p>
<blockquote>
<p><strong>Note:</strong> An inhomogeneous version of MPH* can be as follows: <code>[G1(mph[1]), G2(mph[2]), ..., GD(mph[d])]</code>, where <code>mph</code> is a <code>d</code>-dimensional MPH* vector and <code>G1</code>, <code>G2</code>, …, <code>GD</code> are strictly increasing functions whose domain and range are the positive real line and whose “slope” is defined on the whole domain (Albrecher et al. 2022)<sup id="fnref:23"><a href="#fn:23" class="footnote" rel="footnote" role="doc-noteref">24</a></sup>.</p>
</blockquote>
<p><a id="Notes"></a></p>
<h2 id="notes">Notes</h2>
<p><a id="Appendix"></a></p>
<h2 id="appendix">Appendix</h2>
<p><a id="Exact_Error_Bounded_and_Approximate_Algorithms"></a></p>
<h3 id="exact-error-bounded-and-approximate-algorithms">Exact, Error-Bounded, and Approximate Algorithms</h3>
<p>There are three kinds of randomization algorithms:</p>
<ol>
<li>
<p>An <em>exact algorithm</em> is an algorithm that samples from the exact distribution requested, assuming that computers—</p>
<ul>
<li>can store and operate on real numbers (which have unlimited precision), and</li>
<li>can generate independent uniform random real numbers</li>
</ul>
<p>(Devroye 1986, p. 1-2)<sup id="fnref:16:1"><a href="#fn:16" class="footnote" rel="footnote" role="doc-noteref">16</a></sup>. However, an exact algorithm implemented on real-life computers can incur error due to the use of fixed precision (especially floating-point numbers), such as rounding and cancellations. An exact algorithm can achieve a guaranteed bound on accuracy (and thus be an <em>error-bounded algorithm</em>) using either arbitrary-precision or interval arithmetic (see also Devroye 1986, p. 2)<sup id="fnref:16:2"><a href="#fn:16" class="footnote" rel="footnote" role="doc-noteref">16</a></sup>. All methods given on this page are exact unless otherwise noted. Note that the <code>RNDRANGEMinMaxExc</code> method is exact in theory, but has no required implementation.</p>
</li>
<li>
<p>An <em>error-bounded algorithm</em> is a sampling algorithm with the following requirements:</p>
<ul>
<li>If the ideal distribution is discrete (takes on values that can map to integers and back without loss), the algorithm samples exactly from that distribution. (But see the note later.)</li>
<li>If the ideal distribution is not discrete, the algorithm samples from a distribution that is close to the ideal within a user-specified error tolerance (see later for details). The algorithm can instead sample a number from the distribution only partially, as long as the fully sampled number can be made close to the ideal within any error tolerance desired.</li>
<li>In sampling from a distribution, the algorithm incurs no approximation error not already present in the inputs (except errors needed to round the final result to the user-specified error tolerance).</li>
</ul>
<p>Many error-bounded algorithms use random bits as their only source of randomness. An application should use error-bounded algorithms whenever possible.</p>
<p>Most algorithms on this page, though, are not <em>error-bounded</em> when naïvely implemented in most number formats (including floating-point numbers). (There are number formats such as “constructive reals” or “recursive reals” that allow real numbers to be approximated to a user-specified error (Boehm 2020)<sup id="fnref:24"><a href="#fn:24" class="footnote" rel="footnote" role="doc-noteref">25</a></sup>.)</p>
</li>
<li>
<p>An <em>inexact</em>, <em>approximate</em>, or <em>biased algorithm</em> is any sampling algorithm that is neither exact nor error-bounded. This includes algorithms that sample from a distribution that is close to the desired distribution, but not within a user-specified error tolerance (see also Devroye 1986, p. 2)<sup id="fnref:16:3"><a href="#fn:16" class="footnote" rel="footnote" role="doc-noteref">16</a></sup>. An application should use this kind of algorithm only if it’s willing to trade accuracy for speed.</p>
</li>
</ol>
<p>There are many ways to describe closeness between two distributions. One suggestion by Devroye and Gravel (2020)<sup id="fnref:25"><a href="#fn:25" class="footnote" rel="footnote" role="doc-noteref">26</a></sup> is Wasserstein distance (or “earth-mover distance”), which they proved has a simple definition in terms of the quantile function (Theorem 8). Here, an algorithm has accuracy ε (the user-specified error tolerance) if it samples from a distribution that is close to the ideal distribution by a Wasserstein distance of not more than ε.</p>
<blockquote>
<p><strong>Examples:</strong></p>
<ol>
<li>Sampling from the exponential distribution via <code>-ln(RNDRANGEMinMaxExc(0, 1))</code> is an <em>exact algorithm</em> (in theory), but not an <em>error-bounded</em> one for common floating-point number formats. The same is true of the Box–Muller transformation.</li>
<li>Karney’s algorithm for the normal distribution (Karney 2016)<sup id="fnref:1:2"><a href="#fn:1" class="footnote" rel="footnote" role="doc-noteref">1</a></sup>, as well as Karney’s implementation of von Neumann’s exponential distribution sampler (Karney 2016)<sup id="fnref:1:3"><a href="#fn:1" class="footnote" rel="footnote" role="doc-noteref">1</a></sup> are both <em>error-bounded</em>, because they return a result that can be made to come close to the normal or exponential distribution, respectively, within any error tolerance desired simply by appending more random digits to the end. See also (Oberhoff 2018)<sup id="fnref:26"><a href="#fn:26" class="footnote" rel="footnote" role="doc-noteref">27</a></sup>.</li>
<li>Examples of <em>approximate algorithms</em> include sampling from a Gaussian-like distribution via a sum of <code>RNDRANGEMinMaxExc(0, 1)</code>, or most cases of modulo reduction to produce uniform-like integers at random (see notes in the section “<a href="https://peteroupc.github.io/randomfunc.html#RNDINT_Random_Integers_in_0_N"><strong>RNDINT</strong></a>”). The following approximate algorithm for the Poisson distribution is another example (Giammatteo and Di Mascio (2020)<sup id="fnref:27:1"><a href="#fn:27" class="footnote" rel="footnote" role="doc-noteref">19</a></sup>): <code>floor(1.0/3 + pow(max(0, Normal(0, 1)*pow(mean,1/6.0)*2/3 + pow(mean, 2.0/3)), 3.0/2))</code>, where <code>mean</code> is greater than 50.</li>
</ol>
<p><strong>Note:</strong> A discrete distribution can be sampled in finite time on average if and only if its so-called <em>Shannon entropy</em> is finite (Knuth and Yao 1976)<sup id="fnref:28"><a href="#fn:28" class="footnote" rel="footnote" role="doc-noteref">28</a></sup>. Unfortunately, some discrete distributions have infinite Shannon entropy, such as some members of the zeta Dirichlet family of distributions (Devroye and Gravel 2020)<sup id="fnref:25:1"><a href="#fn:25" class="footnote" rel="footnote" role="doc-noteref">26</a></sup>. Thus, in practice, an approximate or error-bounded sampler is needed for these distributions. Saad et al. (2020)<sup id="fnref:29"><a href="#fn:29" class="footnote" rel="footnote" role="doc-noteref">29</a></sup> discuss how to sample an approximation of a discrete distribution with a user-specified error tolerance, but only if the ideal distribution takes on a finite number of values (and thus has finite Shannon entropy). On the other hand, a distribution has finite Shannon entropy whenever—</p>
<ul>
<li>it takes on only integers 1 or greater and has a finite <em>t</em><sup>th</sup> moment for some <em>t</em> > 0 (“long-run average” of values raised to <em>t</em><sup>th</sup> power) (Baccetti and Visser 2013)<sup id="fnref:30"><a href="#fn:30" class="footnote" rel="footnote" role="doc-noteref">30</a></sup>, or as a special case,</li>
<li>it takes on only integers 1 or greater and has a finite mean (“long-run average”), or</li>
<li>it has the form <em>X</em> + <em>n</em>, where <em>n</em> is a constant and <em>X</em> is a random variate whose distribution has finite Shannon entropy.</li>
</ul>
</blockquote>
<p><a id="License"></a></p>
<h2 id="license">License</h2>
<p>Any copyright to this page is released to the Public Domain. In case this is not possible, this page is also licensed under <a href="https://creativecommons.org/publicdomain/zero/1.0/"><strong>Creative Commons Zero</strong></a>.</p>
<div class="footnotes" role="doc-endnotes">
<ol>
<li id="fn:1">
<p>Karney, C.F.F., 2016. Sampling exactly from the normal distribution. ACM Transactions on Mathematical Software (TOMS), 42(1), pp.1-14. Also: “<a href="https://arxiv.org/abs/1303.6257v2"><strong>Sampling exactly from the normal distribution</strong></a>”, arXiv:1303.6257v2 [physics.comp-ph], 2014. <a href="#fnref:1" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:1:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a> <a href="#fnref:1:2" class="reversefootnote" role="doc-backlink">↩<sup>3</sup></a> <a href="#fnref:1:3" class="reversefootnote" role="doc-backlink">↩<sup>4</sup></a></p>
</li>
<li id="fn:2">
<p>Thomas, D., et al., “Gaussian Random Number Generators”, <em>ACM Computing Surveys</em> 39(4), 2007. <a href="#fnref:2" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:3">
<p>Malik, J.S., Hemani, A., “Gaussian random number generation: A survey on hardware architectures”, <em>ACM Computing Surveys</em> 49(3), 2016. <a href="#fnref:3" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:3:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a></p>
</li>
<li id="fn:4">
<p>Du, Yusong, Baoying Fan, and Baodian Wei. “An improved exact sampling algorithm for the standard normal distribution.” Computational Statistics (2021): 1-17, also arXiv:2008.03855 [cs.DS]. <a href="#fnref:4" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:4:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a></p>
</li>
<li id="fn:5">
<p>Micciancio, D. and Walter, M., “Gaussian sampling over the integers: Efficient, generic, constant-time”, in Annual International Cryptology Conference, August 2017 (pp. 455-485). <a href="#fnref:5" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:6">
<p>Kabal, P., “Generating Gaussian Pseudo-Random Variates”, McGill University, 2000/2019. <a href="#fnref:6" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:7">
<p>Thomas, D.B., 2014, May. FPGA Gaussian random number generators with guaranteed statistical accuracy. In <em>2014 IEEE 22nd Annual International Symposium on Field-Programmable Custom Computing Machines</em> (pp. 149-156). <a href="#fnref:7" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:8">
<p>Luu, T., “Fast and Accurate Parallel Computation of Quantile Functions for Random Number Generation”, Dissertation, University College London, 2016. <a href="#fnref:8" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:9">
<p>Thistleton, W., Marsh, J., et al., “Generalized Box-Müller Method for Generating q-Gaussian Random Deviates”, <em>IEEE Transactions on Information Theory</em> 53(12), 2007. <a href="#fnref:9" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:10">
<p>Marsaglia, G., Tsang, W.W., “A simple method for generating gamma variables”, <em>ACM Transactions on Mathematical Software</em> 26(3), 2000. <a href="#fnref:10" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:11">
<p>Liu, C., Martin, R., Syring, N., “<a href="https://arxiv.org/abs/1302.1884v3"><strong>Simulating from a gamma distribution with small shape parameter</strong></a>”, arXiv:1302.1884v3 [stat.CO], 2015. <a href="#fnref:11" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:12">
<p>A. Stuart, “Gamma-distributed products of independent random variables”, <em>Biometrika</em> 49, 1962. <a href="#fnref:12" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:13">
<p>Luengo, E.A., “<a href="https://dl.acm.org/doi/abs/10.1145/3527157"><strong>Gamma Pseudo Random Number Generators</strong></a>”, <em>ACM Computing Surveys</em>, 2022. <a href="#fnref:13" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:14">
<p>Bini, B., Buttazzo, G.C., “Measuring the Performance of Schedulability Tests”, <em>Real-Time Systems</em> 30, 129-154 (2005) <a href="#fnref:14" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:15">
<p>Mai, J., Craig, J.R., Tolson, B.A., “<a href="https://www.sciencedirect.com/science/article/pii/S1364815221003248"><strong>The pie-sharing problem: Unbiased sampling of N+1 summative weights</strong></a>”, <em>Environmental Modelling & Software</em> 148, February 2022. <a href="#fnref:15" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:16">
<p>Devroye, L., <a href="http://luc.devroye.org/rnbookindex.html"><strong><em>Non-Uniform Random Variate Generation</em></strong></a>, 1986. <a href="#fnref:16" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:16:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a> <a href="#fnref:16:2" class="reversefootnote" role="doc-backlink">↩<sup>3</sup></a> <a href="#fnref:16:3" class="reversefootnote" role="doc-backlink">↩<sup>4</sup></a></p>
</li>
<li id="fn:17">
<p>Tomasz J. Kozubowski, “<a href="https://www.sciencedirect.com/science/article/pii/S0377042799003180"><strong>Computer simulation of geometric stable distributions</strong></a>”, <em>Journal of Computational and Applied Mathematics</em> 116(2), pp. 221-229, 2000. <a href="https://doi.org/10.1016/S0377-0427%2899%2900318-0"><strong>https://doi.org/10.1016/S0377-0427(99)00318-0</strong></a> <a href="#fnref:17" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:17:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a></p>
</li>
<li id="fn:18">
<p>Bladt, Martin. “Phase-type distributions for claim severity regression modeling.” ASTIN Bulletin: The Journal of the IAA (2021): 1-32. <a href="#fnref:18" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:27">
<p>Giammatteo, P., and Di Mascio, T., “Wilson-Hilferty-type approximation for Poisson Random Variable”, Advances in Science, Technology and Engineering Systems Journal 5(2), 2020. <a href="#fnref:27" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:27:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a></p>
</li>
<li id="fn:19">
<p>Fang et al., Symmetric multivariate and related distributions, 1990. <a href="#fnref:19" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:20">
<p>Cario, M. C., B. L. Nelson, “Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix”, 1997. <a href="#fnref:20" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:21">
<p>Knudson, A.D., Kozubowski, T.J., et al., “A flexible multivariate model for high-dimensional correlated count data”, <em>Journal of Statistical Distributions and Applications</em> 8:6, 2021. <a href="#fnref:21" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:22">
<p>Hofert, M., and Maechler, M. “Nested Archimedean Copulas Meet R: The nacopula Package”. <em>Journal of Statistical Software</em> 39(9), 2011, pp. 1-20. <a href="#fnref:22" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:23">
<p>Albrecher, Hansjörg, Mogens Bladt, and Jorge Yslas. “Fitting inhomogeneous phase‐type distributions to data: the univariate and the multivariate case.” Scandinavian Journal of Statistics 49, no. 1 (2022): 44-77. <a href="#fnref:23" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:24">
<p>Boehm, Hans-J. “Towards an API for the real numbers.” In Proceedings of the 41st ACM SIGPLAN Conference on Programming Language Design and Implementation, pp. 562-576. 2020. <a href="#fnref:24" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:25">
<p>Devroye, L., Gravel, C., “<a href="https://arxiv.org/abs/1502.02539v6"><strong>Random variate generation using only finitely many unbiased, independently and identically distributed random bits</strong></a>”, arXiv:1502.02539v6 [cs.IT], 2020. <a href="#fnref:25" class="reversefootnote" role="doc-backlink">↩</a> <a href="#fnref:25:1" class="reversefootnote" role="doc-backlink">↩<sup>2</sup></a></p>
</li>
<li id="fn:26">
<p>Oberhoff, Sebastian, “<a href="https://dc.uwm.edu/etd/1888"><strong>Exact Sampling and Prefix Distributions</strong></a>”, <em>Theses and Dissertations</em>, University of Wisconsin Milwaukee, 2018. <a href="#fnref:26" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:28">
<p>Knuth, Donald E. and Andrew Chi-Chih Yao. “The complexity of nonuniform random number generation”, in <em>Algorithms and Complexity: New Directions and Recent Results</em>, 1976. <a href="#fnref:28" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:29">
<p>Feras A. Saad, Cameron E. Freer, Martin C. Rinard, and Vikash K. Mansinghka, “<a href="https://arxiv.org/abs/2001.04555"><strong>Optimal Approximate Sampling From Discrete Probability Distributions</strong></a>”, arXiv:2001.04555 [cs.DS], also in Proc. ACM Program. Lang. 4, POPL, Article 36 (January 2020), 33 pages. <a href="#fnref:29" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
<li id="fn:30">
<p>Baccetti, Valentina, and Matt Visser. “Infinite Shannon entropy.” Journal of Statistical Mechanics: Theory and Experiment 2013, no. 04 (2013): P04010, also in arXiv:1212.5630. <a href="#fnref:30" class="reversefootnote" role="doc-backlink">↩</a></p>
</li>
</ol>
</div>
</div><nav id="navigation"><ul>
<li><a href="/">Back to start site.</a>
<li><a href="https://github.com/peteroupc/peteroupc.github.io">This site's repository (source code)</a>
<li><a href="https://github.com/peteroupc/peteroupc.github.io/issues">Post an issue or comment</a></ul>
<div class="noprint">
<p>
<a href="//x.com/intent/post?text=More+Random+Sampling+Methods&url=https%3A%2F%2Fpeteroupc.github.io%2Frandomnotes.html">Share via X (Twitter)</a>, <a href="//www.facebook.com/sharer/sharer.php?u=https%3A%2F%2Fpeteroupc.github.io%2Frandomnotes.html">Share via Facebook</a>, <a href="//news.ycombinator.com/submitlink?u=https%3A%2F%2Fpeteroupc.github.io%2Frandomnotes.html&t=More+Random+Sampling+Methods">Share via Hacker News</a>, <a href="//www.reddit.com/submit?url=https%3A%2F%2Fpeteroupc.github.io%2Frandomnotes.html&title=More+Random+Sampling+Methods">Share via Reddit</a>
</p>
</div>
<p style='font-size:120%;font-weight:bold'><a href='https://peteroupc.github.io/randomnotes.pdf'>Download a PDF of this page</a></p><p>Of interest to the computer graphics and music community:</p><ul><li><a href='https://peteroupc.github.io/graphics.html'><b>Graphics for classic-style game development</b></a>.<li><a href='https://peteroupc.github.io/graphics.html#Building_a_Public_Domain_music_synthesis_library_and_instrument_banks'><b>MIDI music synthesis for classic-style games and apps</b></a>.<li><a href='https://peteroupc.github.io/classic-wallpaper'><b>Tileable wallpapers with limited colors and resolution</b></a>.<li><a href='https://peteroupc.github.io/graphicsapi.html'><b>Lean programming interfaces for classic graphics</b></a>.<li><a href='https://peteroupc.github.io/classic-wallpaper/docs/uielements.html'><b>Traditional user-interface graphics</b></a>.<li><a href='https://peteroupc.github.io/classic-wallpaper/docs/pixeltovector.html'><b>Converting images to vector graphics</b></a>.</ul><p>Open questions for math and probability experts:</p><ul><li><a href='https://peteroupc.github.io/requests.html'><b>Requests and open questions</b></a>.<li><a href='https://peteroupc.github.io/bernreq.html'><b>Open questions on the Bernoulli factory problem (the new-coins-from-old problem)</b></a>.<li><a href='https://peteroupc.github.io/requestsother.html'><b>Other open questions on probability</b></a>.<li><a href='https://peteroupc.github.io/sampling.html'><b>The sampling problem</b></a>.</ul></nav></body></html>