[4] | 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) = b*normpdf(a-b./y) ./ (y.^2*normcdf(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(2*pi) * 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 ./ (pi*sqrt(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 = (a*b)/((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 = (n*a + m*b)/(a+b); |
---|
| 396 | % Variance = (a*b)*(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)/(2*a*b-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 < 5*z) & (r+sum_ab*log(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 >= 5*z) | (r+sum_ab*log(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_ab*t/(max_ab*t + 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_ab*t / (max_ab*t + 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).^(2*y+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.5*a*besseli(nu+1,a)/besseli(nu,a); |
---|
| 711 | % chi2 = mu + 0.25*a^2*besseli(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^30*p0*cumprod([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( (2*m+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^30*pm*cumprod([1, aa./(((m+1):last).*(((m+1):last)+nu))] ) ); |
---|
| 784 | P(m-offset:-1:1) = ... |
---|
| 785 | round( 2^30*pm*cumprod((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 | % y*log(p) + (n-y)*log(1-p) ); 0<p<1, n>1 |
---|
| 801 | % |
---|
| 802 | % Mean = n*p; |
---|
| 803 | % Variance = n*p*(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( n*p + sqrt(n*p*(1-p))*randn( sampleSize ) ); |
---|
| 864 | |
---|
| 865 | elseif p<1e-4 & n*p > 1 & n*p < 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 | mode*log(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^30*pmode*cumprod([1, h*(n-(mode+1:last)+1)./(mode+1:last)] ) ); |
---|
| 902 | P(mode-offset:-1:1) = ... |
---|
| 903 | round( 2^30*pmode*cumprod( (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+b*y) ); 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) + 2*log(1+b)) / (2*b*(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) = c*d * 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)*b*gamma((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 = 2*randraw('gamma', 0.5*r, 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*(lambda*y).^(r/4)) .* ... |
---|
| 1193 | % besseli(r/2-1, sqrt(lambda*y)); lambda>=0; r=positive integer; |
---|
| 1194 | % |
---|
| 1195 | % Mean = lambda+r; |
---|
| 1196 | % Variance = 2*(2*lambda+r); |
---|
| 1197 | % Skewness = 2*sqrt(2)*(3*lambda+r)/(2*lambda+r)^(3/2); |
---|
| 1198 | % Kurtosis = 12*(4*lambda+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))/(2*s); 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 - (coeff1*sin(coeff3*out) + coeff2*out + u) ./ ... |
---|
| 1302 | (coeff2*cos(coeff3*out)+coeff2); |
---|
| 1303 | |
---|
| 1304 | indxs = find(abs(outNext - out)>tol); |
---|
| 1305 | outPrev = out(indxs); |
---|
| 1306 | while ~isempty(indxs) |
---|
| 1307 | |
---|
| 1308 | outNext = outPrev - (coeff1*sin(coeff3*outPrev) + coeff2*outPrev + u(indxs)) ./ ... |
---|
| 1309 | (coeff2*cos(coeff3*outPrev)+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 = 12*sqrt(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) * (d1*y./(d1*y+d2)).^(d1/2) .* (1 - d1*y./(d1*y+d2)).^(d2/2) ./ y; |
---|
| 1490 | % cdf(y) = beatinc(d1*y./(d1*y+d2), d1/2, d2/2); |
---|
| 1491 | % |
---|
| 1492 | % Mean = d2/(d2-2), provided d2 > 2; |
---|
| 1493 | % Variance = 2*d2^2*(d1+d2-2)/(d1*(d2-2)^2*(d2-4)), provided d2>4; |
---|
| 1494 | % Skewness = (2*d1+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.5*d1 0.5*d2], sampleSize); |
---|
| 1527 | out = d2*out ./ (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)^2*d2^2/((d2-2)^2*d1^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.0331*x2.^2) & (log(u)>=0.5*x2+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.0331*x2.*x2) | (log(u)<0.5*x2+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)/(2*besselk(lam, sqrt(chi*psi))) .* exp(-1/2*(chi./y + psi*y)); y > 0 |
---|
| 1786 | % |
---|
| 1787 | % Mean = sqrt( chi / psi ) * besselk(lam+1,sqrt(chi*psi),1)/besselk(lam,sqrt(chi*psi),1); |
---|
| 1788 | % Variance = chi/psi * besselk(lam+2,sqrt(chi*psi),1)/besselk(lam,sqrt(chi*psi),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(psi*chi), sqrt(psi*chi)); |
---|
| 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.5*h*log(xm*ym) + 0.5*log(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)/2*log(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)/2*log(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)/2*log(m) + b/4*(m + (1/m)); |
---|
| 1895 | |
---|
| 1896 | r = (6*m + 2*h*m - b*m^2 + b)/(4*m^2); |
---|
| 1897 | s = (1 + h - b*m)/(2*m^2); |
---|
| 1898 | p = (3*s - r^2)/3; |
---|
| 1899 | q = (2*r^3)/27 - (r*s)/27 + b/(-4*m^2); |
---|
| 1900 | eta = sqrt(-(p^3)/27); |
---|
| 1901 | |
---|
| 1902 | y1 = 2*exp(log(eta)/3) * cos(acos(-q/(2*eta))/3) - r/3; |
---|
| 1903 | y2 = 2*exp(log(eta)/3) * cos(acos(-q/(2*eta))/3 + 2/3*pi) - 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)/2*log(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)/2*log(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(2*pi) * alpha^(lam-1/2) * delta^lam * ... |
---|
| 1966 | % besselk(lam, delta*sqrt(alpha^2-beta^2) ) ) * ... |
---|
| 1967 | % (delta^2 + (y-mu).^2).^(1/2*(lam-1/2)) .* ... |
---|
| 1968 | % besselk( lam-1/2, alpha*sqrt(delta^2 + (y-mu).^2) ) .* ... |
---|
| 1969 | % exp( beta*(y-mu) ); |
---|
| 1970 | % |
---|
| 1971 | % Mean = mu + beta*delta^2/(delta*sqrt(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^2*delta^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 + beta*ygig + 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/(4*a) * cos( pi*(y-t)/(2*s) ); t-s <= y <= t+s |
---|
| 2092 | % cdf(y) = 1/2 * ( 1 + sin( pi*(y-t)/(2*s) ) ); 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 + s*2/pi*asin(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)*n*M/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^30*pmode*cumprod([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^30*pmode*cumprod((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/(2*pi*y^3)) * exp(-chi./(2*y).*(y/theta-1).^2); |
---|
| 2340 | % cdf(y) = normcdf(sqrt(chi./y).*(y/theta-1)) + ... |
---|
| 2341 | % exp(2*chi/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/(2*chi)*(sqrt(9*theta^2+4*chi^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/(Y*theta^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.5*theta/chi * ( theta*V - sqrt(4*theta*chi*V + ... |
---|
| 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.5*theta/chi * ( theta*chisq1 - ... |
---|
| 2427 | sqrt(4*theta*chi*chisq1 + theta^2*chisq1.^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/(2*lam)*exp(-abs(y-theta)/lam); |
---|
| 2438 | % cdf = (y<=theta) .* 1/2*exp((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 + lam*log(2*u(u<=0.5)); |
---|
| 2485 | out(u>0.5) = theta - lam*log(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/(k*sqrt(2*pi)) * exp(-1/2*((log(y)-a)/k)^2) |
---|
| 2553 | % cdf = 1/2*(1 + erf((log(y)-a)/(k*sqrt(2)))); |
---|
| 2554 | % |
---|
| 2555 | % Mean = exp( a + k^2/2 ); |
---|
| 2556 | % Variance = exp(2*a+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) + 2*exp(3) + 3*exp(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 = 2*a*sqrt(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) + r*log(p) + y*log(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-6*p+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) + r*log(p) + y*log(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 | r*log(p) + mode*log(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^30*pmode*cumprod([1, (r+ii-1)./ii * q] ) ); |
---|
| 2742 | |
---|
| 2743 | ii = (mode-1):-1:offset; |
---|
| 2744 | P(mode-offset:-1:1) = round( 2^30*pmode*cumprod((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(2*pi) * 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/(sigma*sqrt(2*pi)) * exp(-1/2*((y-mu)/sigma).^2); |
---|
| 2761 | % cdf(y) = 1/2*(1+erf((y-mu)/(sigma*sqrt(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(delta*sqrt(alpha^2-beta^2) - beta*mu) * ... |
---|
| 2883 | % 1/sqrt(1+((y-mu)/delta).^2) .* ... |
---|
| 2884 | % besselk(1, alpha*delta*sqrt(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 = 3*beta/(alpha*sqrt(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 + beta*invGaussY + 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^2*a/((a-2)*(a-1)^2); |
---|
| 2952 | % Skewness = 2*(a+1)*sqrt(a-2)/((a-3)*sqrt(a)); |
---|
| 2953 | % Kurtosis = 6*(a^3+a^2-6*a-2)/(a*(a^2-7*a+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^30*p*cumprod([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^30*pmode*cumprod([1, lam./ii]) ); |
---|
| 3205 | |
---|
| 3206 | ii = (mode-1):-1:offset; |
---|
| 3207 | P(mode-offset:-1:1) = round( 2^30*pmode*cumprod((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/(4*s) * (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 * 2*sin(1/3*asin(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/(pi*R^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 + R*sqrt(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,2*sqrt(lambda1*lambda2)); |
---|
| 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(nu*pi)*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)^2*gamma(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(nu*pi)*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)^2*gamma(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-a*b-a*c-b*c); |
---|
| 3562 | % Skewness = sqrt(2)*(a+b-2*c)*(2*a-b-c)*(a-2*b+c) / (5*(a^2+b^2+c^2-a*b-a*c-b*c)^(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( t*u1 ) + (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(k*cos(y-a)) ./ (2*pi*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(2*pi*u2(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(-2*k*sqrt(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(-2*k*sqrt(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 + (2*round(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/(2*pi*y^3)) * exp(-chi./(2*y).*(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(2*pi) * 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; |
---|