Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

randraw.m @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago
Initial import--MML version from SOLEIL@2013
File size: 205.8 KB

Line
1	function varargout = randraw(distribName, distribParams, varargin)
2	%RANDRAW - EFFICIENT RANDOM VARIATES GENERATOR
3	%
4	% See alphabetical list of the supported distributions below (over 50 distributions)
5	%
6	% 1) randraw
7	% presents general help.
8	% 2) randraw( distribName )
9	% presents help for the specific distribution defined
10	% by usage string distribName (see table below).
11	% 3) Y = randraw( distribName, distribParams, sampleSize );
12	% returns array Y of size = sampleSize of random variates from distribName
13	% distribution with parameters distribParams
14	%
15	% ALPHABETICAL LIST OF THE SUPPORTED DISTRIBUTIONS:
16	% ____________________________________________________________________
17	% \| DISTRIBUTION NAME \| USAGE STRING \|
18	% \|___________________________________________\|________________________\|
19	% \| Alpha \| 'alpha' \|
20	% \| Anglit \| 'anglit' \|
21	% \| Antilognormal \| 'lognorm' \|
22	% \| Arcsin \| 'arcsin' \|
23	% \| Bernoulli \| 'bern' \|
24	% \| Bessel \| 'bessel' \|
25	% \| Beta \| 'beta' \|
26	% \| Binomial \| 'binom' \|
27	% \| Bradford \| 'bradford' \|
28	% \| Burr \| 'burr' \|
29	% \| Cauchy \| 'cauchy' \|
30	% \| Chi \| 'chi' \|
31	% \| Chi-Square (Non-Central) \| 'chisqnc' \|
32	% \| Chi-Square (Central) \| 'chisq' \|
33	% \| Cobb-Douglas \| 'lognorm' \|
34	% \| Cosine \| 'cosine' \|
35	% \| Double-Exponential \| 'laplace' \|
36	% \| Erlang \| 'erlang' \|
37	% \| Exponential \| 'exp' \|
38	% \| Extreme-Value \| 'extrval' \|
39	% \| F (Central) \| 'f' \|
40	% \| F (Non-Central) \| 'fnc' \|
41	% \| Fisher-Tippett \| 'extrval' \|
42	% \| Fisk \| 'fisk' \|
43	% \| Frechet \| 'frechet' \|
44	% \| Furry \| 'furry' \|
45	% \| Gamma \| 'gamma' \|
46	% \| Generalized Inverse Gaussian \| 'gig' \|
47	% \| Generalized Hyperbolic \| 'gh' \|
48	% \| Geometric \| 'geom' \|
49	% \| Gompertz \| 'gompertz' \|
50	% \| Gumbel \| 'gumbel' \|
51	% \| Half-Cosine \| 'hcos' \|
52	% \| Hyperbolic Secant \| 'hsec' \|
53	% \| Hypergeometric \| 'hypergeom' \|
54	% \| Inverse Gaussian \| 'ig' \|
55	% \| Laplace \| 'laplace' \|
56	% \| Logistic \| 'logistic' \|
57	% \| Lognormal \| 'lognorm' \|
58	% \| Lomax \| 'lomax' \|
59	% \| Lorentz \| 'lorentz' \|
60	% \| Maxwell \| 'maxwell' \|
61	% \| Negative Binomial \| 'negbinom' \|
62	% \| Normal \| 'norm' \|
63	% \| Normal-Inverse-Gaussian (NIG) \| 'nig' \|
64	% \| Pareto \| 'pareto' \|
65	% \| Pareto2 \| 'pareto2' \|
66	% \| Pascal \| 'pascal' \|
67	% \| Planck \| 'planck' \|
68	% \| Poisson \| 'po' \|
69	% \| Quadratic \| 'quadr' \|
70	% \| Rademacher \| 'rademacher' \|
71	% \| Rayleigh \| 'rayl' \|
72	% \| Semicircle \| 'semicirc' \|
73	% \| Skellam \| 'skellam' \|
74	% \| Student's-t \| 't' \|
75	% \| Triangular \| 'tri' \|
76	% \| Truncated Normal \| 'normaltrunc' \|
77	% \| Tukey-Lambda \| 'tukeylambda' \|
78	% \| U-shape \| 'u' \|
79	% \| Uniform (continuous) \| 'uniform' \|
80	% \| Von Mises \| 'vonmises' \|
81	% \| Wald \| 'wald' \|
82	% \| Weibull \| 'weibull' \|
83	% \| Wigner Semicircle \| 'wigner' \|
84	% \| Yule \| 'yule' \|
85	% \| Zeta \| 'zeta' \|
86	% \| Zipf \| 'zipf' \|
87	% \|___________________________________________\|________________________\|
88
89	% Version 1.5 - December 2005
90	% 'true' and 'false' functions were replased by ones and zeros to support Matlab releases
91	% below 6.5
92	% Version 1.4 - September 2005 -
93	% Bugs fix:
94	% 1) GAMMA distribution (thanks to Earl Lawrence):
95	% special case for a<1
96	% 2) GIG distribution (thanks to Panagiotis Braimakis):
97	% typo in help
98	% code adjustment to overcome possible computational overflows
99	% 3) CHI SQUARE distribution
100	% typo in help
101	% Version 1.3 - July 2005 -
102	% Bug fix:
103	% Typo in GIG distribution generation:
104	% should be 'out' instead of 'x' in lines 1852 and 1858
105	% Version 1.2 - May 2005 -
106	% Bugs fix:
107	% 1) Poisson distribution did not work for lambda < 21.4. Typo ( ti instead of t )
108	% 2) GIG distribution: support to chi=0 or psi=0 cases
109	% 3) Beta distribution: column sampleSize
110	% 4) Cauchy distribution: typo in example
111	% 5) Chi distribution: typo in example
112	% 6) Non-central F distribution: number of input parameters
113	% 7) INVERSE GAUSSIAN (IG) distribution: typo in example
114	%
115	% Version 1.1 - April 2005 - Bug fix: Generation from binomial distribution using only 'binomial'
116	% usage string was changed to 'binom' ( 'binomial' works too ).
117	% Version 1.0 - March 2005 - Initial version
118	% Alex Bar Guy & Alexander Podgaetsky
119	% alex@wavion.co.il
120	%
121
122	% Reference links:
123	% 1) http://mathworld.wolfram.com/topics/StatisticalDistributions.html
124	% 2) http://en.wikipedia.org/wiki/Category:Probability_distributions
125	% 3) http://www.brighton-webs.co.uk/index.asp
126	% 4) http://www.jstatsoft.org/v11/i03/v11i03.pdf
127	% 5) http://www.quantlet.com/mdstat/scripts/csa/html/node236.html
128
129	% These programs are distributed in the hope that they will be useful,
130	% but WITHOUT ANY WARRANTY; without even the implied warranty of
131	% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
132
133	% Any comments and suggestions please send to:
134	% alex@wavion.co.il
135
136	funcName = mfilename;
137
138	if nargin == 0
139	help(funcName);
140	return;
141	elseif nargin == 1
142	runMode = 'distribHelp';
143	elseif nargin == 2
144	runMode = 'genRun';
145	sampleSize = [1 1];
146	else
147	runMode = 'genRun';
148	sampleSize = [varargin{1:end}];
149	end
150
151	distribNameInner = lower( distribName( ~isspace( distribName ) ) );
152
153	if strcmp(runMode, 'distribHelp')
154	fid = fopen( [ funcName '.m' ], 'r' );
155	printHelpFlag = 0;
156	while 1
157	tline = fgetl( fid );
158	if ~ischar( tline )
159	fprintf( '\n Unknown distribution name ''%s''.\n', distribName );
160	break;
161	end
162	if ~isempty( strfind( tline, [ 'END ', distribNameInner,' HELP' ] ) )
163	printHelpFlag = 0;
164	break;
165	end
166	if printHelpFlag
167	startPosition = strfind( tline, ' % ' ) + 3;
168	printLine = tline( startPosition : end );
169	if ~strcmp( funcName, 'randraw' )
170	indxs = strfind( printLine, 'randraw' );
171	while ~isempty( indxs )
172	headLine = printLine( 1:indxs(1)-1 );
173	tailLine = printLine( indxs(1)+7:end );
174	printLine = [ headLine, funcName, tailLine ];
175	indxs = strfind( printLine, 'randraw' );
176	end
177	end
178	pause(0.02);
179	fprintf( '\n%s', printLine );
180	end
181	if ~isempty( strfind( tline, [ 'START ', distribNameInner,' HELP' ] ) )
182	printHelpFlag = 1;
183	end
184	end
185	fprintf( '\n\n' );
186	fclose( fid );
187	return;
188	end
189
190	if length(sampleSize) == 1
191	sampleSize = [ sampleSize, 1 ];
192	end
193
194	if strcmp(runMode, 'genRun')
195	runExample = 0;
196	plotFlag = 0;
197
198	dbclear if warning;
199	out = [];
200	if prod(sampleSize) > 0
201	switch lower( distribNameInner )
202	case {'alpha'}
203	% START alpha HELP
204	% THE ALPHA DISTRIBUTION
205	%
206	% pdf(y) = bnormpdf(a-b./y) ./ (y.^2normcdf(a)); y>0; a>0; b>0;
207	% cdf(y) = normcdf(a-b./y)/normcdf(a); y>0; a>0; b>0;
208	% where normpdf(x) = 1/sqrt(2pi) exp(-1/2*x.^2); is the standard normal PDF
209	% normcdf(x) = 0.5*(1+erf(y/sqrt(2))); is the standard normal CDF
210	%
211	% PARAMETERS:
212	% a - shape parameter (a>0)
213	% b - shape parameter (b>0)
214	%
215	% SUPPORT:
216	% y, y>0
217	%
218	% CLASS:
219	% Continuous skewed distributions
220	%
221	% USAGE:
222	% randraw('alpha', [], sampleSize) - generate sampleSize number
223	% of variates from Alpha distribution with shape parameters a and b;
224	% randraw('alpha') - help for Alpha distribution;
225	%
226	% EXAMPLES:
227	% 1. y = randraw('alpha', [1 2], [1 1e5]);
228	% 2. y = randraw('alpha', [2 3], 1, 1e5);
229	% 3. y = randraw('alpha', [10 50], 1e5 );
230	% 4. y = randraw('alpha', [20.5 30.5], [1e5 1] );
231	% 5. randraw('alpha');
232	% END alpha HELP
233
234	% References:
235	% 1. doc erf
236
237	checkParamsNum(funcName, 'Alpha', 'alpha', distribParams, [2]);
238	a = distribParams(1);
239	b = distribParams(2);
240	validateParam(funcName, 'Alpha', 'alpha', '[a, b]', 'a', a, {'> 0'});
241	validateParam(funcName, 'Alpha', 'alpha', '[a, b]', 'b', b, {'> 0'});
242
243	out = b ./ ( a - norminv(normcdf(a)*rand(sampleSize)) );
244
245	case {'anglit'}
246	% START anglit HELP
247	% THE ANGLIT DISTRIBUTION
248	%
249	% Standard form of anglit distribution:
250	% pdf(y) = sin(2*y+pi/2); -pi/4 <= y <= pi/4;
251	% cdf(y) = sin(y+pi/4).^2; -pi/4 <= y <= pi/4;
252	%
253	% Mean = Median = Mode = 0;
254	% Variance = (pi/4)^2 - 0.5;
255	%
256	% General form of anglit distribution:
257	% pdf(y) = sin(pi/2*(y-t)/s+pi/2); t-s <= y <= t+s; s>0
258	% cdf(y) = sin(pi/4*(y-t)/s+pi/4).^2; t-s <= y <= t+s; s>0
259	%
260	% Mean = Median = Mode = t;
261	% Variance = ???????;
262	%
263	% PARAMETERS:
264	% t - location
265	% s -scale; s>0
266	%
267	% SUPPORT:
268	% y, -pi/4 <= y <= pi.4 - standard Anglit distribution
269	% or
270	% y, t-s <= y <= t+s - generalized Anglit distribution
271	%
272	% CLASS:
273	% Continuous distributions
274	%
275	% USAGE:
276	% randraw('anglit', [], sampleSize) - generate sampleSize number
277	% of variates from standard Anglit distribution;
278	% randraw('anglit', [t, s], sampleSize) - generate sampleSize number
279	% of variates from generalized Anglit distribution
280	% with location parameter 't' and scale parameter 's';
281	% randraw('anglit') - help for Anglit distribution;
282	%
283	% EXAMPLES:
284	% 1. y = randraw('anglit', [], [1 1e5]);
285	% 2. y = randraw('anglit', [], 1, 1e5);
286	% 3. y = randraw('anglit', [], 1e5 );
287	% 4. y = randraw('anglit', [10 3], [1e5 1] );
288	% 5. randraw('anglit');
289	%
290	% END anglit HELP
291
292	checkParamsNum(funcName, 'Anglit', 'anglit', distribParams, [0, 2]);
293	if numel(distribParams)==2
294	t = distribParams(1);
295	s = distribParams(2);
296	validateParam(funcName, 'Anglit', 'anglit', '[t, s]', 's', s, {'> 0'});
297	else
298	t = 0;
299	s = pi/4;
300	end
301
302	out = t + s * (4/pi*asin(sqrt(rand(sampleSize)))-1);
303
304	case {'arcsin'}
305	% START arcsin HELP
306	% THE ARC-SINE DISTRIBUTION
307	%
308	% pdf(y) = 1 ./ (pisqrt(y.(1-y))); 0<y<1;
309	% cdf(y) = 2*asin(sqrt(y))/pi; 0<y<1;
310	%
311	% Mean = 0.5;
312	% Variance = 0.125;
313	%
314	% PARAMETERS:
315	% None
316	%
317	% SUPPORT:
318	% y, 0<y<1
319	%
320	% CLASS:
321	% Continuous symmetric distributions
322	% NOTES:
323	% The arc-sine distribution is a special case of the beta distribution
324	% with both parameters equal to 1/2. The generalized arc-sine distribution
325	% is the special case of the beta distribution where the two parameters sum
326	% to 1 but are not necessarily equal to 1/2.
327	%
328	% USAGE:
329	% randraw('arcsin', [], sampleSize) - generate sampleSize number
330	% of variates from the Arc-sine distribution;
331	% randraw('arcsin') - help for the Arc-sine distribution;
332	%
333	% EXAMPLES:
334	% 1. y = randraw('arcsin', [], [1 1e5]);
335	% 2. y = randraw('arcsin', [], 1, 1e5);
336	% 3. y = randraw('arcsin', [], 1e5 );
337	% 4. randraw('arcsin');
338	% SEE ALSO:
339	% U distribution
340	% END arcsin HELP
341
342	checkParamsNum(funcName, 'Arcsin', 'arcsin', distribParams, 0);
343	out = sin( rand(sampleSize)*pi/2 ).^2;
344
345	case {'bernoulli', 'bern'}
346	% START bernoulli HELP START bern HELP
347	% THE BERNOULLI DISTRIBUTION
348	%
349	% pdf(y) = p.^y .* (1-p).^(1-y);
350	% cdf(y) = (y==0)(1-p) + (y==1)1;
351	%
352	% PARAMETERS:
353	% p is a probability of success; (0<p<1)
354	%
355	% SUPPORT:
356	% y = [0 1];
357	%
358	% CLASS:
359	% Discrete distributions
360	%
361	% USAGE:
362	% randraw('bern', p, sampleSize) - generate sampleSize number
363	% of variates from the Bernoulli distribution with probability of success p
364	% randraw('bern') - help for the Bernoulli distribution;
365	%
366	% EXAMPLES:
367	% 1. y = randraw('bern', 0.5, [1 1e5]);
368	% 2. y = randraw('bern', 0.1, 1, 1e5);
369	% 3. y = randraw('bern', 0.9, 1e5 );
370	% 4. randraw('bern');
371	% END bernoulli HELP END bern HELP
372
373	checkParamsNum(funcName, 'Bernoulli', 'bernoulli', distribParams, 1);
374	validateParam(funcName, 'Bernoulli', 'bernoulli', 'p', 'p', distribParams(1), {'>=0','<=1'});
375
376	out = double( rand( sampleSize ) < distribParams );
377
378	case {'beta', 'powerfunction', 'powerfunc'}
379	% START beta HELP
380	% THE BETA DISTRIBUTION
381	%
382	% ( sometimes: Power Function distribution )
383	%
384	% Standard form of the Beta distribution:
385	% pdf(y) = y.^(a-1).*(1-y).^(b-1) / beta(a, b);
386	% cdf(y) = betainc(y,a,b), if (y>=0 & y<=1); 0, if x<0; 1, if x>1
387	%
388	% Mean = a/(a+b);
389	% Variance = (ab)/((a+b)^2(a+b+1));
390	%
391	% General form of the Beta distribution:
392	% pdf(y) = (y-m).^(a-1).(n-y).^(b-1) / (beta(a, b)(n-m)^(a+b-1));
393	% cdf(y) = betainc((y-m)/(n-m),a,b), if (y>=m & y<=n); 0, if x<m; 1, if x>n
394	%
395	% Mean = (na + mb)/(a+b);
396	% Variance = (ab)(n-m)^2/((a+b)^2*(a+b+1));
397	%
398	% PARAMETERS:
399	% a>0 - shape parameter
400	% b>0 - shape parameter
401	% m - location
402	% n -scale (upper bound); n>=m
403	%
404	% SUPPORT:
405	% y, 0<=y<=1 - standard beta distribution
406	% or
407	% y, m<=y<=n - generalized beta distribution
408	%
409	% CLASS:
410	% Continuous skewed distributions
411	%
412	% USAGE:
413	% randraw('beta', [a, b], sampleSize) - generate sampleSize number
414	% of variates from standard beta distribution with shape parameters
415	% 'a' and 'b'
416	% randraw('beta', [m, n, a, b], sampleSize) - generate sampleSize number
417	% of variates from generalized beta distribution on the interval [m, n]
418	% with shape parameters 'a' and 'b';
419	% randraw('beta') - help for the Beta distribution;
420	% EXAMPLES:
421	% 1. y = randraw('beta', [0.2 0.9], [1 1e5]);
422	% 2. y = randraw('beta', [0.6 3.2], 1, 1e5);
423	% 3. y = randraw('beta', [-10 20 3.1 6.2], 1e5 );
424	% 4. y = randraw('beta', [3 4 5.3 0.7], [1e5 1] );
425	% 5. randraw('beta');
426	% END beta HELP
427
428	% Refernce:
429	% Dagpunar, John.
430	% Principles of Random Variate Generation.
431	% Oxford University Press, 1988.
432	%
433	% max_ab < 0.5 Joehnk's algorithm
434	% 1 < min_ab Cheng's algortihm BB
435	% min_ab <= 1 <= max_ab Atkinson's switching algorithm
436	% 0.5<= max_ab < 1 Atkinson's switching algorithm
437
438
439	checkParamsNum(funcName, 'Beta', 'beta', distribParams, [2, 4]);
440
441	if numel(distribParams) == 2
442	a = distribParams(1);
443	b = distribParams(2);
444	m = 0;
445	n = 1;
446	validateParam(funcName, 'Beta', 'beta', '[a, b]', 'a', a, {'> 0'});
447	validateParam(funcName, 'Beta', 'beta', '[a, b]', 'b', b, {'> 0'});
448	else
449	m = distribParams(1);
450	n = distribParams(2);
451	a = distribParams(3);
452	b = distribParams(4);
453	validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'n-m', n-m, {'>= 0'});
454	validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'a', a, {'> 0'});
455	validateParam(funcName, 'Beta', 'beta', '[m, n, a, b]', 'b', b, {'> 0'});
456
457	end
458
459	sampleSizeIn = sampleSize;
460	sampleSize = [ 1, prod( sampleSizeIn ) ];
461
462	max_ab = max( a, b );
463	min_ab = min( a, b );
464	if max_ab < 0.5
465	% Use log(u1^a) and log(u2^b), rather than a and b, to avoid
466	% underflow for very small a or b.
467
468	loga = log(rand( sampleSize ))/a;
469	logb = log(rand( sampleSize ))/b;
470	logsum = (loga>logb).*(loga + log(1+ exp(logb-loga))) + ...
471	(loga<=logb).*(logb + log(1+ exp(loga-logb)));
472	out = exp(loga - logsum);
473
474	indxs = find( logsum > 0);
475
476	while ~isempty( indxs )
477	indxsSize = size( indxs );
478	loga = log(rand( indxsSize ))/a;
479	logb = log(rand( indxsSize ))/b;
480	logsum = (loga>logb).*(loga + log(1+ exp(logb-loga))) + ...
481	(loga<=logb).*(logb + log(1+ exp(loga-logb)));
482
483	l = (logsum <= 0);
484	out( indxs( l ) ) = exp(loga(l) - logsum(l));
485	indxs = indxs( ~l );
486	end
487
488	elseif min_ab > 1
489	% Algorithm BB
490
491	sum_ab = a + b;
492	lambda = sqrt((sum_ab-2)/(2ab-sum_ab));
493	c = min_ab+1/lambda;
494
495	u1 = rand( sampleSize );
496	u2 = rand( sampleSize );
497	v = lambda*log(u1./(1-u1));
498	z = u1.u1.u2;
499	clear('u1'); clear('u2');
500	w = min_ab*exp(v);
501	r = c*v-1.38629436112;
502	clear('v');
503	s = min_ab+r-w;
504	if a == min_ab
505	out = w./(max_ab+w);
506	else
507	out = max_ab./(max_ab+w);
508	end
509
510	t = log(z);
511	indxs = find( (s+2.609438 < 5z) & (r+sum_ablog(sum_ab./(max_ab+w)) < t) );
512
513	clear('v');
514	clear('z');
515	clear('w');
516	clear('r');
517	while ~isempty( indxs )
518	indxsSize = size( indxs );
519
520	u1 = rand( indxsSize );
521	u2 = rand( indxsSize );
522	v = lambda*log(u1./(1-u1));
523	z = u1.u1.u2;
524	clear('u1'); clear('u2');
525	w = min_ab*exp(v);
526	r = c*v-1.38629436112;
527	clear('v');
528	s = min_ab+r-w;
529	t = log(z);
530
531	l = (s+2.609438 >= 5z) \| (r+sum_ablog(sum_ab./(max_ab+w)) >= t);
532	if a == min_ab
533	out( indxs( l ) ) = w(l)./(max_ab+w(l));
534	else
535	out( indxs( l ) ) = max_ab./(max_ab+w(l));
536	end
537	indxs = indxs( ~l );
538
539	end
540
541	elseif min_ab < 1 & max_ab > 1
542	% Atkinson's switching method
543
544	t = (1-min_ab)/(1+max_ab - min_ab);
545	r = max_abt/(max_abt + min_ab*(1-t)^max_ab);
546
547	u1 = rand( sampleSize );
548	w = zeros( sampleSize );
549	l = u1 < r;
550	w( l ) = t*(u1( l )/r).^(1/min_ab);
551	l = ~l;
552	w( l ) = 1- (1-t)*((1-u1(l))/(1-r)).^(1/max_ab);
553	if a == min_ab
554	out = w;
555	else
556	out = 1 - w;
557	end
558	u2 = rand( sampleSize );
559
560	indxs1 = find(u1 < r);
561	indxs2 = find(u1 >= r);
562	clear('u1');
563	indxs = [ indxs1( log(u2(indxs1)) >= (max_ab-1)*log(1-w(indxs1)) ), ...
564	indxs2( log(u2(indxs2)) >= (min_ab-1)*log(w(indxs2)/t) ) ];
565
566	clear('u1');
567	clear('u2');
568	while ~isempty( indxs )
569	indxsSize = size( indxs );
570	u1 = rand( indxsSize );
571	w = zeros( indxsSize );
572	l = u1 < r;
573	w( l ) = t*(u1( l )/r).^(1/min_ab);
574	l = ~l;
575	w( l ) = 1- (1-t)*((1-u1(l))/(1-r)).^(1/max_ab);
576
577	u2 = rand( indxsSize );
578
579	indxs1 = find(u1 < r);
580	indxs2 = find(u1 >= r);
581	clear('u1');
582	l = logical( zeros( indxsSize ) );
583	l( [ indxs1( log(u2(indxs1)) < (max_ab-1)*log(1-w(indxs1)) ), ...
584	indxs2( log(u2(indxs2)) < (min_ab-1)*log(w(indxs2)/t) ) ] ) = 1;
585
586	clear('u1');
587	clear('u2');
588	if a == min_ab
589	out( indxs(l) ) = w(l);
590	else
591	out( indxs(l) ) = 1 - w(l);
592	end
593	indxs = indxs( ~l );
594	end
595
596	else
597	% Atkinson's Algorithm
598
599	if min_ab == 1
600	t = 0.5;
601	r = 0.5;
602	else
603	t = 1/(1+sqrt(max_ab(1-max_ab)/(min_ab(1-min_ab))));
604	r = max_abt / (max_abt + min_ab*(1-t));
605	end
606
607	u1 = rand( sampleSize );
608	out = zeros( sampleSize );
609	w = zeros( sampleSize );
610	l1 = u1 < r;
611	w(l1) = t*(u1(l1)/r).^(1/min_ab);
612	l2 = u1 >= r;
613	w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab);
614	if a == min_ab
615	out = w;
616	else
617	out = 1 - w;
618	end
619
620	u2 = rand( sampleSize );
621	indxs1 = find(l1);
622	indxs2 = find(l2);
623	indxs = [ indxs1( log(u2(l1)) >= (max_ab -1)*log((1-w(l1))/(1-t)) ), ...
624	indxs2( log(u2(l2)) >= (min_ab -1) * log(w(l2)/t) ) ];
625	clear('u2');
626
627	while ~isempty( indxs )
628	indxsSize = size( indxs );
629	u1 = rand( indxsSize );
630	w = zeros( indxsSize );
631
632	l1 = u1 < r;
633	w(l1) = t*(u1(l1)/r).^(1/min_ab);
634	l2 = u1 >= r;
635	w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab);
636
637	u2 = rand( indxsSize );
638
639	indxs1 = find(l1);
640	indxs2 = find(l2);
641	clear('u1');
642	l = logical( zeros( indxsSize ) );
643
644	l( [ indxs1(log(u2(l1)) < (max_ab -1)*log((1-w(l1))/(1-t))), ...
645	indxs2(log(u2(l2))< (min_ab -1) * log(w(l2)/t)) ] ) = 1;
646
647	if a == min_ab
648	out(indxs(l)) = w(l);
649	else
650	out(indxs(l)) = 1 - w(l);
651	end
652	indxs = indxs( ~l );
653	end
654
655	end
656
657	out = m + (n-m) * out;
658
659	reshape( out, sampleSizeIn );
660
661	case {'bessel'}
662	% START bessel HELP
663	% THE BESSEL DISTRIBUTION
664	%
665	% Bessel distribution arises in the theory of stochastic processes.
666	% Bessel(nu,a) is a discrete distribution on the non-negative integers with
667	% parameters nu > -1 and a > 0.
668	%
669	% pdf(y) = (a/2).^(2y+nu) ./ (besseli(nu,a).factorial(y).*gamma(y+nu+1));
670	%
671	% PARAMETERS:
672	% nu > -1, a > 0
673	% SUPPORT:
674	% y = 0, 1, 2, 3, ...
675	% CLASS:
676	% Discrete distributions
677	%
678	% USAGE:
679	% randraw('bessel', [nu, a], sampleSize) - generate sampleSize number
680	% of variates from the Bessel distribution with parameters
681	% 'nu' and 'a'
682	% randraw('bessel') - help for the Bessel distribution;
683	% EXAMPLES:
684	% 1. y = randraw('bessel', [2 0.9], [1 1e5]);
685	% 2. y = randraw('bessel', [0.6 3.2], 1, 1e5);
686	% 3. y = randraw('bessel', [-0.2 8.1], 1e5 );
687	% 4. y = randraw('bessel', [4 5.3], [1e5 1] );
688	% 5. randraw('bessel');
689	% END bessel HELP
690
691	% Method:
692	%
693	% We implemented Condensed Table-Lookup method suggested in
694	% George Marsaglia, "Fast Generation Of Discrete Random Variables,"
695	% Journal of Statistical Software, July 2004, Volume 11, Issue 4
696	%
697	% Reference:
698	% L. Devroye, "Simulating Bessel random variables,"
699	% Statistics and Probability Letters, vol. 57, pp. 249-257, 2002.
700	%
701
702	checkParamsNum(funcName, 'Bessel', 'bessel', distribParams, [2]);
703
704	nu = distribParams(1);
705	a = distribParams(2);
706
707	validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'nu', nu, {'> -1'});
708	validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'a', a, {'> 0'});
709
710	% mu = 0.5abesseli(nu+1,a)/besseli(nu,a);
711	% chi2 = mu + 0.25a^2besseli(nu+1,a)/besseli(nu,a)*...
712	% (besseli(nu+2,a)/besseli(nu+1,a)-besseli(nu+1,a)/besseli(nu,a));
713
714	besseliNuA = besseli(nu, a);
715
716	proceed = 1;
717	if ~isfinite( besseliNuA )
718	warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: '];
719	warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns Inf.'];
720	warnStr{3} = ['Unable to proceed, return zeros ...'];
721	warning('%s\n %s\n %s',warnStr{1},warnStr{2},warnStr{3});
722
723	%warning([upper(funcName), ' - Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns Inf. Unable to proceed, return zeros ...']);
724	out = zeros( sampleSize );
725	proceed = 0;
726	end
727	if besseliNuA == 0
728	warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: '];
729	warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns 0.'];
730	warnStr{3} = ['Unable to proceed, return zeros ...'];
731	warning('%s\n %s\n %s',warnStr{1},warnStr{2},warnStr{3});
732	%warning([upper(funcName), '- Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns 0. Unable to proceed, return zeros ...']);
733	out = zeros( sampleSize );
734	proceed = 0;
735	end
736
737	if proceed
738	p0 = exp( nu*log(a/2) - gammaln(nu+1) ) / besseliNuA;
739	if p0 >= 5e-10
740	t = p0;
741
742	aa = (a/2)^2;
743	nu1 = nu+1;
744	i = 1;
745	while t*2147483648 > 1
746	t = t * aa/((i)*(i+nu));
747	i = i + 1;
748	end
749	sizeP = i-1;
750	offset = 0;
751
752	P = round( 2^30p0cumprod([1, aa./((1:sizeP-1).*((1:sizeP-1)+nu))] ) );
753
754	else % if p0 >= 5e-10
755	m = floor(0.5*(sqrt(a^2+nu^2)-nu));
756	pm = exp( (2m+nu)log(a/2) - log(besseliNuA) - ...
757	gammaln(m+1) - gammaln(m+nu+1) );
758
759	aa = (a/2)^2;
760	t = pm;
761	i = m + 1;
762	while t * 2147483648 > 1
763	t = t * aa/((i)*(i+nu));
764	i = i + 1;
765	end
766	last = i-2;
767
768	t = pm;
769	j = -1;
770	for i = m-1:-1:0
771	t = t * (i+1)*(i+1+nu)/aa;
772	if t*2147483648 < 1
773	j=i;
774	break;
775	end
776	end
777
778	offset = j+1;
779	sizeP = last-offset+1;
780
781	P = zeros(1, sizeP);
782	P(m-offset+1:last-offset+1) = ...
783	round( 2^30pmcumprod([1, aa./(((m+1):last).*(((m+1):last)+nu))] ) );
784	P(m-offset:-1:1) = ...
785	round( 2^30pmcumprod((m:-1:(offset+1)).*((m:-1:(offset+1))+nu)/aa) );
786
787	end % if p0 >= 5e-10, else ...
788
789
790	out = randFrom5Tbls( P, offset, sampleSize);
791
792	end % if proceed
793
794	case {'binom', 'binomial'}
795	% START binom HELP START binomial HELP
796	% THE BINOMIAL DISTRIBUTION
797	%
798	% pdf(y) = nchoosek(n,y)p^y(1-p)^(n-y) = ...
799	% exp( gammaln(n+1) - gammaln(n-y+1) - gammaln(y+1) + ...
800	% ylog(p) + (n-y)log(1-p) ); 0<p<1, n>1
801	%
802	% Mean = n*p;
803	% Variance = np(1-p);
804	% Mode = floor( (n+1)*p );
805	%
806	% PARAMETERS:
807	% p - probability of success in a single trial; (0<p<1)
808	% n - total number of trials; (n= 1, 2, 3, 4, ...)
809	%
810	% SUPPORT:
811	% y - number of success, y = 0, 1, 2, 3 ...
812	%
813	% CLASS:
814	% Discrete distributions
815	%
816	% NOTES:
817	% Constructive definition:
818	% We consider a random experiment with n independent trials; in each trial
819	% a certain random event A can occur (the urn model with replacement is
820	% a special case of such an experiment). Let
821	% p = probability of A in a single trial;
822	% n = total number of trials;
823	% y = number of successes (= number of trials where A occurs).
824	%
825	% USAGE:
826	% randraw('binom', [n, p], sampleSize) - generate sampleSize number
827	% of variates from the Binomial distribution with total number of trials
828	% 'n' and probability of success in a single trial 'p'
829	% randraw('binom') - help for the Binomial distribution;
830	%
831	% EXAMPLES:
832	% 1. y = randraw('binom', [10 0.9], [1 1e5]);
833	% 2. y = randraw('binom', [100 0.15], 1, 1e5);
834	% 3. y = randraw('binom', [5 0.5], 1e5 );
835	% 4. y = randraw('binom', [1000 0.02], [1e5 1] );
836	% 5. randraw('binom');
837	% END binom HELP END binomial HELP
838
839	% Method:
840	%
841	% We implemented Condensed Table-Lookup method suggested in
842	% George Marsaglia, "Fast Generation Of Discrete Random Variables,"
843	% Journal of Statistical Software, July 2004, Volume 11, Issue 4
844
845	checkParamsNum(funcName, 'Binomial', 'binomial', distribParams, [2]);
846
847	n = distribParams(1);
848	p = distribParams(2);
849
850	validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'n', n, {'> 0','==integer'});
851	validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'p', p, {'> 0','< 1'});
852
853	% if n large and p near 1, generate j=Binom(n,1-p), return n-j
854
855	switchFlag = 0;
856	if n*p <= 1
857	p = 1-p;
858	switchFlag = 1;
859	end
860
861	if n*p > 1e3 & n > 1e3
862
863	out = round( np + sqrt(np(1-p))randn( sampleSize ) );
864
865	elseif p<1e-4 & np > 1 & np < 100
866
867	out = feval(funcName,'poisson',n*p, sampleSize);
868
869	else
870
871	mode = floor( (n+1)*p );
872	q = 1 - p;
873	h = p/q;
874
875	pmode = exp( gammaln(n+1) - gammaln(n-mode+1) - gammaln(mode+1) + ...
876	modelog(p) + (n-mode)log(1-p) );
877
878	i = mode + 1;
879	t = pmode;
880	while t*2147483648 > 1
881	t = t * h*(n-i+1)/i;
882	i = i + 1;
883	end
884	last = i - 2;
885
886	t = pmode;
887	j = -1;
888	for i=mode-1:-1:0
889	t = t * (i+1)/h/(n-i);
890	if t*2147483648 < 1
891	j=i;
892	break;
893	end
894	end
895	offset=j+1;
896	sizeP = last-offset+1;
897
898	P = zeros(1, sizeP);
899
900	P(mode-offset+1:last-offset+1) = ...
901	round( 2^30pmodecumprod([1, h*(n-(mode+1:last)+1)./(mode+1:last)] ) );
902	P(mode-offset:-1:1) = ...
903	round( 2^30pmodecumprod( (mode:-1:offset+1)./(h*(n-(mode-1:-1:offset)))) );
904
905	out = randFrom5Tbls( P, offset, sampleSize);
906
907	end
908	if switchFlag
909	out = n - out;
910	end
911
912	case {'bradford'}
913	% START bradford HELP
914	% THE BRADFORD DISTRIBUTION
915	%
916	% pdf(y) = b ./ ( log(1+b)(1+by) ); 0<y<1
917	% cdf(y) = log(1+b*x) ./ log(1+b);
918	%
919	% Mean = (b - log(1+b)) / (b*log(1+b));
920	% Variance = (b(log(1+b)-2) + 2log(1+b)) / (2b(log(1+b))^2);
921	%
922	% PARAMETERS:
923	% b>-1,b~=0 - shape parameter;
924	%
925	% SUPPORT:
926	% 0<y<1
927	%
928	% CLASS:
929	% Continuous skewed distributions
930	%
931	% USAGE:
932	% randraw('bradford', [b], sampleSize) - generate sampleSize number
933	% of variates from the Bradford distribution with parameter b;
934	% randraw('bradford') - help for the Bradford distribution;
935	%
936	% EXAMPLES:
937	% 1. y = randraw('bradford', 1, [1 1e5]);
938	% 2. y = randraw('bradford', 2.2, 1, 1e5);
939	% 3. y = randraw('bradford', -0.4, 1e5 );
940	% 4. y = randraw('bradford', 10, [1e5 1] );
941	% 5. randraw('bradford');
942	% END bradford HELP
943
944	checkParamsNum(funcName, 'Bradford', 'bradford', distribParams, [1]);
945	b = distribParams(1);
946	validateParam(funcName, 'Bradford', 'bradford', 'b', 'b', b, {'> -1','~=0'});
947
948	out = ((1+b).^rand( sampleSize ) - 1)/b;
949
950	case {'burr', 'fisk'}
951	% START burr HELP START fisk HELP
952	% THE BURR DISTRIBUTION
953	% pdf(y) = cd y.^(-c-1) .* (1+y.^-c).^(-d-1); y>0
954	% cdf(y) = (1 + y.^-c).^(-d);
955	%
956	% Mean = gamma(1-1/c)*gamma(1/c+d)/gamma(d);
957	% Variance = COEF / gamma(d)^2;
958	% where COEF = gamma(d)gamma(1-2/c)gamma(2/c+d) - gamma(1-1/c)^2*gamma(1/c+d)^2;
959	%
960	% PARAMETERS:
961	% c - shape parameter (c>0)
962	% d - shape parameter (d>0)
963	%
964	% SUPPORT:
965	% y>0
966	%
967	% NOTES:
968	% The Burr distribution with d = 1, is often called the Fisk or
969	% LogLogistic distribution
970	% The Burr distribution is a generalization of the Fisk distribution
971	%
972	% CLASS:
973	% Continuous skewed distributions
974	%
975	% USAGE:
976	% randraw('burr', [c d], sampleSize) - generate sampleSize number
977	% of variates from the Burr distribution with shape parameters
978	% 'c' and 'd';
979	% randraw('burr') - help for the Burr distribution;
980	%
981	% EXAMPLES:
982	% 1. y = randraw('burr', [1 2], [1 1e5]);
983	% 2. y = randraw('burr', [2 3], 1, 1e5);
984	% 3. y = randraw('burr', [1.5 0.2], 1e5 );
985	% 4. y = randraw('burr', [2 2], [1e5 1] );
986	% 5. randraw('burr');
987	% END burr HELP END fisk HELP
988
989
990	checkParamsNum(funcName, 'Burr', 'burr', distribParams, [2]);
991	c = distribParams(1);
992	d = distribParams(2);
993	validateParam(funcName, 'Burr', 'burr', '[c, d]', 'c', c, {'> 0'});
994	validateParam(funcName, 'Burr', 'burr', '[c, d]', 'd', d, {'> 0'});
995
996	out = ( rand(sampleSize).^(-1/d) - 1).^(-1/c);
997
998	case {'cauchy', 'lorentz', 'caushy'}
999	% START cauchy HELP START lorentz HELP START caushy HELP
1000	% THE CAUCHY DISTRIBUTION
1001	% (sometimes: Lorentz or Breit-Wigner distribution)
1002	%
1003	% The standard form of the Caushy distribution:
1004	% pdf = 1 / ( pi*(1+y.^2) );
1005	% cdf = 0.5 + atan(y)/pi;
1006	%
1007	% The general form of the Cauchy distribution:
1008	% pdf = s ./ (pi*(s^2+(y-t).^2)); s>0;
1009	% cdf = 0.5 + atan((y-t)/s)/pi;
1010	%
1011	% The Cauchy distribution does not have a finite mean or
1012	% standard deviation.
1013	% Like the normal distribution, it is symmetric about its median,
1014	% but with longer and flatter tails.
1015	%
1016	% PARAMETERS:
1017	% s>0 - scale parameter;
1018	% t - loacation;
1019	%
1020	% SUPPORT;
1021	% -Inf < y < Inf
1022	%
1023	% CLASS:
1024	% Continuous skewed distributions
1025	%
1026	% USAGE:
1027	% randraw('caushy',[],sampleSize) - generation array of sampleSize
1028	% of variates from standard Cauchy distribution;
1029	% randraw('caushy',[t, s],sampleSize) - generation array of sampleSize
1030	% of variates from general form of Cauchy distribution;
1031	%
1032	% EXAMPLES:
1033	% 1. y = randraw('cauchy', [], [1 1e5]);
1034	% 2. y = randraw('cauchy', [10 1], 1, 1e5);
1035	% 3. y = randraw('cauchy', [-10 1.5], 1e5 );
1036	% 4. y = randraw('cauchy', [5.1 10.3], [1e5 1] );
1037	% 5. randraw('cauchy');
1038	% END cauchy HELP END lorentz HELP END caushy HELP
1039
1040	checkParamsNum(funcName, 'Cauchy', 'cauchy', distribParams, [0, 2]);
1041
1042	if numel(distribParams)==2
1043	t = distribParams(1);
1044	s = distribParams(2);
1045
1046	validateParam(funcName, 'Cauchy', 'cauchy', '[t, s]', 's', s, {'> 0'});
1047	else
1048	t = 0;
1049	s = 1;
1050	end
1051
1052	out = t + s * tan(pi*( rand( sampleSize ) - 0.5));
1053
1054	case {'chi'}
1055	% START chi HELP
1056	% THE CHI DISTRIBUTION
1057	%
1058	% The standard form of the Chi distribution:
1059	%
1060	% pdf(y) = exp(-y.^2/2).y.^(nu-1) / (2^(nu/2-1)gamma(nu/2)); nu>0; y>0
1061	% cdf(y) = gammainc(y.^2/2,nu/2);
1062	%
1063	% Mean = sqrt(2)*gamma((nu+1)/2)/gamma(nu/2);
1064	% Variance = nu - 2*gamma((nu+1)/2)^2/gamma(nu/2)^2;
1065	%
1066	% The general form of the Chi distribution:
1067	%
1068	% pdf(y) = exp(-((y-a)/b).^2/2).((y-a)/b).^(nu-1) / (2^(nu/2-1)b*gamma(nu/2)); nu>0; y>a; b>0
1069	% cdf(y) = gammainc(((y-a)/b).^2/2,nu/2);
1070	%
1071	% Mean = a + sqrt(2)bgamma((nu+1)/2)/gamma(nu/2);
1072	% Variance = b^2 * (nu-2*gamma((nu+1)/2)^2/gamma(nu/2)^2);
1073	%
1074	% PARAMETERS:
1075	% a - location
1076	% b > 0 - scale
1077	% nu > 0 - shape (also, degrees of freedom)
1078	%
1079	% SUPPORT:
1080	% y, y>0 - standard Chi distribution
1081	% or
1082	% y, y>a - generalized Chi distribution
1083	%
1084	% CLASS:
1085	% Continuous skewed distributions
1086	%
1087	% NOTES:
1088	% The chi-distribution includes several distributions as special cases.
1089	% If nu is 1, the chi-distribution reduces to the half-normal distribution.
1090	% If nu is 2, the chi-distribution is a Rayleigh distribution.
1091	% If nu is 3, the chi-distribution is a Maxwell-Boltzmann distribution.
1092	% A generalized Rayleigh distribution is a chi-distribution with a scale parameter equal to 1.
1093	%
1094	% USAGE:
1095	% randraw('chi', nu, sampleSize) - generate sampleSize number
1096	% of variates from the standrad Chi distribution with shape parameter 'nu';
1097	% randraw('chi', [a, b, nu], sampleSize) - generate sampleSize number
1098	% of variates from the generalized Chi distribution with location parameter
1099	% 'a', scale parameter 'b' and shape parameter 'nu';
1100	% randraw('chi') - help for the Chi distribution;
1101	%
1102	% EXAMPLES:
1103	% 1. y = randraw('chi', [2], [1 1e5]);
1104	% 2. y = randraw('chi', [3, 1, 5], 1, 1e5);
1105	% 3. y = randraw('chi', [-10 1, 1], 1e5 );
1106	% 4. y = randraw('chi', [-2.1 3.5 4], [1e5 1] );
1107	% 5. randraw('chi');
1108	% END chi HELP
1109
1110	checkParamsNum(funcName, 'Chi', 'chi', distribParams, [1, 3]);
1111	if numel(distribParams)==3
1112	a = distribParams(1);
1113	b = distribParams(2);
1114	nu = distribParams(3);
1115	validateParam(funcName, 'Chi', 'chi', '[a, b, nu]', 'b', b, {'> 0'});
1116	validateParam(funcName, 'Chi', 'chi', '[a, b, nu]', 'nu', nu, {'> 0'});
1117	else
1118	a = 0;
1119	b = 1;
1120	nu = distribParams(1);
1121	validateParam(funcName, 'Chi', 'chi', 'nu', 'nu', nu, {'> 0'});
1122	end
1123
1124	out = a + b*sqrt( feval(funcName, 'chisq', [nu], sampleSize) );
1125
1126	case {'chisquare', 'chisq', 'chi2'}
1127	% START chisquare HELP START chisq HELP START chi2 HELP
1128	% THE CHI SQUARE DISTRIBUTION (with r degrees of freedom)
1129	%
1130	% pdf(y) = y.^(r/2-1) .* exp(-y/2) / (gamma(r/2)*2^(r/2)); r >=1 (integer); y>0
1131	% cdf(y) = gammainc(y/2, r/2); r >=1 (integer); y>0;
1132	%
1133	% Mean = r;
1134	% Variance = 2*r;
1135	% Skewness = 2*sqrt(2/r);
1136	% Kurtosis = 12/r;
1137	%
1138	% PARAMETERS:
1139	% r - degrees of freedom ( r = 1, 2, 3, ...)
1140	%
1141	% SUPPORT:
1142	% y, y>0
1143	%
1144	% CLASS:
1145	% Continuous skewed distributions
1146	%
1147	% NOTES:
1148	% 1. Chi square distribution with r degrees of freedom is sum of r squared i.i.d Normal
1149	% distributions with zero mean and variance equal to 1;
1150	% 2. It is a special case of the gamma distribution where:
1151	% the scale parameter is 2 and the shape parameter has the value r/2;
1152	%
1153	% USAGE:
1154	% randraw('chisq', r, sampleSize) - generate sampleSize number
1155	% of variates from CHI SQUARE distribution with r degrees of freedom;
1156	% randraw('chisq') - help for CHI SQUARE distribution;
1157	%
1158	% EXAMPLES:
1159	% 1. y = randraw('chisq', [2], [1 1e5]);
1160	% 2. y = randraw('chisq', [3], 1, 1e5);
1161	% 3. y = randraw('chisq', [4], 1e5 );
1162	% 4. y = randraw('chisq', [5], [1e5 1] );
1163	% 5. randraw('chisq');
1164	%
1165	% SEE ALSO:
1166	% GAMMA, NON-CENTRAL CHI SQUARE distributions
1167	% END chisquare HELP END chisq HELP END chi2 HELP
1168
1169
1170	checkParamsNum(funcName, 'Chi Square', 'chisq', distribParams, [1]);
1171	r = distribParams(1);
1172	validateParam(funcName, 'Chi Square', 'chisq', 'r', 'r', r, {'> 0','==integer'});
1173
1174	if r > 1
1175	out = 2randraw('gamma', 0.5r, sampleSize);
1176	else
1177	out = randn(sampleSize).^2;
1178	end
1179
1180	case {'chisqnc','chisqnoncentral', 'chisqnoncentr','chi2noncentral'}
1181	% START chisqnc HELP START chisqnoncentral HELP START chisqnoncentr HELP START chi2noncentral HELP
1182	% THE NON-CENTRAL CHI-SQUARE DISTRIBUTION (with non-centrality parameter lambda and
1183	% r degrees of freedom)
1184	%
1185	% The non-central chi-square distribution with degrees of freedom r and
1186	% non-centrality parameter lambda is the sum of r independent normal
1187	% distributions with standard deviation 1.
1188	% The non-centrality parameter is one half the sum of squares of the normal
1189	% means.
1190	%
1191	%
1192	% pdf(y) = exp(-(y+lambda)/2).y.^((r-1)/2)./(2(lambday).^(r/4)) . ...
1193	% besseli(r/2-1, sqrt(lambda*y)); lambda>=0; r=positive integer;
1194	%
1195	% Mean = lambda+r;
1196	% Variance = 2(2lambda+r);
1197	% Skewness = 2sqrt(2)(3lambda+r)/(2lambda+r)^(3/2);
1198	% Kurtosis = 12(4lambda+r)/(2*lambda+r)^2;
1199	%
1200	% PARAMETERS:
1201	% lambda - non-centrality parameter: lambda>=0
1202	% r - degrees of freedom ( r = 1, 2, 3, ...)
1203	%
1204	% SUPPORT:
1205	% y, y>0
1206	%
1207	% CLASS:
1208	% Continuous skewed distributions
1209	%
1210	% USAGE:
1211	% randraw('chisqnoncentral', [lambda, r], sampleSize) - generate sampleSize number
1212	% of variates from the NON CENTRAL CHI SQUARE distribution with
1213	% non-centrality parameter 'lambda ' and 'r' degrees of freedom;
1214	% randraw('chisqnoncentral') - help for the NON CENTRAL CHI-SQUARE distribution;
1215	%
1216	% EXAMPLES:
1217	% 1. y = randraw('chisqnoncentral', [10 2], [1 1e5]);
1218	% 2. y = randraw('chisqnoncentral', [20 3], 1, 1e5);
1219	% 3. y = randraw('chisqnoncentral', [30 4], 1e5 );
1220	% 4. y = randraw('chisqnoncentral', [40 5], [1e5 1] );
1221	% 5. randraw('chisqnoncentral');
1222	%
1223	% SEE ALSO:
1224	% CHI SQUARE distribution
1225	% END chisqnc HELP START END chisqnoncentral HELP END chisqnoncentr HELP END chi2noncentral HELP
1226
1227	checkParamsNum(funcName, 'Non-Central Chi-Square', 'chisqnoncentral', distribParams, [2]);
1228	lambda = distribParams(1);
1229	r = distribParams(2);
1230	validateParam(funcName, 'Non-Central Chi-Square', 'chisqnoncentral', 'r', 'r', r, {'> 0','==integer'});
1231
1232	normalS = sqrt(lambda) + randn( sampleSize );
1233	out = feval(funcName, 'chisq', r-1, sampleSize);
1234	out = out + normalS.^2;
1235
1236	case {'cosine'}
1237	% START cosine HELP
1238	% THE COSINE DISTRIBUTION
1239	%
1240	% Standard form of the Cosine distribution:
1241	% pdf(y) = (1+cos(y))/(2*pi); -pi <= y <= pi;
1242	% cdf(y) = (pi+y+sin(y))/(2*pi); -pi <= y <= pi;
1243	%
1244	% Mean = Median = Mode = 0;
1245	% Variance = pi^2/3-2;
1246	%
1247	% General form of the Cosine distribution:
1248	% pdf(y) = (1+cos(pi(y-t)/s))/(2s); t-s <= y <= t+s; s>0
1249	% cdf(y) = (pi + pi(y-t)/s + sin(pi(y-t)/s))/(2*pi); t-s <= y <= t+s; s>0
1250	%
1251	% Mean = Median = Mode = t;
1252	% Variance = (pi^2/3-2)*(s/pi)^2;
1253	%
1254	% PARAMETERS:
1255	% t - location
1256	% s -scale; s>0
1257	%
1258	% SUPPORT:
1259	% y, -pi <= y <= pi - standard Cosine distribution
1260	% or
1261	% y, t-s <= y <= t+s - generalized Cosine distribution
1262	%
1263	% CLASS:
1264	% Continuous distributions
1265	%
1266	% USAGE:
1267	% randraw('cosine', [], sampleSize) - generate sampleSize number
1268	% of variates from the standard Cosine distribution;
1269	% randraw('cosine', [t, s], sampleSize) - generate sampleSize number
1270	% of variates from the generalized Cosine distribution
1271	% with location parameter 't' and scale parameter 's';
1272	% randraw('cosine') - help for the Cosine distribution;
1273	%
1274	% EXAMPLES:
1275	% 1. y = randraw('cosine', [], [1 1e5]);
1276	% 2. y = randraw('cosine', [], 1, 1e5);
1277	% 3. y = randraw('cosine', [], 1e5 );
1278	% 4. y = randraw('cosine', [10 3], [1e5 1] );
1279	% 5. randraw('cosine');
1280	% END cosine HELP
1281
1282	checkParamsNum(funcName, 'Cosine', 'cosine', distribParams, [0, 2]);
1283	if numel(distribParams)==2
1284	t = distribParams(1);
1285	s = distribParams(2);
1286	validateParam(funcName, 'Cosine', 'cosine', '[t, s]', 's', s, {'> 0'});
1287	else
1288	t = 0;
1289	s = pi;
1290	end
1291
1292	tol = 1e-9;
1293
1294	coeff1 = 1/(2*pi);
1295	coeff2 = 1/(2*s);
1296	coeff3 = pi/s;
1297
1298
1299	u = 0.5 - rand(sampleSize);
1300	out = -u*s;
1301	outNext = out - (coeff1sin(coeff3out) + coeff2*out + u) ./ ...
1302	(coeff2cos(coeff3out)+coeff2);
1303
1304	indxs = find(abs(outNext - out)>tol);
1305	outPrev = out(indxs);
1306	while ~isempty(indxs)
1307
1308	outNext = outPrev - (coeff1sin(coeff3outPrev) + coeff2*outPrev + u(indxs)) ./ ...
1309	(coeff2cos(coeff3outPrev)+coeff2);
1310	l = (abs(outNext - outPrev)>tol);
1311	out(indxs(~l)) = outNext(~l);
1312	outPrev = outNext(l);
1313	indxs = indxs(l);
1314	end
1315
1316	out = t + out;
1317
1318	case {'erlang'}
1319	% START erlang HELP
1320	% THE ERLANG DISTRIBUTION
1321	%
1322	% pdf = (y/a).^(n-1) .* exp( -y/a ) / (a*gamma(n));
1323	% cdf = gammainc( n, y/sacle );
1324	%
1325	% Mean = a*n;
1326	% Variance = a^2*n;
1327	% Skewness = 2/sqrt(n);
1328	% Kurtosis = 6/n;
1329	% Mode = (a<1)0 + (a>=1)a*(n-1);
1330	%
1331	% PARAMETERS:
1332	% a - scale parameter (a>0)
1333	% n - shape parameter (n = 1, 2, 3, ...)
1334	%
1335	% SUPPORT:
1336	% y, y >= 0
1337	%
1338	% CLASS:
1339	% Continuous skewed distributions
1340	%
1341	% NOTES:
1342	% The Erlang distribution is a special case of the gamma distribution where
1343	% the shape parameter is an integer
1344	%
1345	% USAGE:
1346	% randraw('erlang', [a, n], sampleSize) - generate sampleSize number
1347	% of variates from the Erlang distribution
1348	% with scale parameter 'a' and shape parameter 'n';
1349	% randraw('erlang') - help for the Erlang distribution;
1350	%
1351	% EXAMPLES:
1352	% 1. y = randraw('erlang', [1, 3], [1 1e5]);
1353	% 2. y = randraw('erlang', [0.5, 5], 1, 1e5);
1354	% 3. y = randraw('erlang', [10, 6], 1e5 );
1355	% 4. y = randraw('erlang', [7, 4], [1e5 1] );
1356	% 5. randraw('erlang');
1357	%
1358	% SEE ALSO:
1359	% GAMMA distribution
1360	% END erlang HELP
1361
1362	%
1363	% Inverse CDF transformation method.
1364	%
1365
1366	checkParamsNum(funcName, 'Erlang', 'erlang', distribParams, [2]);
1367	a = distribParams(1);
1368	n = distribParams(2);
1369	validateParam(funcName, 'Erlang', 'erlang', '[a, n]', 'a', a, {'> 0'});
1370	validateParam(funcName, 'Erlang', 'erlang', '[a, n]', 'n', n, {'> 0', '==integer'});
1371
1372	out = feval(funcName, 'gamma', n, sampleSize);
1373	out = a * out;
1374
1375	case {'exp','exponential'}
1376	% START exp HELP START exponential HELP
1377	% THE EXPONENTIAL DISTRIBUTION
1378	%
1379	% pdf = lambda * exp( -lambda*y );
1380	% cdf = 1 - exp(-lambda*y);
1381	%
1382	% Mean = 1/lambda;
1383	% Variance = 1/lambda^2;
1384	% Mode = lambda;
1385	% Median = log(2)/lambda;
1386	% Skewness = 2;
1387	% Kurtosis = 6;
1388	%
1389	% PARAMETERS:
1390	% lambda - inverse scale or rate (lambda>0)
1391	%
1392	% SUPPORT:
1393	% y, y>= 0
1394	%
1395	% CLASS:
1396	% Continuous skewed distributions
1397	%
1398	% NOTES:
1399	% The discrete version of the Exponential distribution is
1400	% the Geometric distribution.
1401	%
1402	% USAGE:
1403	% randraw('exp', lambda, sampleSize) - generate sampleSize number
1404	% of variates from the Exponential distribution
1405	% with parameter 'lambda';
1406	% randraw('exp') - help for the Exponential distribution;
1407	%
1408	% EXAMPLES:
1409	% 1. y = randraw('exp', 1, [1 1e5]);
1410	% 2. y = randraw('exp', 1.5, 1, 1e5);
1411	% 3. y = randraw('exp', 2, 1e5 );
1412	% 4. y = randraw('exp', 3, [1e5 1] );
1413	% 5. randraw('exp');
1414	%
1415	% SEE ALSO:
1416	% GEOMETRIC, GAMMA, POISSON, WEIBULL distributions
1417	% END exp HELP END exponential HELP
1418
1419	checkParamsNum(funcName, 'Exponential', 'exp', distribParams, [1]);
1420	lambda = distribParams(1);
1421	validateParam(funcName, 'Exponential', 'exp', 'lambda', 'lambda', lambda, {'> 0'});
1422
1423	out = -log( rand( sampleSize ) ) / lambda;
1424
1425	case {'extrval', 'extremevalue', 'extrvalue', 'gumbel'}
1426	% START extrval HELP START extremevalue HELP START extrvalue HELP START gumbel HELP
1427	% THE EXTREME VALUE DISTRIBUTION
1428	% Also known as the Fisher-Tippett distribution or log-Weibull distribution or Gumbel
1429	% distribution
1430	%
1431	% pdf(y) = 1/b * exp((mu-y)/b - exp((mu-y)/b)); -Inf<y<Inf; b>0
1432	% cdf(y) = exp(-exp((mu-y)/b)); -Inf<y<Inf; b>0
1433	%
1434	% Mean = mu + b*g; where g=5.772156649015329e-001; is the Euler-Mascheroni constant
1435	% Variance = pi^2*b^2/6;
1436	% Skewness = 12sqrt(6)zeta3/pi^3; where zeta3=1.20205690315732e+000; is Apery's constant
1437	% Kurtosis = 12/5;
1438	%
1439	% PARAMETERS:
1440	% mu - location (-Inf<mu<inf)
1441	% b - scale (b>0)
1442	%
1443	% SUPPORT:
1444	% y, -Inf<y<Inf
1445	%
1446	% CLASS:
1447	% Continuous skewed distributions
1448	%
1449	% NOTES:
1450	%
1451	% USAGE:
1452	% randraw('extrvalue', [mu b], sampleSize) - generate sampleSize number
1453	% of variates from Extreme-Value distribution with location parameter mu
1454	% and scale parameter b;
1455	% randraw('extrvalue') - help for Extreme-Value distribution;
1456	%
1457	% EXAMPLES:
1458	% 1. y = randraw('extrvalue', [0 1], [1 1e5]);
1459	% 2. y = randraw('extrvalue', [10 2], 1, 1e5);
1460	% 3. y = randraw('extrvalue', [15.5 4], 1e5 );
1461	% 4. y = randraw('extrvalue', [100 100], [1e5 1] );
1462	% 5. randraw('extrvalue');
1463	%
1464	% SEE ALSO:
1465	% GUMBEL, WEIBULL distributions
1466	% END extrval HELP END extremevalue HELP END extrvalue HELP END gumbel HELP
1467
1468	checkParamsNum(funcName, 'Extreme Value', 'extrvalue', distribParams, [2]);
1469	mu = distribParams(1);
1470	b = distribParams(2);
1471	validateParam(funcName, 'Extreme Value', 'extrvalue', '[mu, b]', 'b', b, {'> 0'});
1472
1473	out = mu - b * log(-log( rand( sampleSize )));
1474
1475	case {'f','fdistribution', 'fdistrib', 'fdistr', 'fdist', 'fdis' }
1476	% START f HELP START fdistribution HELP START fdistrib HELP START fdistr HELP START fdist HELP START fdis HELP
1477	% THE F-DISTRIBUTION (also Central F-distribution)
1478	%
1479	% In statistics and probability, the F-distribution is a continuous
1480	% probability distribution. It is also known as Snedecor's F distribution or
1481	% the Fisher-Snedecor distribution (after Ronald Fisher and George W. Snedecor).
1482	% A random variate of the F-distribution arises as the ratio of two chi-squared
1483	% variates: (U1/d1)/(U2/d2), where U1 and U2 have chi-square distributions with
1484	% d1 and d2 degrees of freedom respectively, and U1 and U2 are independent.
1485	% The F-distribution arises frequently as the null distribution of a test statistic,
1486	% especially in likelihood-ratio tests, perhaps most notably in the analysis of
1487	% variance;
1488	%
1489	% pdf(y) = 1/beta(d1/2,d2/2) * (d1y./(d1y+d2)).^(d1/2) .* (1 - d1y./(d1y+d2)).^(d2/2) ./ y;
1490	% cdf(y) = beatinc(d1y./(d1y+d2), d1/2, d2/2);
1491	%
1492	% Mean = d2/(d2-2), provided d2 > 2;
1493	% Variance = 2d2^2(d1+d2-2)/(d1(d2-2)^2(d2-4)), provided d2>4;
1494	% Skewness = (2d1+d2-2)sqrt(8(d2-4))/((d2-6)sqrt(d1*(d1+d2-2))), provided d2>6;
1495	% Mode = (d1-2)/d1 * d2/(d2+2), provided d1>2;
1496	%
1497	% PARAMETERS:
1498	% d1 - positive integer
1499	% d2 - positive integer
1500	%
1501	% SUPPORT:
1502	% y, y>0
1503	%
1504	% CLASS:
1505	% Continuous skewed distributions
1506	%
1507	% USAGE:
1508	% randraw('f', [d1 d2], sampleSize) - generate sampleSize number
1509	% of variates from F-distribution with parameters d1 and d2;
1510	% randraw('f') - help for F-distribution;
1511	%
1512	% EXAMPLES:
1513	% 1. y = randraw('f', [2 3], [1 1e5]);
1514	% 2. y = randraw('f', [2 3], 1, 1e5);
1515	% 3. y = randraw('f', [2 3], 1e5 );
1516	% 4. y = randraw('f', [2 3], [1e5 1] );
1517	% 5. randraw('f');
1518	% END f HELP END fdistribution HELP END fdistrib HELP END fdistr HELP END fdist HELP END fdis HELP
1519
1520	checkParamsNum(funcName, 'F', 'f', distribParams, [2]);
1521	d1 = distribParams(1);
1522	d2 = distribParams(2);
1523	validateParam(funcName, 'F', 'f', '[d1, d2]', 'd1', d1, {'> 0','==integer'});
1524	validateParam(funcName, 'F', 'f', '[d1, d2]', 'd2', d2, {'> 0','==integer'});
1525
1526	out = feval(funcName, 'beta', [0.5d1 0.5d2], sampleSize);
1527	out = d2out ./ (d1(1-out));
1528
1529	case {'fnc', 'fnoncentral', 'fnoncentr'}
1530	% START fnc HELP START fnoncentral HELP START fnoncentr HELP
1531	% THE NONCENTRAL F-DISTRIBUTION
1532	%
1533	% The central F distribution is the ratio of 2 central chi-square distributions with
1534	% d1 and d2 degrees of freedom respectively. The noncentral F distribution is the ratio
1535	% of a non-central chi-square distribution with d1 degrees of freedom and non-centrality
1536	% parameter lambda and a central chi-square distribution with degrees of freedom parameter
1537	% d2.
1538	% The non-centrality parameter should be non-negative, and both degrees of freedom parameters
1539	% should be positive.
1540	%
1541	% Mean = (d1+lambda)d2/(d1(d2-2)), provided d2 > 2;
1542	% Variance = ( ( d1+lambda )^2 + 2( d1+lambda )d2^2 )/( (d2-2)(d2-4)d1^2 ) - ...
1543	% (d1+lambda)^2d2^2/((d2-2)^2d1^2); provided d2 > 4;
1544	%
1545	% PARAMETERS:
1546	% lambda - non-centrality parameter (lambda>=0);
1547	% d1 - positive integer
1548	% d2 - positive integer
1549	%
1550	% SUPPORT:
1551	% y, y>0
1552	%
1553	% CLASS:
1554	% Continuous skewed distributions
1555	%
1556	% USAGE:
1557	% randraw('fnoncentral', [lambda d1 d2], sampleSize) - generate sampleSize number
1558	% of variates from noncentral F-distribution with parameters lambda, d1 and d2;
1559	% randraw('fnoncentral') - help for noncentral F-distribution;
1560	%
1561	% EXAMPLES:
1562	% 1. y = randraw('fnoncentral', [1 2 5], [1 1e5]);
1563	% 2. y = randraw('fnoncentral', [2 3 5], 1, 1e5);
1564	% 3. y = randraw('fnoncentral', [6 6 6], 1e5 );
1565	% 4. y = randraw('fnoncentral', [1.1 8 9], [1e5 1] );
1566	% 5. randraw('fnoncentral');
1567	%
1568	% SEE ALSO:
1569	% F distribution
1570	% END fnc HELP END fnoncentral HELP END fnoncentr HELP
1571
1572	checkParamsNum(funcName, 'noncentral F', 'fnoncentral', distribParams, [3]);
1573
1574	lambda = distribParams(1);
1575	d1 = distribParams(2);
1576	d2 = distribParams(3);
1577
1578	validateParam(funcName, 'noncentral F', 'fnoncentral', '[lambda, d1, d2]', 'lambda', lambda, {'>=0'});
1579	validateParam(funcName, 'noncentral F', 'fnoncentral', '[lambda, d1, d2]', 'd1', d1, {'> 0','==integer'});
1580	validateParam(funcName, 'noncentral F', 'fnoncentral', '[lambda, d1, d2]', 'd2', d2, {'> 0','==integer'});
1581
1582	chisq1 = feval(funcName, 'chisqnoncentral', [lambda, d1], sampleSize);
1583	out = feval(funcName, 'chisq', d2, sampleSize);
1584	out = (chisq1/d1) ./ (out/d2);
1585
1586	case {'gamma'}
1587	% START gamma HELP START gama HELP
1588	% THE GAMMA DISTRIBUTION
1589	%
1590	% The standard form of the GAMMA distribution:
1591	%
1592	% pdf(y) = y^(a-1)*exp(-y)/gamma(a); y>=0, a>0
1593	% cdf(y) = gammainc(y, a);
1594	%
1595	% Mean = a;
1596	% Variance = a;
1597	% Skewness = 2/sqrt(a);
1598	% Kurtosis = 6/a;
1599	% Mode = a-1;
1600	%
1601	% The general form of the GAMMA distribution:
1602	%
1603	% pdf(y) = ((y-m)/b).^(a-1) .* exp(-(y-m)/b)/ (b*gamma(a)); y>=m; a>0; b>0
1604	% cdf(y) = gammainc((y-m)/b, a); y>=m; a>0; b>0
1605	%
1606	% Mean = m + a*b;
1607	% Variance = a*b^2;
1608	% Skewness = 2/sqrt(a);
1609	% Kurtosis = 6/a;
1610	% Mode = m + b*(a-1);
1611	%
1612	% PARAMETERS:
1613	% m - location
1614	% b - scale; b>0
1615	% a - shape; a>0
1616	%
1617	% SUPPORT:
1618	% y, y>=0 - standard GAMMA distribution
1619	% or
1620	% y, y>=m - generalized GAMMA distribution
1621	%
1622	% CLASS:
1623	% Continuous skewed distributions
1624	%
1625	% NOTES:
1626	% 1. The GAMMA distribution approaches a NORMAL distribution as a goes to Inf
1627	% 5. GAMMA(m, b, a), where a is an integer, is the Erlang distribution.
1628	% 6. GAMMA(m, b, 1) is the Exponential distribution.
1629	% 7. GAMMA(0, 2, nu/2) is the Chi-square distribution with nu degrees of freedom.
1630	%
1631	% USAGE:
1632	% randraw('gamma', a, sampleSize) - generate sampleSize number
1633	% of variates from standard GAMMA distribution with shape parameter 'a';
1634	% randraw('gamma', [m, b, a], sampleSize) - generate sampleSize number
1635	% of variates from generalized GAMMA distribution
1636	% with location parameter 'm', scale parameter 'b' and shape parameter 'a';
1637	% randraw('gamma') - help for GAMMA distribution;
1638	%
1639	% EXAMPLES:
1640	% 1. y = randraw('gamma', [2], [1 1e5]);
1641	% 2. y = randraw('gamma', [0 10 2], 1, 1e5);
1642	% 3. y = randraw('gamma', [3], 1e5 );
1643	% 4. y = randraw('gamma', [1/3], 1e5 );
1644	% 5. y = randraw('gamma', [1 3 2], [1e5 1] );
1645	% 6. randraw('gamma');
1646	%
1647	% END gamma HELP END gama HELP
1648
1649	% Method:
1650	%
1651	% Reference:
1652	% George Marsaglia and Wai Wan Tsang, "A Simple Method for Generating Gamma
1653	% Variables": ACM Transactions on Mathematical Software, Vol. 26, No. 3,
1654	% September 2000, Pages 363-372
1655
1656	checkParamsNum(funcName, 'Gamma', 'gamma', distribParams, [1, 3]);
1657	if numel(distribParams)==3
1658	m = distribParams(1);
1659	b = distribParams(2);
1660	a = distribParams(3);
1661	validateParam(funcName, 'Gamma', 'gamma', '[m, b, a]', 'a', a, {'> 0'});
1662	validateParam(funcName, 'Gamma', 'gamma', '[m, b, a]', 'b', b, {'> 0'});
1663	else
1664	m = 0;
1665	b = 1;
1666	a = distribParams(1);
1667	validateParam(funcName, 'Gamma', 'gamma', '[m, b, a]', 'a', a, {'> 0'});
1668	end
1669
1670	if a < 1
1671	% If a<1, one can use GAMMA(a)=GAMMA(1+a)*UNIFORM(0,1)^(1/a);
1672	out = m + b(feval(funcName, 'gamma', 1+a, sampleSize)).(rand(sampleSize).^(1/a));
1673
1674	else
1675
1676	d = a - 1/3;
1677	c = 1/sqrt(9*d);
1678
1679	x = randn( sampleSize );
1680	v = 1+c*x;
1681
1682	indxs = find(v <= 0);
1683	while ~isempty(indxs)
1684	indxsSize = size( indxs );
1685	xNew = randn( indxsSize );
1686	vNew = a+c*xNew;
1687
1688	l = (vNew > 0);
1689	v( indxs( l ) ) = vNew(l);
1690	x( indxs( l ) ) = xNew(l);
1691	indxs = indxs( ~l );
1692	end
1693
1694	u = rand( sampleSize );
1695	v = v.^3;
1696	x2 = x.^2;
1697	out = d*v;
1698
1699	indxs = find( (u>=1-0.0331x2.^2) & (log(u)>=0.5x2+d*(1-v+log(v))) );
1700	while ~isempty(indxs)
1701	indxsSize = size( indxs );
1702
1703	x = randn( indxsSize );
1704	v = 1+c*x;
1705	indxs1 = find(v <= 0);
1706	while ~isempty(indxs1)
1707	indxsSize1 = size( indxs1 );
1708	xNew = randn( indxsSize1 );
1709	vNew = a+c*xNew;
1710
1711	l1 = (vNew > 0);
1712	v( indxs1(l1) ) = vNew(l1);
1713	x( indxs1(l1) ) = xNew(l1);
1714	indxs1 = indxs1( ~l1 );
1715	end
1716
1717	u = rand( indxsSize );
1718	v = v .* v .* v;
1719	x2 = x.*x;
1720
1721	l = (u<1-0.0331x2.x2) \| (log(u)<0.5x2+d(1-v+log(v)));
1722	out( indxs( l ) ) = d*v(l);
1723	indxs = indxs( ~l );
1724	end % while ~isempty(indxs)
1725
1726	out = m + b*out;
1727
1728	end % if a < 1, else ...
1729
1730	case {'geometric', 'geom', 'furry'}
1731	% START geometric HELP START geom HELP START furry HELP
1732	% THE GEOMETRIC DISTRIBUTION
1733	%
1734	% pdf(n) = p*(1-p)^(n-1);
1735	%
1736	% Mean = 1/p;
1737	% Variance = (1-p)/p^2;
1738	% Mode = 1;
1739	%
1740	% PARAMETERS:
1741	% p - probability of success (0<p<1)
1742	%
1743	% SUPPORT:
1744	% n, n = 1, 2, 3, 4, ...
1745	%
1746	% CLASS:
1747	% Discrete distributions
1748	%
1749	% NOTES:
1750	% 1. The Geometric distribution is the discrete version of
1751	% the Exponential distribution.
1752	% 2. The Geometric distribution is sometimes called
1753	% the Furry distribution.
1754	% 3. In a series of Bernoulli trials, with Prob(success) = p,
1755	% the number of trials required to realize the first
1756	% success is ~Geometric(p).
1757	% 4. For the k'th success, see the Negative Binomial distribution.
1758	%
1759	% USAGE:
1760	% randraw('geom', p, sampleSize) - generate sampleSize number
1761	% of variates from the Geometric distribution with
1762	% probability of success 'p'
1763	% randraw('geom') - help for the Geometric distribution;
1764	%
1765	% EXAMPLES:
1766	% 1. y = randraw('geom', 0.1, [1 1e5]);
1767	% 2. y = randraw('geom', 0.22, 1, 1e5);
1768	% 3. y = randraw('geom', 0.5, 1e5 );
1769	% 4. y = randraw('geom', 0.99, [1e5 1] );
1770	% 5. randraw('geom');
1771	% END geometric HELP END geom HELP END furry HELP
1772
1773	checkParamsNum(funcName, 'Geometric', 'geom', distribParams, [1]);
1774	p = distribParams(1);
1775	validateParam(funcName, 'Geometric', 'geom', 'p', 'p', p, {'> 0'});
1776	validateParam(funcName, 'Geometric', 'geom', 'p', 'p', p, {'< 1'});
1777
1778	out = ceil( log( rand( sampleSize ) ) / log( 1 - p ) );
1779
1780	case {'gig'}
1781	% START gig HELP
1782	% THE GENERALIZED INVERSE GAUSSIAN DISTRIBUTION
1783	% GIG(lam, chi, psi)
1784	%
1785	% pdf = (psi/chi)^(lam/2)y.^(lam-1)/(2besselk(lam, sqrt(chipsi))) . exp(-1/2(chi./y + psiy)); y > 0
1786	%
1787	% Mean = sqrt( chi / psi ) * besselk(lam+1,sqrt(chipsi),1)/besselk(lam,sqrt(chipsi),1);
1788	% Variance = chi/psi * besselk(lam+2,sqrt(chipsi),1)/besselk(lam,sqrt(chipsi),1) - Mean^2;
1789	%
1790	% PARAMETERS:
1791	% chi>0, psi>=0 if lam<0;
1792	% chi>0, psi>0 if lam=0;
1793	% chi>=0, psi>0 if lam>0;
1794	%
1795	% SUPPORT:
1796	% y, y >= 0
1797	%
1798	% CLASS:
1799	% Continuous skewed distributions
1800	%
1801	% NOTES:
1802	% 1) GIG(lam, chi, psi) = 1/c * GIG(lam, chi*c, psi/c), for all c>=0
1803	% 2) GIG(lam, chi, psi) = sqrt(chi/psi) * GIG(lam, sqrt(psichi), sqrt(psichi));
1804	% 3) GIG(lam, chi, psi) = 1 / GIG(-lam, psi, chi);
1805	%
1806	% Special cases of GIG distribution are the gamma distribution (chi=0), the
1807	% reciprocal gamma distribution (psi=0), the inverse Gaussian distribution
1808	% (lam = -1/2), and the inverse Gaussian or random walk distribution (lam=1/2).
1809	%
1810	% USAGE:
1811	% randraw('gig', [lam, chi, psi], sampleSize) - generate sampleSize number
1812	% of variates from the Generalized Inverse Gaussian distribution with
1813	% parameters 'lam', 'chi' and 'psi'
1814	% randraw('gig') - help for the Generalized Inverse Gaussian distribution;
1815	%
1816	% EXAMPLES:
1817	% 1. y = randraw('gig', [-1, 2, 0], [1 1e5]);
1818	% 2. y = randraw('gig', [2, 3, 4], 1, 1e5);
1819	% 3. y = randraw('gig', [0, 1.1, 2.2], 1e5 );
1820	% 4. y = randraw('gig', [2.5, 3.5, 4.5], [1e5 1] );
1821	% 5. y = randraw('gig', [0.5, 0.6, 0.7], [1e5 1] );
1822	% 6. randraw('gig');
1823	% END gig HELP
1824
1825	% Reference:
1826	% 1. Dagpunar, J.S., "Principles of random variate generation,"
1827	% Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9
1828	% 2. Dagpunar, J.S., "An easily implemented generalized inverse Gaussian generator,"
1829	% Commun. Statist. Simul. 18(2), 1989, pp 703-710.
1830
1831	checkParamsNum(funcName, 'Generalized Inverse Gaussian', 'gig', distribParams, [3]);
1832	lam = distribParams(1);
1833	chi = distribParams(2);
1834	psi = distribParams(3);
1835
1836	if lam < 0,
1837	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'chi', chi, {'> 0'});
1838	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'psi', psi, {'>=0'});
1839	elseif lam > 0,
1840	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'chi', chi, {'>=0'});
1841	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'psi', psi, {'> 0'});
1842	else % lam==0
1843	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'chi', chi, {'> 0'});
1844	validateParam(funcName, 'Generalized Inverse Gaussian', 'gig', '[lambda, chi, psi]', 'psi', psi, {'> 0'});
1845	end
1846
1847	if chi == 0,
1848	% Gamma distribution: Gamma(m=0, b=2/psi, lam)
1849	out = feval(funcName, 'gamma', [0, 2/psi, lam], sampleSize);
1850	varargout{1} = out;
1851	return;
1852	end
1853
1854	if psi == 0,
1855	% Reciprocal Gamma distribution: Gamma(m=0, b=2/chi, -lam)
1856	out = feval(funcName, 'gamma', [0, 2/chi, -lam], sampleSize);
1857	varargout{1} = 1./out;
1858	return;
1859	end
1860
1861	h = lam;
1862	b = sqrt( chi * psi );
1863
1864	if h<=1 & b<=1
1865	% without shifting by m
1866
1867	ym = (-h-1 + sqrt((h+1)^2 + b^2))/b;
1868	xm = ( h-1 + sqrt((h-1)^2 + b^2))/b;
1869	% a = vplus/uplus
1870	a = exp(-0.5hlog(xmym) + 0.5log(xm/ym) + b/4*(xm + 1/xm - ym - 1.0/ym));
1871	% c = 1/log(sqrt(hx(xm)))
1872	c = -(h-1)/2* log(xm) + b/4*(xm + 1/xm);
1873	% vminus = 0
1874
1875	u = rand( sampleSize );
1876	v = rand( sampleSize );
1877	out = a * (v./u);
1878	indxs = find( log(u) > (h-1)/2log(out) - b/4(out + 1./out) + c );
1879	while ~isempty( indxs )
1880	indxsSize = size( indxs );
1881	u = rand( indxsSize );
1882	v = rand( indxsSize );
1883	outNew = a * (v./u);
1884	l = log(u) <= (h-1)/2log(outNew) - b/4(outNew + 1./outNew) + c;
1885	out( indxs( l ) ) = outNew(l);
1886	indxs = indxs( ~l );
1887	end
1888
1889	else % if h<=1 & b<=1
1890	% with shifting by m
1891
1892	% Mode of the reparameterized distribution GIG(lam, b, b)
1893	m = ( h-1+sqrt((h-1)^2+b^2) ) / b; % Mode
1894	log_1_over_pm = -(h-1)/2log(m) + b/4(m + (1/m));
1895
1896	r = (6m + 2hm - bm^2 + b)/(4*m^2);
1897	s = (1 + h - bm)/(2m^2);
1898	p = (3*s - r^2)/3;
1899	q = (2r^3)/27 - (rs)/27 + b/(-4*m^2);
1900	eta = sqrt(-(p^3)/27);
1901
1902	y1 = 2exp(log(eta)/3) cos(acos(-q/(2*eta))/3) - r/3;
1903	y2 = 2exp(log(eta)/3) cos(acos(-q/(2eta))/3 + 2/3pi) - r/3;
1904
1905
1906	vplus = exp( log_1_over_pm + log(1/y1) + (h-1)/2*log(1/y1 + m) - ...
1907	b/4*(1/y1 + m + 1/(1/y1 + m)) );
1908	vminus = -exp( log_1_over_pm + log(-1/y2) + (h-1)/2*log(1/y2 + m) - ...
1909	b/4*(1/y2 + m + 1/(1/y2 + m)) );
1910
1911	u = rand( sampleSize );
1912	v = vminus + (vplus - vminus) * rand( sampleSize );
1913	z = v ./ u;
1914	clear('v');
1915	indxs = find( z < -m );
1916
1917	while ~isempty(indxs),
1918	indxsSize = size( indxs );
1919	uNew = rand( indxsSize );
1920	vNew = vminus + (vplus - vminus) * rand( indxsSize );
1921	zNew = vNew ./ uNew;
1922	l = (zNew >= -m);
1923	z( indxs( l ) ) = zNew(l);
1924	u( indxs( l ) ) = uNew(l);
1925	indxs = indxs( ~l );
1926	end
1927
1928	out = z + m;
1929	indxs = find( log(u) > (log_1_over_pm + (h-1)/2log(out) - b/4(out + 1./out)) );
1930
1931	while ~isempty(indxs),
1932	indxsSize = size( indxs );
1933	u = rand( indxsSize );
1934	v = vminus + (vplus - vminus) * rand( indxsSize );
1935	z = v ./ u;
1936	clear('v');
1937	indxs1 = find( z < -m );
1938	while ~isempty(indxs1),
1939	indxsSize1 = size( indxs1 );
1940	uNew = rand( indxsSize1 );
1941	vNew = vminus + (vplus - vminus) * rand( indxsSize1 );
1942	zNew = vNew ./ uNew;
1943	l = (zNew >= -m);
1944	z( indxs1( l ) ) = zNew(l);
1945	u( indxs1( l ) ) = uNew(l);
1946	indxs1 = indxs1( ~l );
1947	end
1948
1949	outNew = z + m;
1950	l = ( log(u) <= (log_1_over_pm + (h-1)/2log(outNew) - b/4(outNew + 1./outNew)) );
1951	out( indxs(l) ) = outNew( l );
1952	indxs = indxs( ~l );
1953
1954	end
1955
1956	end %% if h<=1 & b<=1, else ...
1957
1958	out = sqrt( chi / psi ) * out;
1959
1960	case {'gh'}
1961	% START gh HELP
1962	% THE GENERALIZED HYPERBOLIC DISTRIBUTION
1963	% GH(lam, alpha, beta, mu, delta)
1964	%
1965	% pdf = (alpha^2-beta^2)^(lam/2) / (sqrt(2pi) alpha^(lam-1/2) * delta^lam * ...
1966	% besselk(lam, deltasqrt(alpha^2-beta^2) ) ) ...
1967	% (delta^2 + (y-mu).^2).^(1/2(lam-1/2)) . ...
1968	% besselk( lam-1/2, alphasqrt(delta^2 + (y-mu).^2) ) . ...
1969	% exp( beta*(y-mu) );
1970	%
1971	% Mean = mu + betadelta^2/(deltasqrt(alpha^2-beta^2)) * besselk(lam+1, delta*sqrt(alpha^2-beta^2) ) / ...
1972	% besselk(lam, delta*sqrt(alpha^2-beta^2) );
1973	% Variance = delta^2 * ( besselk(lam+1, zeta)/(zeta*besselk(lam, zeta)) + ...
1974	% beta^2delta^2/zeta^2 (besselk(lam+2, zeta)/besselk(lam, zeta) - ...
1975	% (besselk(lam+1, zeta)/besselk(lam, zeta))^2) );
1976	% where zeta = delta*sqrt(alpha^2-beta^2);
1977	%
1978	% PARAMETERS:
1979	% lam, -Inf < lam < Inf;
1980	% alpha - shape parameter (alpha>0) (steepness)
1981	% beta - 0 <= abs(beta) < alpha (skewness)
1982	% mu - location parameter (-Inf < mu < Inf)
1983	% delta - scale parameter (delta > 0)
1984	%
1985	% SUPPORT:
1986	% y, -Inf < y < Inf;
1987	%
1988	% CLASS:
1989	% Continuous skewed distributions
1990	%
1991	% USAGE:
1992	% randraw('gh', [lam, alpha, beta, mu, delta], sampleSize) - generate sampleSize number
1993	% of variates from the Generalized Hyperbolic distribution with
1994	% parameters [lam, alpha, beta, mu, delta]
1995	% randraw('gh') - help for the Generalized Hyperbolic distribution;
1996	%
1997	% EXAMPLES:
1998	% 1. y = randraw('gh', [3, 4, 3, 7.5, 1.5], [1 1e5]);
1999	% 2. y = randraw('gh', [3, 4, 3, 8.5, 1.5], 1, 1e5);
2000	% 3. y = randraw('gh', [3, 4, 3, 9.5, 1.5], 1e5 );
2001	% 4. y = randraw('gh', [3, 4, 3, 12.5, 1.5], [1e5 1] );
2002	% 5. randraw('gh');
2003	% END gh HELP
2004
2005	checkParamsNum(funcName, 'GH', 'gh', distribParams, [5]);
2006
2007	lam = distribParams(1);
2008	alpha = distribParams(2);
2009	beta = distribParams(3);
2010	mu = distribParams(4);
2011	delta = distribParams(5);
2012
2013
2014	validateParam(funcName, 'GH', 'gh', '[lam, alpha, beta, mu, delta]', 'alpha', alpha, {'> 0'});
2015	validateParam(funcName, 'GH', 'gh', '[lam, alpha, beta, mu, delta]', 'delta', delta, {'> 0'});
2016	validateParam(funcName, 'GH', 'gh', '[lam, alpha, beta, mu, delta]', 'alpha-abs(beta)', alpha-abs(beta), {'> 0'});
2017
2018	ygig = feval(funcName, 'gig', [lam, delta^2, alpha^2-beta^2], sampleSize);
2019
2020	out = mu + betaygig + sqrt(ygig).randn( sampleSize );
2021
2022	case {'gompertz'}
2023	% START gompertz HELP
2024	% THE GOMPERTZ DISTRIBUTION
2025	%
2026	% pdf(y) = b * c.^y * exp(-b*(c.^y-1)/log(c)); y>=0; b>0; c>1
2027	% cdf(y) = 1 - exp(-b*(c.^y-1)/log(c)); y>=0; b>0; c>1
2028	%
2029	% Mean = exp(b/log(c)) * (-1/log(c)) * (-expint(b/log(c)));
2030	% Variance =
2031	%
2032	% PARAMETERS:
2033	% b - shape parameter (b>0)
2034	% c - shape parameter (c>1)
2035	%
2036	% SUPPORT:
2037	% y, y>=0
2038	%
2039	% CLASS:
2040	% Continuous skewed distributions
2041	%
2042	% NOTES:
2043	% There are several forms given for the Gompertz distribution in the literature.
2044	% In particular, one common form uses the parameter alpha where alpha=log(c).
2045	%
2046	% The Gompertz distribution is frequently used by actuaries as a distribution
2047	% of length of life.
2048	%
2049	% USAGE:
2050	% randraw('gompertz', [b, c], sampleSize) - generate sampleSize number
2051	% of variates from Gompertz distribution with shape parameters 'b'
2052	% and 'c';
2053	% randraw('gompertz') - help for Gompertz distribution;
2054	%
2055	% EXAMPLES:
2056	% 1. y = randraw('gompertz', [1 2], [1 1e5]);
2057	% 2. y = randraw('gompertz', [2 3], 1, 1e5);
2058	% 3. y = randraw('gompertz', [1.1 5], 1e5 );
2059	% 4. y = randraw('gompertz', [2.3 4.8], [1e5 1] );
2060	% 5. randraw('gompertz');
2061	% END gompertz HELP
2062
2063	% Method:
2064	%
2065	% Inverse CDF transformation method.
2066	%
2067	% Reference:
2068	%
2069	% Dennis Kunimura, "The Compertz Distribution - Estimation of Parameters,"
2070	% Actuarial Research Clearing House, 1998, Vol.2
2071	checkParamsNum(funcName, 'Gompertz', 'gompertz', distribParams, [2]);
2072	b = distribParams(1);
2073	c = distribParams(2);
2074	validateParam(funcName, 'Gompertz', 'gompertz', '[b, c]', 'b', b, {'> 0'});
2075	validateParam(funcName, 'Gompertz', 'gompertz', '[b, c]', 'c', c, {'> 1'});
2076
2077	out = log( 1 - log(1-rand(sampleSize))*log(c)/b ) / log(c);
2078
2079	case {'halfcosine', 'hcosine', 'hcos'}
2080	% START halfcosine HELP
2081	% THE HALF-COSINE DISTRIBUTION
2082	%
2083	% Standard Half-cosine distribution:
2084	% pdf(y) = 1/4 * cos( y/2 ); -pi <= y <= pi
2085	% cdf(y) = 1/2 * ( 1 + sin( y/2 ) ); -pi <= y <= pi
2086	%
2087	% Mean = 0;
2088	% Variance = pi^2-8;
2089	%
2090	% General half-cosine distribution:
2091	% pdf(y) = pi/(4a) cos( pi(y-t)/(2s) ); t-s <= y <= t+s
2092	% cdf(y) = 1/2 * ( 1 + sin( pi(y-t)/(2s) ) ); t-s <= y <= t+s
2093	%
2094	% Mean = t;
2095	% Variance = (1-8/pi^2)*s^2;
2096	%
2097	% PARAMETERS:
2098	% t - location
2099	% s -scale; s>0
2100	%
2101	% SUPPORT:
2102	% y, -pi <= y <= pi - standard Half-cosine distribution
2103	% or
2104	% y, t-s <= y <= t+s - generalized Half-cosine distribution
2105	%
2106	% CLASS:
2107	% Continuous distributions
2108	%
2109	% USAGE:
2110	% randraw('halfcosine', [], sampleSize) - generate sampleSize number
2111	% of variates from standard Half-cosine distribution;
2112	% randraw('halfcosine', [t, s], sampleSize) - generate sampleSize number
2113	% of variates from generalized Half-cosine distribution
2114	% with location parameter 't' and scale parameter 's';
2115	% randraw('halfcosine') - help for Half-cosine distribution;
2116	%
2117	% EXAMPLES:
2118	% 1. y = randraw('halfcosine', [], [1 1e5]);
2119	% 2. y = randraw('halfcosine', [], 1, 1e5);
2120	% 3. y = randraw('halfcosine', [], 1e5 );
2121	% 4. y = randraw('halfcosine', [10 3], [1e5 1] );
2122	% 5. randraw('halfcosine');
2123	%
2124	% END halfcosine HELP
2125
2126	checkParamsNum(funcName, 'Half-Cosine', 'halfcosine', distribParams, [0, 2]);
2127	if numel(distribParams)==2
2128	t = distribParams(1);
2129	s = distribParams(2);
2130	validateParam(funcName, 'Half-Cosine', 'halfcosine', '[t, s]', 's', s, {'> 0'});
2131	else
2132	t = 0;
2133	s = pi;
2134	end
2135
2136	out = t + s2/piasin(2*rand(sampleSize)-1);
2137
2138	case {'hyperbolicsecant', 'hsecant', 'hsec'}
2139	% START hsecant HELP START hsec HELP START hyperbolicsecant HELP
2140	% THE HYPERBOLIC SECANT DISTRIBUTION
2141	%
2142	% Standard form of the Hyperbolic Secant distribution
2143	% pdf(y) = sech(y)/pi;
2144	% cdf(y) = 2*atan(exp(y))/pi;
2145	%
2146	% Mean = Median = Mode = 0;
2147	% Variance = pi^2/4;
2148	% Skewness = 0;
2149	% Kurtosis = 2;
2150	%
2151	% General form of the Hyperbolic Secant distribution
2152	% pdf(y) = sech((y-a)/b)/(b*pi);
2153	% cdf(y) = 2*atan(exp((y-a)/b))/pi;
2154	%
2155	% Mean = Median = Mode = a;
2156	% Variance = (pi*b)^2/4;
2157	% Skewness = 0;
2158	% Kurtosis = 2;
2159	%
2160	% PARAMETERS:
2161	% a - location;
2162	% b - scale; (b>0)
2163	%
2164	% SUPPORT:
2165	% y, -Inf < y < Inf
2166	%
2167	% CLASS:
2168	% Continuous symmetric distributions
2169	%
2170	% NOTES:
2171	% 1. The Hyperbolic Secant is related to the Logistic distribution.
2172	% 2. If Z ~Hyperbolic Secant, then W = exp(Z) ~Half Cauchy.
2173	% 3. The Hyperbolic Secant distribution is used in lifetime analysis.
2174	%
2175	% USAGE:
2176	% randraw('hsecant', [], sampleSize) - generate sampleSize number
2177	% of variates from standard Hyperbolic Secant distribution;
2178	% randraw('hsecant', [a, b], sampleSize) - generate sampleSize number
2179	% of variates from generalized Hyperbolic Secant distribution
2180	% with location parameter 'a' and scale parameter 'b';
2181	% randraw('hsecant') - help for Hyperbolic Secant distribution;
2182	%
2183	% EXAMPLES:
2184	% 1. y = randraw('hsecant', [], [1 1e5]);
2185	% 2. y = randraw('hsecant', [-1 4], 1, 1e5);
2186	% 3. y = randraw('hsecant', [], 1e5 );
2187	% 4. y = randraw('hsecant', [10.1 3.2], [1e5 1] );
2188	% 5. randraw('hsecant');
2189	% END hsecant HELP END hsec HELP END hyperbolicsecant HELP
2190
2191	checkParamsNum(funcName, 'Hyperbolic Secant', 'hsecant', distribParams, [0, 2]);
2192	if numel(distribParams)==2
2193	a = distribParams(1);
2194	b = distribParams(2);
2195	validateParam(funcName, 'Hyperbolic Secant', 'hsecant', '[a, b]', 'b', b, {'> 0'});
2196	else
2197	a = 0;
2198	b = 1;
2199	end
2200
2201	out = a + b* log(tan(pi/2*rand(sampleSize)));
2202
2203	case {'hypergeom', 'hypergeometric'}
2204	% START hypergeom HELP START hypergeometric HELP
2205	% THE HYPERGEOMETRIC DISTRIBUTION
2206	% If Y is the number of SUCCESSES in a completely random sample of size n drawn from
2207	% a population consisting of M SUCCESSES and (N-M) FAILURES, then Y distributed according
2208	% to Hypergeometric distribution
2209	%
2210	% pdf(y) = nchoosek(M,y)*nchoosek(N-M,n-y) / nchoosek(N,n) = ...
2211	% exp( gammaln(M+1) - gammaln(M-y+1) - gammaln(y+1) + ...
2212	% gammaln(N-M+1) - gammaln(N-M-n+y+1) - gammaln(n-y+1) + ...
2213	% gammaln(n+1) + gammaln(N-n+1) - gammaln(N+1) );
2214	%
2215	% max(0, n-N+M) <= y <= min(n,M)
2216	%
2217	% Mean = n*(M/N);
2218	% Variance = (N-n)/(N-1)nM/N*(1-M/N);
2219	% Mode = floor( (M+1)*(n+1)/(N+2) );
2220	%
2221	% PARAMETERS:
2222	% N, N = 2, 3, 4, ...
2223	% M, 0 < M < N
2224	% n, 0 < n < N
2225	%
2226	% SUPPORT:
2227	% y, y is integer and max(0, n-N+M) <= y <= min(n,M),
2228	%
2229	% CLASS:
2230	% Discrete distributions
2231	%
2232	% NOTES:
2233	% 1. In the urn model:
2234	% From an urn with white and black balls a random sample is drawn without
2235	% replacement, then
2236	% N = total number of balls in the urn;
2237	% M = number of white balls in the urn;
2238	% n = sample size (number of balls drawn without replacement);
2239	% Y = number of white balls in the sample.
2240	%
2241	% 2. When the population size is large (i.e. N is large) the hypergeometric
2242	% distribution can be approximated reasonably well with a binomial
2243	% distribution with parameters n (number of trials) and p = M / N
2244	% (probability of success in a single trial).
2245	%
2246	% USAGE:
2247	% randraw('hypergeom', [N, M, n], sampleSize) - generate sampleSize number
2248	% of variates from the Hypergeometric distribution
2249	% with parameters N, M and n;
2250	% randraw('hypergeom') - help for the Hypergeometric distribution;
2251	%
2252	% EXAMPLES:
2253	% 1. y = randraw('hypergeom', [20 13 17], [1 1e5]);
2254	% 2. y = randraw('hypergeom', [30 22 14], 1, 1e5);
2255	% 3. y = randraw('hypergeom', [50 3 22], 1e5 );
2256	% 4. y = randraw('hypergeom', [33 32 10], [1e5 1] );
2257	% 5. randraw('hypergeom');
2258	%
2259	% SEE ALSO:
2260	% BINOMIAL distribution
2261	% END hypergeom HELP END hypergeometric HELP
2262
2263	checkParamsNum(funcName, 'Hypergeometric', 'hypergeom', distribParams, [3]);
2264	N = distribParams(1);
2265	M = distribParams(2);
2266	n = distribParams(3);
2267	validateParam(funcName, 'Hypergeometric', 'hypergeom', '[N, M, n]', 'N', N, {'> 1', '==integer'});
2268	validateParam(funcName, 'Hypergeometric', 'hypergeom', '[N, M, n]', 'M', M, {'> 0', '==integer'});
2269	validateParam(funcName, 'Hypergeometric', 'hypergeom', '[N, M, n]', 'n', n, {'> 0', '==integer'});
2270	validateParam(funcName, 'Hypergeometric', 'hypergeom', '[N, M, n]', 'N-M', N-M, {'> 0'});
2271	validateParam(funcName, 'Hypergeometric', 'hypergeom', '[N, M, n]', 'N-n', N-n, {'> 0'});
2272
2273	mode = floor( (M+1)/(N+2)*(n+1) );
2274	if ~isfinite(mode)
2275	warning('Numeric Overflow !');
2276	varargout{1} = [];
2277	return;
2278	end
2279	pmode = exp( gammaln(M+1) - gammaln(M-mode+1) - gammaln(mode+1) + ...
2280	gammaln(N-M+1) - gammaln(N-M-n+mode+1) - gammaln(n-mode+1) + ...
2281	gammaln(n+1) + gammaln(N-n+1) - gammaln(N+1) );
2282
2283	if pmode < 5e-10
2284	varargout{1} = repmat(mode, sampleSize);
2285	return;
2286	end
2287
2288	if ~isfinite(pmode)
2289	varargout{1} = feval(funcName, 'binomial', [n M/N], sampleSize);
2290	end
2291	if pmode==1,
2292	varargout{1} = repmat(mode, sampleSize);
2293	return;
2294	end
2295	% nchoosek(M,y)*nchoosek(N-M,n-y) / nchoosek(N,n)
2296	t=pmode;
2297	ii = mode+1;
2298	while t*2147483648 > 1
2299	t = t * (M-ii+1)/(ii) *(n-ii+1)/(N-M-n+ii);
2300	ii = ii + 1;
2301	end
2302	last=ii-2;
2303
2304	t=pmode;
2305	j=-1;
2306	for ii=mode-1:-1:0
2307	t = t * (ii+1)/(M-ii) *(N-M-n+ii+1)/(n-ii);
2308	if t*2147483648 < 1
2309	j=ii;
2310	break;
2311	end
2312	end
2313	offset=j+1;
2314	sizeP=last-offset+1;
2315
2316	P = zeros(1, sizeP);
2317
2318	ii = (mode+1):last;
2319	P(mode-offset+1:last-offset+1) = round( 2^30pmodecumprod([1, (M-ii+1)./(ii).*(n-ii+1)./(N-M-n+ii)] ) );
2320
2321	ii = (mode-1):-1:offset;
2322	P(mode-offset:-1:1) = round( 2^30pmodecumprod((ii+1)./(M-ii).*(N-M-n+ii+1)./(n-ii)) );
2323
2324	out = randFrom5Tbls( P, offset, sampleSize);
2325
2326
2327	case {'ig', 'inversegauss', 'invgauss'}
2328	%
2329	% START ig HELP START inversegauss HELP START invgauss HELP
2330	% THE INVERSE GAUSSIAN DISTRIBUTION
2331	%
2332	% The Inverse Gaussian distribution is left skewed distribution whose
2333	% location is set by the mean with the profile determined by the
2334	% scale factor. The random variable can take a value between zero and
2335	% infinity. The skewness increases rapidly with decreasing values of
2336	% the scale parameter.
2337	%
2338	%
2339	% pdf(y) = sqrt(chi/(2piy^3)) * exp(-chi./(2y).(y/theta-1).^2);
2340	% cdf(y) = normcdf(sqrt(chi./y).*(y/theta-1)) + ...
2341	% exp(2chi/theta)normcdf(sqrt(chi./y).*(-y/theta-1));
2342	%
2343	% where normcdf(x) = 0.5*(1+erf(y/sqrt(2))); is the standard normal CDF
2344	%
2345	% Mean = theta;
2346	% Variance = theta^3/chi;
2347	% Skewness = sqrt(9*theta/chi);
2348	% Kurtosis = 15*mean/scale;
2349	% Mode = theta/(2chi)(sqrt(9theta^2+4chi^2)-3*theta);
2350	%
2351	% PARAMETERS:
2352	% theta - location; (theta>0)
2353	% chi - scale; (chi>0)
2354	%
2355	% SUPPORT:
2356	% y, y>0
2357	%
2358	% CLASS:
2359	% Continuous skewed distribution
2360	%
2361	% NOTES:
2362	% 1. There are several alternate forms for the PDF,
2363	% some of which have more than two parameters
2364	% 2. The Inverse Gaussian distribution is often called the Inverse Normal
2365	% 3. Wald distribution is a special case of The Inverse Gaussian distribution
2366	% where the mean is a constant with the value one.
2367	% 4. The Inverse Gaussian distribution is a special case of The Generalized
2368	% Hyperbolic Distribution
2369	%
2370	% USAGE:
2371	% randraw('ig', [theta, chi], sampleSize) - generate sampleSize number
2372	% of variates from the Inverse Gaussian distribution with
2373	% parameters theta and chi;
2374	% randraw('ig') - help for the Inverse Gaussian distribution;
2375	%
2376	% EXAMPLES:
2377	% 1. y = randraw('ig', [0.1, 1], [1 1e5]);
2378	% 2. y = randraw('ig', [3.2, 10], 1, 1e5);
2379	% 3. y = randraw('ig', [100.2, 6], 1e5 );
2380	% 4. y = randraw('ig', [10, 10.5], [1e5 1] );
2381	% 5. randraw('ig');
2382	%
2383	% SEE ALSO:
2384	% WALD distribution
2385	% END ig HELP END inversegauss HELP END invgauss HELP
2386
2387	% Method:
2388	%
2389	% There is an efficient procedure that utilizes a transformation
2390	% yielding two roots.
2391	% If Y is Inverse Gauss random variable, then following to [1]
2392	% we can write:
2393	% V = chi(Y-theta)^2/(Ytheta^2) ~ Chi-Square(1),
2394	%
2395	% i.e. V is distributed as a chi-square random variable with
2396	% one degree of freedom.
2397	% So it can be simply generated by taking a square of a
2398	% standard normal random number.
2399	% Solving this equation for Y yields two roots:
2400	%
2401	% y1 = theta + 0.5theta/chi ( thetaV - sqrt(4thetachiV + ...
2402	% theta^2*V.^2) );
2403	% and
2404	% y2 = theta^2/y1;
2405	%
2406	% In [2] showed that Y can be simulated by choosing y1 with probability
2407	% theta/(theta+y1) and y2 with probability 1-theta/(theta+y1)
2408	%
2409	% References:
2410	% [1] Shuster, J. (1968). On the Inverse Gaussian Distribution Function,
2411	% Journal of the American Statistical Association 63: 1514-1516.
2412	%
2413	% [2] Michael, J.R., Schucany, W.R. and Haas, R.W. (1976).
2414	% Generating Random Variates Using Transformations with Multiple Roots,
2415	% The American Statistician 30: 88-90.
2416	%
2417	%
2418
2419	checkParamsNum(funcName, 'Inverse Gaussian', 'ig', distribParams, [2]);
2420	theta = distribParams(1);
2421	chi = distribParams(2);
2422	validateParam(funcName, 'Inverse Gaussian', 'ig', '[theta, chi]', 'theta', theta, {'> 0'});
2423	validateParam(funcName, 'Inverse Gaussian', 'ig', '[theta, chi]', 'chi', chi, {'> 0'});
2424
2425	chisq1 = randn(sampleSize).^2;
2426	out = theta + 0.5theta/chi ( theta*chisq1 - ...
2427	sqrt(4thetachichisq1 + theta^2chisq1.^2) );
2428
2429	l = rand(sampleSize) >= theta./(theta+out);
2430	out( l ) = theta^2./out( l );
2431
2432	case {'laplace' 'doubleexponential', 'doubleexp', 'bilateralexponential', 'bilateralexp'}
2433	% START laplace HELP START doubleexponential HELP START doubleexp HELP START bilateralexponential HELP START bilateralexp HELP
2434	% THE LAPLACE DISTRIBUTION
2435	% (sometimes: double-exponential or bilateral exponential distribution)
2436	%
2437	% pdf = 1/(2lam)exp(-abs(y-theta)/lam);
2438	% cdf = (y<=theta) .* 1/2exp((y-theta)/lam) + (y>theta) . (1 - 1/2*exp((theta-y)/lam));
2439	%
2440	% Mean = Median = Mode = theta;
2441	% Variance = 2*lam^2;
2442	% Skewness = 0;
2443	% Kurtosis = 3;
2444	%
2445	% PARAMETERS:
2446	% theta - location
2447	% lam - scale (lam>0)
2448	%
2449	% SUPPORT:
2450	% y , -Inf < y < Inf
2451	%
2452	% CLASS:
2453	% Continuous symmetric distribution
2454	%
2455	% USAGE:
2456	% randraw('laplace', [], sampleSize) - generate sampleSize number
2457	% of variates from the Laplace distribution with
2458	% loaction parameter theta=0 and scale parameter lam=1;
2459	% randraw('laplace', [theta, lam], sampleSize) - generate sampleSize number
2460	% of variates from the Laplace distribution with
2461	% loaction parameter theta and scale parameter lam;
2462	% randraw('laplace') - help for the Laplace distribution;
2463	%
2464	% EXAMPLES:
2465	% 1. y = randraw('laplace', [0, 1], [1 1e5]);
2466	% 2. y = randraw('laplace', [3.2, 10], 1, 1e5);
2467	% 3. y = randraw('laplace', [100.2, 6], 1e5 );
2468	% 4. y = randraw('laplace', [10, 10.5], [1e5 1] );
2469	% 5. randraw('laplace');
2470	% END laplace HELP END doubleexponential HELP END doubleexp HELP END bilateralexponential HELP END bilateralexp HELP
2471
2472	checkParamsNum(funcName, 'Laplace', 'laplace', distribParams, [0, 2]);
2473	if numel(distribParams)==2
2474	theta = distribParams(1);
2475	lam = distribParams(2);
2476	validateParam(funcName, 'Laplace', 'laplace', '[theta, lam]', 'lam', lam, {'> 0'});
2477	else
2478	theta = 0;
2479	lam = 1;
2480	end
2481
2482	u = rand( sampleSize );
2483	out = zeros( sampleSize );
2484	out(u<=0.5) = theta + lamlog(2u(u<=0.5));
2485	out(u>0.5) = theta - lamlog(2(1-u(u>0.5)));
2486
2487	case {'logistic'}
2488	% START logistic HELP
2489	% THE LOGISTIC DISTRIBUTION
2490	% The logistic distribution is a symmetrical bell shaped distribution.
2491	% One of its applications is an alternative to the Normal distribution
2492	% when a higher proportion of the population being modeled is
2493	% distributed in the tails.
2494	%
2495	% pdf(y) = exp((y-a)/k)./(k*(1+exp((y-a)/k)).^2);
2496	% cdf(y) = 1 ./ (1+exp(-(y-a)/k))
2497	%
2498	% Mean = a;
2499	% Variance = k^2*pi^2/3;
2500	% Skewness = 0;
2501	% Kurtosis = 1.2;
2502	%
2503	% PARAMETERS:
2504	% a - location;
2505	% k - scale (k>0);
2506	%
2507	% SUPPORT:
2508	% y, -Inf < y < Inf
2509	%
2510	% CLASS:
2511	% Continuous symmetric distribution
2512	%
2513	% USAGE:
2514	% randraw('logistic', [], sampleSize) - generate sampleSize number
2515	% of variates from the standard Logistic distribution with
2516	% loaction parameter a=0 and scale parameter k=1;
2517	% randraw('logistic', [a, k], sampleSize) - generate sampleSize number
2518	% of variates from the Logistic distribution with
2519	% loaction parameter 'a' and scale parameter 'k';
2520	% randraw('logistic') - help for the Logistic distribution;
2521	%
2522	% EXAMPLES:
2523	% 1. y = randraw('logistic', [], [1 1e5]);
2524	% 2. y = randraw('logistic', [0, 4], 1, 1e5);
2525	% 3. y = randraw('logistic', [-1, 10.2], 1e5 );
2526	% 4. y = randraw('logistic', [3.2, 0.3], [1e5 1] );
2527	% 5. randraw('logistic');
2528	% END logistic HELP
2529
2530	% Method:
2531	%
2532	% Inverse CDF transformation method.
2533
2534	checkParamsNum(funcName, 'Logistic', 'logistic', distribParams, [0, 2]);
2535	if numel(distribParams)==2
2536	a = distribParams(1);
2537	k = distribParams(2);
2538	validateParam(funcName, 'Laplace', 'laplace', '[a, k]', 'k', k, {'> 0'});
2539	else
2540	a = 0;
2541	k = 1;
2542	end
2543
2544	u1 = rand( sampleSize );
2545	out = a - k*log( 1./u1 - 1 );
2546
2547	case { 'lognorm', 'lognormal', 'cobbdouglas', 'antilognormal' }
2548	% START lognorm HELP START lognormal HELP START cobbdouglas HELP START antilognormal HELP
2549	% THE LOG-NORMAL DISTRIBUTION
2550	% (sometimes: Cobb-Douglas or antilognormal distribution)
2551	%
2552	% pdf = 1/(ksqrt(2pi)) * exp(-1/2*((log(y)-a)/k)^2)
2553	% cdf = 1/2(1 + erf((log(y)-a)/(ksqrt(2))));
2554	%
2555	% Mean = exp( a + k^2/2 );
2556	% Variance = exp(2a+k^2)( exp(k^2)-1 );
2557	% Skewness = (exp(1)+2)*sqrt(exp(1)-1), for a=0 and k=1;
2558	% Kurtosis = exp(4) + 2exp(3) + 3exp(2) - 6; for a=0 and k=1;
2559	% Mode = exp(a-k^2);
2560	%
2561	% PARAMETERS:
2562	% a - location
2563	% k - scale (k>0)
2564	%
2565	% SUPPORT:
2566	% y, y>0
2567	%
2568	% CLASS:
2569	% Continuous skewed distribution
2570	%
2571	% NOTES:
2572	% 1) The LogNormal distribution is always right-skewed
2573	% 2) Parameters a and k are the mean and standard deviation
2574	% of y in (natural) log space.
2575	%
2576	% USAGE:
2577	% randraw('lognorm', [], sampleSize) - generate sampleSize number
2578	% of variates from the standard Lognormal distribution with
2579	% loaction parameter a=0 and scale parameter k=1;
2580	% randraw('lognorm', [a, k], sampleSize) - generate sampleSize number
2581	% of variates from the Lognormal distribution with
2582	% loaction parameter 'a' and scale parameter 'k';
2583	% randraw('lognorm') - help for the Lognormal distribution;
2584	%
2585	% EXAMPLES:
2586	% 1. y = randraw('lognorm', [], [1 1e5]);
2587	% 2. y = randraw('lognorm', [0, 4], 1, 1e5);
2588	% 3. y = randraw('lognorm', [-1, 10.2], 1e5 );
2589	% 4. y = randraw('lognorm', [3.2, 0.3], [1e5 1] );
2590	% 5. randraw('lognorm');
2591	%END lognorm HELP END lognormal HELP END cobbdouglas HELP END antilognormal HELP
2592
2593	checkParamsNum(funcName, 'Lognormal', 'lognorm', distribParams, [0, 2]);
2594	if numel(distribParams)==2
2595	a = distribParams(1);
2596	k = distribParams(2);
2597	validateParam(funcName, 'Lognormal', 'lognorm', '[a, k]', 'k', k, {'> 0'});
2598	else
2599	a = 0;
2600	k = 1;
2601	end
2602
2603	out = exp( a + k * randn( sampleSize ) );
2604
2605	case {'maxwell'}
2606	% START maxwell HELP
2607	% THE MAXWELL DISTRIBUTION
2608	%
2609	% pdf(y) = 1/a^3 * sqrt(2/pi) * y.^2 * exp(-y.^2/(2*a^2)); a > 0, y >= 0
2610	% cdf(y) = gammainc(3/2, y.^2/(2*a^2));
2611	%
2612	% Mean = 2asqrt(2/pi);
2613	% Variance = a^2*(3-8/pi);
2614	% Skewness = 2(16/pi-5)sqrt(2/pi) / (3-8/pi)^(3/2) = 0.48569282804959
2615	% Kurtosis = (15-8/pi)/(3-8/pi)^2 - 3 ???
2616	%
2617	% PARAMETERS:
2618	% a - scale parameter (a > 0)
2619	%
2620	% SUPPORT:
2621	% y, y >= 0
2622	%
2623	% CLASS:
2624	% Continuous skewed distribution
2625	%
2626	% NOTES:
2627	% The distribution of speeds of molecules in thermal equilibrium as given by
2628	% statistical mechanics and named after the famous scottish physicist James
2629	% Clerk Maxwell (1831-1879).
2630	%
2631	% USAGE:
2632	% randraw('maxwell', a, sampleSize) - generate sampleSize number
2633	% of variates from the Maxwell distribution with
2634	% scale parameter 'a';
2635	% randraw('maxwell') - help for the Maxwell distribution;
2636	%
2637	% EXAMPLES:
2638	% 1. y = randraw('maxwell', 1.1, [1 1e5]);
2639	% 2. y = randraw('maxwell', 0.5, 1, 1e5);
2640	% 3. y = randraw('maxwell', 10, 1e5 );
2641	% 4. y = randraw('maxwell', 5.5, [1e5 1] );
2642	% 5. randraw('maxwell');
2643	% END maxwell HELP
2644
2645	checkParamsNum(funcName, 'Maxwell', 'maxwell', distribParams, [1]);
2646	a = distribParams(1);
2647	validateParam(funcName, 'Maxwell', 'maxwell', 'a', 'a', a, {'> 0'});
2648
2649	out = sqrt( randn(sampleSize).^2 + randn(sampleSize).^2 + ...
2650	randn(sampleSize).^2 ) * a;
2651
2652	case {'negativebinomial', 'negbinomial', 'negbinom', 'pascal'}
2653	% START negbinom HELP START pascal HELP
2654	% THE NEGATIVE BINOMIAL DISTRIBUTION
2655	% ( sometimes: Pascal distribution )
2656	%
2657	% Negative Binomial (also known as the Pascal or Polya) distribution
2658	% gives the probability of r-1 successes and y failures in y+r-1 trials
2659	% and success on the (y+r)'th trial (if r is positive integer )
2660	%
2661	% pdf(y) = gamma(r+y)./(gamma(y+1)gamma(r)) p^r * (1-p)^y = ...
2662	% exp( gammaln(r+y) - gammaln(y+1) - gammaln(r) + rlog(p) + ylog(1-p) );
2663	% y>=0
2664	%
2665	% Mode = (r>1)floor( (r-1)(1-p)/p ) + (r<=1)*0;
2666	% Mean = r*(1-p)/p;
2667	% Variance = r*(1-p)/p^2;
2668	% Skewness = (2-p) / sqrt(r*(1-p));
2669	% Kurtosis = (p^2-6p+6)/(r(1-p));
2670	%
2671	% PARAMETERS:
2672	% r>0;
2673	% p - probability of success in a single trial ( 0< p <1 );
2674	%
2675	% SUPPORT:
2676	% y = 0, 1, 2, 3, ...
2677	%
2678	% CLASS:
2679	% Discrete distribution
2680	%
2681	% USAGE:
2682	% randraw('negbinom', [r, p], sampleSize) - generate sampleSize number
2683	% of variates from the Negative Binomial distribution with
2684	% parameters 'r' and 'p';
2685	% randraw('negbinom') - help for the Negative Binomial distribution;
2686	%
2687	% EXAMPLES:
2688	% 1. y = randraw('negbinom', [10 0.2], [1 1e5]);
2689	% 2. y = randraw('negbinom', [100, 0.9], 1, 1e5);
2690	% 3. y = randraw('negbinom', [20 0.1], 1e5 );
2691	% 4. y = randraw('negbinom', [30 0.99], [1e5 1] );
2692	% 5. randraw('negbinom');
2693	% END negbinom HELP END pascal HELP
2694
2695	% Method:
2696	%
2697	% We implemented Condensed Table-Lookup method suggested in
2698	% George Marsaglia, "Fast Generation Of Discrete Random Variables,"
2699	% Journal of Statistical Software, July 2004, Volume 11, Issue 4
2700
2701	% pdf = exp( gammaln(r+y) - gammaln(y+1) - gammaln(r) + rlog(p) + ylog(1-p) );
2702
2703	checkParamsNum(funcName, 'Negative Binomial', 'negbinom', distribParams, [2]);
2704	r = distribParams(1);
2705	validateParam(funcName, 'Negative Binomial', 'negbinom', '[r, p]', 'r', r, {'> 0'});
2706	p = distribParams(2);
2707	validateParam(funcName, 'Negative Binomial', 'negbinom', '[r, p]', 'p', p, {'> 0','< 1'});
2708
2709	q = 1 - p;
2710
2711	if r*q/p^2 > 1e8
2712	out = 1/p * feval(funcName, 'gamma', r*(1-p), sampleSize);
2713	else
2714	mode = (r>1)floor( (r-1)(1-p)/p ) + (r<=1)*0;
2715	pmode = exp( gammaln(r+mode) - gammaln(mode+1) - gammaln(r) + ...
2716	rlog(p) + modelog(1-p) );
2717
2718	t=pmode;
2719	ii = mode+1;
2720	while t*2147483648 > 1
2721	t = t * (r+ii-1)/ii * q;
2722	ii = ii + 1;
2723	end
2724	last=ii-2;
2725
2726	t=pmode;
2727	j=-1;
2728	for ii=mode-1:-1:0
2729	t = t * (ii+1)/(r+ii)/q;
2730	if t*2147483648 < 1
2731	j=ii;
2732	break;
2733	end
2734	end
2735	offset=j+1;
2736	sizeP=last-offset+1;
2737
2738	P = zeros(1, sizeP);
2739
2740	ii = (mode+1):last;
2741	P(mode-offset+1:last-offset+1) = round( 2^30pmodecumprod([1, (r+ii-1)./ii * q] ) );
2742
2743	ii = (mode-1):-1:offset;
2744	P(mode-offset:-1:1) = round( 2^30pmodecumprod((ii+1)./(r+ii)/q) );
2745
2746	out = randFrom5Tbls( P, offset, sampleSize);
2747	end
2748
2749	case {'normal', 'gaussian', 'gauss', 'norm'} % Gaussian distribution
2750	% START normal HELP START gaussian HELP START gauss HELP START norm HELP
2751	% THE NORMAL DISTRIBUTION
2752	% Standard form of the Normal distribution:
2753	% pdf(y) = 1/sqrt(2pi) exp(-1/2*y.^2);
2754	% cdf(y) = 0.5*(1+erf(y/sqrt(2)));
2755	%
2756	% Mean = 0;
2757	% Variance = 1;
2758	%
2759	% General form of the Normal distribution:
2760	% pdf(y) = 1/(sigmasqrt(2pi)) * exp(-1/2*((y-mu)/sigma).^2);
2761	% cdf(y) = 1/2(1+erf((y-mu)/(sigmasqrt(2))));
2762	%
2763	% Mean = mu;
2764	% Variance = sigma^2;
2765	%
2766	% PARAMETERS:
2767	% mu - location (mean)
2768	% sigma>0 - scale (std)
2769	%
2770	% SUPPORT:
2771	% y, -Inf < y < Inf
2772	%
2773	% CLASS:
2774	% Continuous symmetric distributions
2775	%
2776	% USAGE:
2777	% randraw('norm', [], sampleSize) - generate sampleSize number
2778	% of variates from the standard Normal distribution;
2779	% randraw('norm', [a, k], sampleSize) - generate sampleSize number
2780	% of variates from the Normal distribution with
2781	% mean 'mu' and std 'sigma';
2782	% randraw('norm') - help for the Lognormal distribution;
2783	%
2784	% EXAMPLES:
2785	% 1. y = randraw('norm', [], [1 1e5]);
2786	% 2. y = randraw('norm', [0, 4], 1, 1e5);
2787	% 3. y = randraw('norm', [-1, 10.2], 1e5 );
2788	% 4. y = randraw('norm', [3.2, 0.3], [1e5 1] );
2789	% 5. randraw('norm');
2790	% END normal HELP END gaussian HELP END gauss HELP END norm HELP
2791
2792	checkParamsNum(funcName, 'Normal', 'normal', distribParams, [0, 2]);
2793
2794	if numel(distribParams)==2
2795	mu = distribParams(1);
2796	sigma = distribParams(2);
2797	validateParam(funcName, 'Normal', 'normal', '[mu, sigma]', 'sigma', sigma, {'> 0'});
2798	else
2799	mu = 0;
2800	sigma = 1;
2801	end
2802
2803	out = mu + sigma * randn( sampleSize );
2804
2805	case {'normaltrunc', 'normaltruncated', 'gausstrunc'}
2806	% START normaltrunc HELP START normaltruncated HELP START gausstrunc HELP
2807	% THE TRUNCATED NORMAL DISTRIBUTION
2808	%
2809	% pdf(y) = normpdf((y-mu)/sigma) / (sigma*(normcdf((b-mu)/sigma)-normcdf((a-mu)/sigma))); a<=y<=b;
2810	% cdf(y) = (normcdf((y-mu)/sigma)-normcdf((a-mu)/sigma)) / (normcdf((b-mu)/sigma)-normcdf((a-mu)/sigma)); a<=y<=b;
2811	% where mu and sigma are the mean and standard deviation of the parent normal
2812	% distribution and a and b are the lower and upper truncation points.
2813	% normpdf and normcdf are the PDF and CDF for the standard normal distribution respectvely
2814	% ( run randraw('normal') for help).
2815	%
2816	%
2817	% PARAMETERS:
2818	% a - lower truncation point;
2819	% b - upper truncation point; (b>=a)
2820	% mu - Mean of the parent normal distribution
2821	% sigma - standard deviation of the parent normal distribution (sigma>0)
2822	%
2823	%
2824	% SUPPORT:
2825	% y, a <= y <= b
2826	%
2827	% CLASS:
2828	% Continuous distributions
2829	%
2830	% USAGE:
2831	% randraw('normaltrunc', [a, b, mu, sigma], sampleSize) - generate sampleSize number
2832	% of variates from Truncated Normal distribution on the interval (a, b) with
2833	% parameters 'mu' and 'sigma';
2834	% randraw('normaltrunc') - help for Truncated Normal distribution;
2835	%
2836	% EXAMPLES:
2837	% 1. y = randraw('normaltrunc', [0, 1, 0, 1], [1 1e5]);
2838	% 2. y = randraw('normaltrunc', [0, 1, 10, 3], 1, 1e5);
2839	% 3. y = randraw('normaltrunc', [-10, 10, 0, 1], 1e5 );
2840	% 4. y = randraw('normaltrunc', [-13.1, 15.2, 20.1, 3.3], [1e5 1] );
2841	% 5. randraw('normaltrunc');
2842	% END normaltrunc HELP END normaltruncated HELP END gausstrunc HELP
2843
2844	checkParamsNum(funcName, 'Truncated Normal', 'normaltrunc', distribParams, [4]);
2845
2846	a = distribParams(1);
2847	b = distribParams(2);
2848	mu = distribParams(3);
2849	sigma = distribParams(4);
2850	validateParam(funcName, 'Truncated Normal', 'normaltrunc', '[a, b, mu, sigma]', 'b-a', b-a, {'>=0'});
2851	validateParam(funcName, 'Truncated Normal', 'normaltrunc', '[a, b, mu, sigma]', 'sigma', sigma, {'> 0'});
2852
2853	PHIl = normcdf((a-mu)/sigma);
2854	PHIr = normcdf((b-mu)/sigma);
2855
2856	out = mu + sigma( sqrt(2)erfinv(2(PHIl+(PHIr-PHIl)rand(sampleSize))-1) );
2857
2858	case {'nig'}
2859	% START nig HELP
2860	% THE NORMAL INVERSE GAUSSIAN (NIG) DISTRIBUTION
2861	% NIG(alpha, beta, mu, delta)
2862	%
2863	% Heavy-tailed distributions such as the normal inverse Gauss distribution
2864	% (NIG) play a prominent role in the statistical analysis of economic time-series.
2865	% A number of empirical studies have shown that the marginal distribution of the
2866	% daily returns of liquid shares are NIG.
2867	%
2868	% The NIG density is a variance-mean mixture of a Gaussian density with
2869	% an inverse Gaussian.
2870	% The shape of the NIG density is specified by the four-dimensional
2871	% parameter vector [alpha, beta , mu, delta]. The rich parametrization
2872	% makes the NIG density a suitable model for a variety of unimodal positive
2873	% kurtotic data. The alpha-parameter controls the steepness or pointiness
2874	% of the density, which increases monotonically with increasing alpha.
2875	% A large alpha implies light tails, a small value implies heavy tails.
2876	% The beta-parameter controls the skewness. For beta<0, the density is skewed to
2877	% the left, for beta>0 the density is skewed to the right, while
2878	% beta=0 implies a symmetric density around mu, which is a centrality parameter.
2879	% The delta-parameter is a scale-like parameter.
2880	%
2881	% pdf(y; alpha, beta, mu, delta) = ...
2882	% alpha/pi * exp(deltasqrt(alpha^2-beta^2) - betamu) * ...
2883	% 1/sqrt(1+((y-mu)/delta).^2) .* ...
2884	% besselk(1, alphadeltasqrt(1+((y-mu)/delta).^2)) .*...
2885	% exp(beta*y);
2886	%
2887	% Mean = mu+beta*delta/sqrt(alpha^2-beta^2);
2888	% Variance = delta * (alpha^2 / sqrt(alpha^2 - beta^2)^3);
2889	% Skewness = 3beta/(alphasqrt(delta*sqrt(alpha^2 - beta^2)));
2890	% Kurtosis = 3(1 + 4(beta/alpha)^2)/(delta*sqrt(alpha^2 - beta^2));
2891	%
2892	% PARAMETERS:
2893	% alpha, alpha > 0
2894	% beta, 0 <= abs(beta) < alpha
2895	% mu, -Inf < mu < Inf
2896	% delta, delta > 0
2897	%
2898	% SUPPORT:
2899	% y, -Inf < y < Inf
2900	%
2901	% CLASS:
2902	% Continuous skewed distribution
2903	%
2904	% USAGE:
2905	% randraw('nig', [alpha, beta, mu, delta], sampleSize) - generate sampleSize
2906	% number of variates from NIG distribution with parameters
2907	% alpha, beta, mu and delta
2908	% randraw('nig') - help for NIG distribution;
2909	%
2910	% EXAMPLES:
2911	% 1. y = randraw('nig', [2, 1, 4, 2], [1 1e5]);
2912	% 2. y = randraw('nig', [2, 1, 4, 2], 1, 1e5);
2913	% 3. y = randraw('nig', [2, 1, 4, 2], 1e5 );
2914	% 4. y = randraw('nig', [2, 1, 4, 2], [1e5 1] );
2915	% 5. randraw('nig');
2916	%
2917	% SEE ALSO:
2918	% INVERSE GAUSSIAN distribution
2919	% END nig HELP
2920
2921	% REFERENCES:
2922	% 1. http://www.quantlet.com/mdstat/scripts/csa/html/node236.html
2923	% 2. http://www.anst.uu.se/larsfors/APRv1_5.pdf
2924	% 3. http://ica2001.ucsd.edu/index_files/pdfs/048-jenssen.pdf
2925	% 4. http://www.freidok.uni-freiburg.de/volltexte/15/pdf/15_1.pdf
2926
2927
2928	checkParamsNum(funcName, 'NIG', 'nig', distribParams, [4]);
2929
2930	alpha = distribParams(1);
2931	beta = distribParams(2);
2932	mu = distribParams(3);
2933	delta = distribParams(4);
2934
2935	validateParam(funcName, 'NIG', 'nig', '[alpha, beta, mu, delta]', 'alpha', alpha, {'> 0'});
2936	validateParam(funcName, 'NIG', 'nig', '[alpha, beta, mu, delta]', 'delta', delta, {'> 0'});
2937	validateParam(funcName, 'NIG', 'nig', '[alpha, beta, mu, delta]', 'alpha-abs(beta)', alpha-abs(beta), {'> 0'});
2938
2939	invGaussY = feval(funcName, 'ig', [delta/sqrt(alpha^2-beta^2), delta^2], sampleSize);
2940	out = mu + betainvGaussY + sqrt(invGaussY).randn(sampleSize);
2941
2942	case {'pareto'}
2943	% START pareto HELP
2944	% Pareto or "power law" distribution, used in the analysis of financial data
2945	% and critical behavior
2946	%
2947	% pdf = a*k^a ./ y.^(a+1);
2948	% cdf = 1 - (k./y).^a;
2949	%
2950	% Mean = k*a/(a-1);
2951	% Variance = k^2a/((a-2)(a-1)^2);
2952	% Skewness = 2(a+1)sqrt(a-2)/((a-3)*sqrt(a));
2953	% Kurtosis = 6(a^3+a^2-6a-2)/(a(a^2-7a+12));
2954	%
2955	% PARAMETERS:
2956	% k - location parameter (k>0)
2957	% a - shape parameter (a>0)
2958	%
2959	% SUPPORT:
2960	% y, y > k
2961	%
2962	% CLASS:
2963	% Continuous skewed distributions
2964	%
2965	% USAGE:
2966	% randraw('pareto', [k, a], sampleSize) - generate sampleSize
2967	% number of variates from the Pareto distribution with parameters
2968	% 'k' and 'a'
2969	% randraw('pareto') - help for the Pareto distribution;
2970	%
2971	% EXAMPLES:
2972	% 1. y = randraw('pareto', [1, 2], [1 1e5]);
2973	% 2. y = randraw('pareto', [3, 8], 1, 1e5);
2974	% 3. y = randraw('pareto', [0.5, 2.4], 1e5 );
2975	% 4. y = randraw('pareto', [3.5, 4.5], [1e5 1] );
2976	% 5. randraw('pareto');
2977	% END pareto HELP
2978
2979	checkParamsNum(funcName, 'Pareto', 'pareto', distribParams, [2]);
2980
2981	k = distribParams(1);
2982	a = distribParams(2);
2983	validateParam(funcName, 'Pareto', 'pareto', '[k, a]', 'k', k, {'> 0'});
2984	validateParam(funcName, 'Pareto', 'pareto', '[k, a]', 'a', a, {'> 0'});
2985
2986	out = k*rand( sampleSize ).^(-1/a);
2987
2988	case {'pareto2', 'lomax'}
2989	% START pareto2 HELP
2990	% THE PARETO DISTRIBUTION OF THE SECOND TYPE
2991	% (sometimes Lomax distribution )
2992	%
2993	% pdf = b*k^b ./ (k+y).^(b+1); b>0, y>0;
2994	% cdf = 1 - k^b./(k+y).^b;
2995	%
2996	% PARAMETERS:
2997	% k - location parameters (k>0)
2998	% b - shape parameters (b>0)
2999	%
3000	% SUPPORT:
3001	% y, y>0
3002	%
3003	% CLASS;
3004	% Continuous skewed distribution
3005	%
3006	% USAGE:
3007	% randraw('pareto2', [k, b], sampleSize) - generate sampleSize
3008	% number of variates from the Pareto Second Type distribution with parameters
3009	% 'k' and 'b'
3010	% randraw('pareto2') - help for the Pareto Second Type distribution;
3011	%
3012	% EXAMPLES:
3013	% 1. y = randraw('pareto2', [1, 2], [1 1e5]);
3014	% 2. y = randraw('pareto2', [3, 8], 1, 1e5);
3015	% 3. y = randraw('pareto2', [0.5, 2.4], 1e5 );
3016	% 4. y = randraw('pareto2', [3.5, 4.5], [1e5 1] );
3017	% 5. randraw('pareto2');
3018	% END pareto2 HELP
3019
3020	checkParamsNum(funcName, 'Pareto Second Type', 'pareto2', distribParams, [2]);
3021
3022	k = distribParams(1);
3023	b = distribParams(2);
3024	validateParam(funcName, 'Pareto Second Type', 'pareto2', '[k, b]', 'k', k, {'> 0'});
3025	validateParam(funcName, 'Pareto Second Type', 'pareto2', '[k, b]', 'b', b, {'> 0'});
3026
3027	out = k*(1 - rand( sampleSize )).^(-1/b) - k;
3028
3029	case 'planck'
3030	% START planck HELP
3031	% THE PLANCK DISTRIBUTION
3032	% The Planck distribution widely used in Physics.
3033	%
3034	% The Planck distribution ia a two parameter distribution:
3035	% pdf(y) = b^(a+1)/(gamma(a+1)zeta(a+1)) y^a/(exp(b*y)-1);
3036	% where zeta(c) is the Riemann zeta function defined as
3037	% zeta(c) = sum from k=1 to Inf of 1/k^c.
3038	%
3039	% PARAMETERS:
3040	% a > 0 is a shape parameter
3041	% b > 0 is a scale parameter
3042	%
3043	% SUPPORT:
3044	% y, y>0
3045	%
3046	% CLASS:
3047	% Continuous skewed distributions
3048	%
3049	% USAGE:
3050	% randraw('planck', [a, b], sampleSize) - generate sampleSize
3051	% number of variates from the Planck distribution with parameters
3052	% 'a' and 'b'
3053	% randraw('planck') - help for the Planck distribution;
3054	%
3055	% EXAMPLES:
3056	% 1. y = randraw('planck', [1, 2], [1 1e5]);
3057	% 2. y = randraw('planck', [3, 8], 1, 1e5);
3058	% 3. y = randraw('planck', [0.5, 2.4], 1e5 );
3059	% 4. y = randraw('planck', [3.5, 4.5], [1e5 1] );
3060	% 5. randraw('planck');
3061	% END planck HELP
3062
3063	% Reference:
3064	% Luc Devroye, "Non-Uniform Random Variate Generation,"
3065	% Springer 1986, 850p. 3-540-96305-7
3066
3067	checkParamsNum(funcName, 'Planck', 'planck', distribParams, [2]);
3068
3069	a = distribParams(1);
3070	b = distribParams(2);
3071	validateParam(funcName, 'Planck', 'planck', '[a, b]', 'a', a, {'> 0'});
3072	validateParam(funcName, 'Planck', 'planck', '[a, b]', 'b', b, {'> 0'});
3073
3074	zetav = feval(funcName, 'zeta', a+1, sampleSize);
3075	out = feval(funcName, 'gamma', a+1, sampleSize) ./ ...
3076	(b * zetav);
3077
3078	case {'poisson', 'po'}
3079	% START po HELP START poisson HELP
3080	% THE POISSON DISTRIBUTION
3081	% ~ Poisson(lambda)
3082	%
3083	% pdf(y) = exp(-lambda)*lambda^y/factorial(y) = ...
3084	% exp( -lambda + y*log(lambda) - gammaln(y+1) ); lambda>0
3085	%
3086	% Mean = lambda;
3087	% Variance = lambda
3088	% Mode = floor(lambda);
3089	%
3090	% PARAMETERS:
3091	% lambda, lambda > 0
3092	%
3093	% SUPPORT:
3094	% y, y = 0, 1, 2, 3, 4, ...
3095	%
3096	% CLASS:
3097	% Discrete distributions
3098	%
3099	% NOTES:
3100	% 1. If lambda is an integer, Mode also equals (lambda+1).
3101	% 2. The Poisson distribution is commonly used as an approximation
3102	% to the Binomial distribution when probability of success is very small.
3103	% 3. In queueing theory, when interarrival times are ~Exponential, the number of arrivals in a
3104	% fixed interval are ~Poisson.
3105	% 4. Errors in observations with integer values (i.e., miscounting) are ~Poisson.
3106	%
3107	% USAGE:
3108	% randraw('po', lambda, sampleSize) - generate sampleSize
3109	% number of variates from the Poisson distribution with
3110	% parameter lambda;
3111	% randraw('po') - help for the Poisson distribution;
3112	%
3113	% EXAMPLES:
3114	% 1. y = randraw('po', [2], [1 1e5]);
3115	% 2. y = randraw('po', [3], 1, 1e5);
3116	% 3. y = randraw('po', [10.4], 1e5 );
3117	% 4. y = randraw('po', [100.25], [1e5 1] );
3118	% 5. randraw('po');
3119	%
3120	% SEE ALSO:
3121	% BINOMIAL distribution
3122	%
3123	% END po HELP END poisson HELP
3124
3125	% Method:
3126	% 1) If lambda > 1000 we use normal approximation to poisson distribution
3127	% mean=lambda, variance=lambda^2
3128	% 2) If lambda < 1000 we use the following reference:
3129	% George Marsaglia, Fast Generation Of Discrete Random Variables,
3130	% Journal of Statistical Software, July 2004, Volume 11, Issue 4
3131	%
3132	% (Note: this method does not support lambda<5e-10.
3133	% So for lambda<5e-10 we return 0 )
3134
3135	if numel(distribParams) ~= 1
3136	error('Poisson Distribution: Wrong numebr of parameters (run randraw(''poisson'') for help) ');
3137	end
3138	lam = distribParams(1);
3139	if lam < 0
3140	error('Poisson Distribution: Parameter ''lambda'' should be positive (run randraw(''poisson'') for help) ');
3141	end
3142
3143	if lam > 1e3
3144	% For sufficiently large values of lambda (lambda > 1000 say),
3145	% the normal distribution with mean lambda and variance lambda is
3146	% an excellent approximation to the Poisson distribution.
3147
3148	out = round( lam + sqrt(lam) * randn( sampleSize ) );
3149
3150	else % lam <= 1e3
3151
3152	if lam<21.4
3153	if lam<5e-10
3154	varargout{1} = zeros( sampleSize );
3155	return;
3156	end
3157	t=exp(-lam);
3158	p = t;
3159	ii = 1;
3160	while t*2147483648 > 1
3161	t = t * (lam/ii);
3162	ii = ii + 1;
3163	end
3164	sizeP = ii-1;
3165	offset = 0;
3166	%/* Given size, fill P array (30-bit integers) */
3167
3168	P = round( 2^30pcumprod([1, lam./(1:sizeP-1)] ) );
3169	else %lam>21.4
3170
3171	% maximum lam = 1940;
3172
3173	mode = floor(lam);
3174
3175	loglam = log(lam);
3176	log2147483648 = log(2147483648);
3177	tmode = -lam + mode*loglam - gammaln(mode+1);
3178	pmode = exp( tmode );
3179
3180	t = tmode;
3181	ii = mode + 1;
3182	while t + log2147483648 > 0
3183	t = t + loglam - log(ii);
3184	ii = ii + 1;
3185	end
3186	last = ii-2;
3187
3188	t = tmode;
3189	j=-1;
3190	for ii=mode-1:-1:0
3191	t = t - loglam + log(ii+1);
3192	if t + log2147483648 < 0
3193	j=ii;
3194	break;
3195	end
3196	end
3197
3198	offset = j+1;
3199	sizeP = last-offset+1;
3200
3201	P = zeros(1, sizeP);
3202
3203	ii = (mode+1):last;
3204	P(mode-offset+1:last-offset+1) = round( 2^30pmodecumprod([1, lam./ii]) );
3205
3206	ii = (mode-1):-1:offset;
3207	P(mode-offset:-1:1) = round( 2^30pmodecumprod((ii+1)/lam) );
3208
3209	end
3210
3211	out = randFrom5Tbls( P, offset, sampleSize);
3212
3213	end
3214
3215	case {'quadratic', 'quad', 'quadr'}
3216	% START quadratic HELP START quad HELP START quadr HELP
3217	% THE QUADRATIC DISTRIBUTION
3218	%
3219	% Standard form of the quadratic distribution:
3220	%
3221	% pdf(y) = 3/4*(1-y.^2); -1<=y<=1;
3222	%
3223	% Mean = 0;
3224	% Variance = 1/5;
3225	%
3226	% General form of the quadratic distribution:
3227	%
3228	% pdf(y) = 3/(4s) (1-((y-t)/s).^2); t-s<=y<=t+s; s>0
3229	% cdf(y) = 1/2 + 3/4(y-t)/s - 1/4((y-t)/s).^3; ; t-s<=y<=t+s; s>0
3230	%
3231	% Mean = t;
3232	% Variance = s^2/5;
3233	%
3234	% PARAMETERS:
3235	% t - location
3236	% s -scale; s>0
3237	%
3238	% SUPPORT:
3239	% y, -1 <= y <= 1 - standard Quadratic distribution
3240	% or
3241	% y, t-s <= y <= t+s - generalized Quadratic distribution
3242	%
3243	% CLASS:
3244	% Continuous distributions
3245	%
3246	% USAGE:
3247	% randraw('quadr', [], sampleSize) - generate sampleSize number
3248	% of variates from standard Quadratic distribution;
3249	% randraw('quadr', [t, s], sampleSize) - generate sampleSize number
3250	% of variates from generalized Quadratic distribution
3251	% with location parameter 't' and scale parameter 's';
3252	% randraw('quadr') - help for Quadratic distribution;
3253	%
3254	% EXAMPLES:
3255	% 1. y = randraw('quadr', [], [1 1e5]);
3256	% 2. y = randraw('quadr', [], 1, 1e5);
3257	% 3. y = randraw('quadr', [], 1e5 );
3258	% 4. y = randraw('quadr', [10 3], [1e5 1] );
3259	% 5. randraw('quadr');
3260	% END quadratic HELP END quad HELP END quadr HELP
3261
3262
3263	% Method:
3264	%
3265	% Inverse CDF transformation method.
3266	% We use Vi`ete formula to solve cubic equation.
3267
3268	checkParamsNum(funcName, 'Quadratic', 'quadratic', distribParams, [0, 2]);
3269	if numel(distribParams)==2
3270	t = distribParams(1);
3271	s = distribParams(2);
3272	validateParam(funcName, 'Quadratic', 'quadratic', '[t, s]', 's', s, {'> 0'});
3273	else
3274	t = 0;
3275	s = 1;
3276	end
3277
3278	out = t + s * 2sin(1/3asin(2*rand( sampleSize )-1));
3279
3280	case {'rademacher'}
3281	% START rademacher HELP
3282	% THE RADEMACHER DISTRIBUTION
3283	% The Rademacher distribution takes value 1 with probability 1/2
3284	% and value -1 with probability 1/2 (it is simply a random sign)
3285	% ( Hans Rademacher (1892-1969) )
3286	%
3287	% PARAMETERS:
3288	% None
3289	%
3290	% SUPPORT:
3291	% y = -1 or 1
3292	%
3293	% CLASS:
3294	% Descrete distributions
3295	%
3296	% USAGE:
3297	% randraw('rademacher', [], sampleSize) - generate sampleSize number
3298	% of variates from the Rademacher distribution;
3299	% randraw('rademacher') - help for the Rademacher distribution;
3300	%
3301	% EXAMPLES:
3302	% 1. y = randraw('rademacher', [], [1 1e5]);
3303	% 2. y = randraw('rademacher', [], 1, 1e5);
3304	% 3. y = randraw('rademacher', [], 1e5 );
3305	% 4. y = randraw('rademacher', [], [1e5 1] );
3306	% 5. randraw('rademacher');
3307	% END rademacher HELP
3308
3309	checkParamsNum(funcName, 'Rademacher', 'rademacher', distribParams, [0]);
3310
3311	out = 2*round(rand(sampleSize)) - 1;
3312
3313	case {'rayl', 'rayleigh'}
3314	% START rayl HELP START rayleigh HELP
3315	% THE RAYLEIGH DISTRIBUTION
3316	%
3317	% pdf = y./sigma^2 * exp(-y.^2/(2*sigma^2)); y >= 0
3318	% cdf = 1 - exp(-y.^2/(2*sigma^2));
3319	%
3320	% Mean = sqrt(pi/2)*sigma;
3321	% Variance = (4-pi)/2*sigma^2;
3322	%
3323	% PARAMETERS:
3324	% sigma - scale parameter (-Inf < sigma < Inf)
3325	%
3326	% SUPPORT:
3327	% y, y >= 0
3328	%
3329	% CLASS:
3330	% Continuous skewed distributions
3331	%
3332	% USAGE:
3333	% randraw('rayl', sigma, sampleSize) - generate sampleSize number
3334	% of variates from the Rayleigh distribution with scale parameter 'sigma';
3335	% randraw('rayl') - help for the Rayleigh distribution;
3336	%
3337	% EXAMPLES:
3338	% 1. y = randraw('rayl', 1, [1 1e5]);
3339	% 2. y = randraw('rayl', 2.5, 1, 1e5);
3340	% 3. y = randraw('rayl', 3, 1e5 );
3341	% 4. y = randraw('rayl', 4, [1e5 1] );
3342	% 5. randraw('rayl');
3343	%
3344	% SEE ALSO:
3345	% CHI, MAXWELL, WEIBULL distributions
3346	% END rayl HELP END rayleigh HELP
3347
3348	checkParamsNum(funcName, 'Rayleigh', 'rayl', distribParams, [1]);
3349	sigma = distribParams(1);
3350
3351	out = sqrt(-2 * sigma^2 * log(rand( sampleSize ) ));
3352
3353	case {'semicirc', 'semicircle', 'wigner'}
3354	% START semicirc HELP START semicircle HELP START wigner HELP
3355	% THE SEMICIRCLE DISTRIBUTION
3356	% ( Wigner semicircle distribution)
3357	%
3358	% The Wigner semicircle distribution, named after the physicist Eugene Wigner,
3359	% is the probability distribution supported on the interval [m?R, m+R] the graph
3360	% of whose probability density function is a semicircle of radius R centered at
3361	% (m, 0) and then suitably normalized (so that it is really a semi-ellipse).
3362	%
3363	% pdf = 2/(piR^2) sqrt(1-(y-m).^2);
3364	%
3365	% Mean = m;
3366	% Variance = R^2/4;
3367	%
3368	% PARAMETERS:
3369	% m - location;
3370	% R - semicircle radius; (R>0)
3371	%
3372	% SUPPORT:
3373	% y, m-R <= y <= m+R
3374	%
3375	% CLASS:
3376	% Continuous symmetric distributions
3377	%
3378	% USAGE:
3379	% randraw('semicirc', [m, R], sampleSize) - generate sampleSize number
3380	% of variates from the Semicircle distribution on the interval [m-R, m+R];
3381	% randraw('semicirc') - help for the Semicircle distribution;
3382	%
3383	% EXAMPLES:
3384	% 1. y = randraw('semicirc', [0, 1], [1 1e5]);
3385	% 2. y = randraw('semicirc', [-1.5, 5], 1, 1e5);
3386	% 3. y = randraw('semicirc', [2, 10], 1e5 );
3387	% 4. y = randraw('semicirc', [11, 1], [1e5 1] );
3388	% 5. randraw('semicirc');
3389	% END semicirc HELP END semicircle HELP END wigner HELP
3390
3391	checkParamsNum(funcName, 'Semicircle', 'semicirc', distribParams, [2]);
3392	m = distribParams(1);
3393	R = distribParams(2);
3394	validateParam(funcName, 'Semicircle', 'semicirc', '[m, R]', 'R', R, {'> 0'});
3395
3396	out = m + Rsqrt(rand(sampleSize)) . cos( rand(sampleSize)*pi );
3397
3398	case {'skellam'}
3399	% START skellam HELP
3400	% THE SKELLAM DISTRIBUTION
3401	%
3402	% The Skellam distribution is the probability distribution of the difference N1 - N2
3403	% of two uncorrelated random variables N1 and N2 having Poisson distributions
3404	% with different expected values lambda1 and lambda2.
3405	% The Skellam probability density function is:
3406	%
3407	% pdf(n) = exp(-lambda1+lambda2)(lambda1/lambda2)^(n/2)besseli(n,2sqrt(lambda1lambda2));
3408	% where besseli is the modified Bessel function of the first kind.
3409	%
3410	% Mean = lambda1 - lambda2;
3411	% Variance = lambda1 + lambda2;
3412	% Kurtosis = 1/(lambda1 + lambda2);
3413	% Skewness = (lambda1 - lambda2)/(lambda1 + lambda2)^(3/2);
3414	%
3415	% PARAMETERS:
3416	% lambda1 >= 0;
3417	% lambda2 >= 0;
3418	%
3419	% SUPPORT:
3420	% n = ..., -2, -1, 0, 1, 2, 3 ...
3421	%
3422	% CLASS:
3423	% Discrete distributions
3424	%
3425	% USAGE:
3426	% randraw('skellam', [lambda1, lambda2], sampleSize) - generate sampleSize number
3427	% of variates from the Skellam distribution with parameters
3428	% lambda1 and lambda2;
3429	% randraw('skellam') - help for the Skellam distribution;
3430	%
3431	% EXAMPLES:
3432	% 1. y = randraw('skellam', [1, 2], [1 1e5]);
3433	% 2. y = randraw('skellam', [3, 3], 1, 1e5);
3434	% 3. y = randraw('skellam', [5, 6], 1e5 );
3435	% 4. y = randraw('skellam', [1.5, 5.6], [1e5 1] );
3436	% 5. randraw('skellam');
3437	%
3438	% SEE ALSO:
3439	% Poisson Distribution ( randraw('po') );
3440	% END skellam HELP
3441
3442
3443	checkParamsNum(funcName, 'Skellam', 'skellam', distribParams, [2]);
3444	lambda1 = distribParams(1);
3445	lambda2 = distribParams(2);
3446	validateParam(funcName, 'Skellam', 'skellam', '[lambda1, lambda2]', 'lambda1', lambda1, {'> 0'});
3447	validateParam(funcName, 'Skellam', 'skellam', '[lambda1, lambda2]', 'lambda2', lambda2, {'> 0'});
3448
3449	poiss1 = feval(funcName,'poisson',lambda1, sampleSize);
3450	out = feval(funcName,'poisson',lambda2, sampleSize);
3451
3452	out = poiss1 - out;
3453
3454	case {'studentst', 't'}
3455	% START studentst HELP START t HELP
3456	% THE STUDENT'S-T DISTRIBUTION
3457	% ( t-distribution )
3458	%
3459	% Standard form of Student's-t distribution:
3460	%
3461	% pdf = gamma((nu+1)/2)/(sqrt(nupi)gamma(nu/2)) *(1+y.^2/nu)^(-(nu+1)/2);
3462	% or alternatively:
3463	% pdf = (1+y.^2/nu)^(-(nu+1)/2) / (sqrt(nu)*beta(1/2,nu/2));
3464	%
3465	% cdf = 1/2 + ( -(y<0) + (y>=0) ) .* ...
3466	% 1/2*betainc( y.^2./(nu+y.^2), 1/2, nu/2 );
3467	%
3468	% Mean = 0;
3469	% Variance = nu/(nu-2);
3470	% Skewness = 0;
3471	% Kurtosis = 3( (nu-2)^2gamma(nu/2-2)/(4*gamma(nu/2)) - 1 );
3472	%
3473	% General form of Student's-t distribution:
3474	%
3475	% pdf = gamma((nu+1)/2)/(sqrt(nupi)gamma(nu/2)) *(1+((y-chi)/eta).^2/nu)^(-(nu+1)/2);
3476	% or alternatively:
3477	% pdf = (1+((y-chi)/eta).^2/nu)^(-(nu+1)/2) / (sqrt(nu)*beta(1/2,nu/2));
3478	%
3479	% cdf = 1/2 + ( -(y<chi) + (y>=chi) ) .* ...
3480	% 1/2*betainc( ((y-chi)/eta).^2./(nu+((y-chi)/eta).^2), 1/2, nu/2 );
3481	%
3482	% Mean = chi;
3483	% Variance = nu/(nu-2)*eta^2; (nu>2)
3484	% Skewness = 0;
3485	% Kurtosis = 3( (nu-2)^2gamma(nu/2-2)/(4*gamma(nu/2)) - 1 );
3486	%
3487	% PARAMETERS:
3488	% nu - degrees of freedom (nu = 1, 2, 3, ...)
3489	% chi - location parameter
3490	% eta - scale parameter ( eta > 0 )
3491	%
3492	% SUPPORT:
3493	% y , -Inf < y < Inf
3494	%
3495	% CLASS:
3496	% Continuous symmetric distributions
3497	%
3498	% USAGE:
3499	% randraw('t', nu, sampleSize) - generate sampleSize number
3500	% of variates from the standard Student's t-distribution with degrees of
3501	% freedom 'nu'
3502	% randraw('t', [nu, chi, eta], sampleSize) - generate sampleSize number
3503	% of variates from the generalized Student's t-distribution with degrees of
3504	% freedom 'nu', location 'chi' and scale parameter 'eta'
3505	% randraw('t') - help for the Student's t-distribution;
3506	%
3507	% EXAMPLES:
3508	% 1. y = randraw('t', 3, [1 1e5]);
3509	% 2. y = randraw('t', [4, -10, 3], 1, 1e5);
3510	% 3. y = randraw('t', [5, 6.5, 10.5], 1e5 );
3511	% 4. y = randraw('t', [6, 7, 8], [1e5 1] );
3512	% 5. randraw('t');
3513	% END studentst HELP END t HELP
3514
3515	% Method:
3516	%
3517	% If nu<=100 we utilize the following transformation:
3518	%
3519	% Y = X/sqrt(Z/nu) ~ Student's-t(nu),
3520	% where X~Normal(0,1) and Z~Chi^2(nu);
3521	%
3522	% Else, we use Normal(0, 1) instead of Student's-t
3523	%
3524
3525	checkParamsNum(funcName, 'Student''s'' t', 't', distribParams, [1 3]);
3526	if numel(distribParams)==3
3527	nu = distribParams(1);
3528	chi = distribParams(2);
3529	eta = distribParams(3);
3530	validateParam(funcName, 'Student''s'' t', 't', '[nu, chi, eta]', 'nu', nu, {'> 0','==integer'});
3531	validateParam(funcName, 'Student''s'' t', 't', '[nu, chi, eta]', 'eta', eta, {'> 0'});
3532	else
3533	nu = distribParams(1);
3534	validateParam(funcName, 'Student''s'' t', 't', 'nu', 'nu', nu, {'> 0','==integer'});
3535	chi = 0;
3536	eta = 1;
3537	end
3538
3539
3540	if nu <= 100
3541	chisq = feval(funcName, 'chisq', nu, sampleSize);
3542	out = chi + eta * sqrt(nu)*randn( sampleSize ) ./ sqrt( chisq );
3543	else
3544	out = chi + eta * randn( sampleSize );
3545	end
3546
3547	case {'tri', 'triangular'}
3548	% START tri HELP START triangular HELP
3549	% THE TRIANGULAR DISTRIBUTION
3550	%
3551	% pdf = (a <= y & y <= c) .* ( 2(y-a)/((b-a)(c-a)) ) + ...
3552	% (c < y & y <= b) .* ( 2(b-y)/((b-a)(b-c)) ) + ...
3553	% (y < a \| y > b) .* 0;
3554	%
3555	% cdf = ( y < a ) .* 0 + ...
3556	% (a <= y & y <= c) .* ( (y-a).^2/((b-a)*(c-a)) ) + ...
3557	% (c < y & y <= b) .* ( 1- (b-y).^2/((b-a)*(b-c)) ) + ...
3558	% (y > b) .* 1;
3559	%
3560	% Mean = 1/3*(a+b+c);
3561	% Variance = 1/18(a^2+b^2+c^2-ab-ac-bc);
3562	% Skewness = sqrt(2)(a+b-2c)(2a-b-c)(a-2b+c) / (5(a^2+b^2+c^2-ab-ac-bc)^(3/2));
3563	% Kurtosis = -3/5;
3564	%
3565	% PARAMETERS:
3566	% a - lower bound
3567	% c - mode (c>a)
3568	% b - upper bound (b>c>a)
3569	%
3570	% SUPPORT:
3571	% y, a <= y <= c
3572	%
3573	% CLASS:
3574	% Continuous skewed distributions
3575	%
3576	% USAGE:
3577	% randraw('tri', [a, c, b], sampleSize) - generate sampleSize number
3578	% of variates from the Triangular distribution with parameters
3579	% 'a', 'c' and 'b';
3580	% randraw('tri') - help for the Triangular distribution;
3581	%
3582	% EXAMPLES:
3583	% 1. y = randraw('tri', [0, 1, 2], [1 1e5]);
3584	% 2. y = randraw('tri', [1, 10, 20], 1, 1e5);
3585	% 3. y = randraw('tri', [0.5, 5, 10.5], 1e5 );
3586	% 4. y = randraw('tri', [2.4, 3.4, 5.4], [1e5 1] );
3587	% 5. randraw('tri');
3588	% END tri HELP END triangular HELP
3589
3590	checkParamsNum(funcName, 'Triangular', 'tri', distribParams, [3]);
3591	a = distribParams(1);
3592	c = distribParams(2);
3593	b = distribParams(3);
3594	validateParam(funcName, 'Triangular', 'tri', '[a, b, c]', 'c-a', c-a, {'> 0'});
3595	validateParam(funcName, 'Triangular', 'tri', '[a, b, c]', 'b-c', b-c, {'> 0'});
3596	validateParam(funcName, 'Triangular', 'tri', '[a, b, c]', 'b-a', b-a, {'> 0'});
3597
3598	t = (c-a) / (b-a);
3599	u1 = rand( sampleSize );
3600	out = a + (b-a) * ...
3601	((u1 <= t) .* sqrt( tu1 ) + (u1 > t) . ( 1 - sqrt((1-t)*(1-u1)) ));
3602
3603	case {'tukeylambda'}
3604	% START tukeylambda HELP
3605	% THE TUKEY-LAMBDA DISTRIBUTION
3606	%
3607	% The Tukey-Lambda Distribution with shape parameter lambda.
3608	%
3609	% The Tukey-Lambda distribution does not have a simple closed form
3610	% for either the probability density function or the cumulative
3611	% distribution function. The cumulative distribution function is
3612	% calculated numerically. Some special cases are:
3613	% lambda = -1 - approximately Cauchy;
3614	% lambda = 0 - exactly logistic;
3615	% lambda = 0.14 - approximately normal;
3616	% lambda = 0.5 - U-shaped;
3617	% lambda= 1 - exactly uniform.
3618	%
3619	% PARAMETERS:
3620	% lambda - shape parameter
3621	%
3622	% SUPPORT:
3623	% y, -1/lambdal <= y <= 1/lambda
3624	%
3625	% CLASS:
3626	% Continuous symmetric distributions
3627	%
3628	% USAGE:
3629	% randraw('tukeylambda', lambda, sampleSize) - generate sampleSize number
3630	% of variates from the Tukey-Lambda distribution with shale parameter
3631	% 'lambda';
3632	% randraw('tukeylambda') - help for the Tukey-Lambda distribution;
3633	%
3634	% EXAMPLES:
3635	% 1. y = randraw('tukeylambda', -1, [1 1e5]);
3636	% 2. y = randraw('tukeylambda', 0, 1, 1e5);
3637	% 3. y = randraw('tukeylambda', 0.14, 1e5 );
3638	% 4. y = randraw('tukeylambda', 0.5, [1e5 1] );
3639	% 5. randraw('tukeylambda');
3640	% END tukeylambda HELP
3641
3642	checkParamsNum(funcName, 'Tukey-Lambda', 'tukeylambda', distribParams, [1]);
3643	lambda = distribParams(1);
3644	if lambda ~= 0
3645	u = rand( sampleSize );
3646	out = (u.^lambda - (1-u).^lambda) / lambda;
3647	else
3648	out = feval(funcName,'logistic', [0 1], sampleSize);
3649	end
3650
3651	case {'u', 'ushape'}
3652	% START u HELP START ushape HELP
3653	% THE U DISTRIBUTION
3654	% ( U-shape distribution )
3655	%
3656	% Standard form of the U distribution:
3657	%
3658	% pdf(y) = 1./(pi*sqrt(1-y.^2)); -1<=y<=1;
3659	% cdf(y) = 1/2 + 1/pi*asin(y); -1<=y<=1;
3660	%
3661	% Mean = 0;
3662	% Variance = 1/2;
3663	%
3664	% General form of the U distribution:
3665	%
3666	% pdf(y) = 1./(pi*sqrt(s^2-(y-t).^2)); t-s<=y<=t+s; s>0
3667	% cdf(y) = 1/2 + 1/pi*asin((y-t)/a); -1<=y<=1;
3668	%
3669	% Mean = t;
3670	% Variance = s^2/2;
3671	%
3672	% PARAMETERS:
3673	% t - location
3674	% s -scale; s>0
3675	%
3676	% SUPPORT:
3677	% y, -1 <= y <= 1 - standard U distribution
3678	% or
3679	% y, t-s <= y <= t+s - generalized U distribution
3680	%
3681	% CLASS:
3682	% Continuous symmetric distributions
3683	%
3684	% USAGE:
3685	% randraw('u', [], sampleSize) - generate sampleSize number
3686	% of variates from standard U distribution;
3687	% randraw('u', [t, s], sampleSize) - generate sampleSize number
3688	% of variates from generalized U distribution
3689	% with location parameter 't' and scale parameter 's';
3690	% randraw('u') - help for U distribution;
3691	%
3692	% EXAMPLES:
3693	% 1. y = randraw('u', [], [1 1e5]);
3694	% 2. y = randraw('u', [], 1, 1e5);
3695	% 3. y = randraw('u', [], 1e5 );
3696	% 4. y = randraw('u', [10 3], [1e5 1] );
3697	% 5. randraw('u');
3698	%
3699	% SEE ALSO:
3700	% ARCSIN distribution
3701	% END u HELP END ushape HELP
3702
3703	checkParamsNum(funcName, 'U', 'u', distribParams, [0, 2]);
3704	if numel(distribParams)==2
3705	t = distribParams(1);
3706	s = distribParams(2);
3707	validateParam(funcName, 'U', 'u', '[t, s]', 's', s, {'> 0'});
3708	else
3709	t = 0;
3710	s = 1;
3711	end
3712
3713	out = t + s * sin(pi*(rand(sampleSize)-0.5));
3714
3715	case {'uniform', 'unif'}
3716	% START uniform HELP START unif HELP
3717	% THE UNIFORM DISTRIBUTION
3718	%
3719	% pdf = 1/(b-a);
3720	% cdf = (y-a)/(b-a);
3721	%
3722	% Mean = (a+b)/2;
3723	% Variance = (b-a)^2/12;
3724	%
3725	% PARAMETERS:
3726	% a is location of y (lower bound);
3727	% b is scale of y (upper bound) (b > a);
3728	%
3729	% SUPPORT:
3730	% y, a < y < b
3731	%
3732	% CLASS:
3733	% Continuous symmetric distributions
3734	%
3735	% USAGE:
3736	% randraw('uniform', [], sampleSize) - generate sampleSize number
3737	% of variates from standard Uniform distribution (a=0, b=1);
3738	% randraw('uniform', [a, b], sampleSize) - generate sampleSize number
3739	% of variates from the Uniform distribution
3740	% with parameters 'a' and 'b';
3741	% randraw('uniform') - help for the Uniform distribution;
3742	%
3743	% EXAMPLES:
3744	% 1. y = randraw('uniform', [], [1 1e5]);
3745	% 2. y = randraw('uniform', [2, 3], 1, 1e5);
3746	% 3. y = randraw('uniform', [5, 6], 1e5 );
3747	% 4. y = randraw('uniform', [10.5, 11.5], [1e5 1] );
3748	% 5. randraw('uniform');
3749	% END uniform HELP END unif HELP
3750
3751	checkParamsNum(funcName, 'Uniform', 'uniform', distribParams, [0, 2]);
3752	if numel(distribParams)==2
3753	a = distribParams(1);
3754	b = distribParams(2);
3755	validateParam(funcName, 'Uniform', 'uniform', '[a, b]', 'b-a', b-a, {'> 0'});
3756	else
3757	a = 0;
3758	b = 1;
3759	end
3760
3761	out = a + (b-a)*rand( sampleSize );
3762
3763	case {'vonmises'}
3764	% START vonmises HELP
3765	% THE VON MISES DISTRIBUTION
3766	% A continuous distribution defined on the range [0, 2*pi)
3767	% with probability density function:
3768	%
3769	% pdf(y) = exp(kcos(y-a)) ./ (2pi*besseli(0,k));
3770	%
3771	% where besseli(0,x) is a modified Bessel function of the
3772	% first kind of order 0.
3773	% Here, 'a' (a>=0 and a<2*pi) is the mean direction and k > 0 is a
3774	% concentration parameter
3775	% The von Mises distribution is the circular analog of the normal
3776	% distribution on a line.
3777	%
3778	% Mean = a;
3779	%
3780	% PARAMETERS:
3781	% a - location parameter, (a>=0 and a<2*pi)
3782	% k - shape parameter, (k>0)
3783	%
3784	% SUPPORT:
3785	% y, -pi < y < pi
3786	%
3787	% CLASS:
3788	% Continuous distribution
3789	%
3790	% USAGE:
3791	% randraw('vonmises', [a, k], sampleSize) - generate sampleSize number
3792	% of variates from the von Mises distribution with parameters 'a' and 'k';
3793	% randraw('vonmises') - help for the von Mises distribution;
3794	%
3795	% EXAMPLES:
3796	% 1. y = randraw('vonmises', [pi/2, 3], [1 1e5]);
3797	% 2. y = randraw('vonmises', [2*pi/3, 2], 1, 1e5);
3798	% 3. y = randraw('vonmises', [pi/4, 10], 1e5 );
3799	% 4. y = randraw('vonmises', [pi, 2.2], [1e5 1] );
3800	% 5. randraw('vonmises');
3801	% END vonmises HELP
3802
3803	% Method:
3804	% 1) For large k (say, k>700) von Mises distribution tends
3805	% to a Normal Distribution with variance 1/?
3806	% 2) For a small k we implement method suggested in:
3807	% L. Yuan and J.D. Kalbleisch, "On the Bessel distribution and
3808	% related problems," Annals of the Institute of Statistical
3809	% Mathematics, vol. 52, pp. 438-447, 2000
3810	% and described in:
3811	% Luc Devroye, "Simulating Bessel Random Variables,"
3812	% Statistics and Probability Letters, vol. 57, pp. 249-257, 2002.
3813	%
3814
3815
3816	a = distribParams(1);
3817	k = distribParams(2);
3818
3819	if k > 700
3820	% for large k it tends to a Normal Distribution with variance 1/k
3821	out = a + sqrt(1/k)*randn( sampleSize );
3822	else
3823	% Generate X <- Bessel(0,k)
3824	x = feval(funcName,'bessel', [0 k], sampleSize);
3825
3826	% Generate B <- beta(X+1/2, 1/2);
3827	u2 = rand( sampleSize );
3828	l = (x>0);
3829	d1 = (cos(2piu2(l))).^2;
3830	d2 = (cos(pi*u2(~l))).^2;
3831	clear('u2');
3832	t = 2./(2*(x(l)+0.5)-1);
3833	clear('x');
3834	b = zeros( sampleSize );
3835	u1 = rand(sum(l(:)),1);
3836	b(l) = 1 - (1-u1.^t(:)) .* d1(:);
3837	clear('t');
3838	b(~l) = d2;
3839	clear('d1');
3840	clear('d2');
3841
3842	% if U < 1/(1+exp(-2ksqrt(B))),
3843	% then return S*acos(sqrt(B))
3844	% else return S*acos(-sqrt(B))
3845	% where S is a random sign
3846	l = rand(sampleSize) < 1./(1+exp(-2ksqrt(b)));
3847	out = zeros( sampleSize );
3848	out(l) = acos(sqrt(b(l)));
3849	out(~l) = acos(-sqrt(b(~l)));
3850	clear('b');
3851	clear('l');
3852
3853	out = a + (2round(rand(sampleSize))-1) . out;
3854	end
3855
3856	case {'wald'}
3857	% START wald HELP
3858	% THE WALD DISTRIBUTION
3859	%
3860	% The Wald distribution is as special case of the Inverse Gaussian Distribution
3861	% where the mean is a constant with the value one.
3862	%
3863	% pdf = sqrt(chi/(2piy^3)) * exp(-chi./(2y).(y-1).^2);
3864	%
3865	% Mean = 1;
3866	% Variance = 1/chi;
3867	% Skewness = sqrt(9/chi);
3868	% Kurtosis = 3+ 15/scale;
3869	%
3870	% PARAMETERS:
3871	% chi - scale parameter; (chi>0)
3872	%
3873	% SUPPORT:
3874	% y, y>0
3875	%
3876	% CLASS:
3877	% Continuous skewed distributions
3878	%
3879	% USAGE:
3880	% randraw('wald', chi, sampleSize) - generate sampleSize number
3881	% of variates from the Wald distribution with scale parameter 'chi';
3882	% randraw('wald') - help for the Wald distribution;
3883	%
3884	% EXAMPLES:
3885	% 1. y = randraw('wald', 0.5, [1 1e5]);
3886	% 2. y = randraw('wald', 1, 1, 1e5);
3887	% 3. y = randraw('wald', 1.5, 1e5 );
3888	% 4. y = randraw('wald', 2, [1e5 1] );
3889	% 5. randraw('wald');
3890	% END wald HELP
3891
3892	checkParamsNum(funcName, 'Wald', 'wald', distribParams, [1]);
3893	chi = distribParams(1);
3894	validateParam(funcName, 'Wald', 'wald', 'chi', 'chi', chi, {'> 0'});
3895
3896	out = feval(funcName, 'ig', [1 chi], sampleSize);
3897
3898	case {'weibull', 'frechet', 'wbl'}
3899	% START weibull HELP START frechet HELP START wbl HELP
3900	% THE WEIBULL DISTRIBUTION
3901	% ( sometimes: Frechet distribution )
3902	%
3903	% cdf = 1 - exp(-((y-theta)/beta).^alpha)
3904	% pdf = alpha / beta * ((y-theta)/beta).^(alpha-1) .* exp(-((y-theta)/beta).^alpha);
3905	%
3906	% Mean = theta + beta*gamma((alpha+1)/alpha);
3907	% Variance = beta^2 * ( gamma((alpha+2)/alpha) - gamma((alpha+1)/alpha)^2 );
3908	% where gamma is the gamma function
3909	%
3910	% PARAMETERS:
3911	% theta - location parameter;
3912	% alpha - shape parameter ( alpha>0 );
3913	% beta - scale parameter ( beta>0 );
3914	%
3915	% SUPPORT:
3916	% y, y > theta
3917	%
3918	% CLASS:
3919	% Continuous skewed distributions
3920	%
3921	% NOTES:
3922	% If alpha=1 , it is the exponential distribution
3923	% If beta=1; alpha=2, theta=0; it is the standard Rayleigh distribution (sigma=1)
3924	%
3925	% USAGE:
3926	% randraw('weibull', [theta, alpha, beta], sampleSize) - generate sampleSize number
3927	% of variates from the Weibull distribution with location parameter 'theta',
3928	% shape parameter 'alpha' and scale parameter 'beta';
3929	% randraw('weibull') - help for the Weibull distribution;
3930	%
3931	% EXAMPLES:
3932	% 1. y = randraw('weibull', [-10, 2, 3], [1 1e5]);
3933	% 2. y = randraw('weibull', [0, 2, 1], 1, 1e5);
3934	% 3. y = randraw('weibull', [0, 1, 2], 1e5 );
3935	% 4. y = randraw('weibull', [2.1, 5.4, 10.2], [1e5 1] );
3936	% 5. randraw('weibull');
3937	% END weibull HELP END frechet HELP END wbl HELP
3938
3939	checkParamsNum(funcName, 'Weibull', 'weibull', distribParams, [3]);
3940	theta = distribParams(1);
3941	alpha = distribParams(2);
3942	beta = distribParams(3);
3943	validateParam(funcName, 'Weibull', 'weibull', '[theta, alpha, beta]', 'alpha', alpha, {'> 0'});
3944	validateParam(funcName, 'Weibull', 'weibull', '[theta, alpha, beta]', 'beta', beta, {'> 0'});
3945
3946	out = theta + beta * (-log(rand( sampleSize ))).^(1/alpha);
3947
3948	case {'yule', 'yulesimon'}
3949	% START yule HELP START yulesimon HELP
3950	% THE YULE-SIMON DISTRIBUTION
3951	%
3952	% pmf(y) = (p-1)*beta(y, p); p>1; y = 1, 2, 3, 4, ...
3953	% cdf(y) = 1-(y+1).*beta(y+1,p);
3954	%
3955	% Mean = (p-1)/(p-2); for p>2;
3956	% Variance = (p-1)^2/((p-2)^2*(p-3)); for p>3;
3957	%
3958	% PARAMETERS:
3959	% p, p>1
3960	%
3961	% SUPPORT:
3962	% y, y = 1, 2, 3, 4, ...
3963	%
3964	% CLASS:
3965	% Discrete distributions
3966	%
3967	% NOTES:
3968	% 1. It is named after George Udny Yule and Herbert Simon.
3969	% Simon originally called it the Yule distribution.
3970	% 2. The probability mass function pmf(y) has the property that
3971	% for sufficiently large y we have
3972	% pmf(y) = (p-1)*gamma(p)./y.^p;
3973	% This means that the tail of the Yule-Simon distribution is a
3974	% realization of Zipf's law: pmf(y) can be used to model,
3975	% for example, the relative frequency of the y'th most frequent
3976	% word in a large collection of text, which according to Zipf's law
3977	% is inversely proportional to a (typically small) power of y.
3978	%
3979	% USAGE:
3980	% randraw('yule', p, sampleSize) - generate sampleSize number
3981	% of variates from the Yule-Simon distribution with parameter 'p';
3982	% randraw('yule') - help for the Yule-Simon distribution;
3983	%
3984	% EXAMPLES:
3985	% 1. y = randraw('yule', 2, [1 1e5]);
3986	% 2. y = randraw('yule', 3.2, 1, 1e5);
3987	% 3. y = randraw('yule', 100.5, 1e5 );
3988	% 4. y = randraw('yule', 33, [1e5 1] );
3989	% 5. randraw('yule');
3990	% END yule HELP END yulesimon HELP
3991
3992	% Rference:
3993	% Luc Devroye,
3994	% "Non-Uniform Random Variate Generation,"
3995	% Springer Verlag, 1986, pages 550-551.
3996
3997	checkParamsNum(funcName, 'Yule-Simon', 'yule', distribParams, [1]);
3998	p = distribParams(1);
3999	validateParam(funcName, 'Yule-Simon', 'yule', 'p', 'p', p, {'> 1'});
4000
4001	out = ceil( log(rand(sampleSize)) ./ log(1-exp(log(rand(sampleSize))/(p-1))) );
4002
4003	case {'zeta', 'zipf'}
4004	% START zeta HELP START zipf HELP
4005	% ZETA DISTRIBUTION
4006	% (sometimes: Zipf distribution)
4007	%
4008	% pmf(n) = 1. / (n^a * zeta(a)),
4009	% where zeta(s) is a Riemann Zeta function:
4010	% sum from i=1 to Inf of 1/i^a
4011	% a>1
4012	%
4013	% Mean = zeta(a-1)/zeta(a), for a>2
4014	% Variance = zeta(a-2)/zeta(a) - (zeta(a-1)/zeta(a))^2;
4015	%
4016	% PARAMETERS:
4017	% a > 1
4018	%
4019	% SUPPORT:
4020	% n = 1, 2, 3, ... (positive integers)
4021	%
4022	% CLASS:
4023	% Discrete distributions
4024	%
4025	% NOTES:
4026	% The zeta distribution is a long-tailed distribution that is useful for
4027	% size-frequency data. It is sometimes used in insurance as a model for
4028	% the number of policies held by a single person in an insurance portfolio.
4029	% It is also used for the analysis of the frequency of words in long
4030	% sequences of text. When used in linguistics the zeta distribution is
4031	% known as the Zipf distribution.
4032	%
4033	% USAGE:
4034	% randraw('zeta', a, sampleSize) - generate sampleSize number
4035	% of variates from standard Zeta distribution with parameter a;
4036	% randraw('zeta') - help for Zeta distribution;
4037	%
4038	% EXAMPLES:
4039	% 1. y = randraw('zeta', 2, [1 1e5]);
4040	% 2. y = randraw('zeta', 3.5, 1, 1e5);
4041	% 3. y = randraw('zeta', 10, 1e5 );
4042	% 4. y = randraw('zeta', 100, [1e5 1] );
4043	% 5. randraw('zeta');
4044	% END zeta HELP END zipf HELP
4045
4046	% Reference:
4047	% Luc Devroye,
4048	% "Non-Uniform Random Variate Generation,"
4049	% Springer Verlag, 1986, pages 550-551.
4050
4051	a = distribParams(1);
4052
4053	b = 2^(a-1);
4054	am1 = a - 1;
4055	bm1 = b - 1;
4056
4057	u1 = rand( sampleSize );
4058	out = floor( u1.^(-1/am1) );
4059	clear('u1');
4060	u2 = rand( sampleSize );
4061	t = ( 1 + 1./out ).^(a-1);
4062
4063	indxs = find( u2.out.(t-1)/bm1 > t/b );
4064
4065	while ~isempty(indxs)
4066	indxsSize = size( indxs );
4067	u1 = rand( indxsSize );
4068	outNew = floor( u1.^(-1/am1) );
4069	clear('u1');
4070	u2 = rand( indxsSize );
4071	t = ( 1 + 1./outNew ).^(a-1);
4072
4073	l = u2.outNew.(t-1)/bm1 <= t/b;
4074
4075	out( indxs(l) ) = outNew(l);
4076	indxs = indxs(~l);
4077	end
4078
4079	otherwise
4080	fprintf('\n RANDRAW: Unknown distribution name: %s \n', distribName);
4081
4082	end % switch lower( distribNameInner )
4083
4084	end % if prod(sampleSize)>0
4085
4086	varargout{1} = out;
4087
4088	return;
4089
4090	end % if strcmp(runMode, 'genRun')
4091
4092	return;
4093
4094
4095	function checkParamsNum(funcName, distribName, runDistribName, distribParams, correctNum)
4096	if ~any( numel(distribParams) == correctNum )
4097	error('%s Variates Generation:\n %s%s%s%s%s', ...
4098	distribName, ...
4099	'Wrong numebr of parameters (run ',...
4100	funcName, ...
4101	'(''', ...
4102	runDistribName, ...
4103	''') for help) ');
4104	end
4105	return;
4106
4107
4108	function validateParam(funcName, distribName, runDistribName, distribParamsName, paramName, param, conditionStr)
4109	condLogical = 1;
4110	eqCondStr = [];
4111	for nn = 1:length(conditionStr)
4112	if nn==1
4113	eqCondStr = [eqCondStr conditionStr{nn}];
4114	else
4115	eqCondStr = [eqCondStr ' and ' conditionStr{nn}];
4116	end
4117	eqCond = conditionStr{nn}(1:2);
4118	eqCond = eqCond(~isspace(eqCond));
4119	switch eqCond
4120	case{'<'}
4121	condLogical = condLogical & (param<str2num(conditionStr{nn}(3:end)));
4122	case{'<='}
4123	condLogical = condLogical & (param<=str2num(conditionStr{nn}(3:end)));
4124	case{'>'}
4125	condLogical = condLogical & (param>str2num(conditionStr{nn}(3:end)));
4126	case{'>='}
4127	condLogical = condLogical & (param>=str2num(conditionStr{nn}(3:end)));
4128	case{'~='}
4129	condLogical = condLogical & (param~=str2num(conditionStr{nn}(3:end)));
4130	case{'=='}
4131	if strcmp(conditionStr{nn}(3:end),'integer')
4132	condLogical = condLogical & (param==floor(param));
4133	else
4134	condLogical = condLogical & (param==str2num(conditionStr{nn}(3:end)));
4135	end
4136	end
4137	end
4138
4139	if ~condLogical
4140	error('%s Variates Generation: %s(''%s'',%s, SampleSize);\n Parameter %s should be %s\n (run %s(''%s'') for help)', ...
4141	distribName, ...
4142	funcName, ...
4143	runDistribName, ...
4144	distribParamsName, ...
4145	paramName, ...
4146	eqCondStr, ...
4147	funcName, ...
4148	runDistribName);
4149	end
4150	return;
4151
4152	function cdf = normcdf(y)
4153	cdf = 0.5*(1+erf(y/sqrt(2)));
4154	return;
4155
4156	function pdf = normpdf(y)
4157	pdf = 1/sqrt(2pi) exp(-1/2*y.^2);
4158	return;
4159
4160	function cdfinv = norminv(y)
4161	cdfinv = sqrt(2) * erfinv(2*y - 1);
4162	return;
4163
4164	function out = randFrom5Tbls( P, offset, sampleSize)
4165	sizeP = length(P);
4166
4167	if sizeP == 0
4168	out = [];
4169	return;
4170	end
4171
4172	a = mod(floor([0 P]/16777216), 64);
4173	na = cumsum( a );
4174	b = mod(floor([0 P]/262144), 64);
4175	nb = cumsum( b );
4176	c = mod(floor([0 P]/4096), 64);
4177	nc = cumsum( c );
4178	d = mod(floor([0 P]/64), 64);
4179	nd = cumsum( d );
4180	e = mod([0 P], 64);
4181	ne = cumsum( e );
4182
4183	AA = zeros(1, na(end));
4184	BB = zeros(1, nb(end));
4185	CC = zeros(1, nc(end));
4186	DD = zeros(1, nd(end));
4187	EE = zeros(1, ne(end));
4188
4189	t1 = na(end)*16777216;
4190	t2 = t1 + nb(end)*262144;
4191	t3 = t2 + nc(end)*4096;
4192	t4 = t3 + nd(end)*64;
4193
4194	k = (1:sizeP)+offset-1;
4195	for ii = 1:sizeP
4196	AA(na(ii)+(0:a(ii+1))+1) = k(ii);
4197	BB(nb(ii)+(0:b(ii+1))+1) = k(ii);
4198	CC(nc(ii)+(0:c(ii+1))+1) = k(ii);
4199	DD(nd(ii)+(0:d(ii+1))+1) = k(ii);
4200	EE(ne(ii)+(0:e(ii+1))+1) = k(ii);
4201	end
4202
4203	%jj = round(1073741823*rand(sampleSize));
4204	jj = round(min(sum(P),1073741823) *rand(sampleSize));
4205	out = zeros(sampleSize);
4206	N = prod(sampleSize);
4207	for ii = 1:N
4208	if jj(ii) < t1
4209	out(ii) = AA( floor(jj(ii)/16777216)+1 );
4210	elseif jj(ii) < t2
4211	out(ii) = BB(floor((jj(ii)-t1)/262144)+1);
4212	elseif jj(ii) < t3
4213	out(ii) = CC(floor((jj(ii)-t2)/4096)+1);
4214	elseif jj(ii) < t4
4215	out(ii) = DD(floor((jj(ii)-t3)/64)+1);
4216	else
4217	out(ii) = EE(floor(jj(ii)-t4) + 1);
4218	end
4219	end
4220
4221	return;

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: MML/trunk/applications/common/randraw.m @ 4

Download in other formats: