Generating Random Numbers: what changed

Revision 1 to 2. Removed lines have a minus sign; added lines have a plus sign.

Old	New		Line
1		-	~~<figure>~~
2		-	~~<img src="fpc_source_logo.png" title="fpc_source_logo.png" />~~
3		-	~~<figcaption>fpc_source_logo.png</figcaption>~~
4		-	~~</figure>~~
	1	+	![](/wiki-images/8/84/Fpc_source_logo.png)
5	2
6	3		Random numbers are important resources for a diverse field of applications, including scientific analysis, technology, medicine, economy, education, game development and visualization. They play a key role in numeric simulation.
7	4
8	5		Algorithm-generated random numbers are pseudo-random numbers. They belong to a (large) set of repeating numbers, whose sequence is impossible or at least difficult to predict. Unlike Delphi, that uses a linear congruential generator.(See [[Delphi compatible LCG Random]]), Free Pascal uses a MersenneTwister algorithm for its standard [[doc:rtl/system/random.html\|random]] function as defined in [[RTL]]. Before its first use, FPC's random number generator has to be initialized with a single call of the [[doc:rtl/system/randomize.html\|randomize]] function, which sets the seed of the generator. Preferably this is done in the launch phase of the program.
9	6
10	7		Alternatively, on Unix- and Linux-based systems, the virtual devices [[Dev random\|/dev/random]] and `/dev/urandom` are available. They generate (pseudo) random numbers based on hardware.
11	8
12	9		A third option is to use random numbers from external sources, either from specialised hardware devices or from public sources, e.g. based on radioactive decay data.
13	10
14	11		## Uniform Distribution
15	12
16	13		The continuous uniform distribution (also referred to as rectangular distribution) represents a family of symmetric probability distributions. Here, for each member of the family all intervals of the same length on the distribution's support are equally probable.
17	14
18	15		The standard [[RTL]] function [[doc:rtl/system/random.html\|random]] generates random numbers that fulfill a uniform distribution. If called without parameter `random` delivers a floating point pseudorandom number in the interval \[0, 1), i.e. 0 \<= result \< 1. if `random` is called with a longint argument L it delivers a longint random in the interval \[0, L).
19	16
20	17		A further set of uniformly distributed random number generators is presented in [[Marsaglia&\|Marsaglia's pseudo random number generators.]]\
21		-	\
	18	+
22	19		Uniformly distributed random numbers are not useful for every application. In order to create random numbers of other distributions special algorithms are necessary:
23	20
24	21		## Normal (Gaussian) Distribution
25	22
26	23		One of the more common algorithms to produce normally distributed random numbers from uniformly distributed random numbers is the [Box-Müller approach](http://en.wikipedia.org/wiki/Box–Muller_transform). The following function calculates Gaussian-distributed random numbers:
27	24
28	25		``` pascal
29	26		function rnorm (mean, sd: real): real;
30	27		{Calculates Gaussian random numbers according to the Box-Müller approach}
31	28		var
32	29		u1, u2: real;
33	30		begin
34	31		u1 := random;
35	32		u2 := random;
36	33		rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
37	34		end;
38	35		```
39	36
40	37		The same algorithm is used by the [[doc:rtl/math/randg.html\|randg]] randg function from the [[RTL]] [[doc:rtl/math/index.html\|math]] unit:
41	38
42	39		``` pascal
43	40		function randg(mean,stddev: float): float;
44	41		```
45	42
46	43		## Exponential Distribution
47	44
48	45		An exponential distribution occurs frequently in real-world problems. A classical example is the distribution of waiting times between independent Poisson-random events, e.g. the radioactive decay of nuclei \[Press et al. 1989\].
49	46
50	47		The following function delivers a single real random number out of an exponential distribution. Rate is the inverse of the mean, and the constant RESOLUTION determines the granularity of generated random numbers.
51	48
52	49		``` pascal
53	50		function randomExp(a, rate: real): real;
54	51		const
55	52		RESOLUTION = 1000;
56	53		var
57	54		unif: real;
58	55		begin
59	56		if rate = 0 then
60	57		randomExp := NaN
61	58		else
62	59		begin
63	60		repeat
64	61		unif := random(RESOLUTION) / RESOLUTION;
65	62		until unif <> 0;
66	63		randomExp := a - rate * ln(unif);
67	64		end;
68	65		end;
69	66		```
70	67
71	68		## Gamma Distribution
72	69
73	70		The gamma distribution is a two-parameter family of continuous random distributions. It is a generalization of both the exponential distribution and the Erlang distribution. Possible applications of the gamma distribution include modelling and simulation of waiting lines, or queues, and actuarial science.
74	71
75	72		The following function delivers a single real random number out of a gamma distribution. The shape of the distribution is defined by the parameters a, b and c. The function makes use of the function randomExp as defined above.
76	73
77	74		``` pascal
78	75		function randomGamma(a, b, c: real): real;
79	76		const
80	77		RESOLUTION = 1000;
81	78		T = 4.5;
82	79		D = 1 + ln(T);
83	80		var
84	81		unif: real;
85	82		A2, B2, C2, Q, p, y: real;
86	83		p1, p2, v, w, z: real;
87	84		found: boolean;
88	85		begin
89	86		A2 := 1 / sqrt(2 * c - 1);
90	87		B2 := c - ln(4);
91	88		Q := c + 1 / A2;
92	89		C2 := 1 + c / exp(1);
93	90		found := False;
94	91		if c < 1 then
95	92		begin
96	93		repeat
97	94		repeat
98	95		unif := random(RESOLUTION) / RESOLUTION;
99	96		until unif > 0;
100	97		p := C2 * unif;
101	98		if p > 1 then
102	99		begin
103	100		repeat
104	101		unif := random(RESOLUTION) / RESOLUTION;
105	102		until unif > 0;
106	103		y := -ln((C2 - p) / c);
107	104		if unif <= power(y, c - 1) then
108	105		begin
109	106		randomGamma := a + b * y;
110	107		found := True;
111	108		end;
112	109		end
113	110		else
114	111		begin
115	112		y := power(p, 1 / c);
116	113		if unif <= exp(-y) then
117	114		begin
118	115		randomGamma := a + b * y;
119	116		found := True;
120	117		end;
121	118		end;
122	119		until found;
123	120		end
124	121		else if c = 1 then
125	122		{ Gamma distribution becomes exponential distribution, if c = 1 }
126	123		begin
127	124		randomGamma := randomExp(a, b);
128	125		end
129	126		else
130	127		begin
131	128		repeat
132	129		repeat
133	130		p1 := random(RESOLUTION) / RESOLUTION;
134	131		until p1 > 0;
135	132		repeat
136	133		p2 := random(RESOLUTION) / RESOLUTION;
137	134		until p2 > 0;
138	135		v := A2 * ln(p1 / (1 - p1));
139	136		y := c * exp(v);
140	137		z := p1 * p1 * p2;
141	138		w := B2 + Q * v - y;
142	139		if (w + D - T * z >= 0) or (w >= ln(z)) then
143	140		begin
144	141		randomGamma := a + b * y;
145	142		found := True;
146	143		end;
147	144		until found;
148	145		end;
149	146		end;
150	147		```
151	148
152	149		## Erlang Distribution
153	150
154	151		The Erlang distribution is a two parameter family of continuous probability distributions. It is a generalization of the exponential distribution and a special case of the gamma distribution, where c is an integer. The Erlang distribution has been first described by Agner Krarup Erlang in order to model the time interval between telephone calls. It is used for queuing theory and for simulating waiting lines.
155	152
156	153		``` pascal
157	154		function randomErlang(mean: real; k: integer): real;
158	155		const
159	156		RESOLUTION = 1000;
160	157		var
161	158		i: integer;
162	159		unif, prod: real;
163	160		begin
164	161		if (mean <= 0) or (k < 1) then
165	162		randomErlang := NaN
166	163		else
167	164		begin
168	165		prod := 1;
169	166		for i := 1 to k do
170	167		begin
171	168		repeat
172	169		unif := random(RESOLUTION) / RESOLUTION;
173	170		until unif <> 0;
174	171		prod := prod * unif;
175	172		end;
176	173		randomErlang := -mean * ln(prod);
177	174		end;
178	175		end;
179	176		```
180	177
181	178		## Poisson Distribution
182	179
183	180		The Poisson distribution applies to integer values. It represents the probability of k successes, when the probability of a success in each trial is small and the rate of occurrence (the mean value) is constant.
184	181
185	182		``` pascal
186	183		function randomPoisson(mean: integer): integer;
187	184		{ Generator for Poisson distribution (Donald Knuth's algorithm) }
188	185		const
189	186		RESOLUTION = 1000;
190	187		var
191	188		k: integer;
192	189		b, l: real;
193	190		begin
194	191		assert(mean > 0, 'mean < 1');
195	192		k := 0;
196	193		b := 1;
197	194		l := exp(-mean);
198	195		while b > l do
199	196		begin
200	197		k := k + 1;
201	198		b := b * random(RESOLUTION) / RESOLUTION;
202	199		end;
203	200		randomPoisson := k - 1;
204	201		end;
205	202		```
206	203
207	204		## t Distribution
208	205
209	206		The t distribution (also referred to a Student's t distribution, since it was published by William Sealy Gosset in 1908 under the pseudonym Student) is a continuous probability distribution. Its shape is defined by one parameter, the degrees of freedom (df). In statistics, many estimators are t distributed. Therefore, Student's t-distribution plays a major role in a number of widely used statistical analyses, including Student's t-test for assessing the statistical significance of the difference between two sample means, the construction of confidence intervals for the difference between two population means, and in linear regression analysis. The t-distribution also arises in Bayesian analysis of data from a normal family.
210	207
211	208		The following algorithm depends on the [[RTL]] function [[doc:rtl/system/random.html\|random]] and on the randomChisq function
212	209
213	210		``` pascal
214	211		function randomT(df: integer): real;
215	212		{ Generator for Student's t distribution }
216	213		begin
217	214		if df < 1 then
218	215		randomT := NaN
219	216		else begin
220	217		randomT := randg(0, 1) / sqrt(randomChisq(df) / df);
221	218		end;
222	219		end;
223	220		```
224	221
225	222		## Chi Squared Distribution
226	223
227	224		The chi squared distribution is a continuous distribution of random numbers with df degrees of freedom. It is the distribution of a sum of the squares of df independent standard normal random variables. The chi squared distribution has numerous applications in inferential statistics, e.g. in estimating variances and for chi-squared tests. It is a special gamma distribution with c = df/ 2 and b = 2. Therefore the following function depends on the function randomGamma.
228	225
229	226		``` pascal
230	227		function randomChisq(df: integer): real;
231	228		begin
232	229		if df < 1 then
233	230		randomChisq := NaN
234	231		else
235	232		randomChisq := randomGamma(0, 2, 0.5 * df);
236	233		end;
237	234		```
238	235
239	236		## F Distribution
240	237
241	238		The F distribution, also referred to as Fisher-Snedecor distribution, is a continuous probability distribution. It is used for F Test and ANOVA. It has two degrees of freedom that serve as shape parameters v and w and that are positive integers. The following function randomF makes use of randomChisq.
242	239
243	240		``` pascal
244	241		function randomF(v, w: integer): real;
245	242		begin
246	243		if (v < 1) or (w < 1) then
247	244		randomF := NaN
248	245		else
249	246		randomF := randomChisq(v) / v / (randomChisq(w) / w);
250	247		end;
251	248		```
252	249
253	250		## See also
254	251
255	252		- [[Dev random]]
256	253		- [[Functions for descriptive statistics]]
257	254		- [[Samples and permutations]]
258	255		- [[Marsaglia&\|Marsaglia's pseudo random number generators]]
259	256		- [[A simple implementation of the Mersenne twister]]
260	257		- [[Delphi compatible LCG Random]]
261	258		- [[EverettRandom]]
262	259
263	260		## References
264	261
265	262		1. [G. E. P. Box and Mervin E. Muller, A Note on the Generation of Random Normal Deviates, The Annals of Mathematical Statistics (1958), Vol. 29, No. 2 pp. 610–611](http://projecteuclid.org/DPubS/Repository/1.0/Disseminate?view=body&id=pdf_1&handle=euclid.aoms/1177706645)
266	263		2. Dietrich, J. W. (2002). [Der Hypophysen-Schilddrüsen-Regelkreis](http://openlibrary.org/books/OL24586469M/Der_Hypophysen-Schilddrüsen-Regelkreis). Berlin, Germany: Logos-Verlag Berlin. ISBN 978-3-89722-850-4. OCLC 50451543.
267	264		3. Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling (1989). [Numerical Recipes in Pascal](http://openlibrary.org/works/OL16807779W/). The Art of Scientific Computing, Cambridge University Press, ISBN 0-521-37516-9.
268	265		4. Richard Saucier, [Computer Generation of Statistical Distributions](http://ftp.arl.mil/random/random.pdf), ARL-TR-2168, US Army Research Laboratory, Aberdeen Proving Ground, MD, 21005-5068, March 2000.
269	266		5. R.U. Seydel, Generating Random Numbers with Specified Distributions. In: Tools for Computational Finance, Universitext, [DOI 10.1007/978-1-4471-2993-6_2](http://dx.doi.org/10.1007/978-1-4471-2993-6_2), © Springer-Verlag London Limited 2012
270	267		6. Christian Walck, Hand-book on STATISTICAL DISTRIBUTIONS for experimentalists, Internal Report SUF–PFY/96–01, University of Stockholm 2007
271		-
272		-	\
273	268
274	269		---
275	270		From the Free Pascal wiki, page "Generating Random Numbers", as it stood on 2020-11-29 after a change by Jwdietrich, out of the wiki's dump of 2024-05-26. Copyright stays with the author who wrote it, according to the old wiki's [[Licensing]] page. If you have concerns about this material, [message the admin](/wiki-contact).

Back to history