share/distrib/rtest_distrib.mac

   1 (kill (all),
   2  load(distrib),
   3  rnd_state: make_random_state(false),
   4  set_random_state (make_random_state (123456789)),
   5  save_tolerance : float_approx_equal_tolerance,
   6  float_approx_equal_tolerance : 1e-12,
   7  0);
   8 0;
   9
  10
  11
  12
  13 /* normal distribution */
  14
  15 pdf_normal(4/3,n,1);
  16 %e^-((4/3-n)^2/2)/(sqrt(2)*sqrt(%pi));
  17
  18 (assume(s>0), cdf_normal(x,m,s));
  19 erf((x-m)/(sqrt(2)*s))/2+1/2;
  20
  21 quantile_normal(1/7,2,1);
  22 sqrt(2)*inverse_erf(-5/7)+2;
  23
  24 random_normal(2,1/3,3);
  25 [1.967829590026135, 2.120195623477506, 2.230155261001342];
  26
  27
  28
  29
  30 /* student distribution */
  31
  32 pdf_student_t(2,5/9);
  33 3*5^(5/18)*gamma(7/9)/(41^(7/9)*sqrt(%pi)*gamma(5/18));
  34
  35 cdf_student_t(2,5/9);
  36 1-beta_incomplete_regularized(5/18,1/2,5/41)/2;
  37
  38 cdf_student_t(1/10,5);
  39 1-beta_incomplete_regularized(5/2,1/2,500/501)/2;
  40
  41 quantile_student_t(0.01,26);
  42 -2.478629732546409;
  43
  44 random_student_t(102,3);
  45 [0.891862660419367,0.217299210064025,.5425534835346298];
  46
  47
  48
  49
  50 /* chi square distribution is based on gamma */
  51
  52 random_chi2(0.01,3);
  53 [5.143726638818105E-81,3.3017730563278E-188,1.338137218695894E-39];
  54
  55 random_chi2(50,3);
  56 [47.84885101221457,49.4165250891075,52.77360025014328];
  57
  58 random_chi2(500,3);
  59 [450.450023060629,491.1645455437304,448.2236716394652];
  60
  61
  62
  63
  64 /* noncentral chi square distribution */
  65
  66 pdf_noncentral_chi2(13/10,2,8);
  67 %e^-(93/20)*bessel_i(0,2*sqrt(13)/sqrt(5))/2;
  68
  69 cdf_noncentral_chi2(13.1,2,8);
  70 0.7364674817837453;
  71
  72 quantile_noncentral_chi2(0.9,52,1.9);
  73 67.80391983267481;
  74
  75 random_noncentral_chi2(12,6.5,3);
  76 [15.00674323768026,24.34835852498128,17.89531280985209];
  77
  78
  79
  80
  81 /* F distribution */
  82
  83 float(pdf_f(2,3,8/3));
  84 0.1303747102124378;
  85
  86 cdf_f(2,30,5);
  87 1-542856937402935361/(665416609183179841*sqrt(13));
  88
  89 quantile_f(2/3,2,20.1);
  90 1.160908609712148;
  91
  92 random_f(10,31,3);
  93 [1.292249504147881,.7833929757718893,1.312838168471361];
  94
  95
  96
  97
  98 /* exponential distribution is based on gamma */
  99
 100 random_exp(53.21,3);
 101 [.02160755103125626,.03258741884383696,.02294997277868339];
 102
 103 random_exp(3.21,3);
 104 [.4463703709164452,.04871274913659535,.09773054426322608];
 105
 106 random_exp(302.5,3);
 107 [.006696757079634556,.001445392494195325,.002343015491586471];
 108
 109
 110
 111
 112 /* lognormal distribution */
 113
 114 pdf_lognormal(55,-56,3);
 115 %e^-((log(55)+56)^2/18)/(165*sqrt(2)*sqrt(%pi));
 116
 117 cdf_lognormal(3.5,2,1.2);
 118 1/2-erf(0.62269752625386/sqrt(2))/2;
 119
 120 quantile_lognormal(2/5,2.2,0.2);
 121 %e^(0.2*sqrt(2)*inverse_erf(-1/5)+2.2);
 122
 123 random_lognormal(10,15.23,3);
 124 [2.353791985229249E-11,0.187574637106219,2.9149310594062033E-9];
 125
 126
 127
 128
 129 /* gamma distribution */
 130
 131 pdf_gamma(5,2.1,3.5);
 132 0.09686570955466306;
 133
 134 cdf_gamma(5,2,3);
 135 1-gamma_incomplete_regularized(2,5/3);
 136
 137 cdf_gamma(5,2.1,30.5);
 138 0.009139841119737905;
 139
 140 cdf_gamma(5,21,30.5);
 141 0.0;
 142
 143 quantile_gamma(1/6,15,0.2);
 144 2.254909620526077;
 145
 146 random_gamma(1/2,12,3);
 147 [26.21845843519933,1.74922521393547,12.67589121569082];
 148
 149
 150
 151
 152 /* beta distribution */
 153
 154 pdf_beta(1/8,5/8,666/57);
 155 282475249*7^(13/19)/(1073741824*8^(47/152)*beta(5/8,222/19));
 156
 157 cdf_beta(1/2,3,8);
 158 121/128;
 159
 160 quantile_beta(1/8,2,3.2);
 161 0.1534071622832576;
 162
 163 random_beta(1,18,3);
 164 [.01590083925830285,.07944026842342344,.01991813424220185];
 165
 166 random_beta(10,1,3);
 167 [.9768714741444106,.8703877784295725,0.985410331444349];
 168
 169
 170
 171
 172 /* logistic distribution */
 173
 174 random_logistic(-5,6,3);
 175 [-.5460542049418455,-.1005197360489252,-13.85839539776822];
 176
 177 random_logistic(56,60,3);
 178 [16.18610806632373,-24.27721421793832,210.1491259203347];
 179
 180
 181
 182
 183 /* pareto distribution */
 184
 185 random_pareto(56,60,3);
 186 [60.3980718433857,60.02319863764327,61.15889083104624];
 187
 188 random_pareto(6,60.2,3);
 189 [82.00840337695395,68.26980986249082,62.88591356337489];
 190
 191
 192
 193
 194 /* weibull distribution */
 195
 196 kurtosis_weibull(4,5);
 197 (-3*gamma(1/4)*gamma(3/4)/4-3*gamma(1/4)^4/256
 198   +3*sqrt(%pi)*gamma(1/4)^2/16+1) /(sqrt(%pi)/2-gamma(1/4)^2/16)^2 -3;
 199
 200 random_weibull(0.5,5,3);
 201 [2.683098300894825,.3090168279401363,.1236639157364477];
 202
 203 random_weibull(15,0.05,3);
 204 [.04977109163672318,.04870054424691551,.04386166879841022];
 205
 206
 207
 208
 209 /* rayleigh distribution is based on weibull */
 210
 211 pdf_rayleigh(2.5,6);
 212 3.459505910082928E-96;
 213
 214 pdf_rayleigh(5/2,6);
 215 180*%e^-225;
 216
 217 cdf_rayleigh(2.5,6);
 218 1.0;
 219
 220 cdf_rayleigh(5/2,6);
 221 1-%e^-225;
 222
 223
 224
 225
 226 /* laplace distribution */
 227
 228 pdf_laplace(2.5,-5,3);
 229 0.01368083310398313;
 230
 231 cdf_laplace(2/5,-5,30);
 232 (2-%e^-(9/50))/2;
 233
 234 random_laplace(-12,17,3);
 235 [12.6464246671848,-28.88305575541671,-25.390527331018];
 236
 237 random_laplace(-0.02,170,3);
 238 [.1421750180506419,-142.1776353512574,329.313331828124];
 239
 240
 241
 242
 243 /* cauchy distribution */
 244
 245 quantile_cauchy(5/8,15,21);
 246 21*tan(%pi/8)+15;
 247
 248 random_cauchy(1,3,3);
 249 [1.44255833281213,-.4527321151148478,2.236170110804809];
 250
 251 random_cauchy(-65,38,3);
 252 [-66.15930541384257,-53.93341688372076,4099.33722348125];
 253
 254
 255
 256
 257 /* gumbel distribution */
 258
 259 pdf_gumbel(5,-5,3);
 260 %e^(-%e^-(10/3)-10/3)/3;
 261
 262 random_gumbel(-8,6,3);
 263 [-8.522847585460354,6.74085360167784,7.34394651604513];
 264
 265 random_gumbel(80,1/6,3);
 266 [80.07101974101381,80.43036203957068,79.73692916250752];
 267
 268
 269
 270
 271 /* binomial distribution */
 272
 273 pdf_binomial(x,n,1);
 274 kron_delta(n,x);
 275
 276 pdf_binomial(x,n,0);
 277 kron_delta(0,x);
 278
 279 pdf_binomial(0,0,0);
 280 1;
 281
 282 pdf_binomial(0,0,1);
 283 1;
 284
 285 cdf_binomial(3.2,6,1/8);
 286 130683/131072;
 287
 288 cdf_binomial(4,6,1/8);
 289 262101/262144;
 290
 291 quantile_binomial(0.74,23,1/5);
 292 6;
 293
 294 random_binomial(50,1/8,3);
 295 [7,4,6];
 296
 297
 298
 299
 300 /* poisson distribution */
 301
 302 cdf_poisson(13,19/2);
 303 gamma_incomplete_regularized(14,19/2);
 304
 305 cdf_poisson(0.1,0.5);
 306 0.6065306597126336;
 307
 308 quantile_poisson(1/9,10);
 309 6;
 310
 311 random_poisson(110,3);
 312 [105,132,104];
 313
 314 random_poisson(0);
 315 0;
 316
 317
 318 /* geometric distribution */
 319
 320 cdf_geometric(5,1/8);
 321 144495/262144;
 322
 323 random_geometric(1/89,3);
 324 [30,38,58];
 325
 326
 327
 328
 329 /* hypergeometric distribution */
 330
 331 pdf_hypergeometric(4,8,7,5);
 332 70/429;
 333
 334 cdf_hypergeometric(10,18,70,25);
 335 19928255517766/19951839306549;
 336
 337 quantile_hypergeometric(1/21,15,20,30);
 338 11;
 339
 340 random_hypergeometric(19,22,31,3);
 341 [14,14,13];
 342
 343
 344
 345
 346 /* negative binomial distribution */
 347
 348 pdf_negative_binomial(4,7,1/8);
 349 252105/4294967296;
 350
 351 cdf_negative_binomial(4,7,1/8);
 352 425849/4294967296;
 353
 354 quantile_negative_binomial(7/10,10,1/7);
 355 69;
 356
 357 quantile_negative_binomial(1/8,9,6/17);
 358 9;
 359
 360 quantile_negative_binomial(7/1000,2,1/17);
 361 1.0;
 362
 363 random_negative_binomial(6,3/5,3);
 364 [5,8,9];
 365
 366
 367 /* inverse gamma distribution */
 368
 369 (mycontext: newcontext (), 0);
 370 0;
 371
 372 pdf_inverse_gamma (0, x, y);
 373 0;
 374
 375 pdf_inverse_gamma (u, x, y);
 376 if u > 0 then ''(y^x/gamma(x)/u^(1 + x)*exp(-y/u)) else 0;
 377
 378 (assume (u > 0),
 379  pdf_inverse_gamma (u, x, y));
 380 ''(y^x/gamma(x)/u^(1 + x)*exp(-y/u));
 381
 382 (assume (x > 0),
 383  limit (pdf_inverse_gamma (u, x, y), u, inf));
 384 0;
 385
 386 pdf_inverse_gamma (h, 7/2, 20);
 387 if h > 0 then ''(20^(7/2)/gamma(7/2)/h^(9/2)*exp(-20/h)) else 0;
 388
 389 pdf_inverse_gamma (12/7, 7/2, 20);
 390 ''(20^(7/2)/gamma(7/2)/(12/7)^(9/2)*exp(-20/(12/7)));
 391
 392 cdf_inverse_gamma (0, f, g);
 393 0;
 394
 395 cdf_inverse_gamma (w, f, g);
 396 if w > 0 then gamma_incomplete(f, g/w)/gamma(f) else 0;
 397
 398 (assume (w > 0),
 399  cdf_inverse_gamma (w, f, g));
 400 gamma_incomplete(f, g/w)/gamma(f);
 401
 402 cdf_inverse_gamma (5, 3/4, d);
 403 gamma_incomplete(3/4, d/5)/gamma(3/4);
 404
 405 cdf_inverse_gamma (5, 3/4, 11/3);
 406 ''(gamma_incomplete(3/4, (11/3)/5)/gamma(3/4));
 407
 408 (assume (f > 0),
 409  limit (cdf_inverse_gamma (w, f, g), w, inf));
 410 1;
 411
 412 killcontext (mycontext);
 413 done;
 414
 415 quantile_inverse_gamma (v, s, t);
 416 if equal(v, 0) then 0
 417     elseif equal(v, 1) then inf
 418     else 1/quantile_gamma(1 - v, s, 1/t);
 419
 420 quantile_inverse_gamma(u, 1+z, 2-y);
 421 if equal(u,0) then 0 elseif equal(u,1) then inf else 1/quantile_gamma(1-u,z+1,1/(2-y));
 422
 423 quantile_inverse_gamma(1-u^2, 1+a*z, 2 - y/x);
 424 if equal(1-u^2,0) then 0 elseif equal(1-u^2,1) then inf else 1/quantile_gamma(u^2,a*z+1,1/(2-y/x));
 425
 426 quantile_inverse_gamma (0, s, t);
 427 0;
 428
 429 quantile_inverse_gamma (1, s, t);
 430 inf;
 431
 432 (foo(q) := find_root (cdf_inverse_gamma (p, 4.4, 2) = q, p, 0.001, 10*2/3.4),
 433  qq: [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
 434  map (lambda ([q], quantile_inverse_gamma (q, 4.4, 2)), qq));
 435 ''(map (foo, qq));
 436
 437 mean_inverse_gamma (h, i);
 438 i/(h - 1);
 439
 440 mean_inverse_gamma (12, 17);
 441 17/11;
 442
 443 mode_inverse_gamma (u, v);
 444 v/(1 + u);
 445
 446 mode_inverse_gamma (0.5, 1.5);
 447 1.0;
 448
 449 var_inverse_gamma (p, q);
 450 q^2/(p - 1)^2/(p - 2);
 451
 452 var_inverse_gamma (2.5, 0.25);
 453 ''(0.25^2/1.5^2/0.5);
 454
 455 skewness_inverse_gamma (y, z);
 456 4*sqrt(y - 2)/(y - 3);
 457
 458 skewness_inverse_gamma (14.75, 1e6);
 459 ''(4*sqrt(12.75)/11.75);
 460
 461 kurtosis_inverse_gamma (l, m);
 462 6*(5*l - 11)/(l - 3)/(l - 4);
 463
 464 kurtosis_inverse_gamma (4.75, 1e-6);
 465 ''(6*(5*4.75 - 11)/1.75/0.75);
 466
 467 /* simple probabilitistic test for random number generator;
 468  * construct interval with tail mass approximately 1/10^6.
 469  * this is a pretty imprecise test, a test based on the
 470  * Kolgomorov-Smirnov statistic would be more precise.
 471  */
 472 (foo (a, b, n, tail_mass) :=
 473     block ([x, mean, sd, sample_mean, q, low, high],
 474            x: random_inverse_gamma (a, b, n),
 475            mean: mean_inverse_gamma (a, b),
 476            sd: std_inverse_gamma (a, b),
 477            sample_mean: apply ("+", x)/n,
 478            q: quantile_normal (tail_mass/2, 0, 1),
 479            /* bear in mind q < 0 */
 480            low: mean + q*sd/sqrt(n),
 481            high: mean - q*sd/sqrt(n),
 482            if sample_mean > low and sample_mean < high then true
 483                else 'failed_random_inverse_gamma (a, b, n, tail_mass, mean, sd, sample_mean, q, low, high)),
 484  foo (5, 4.25, 1000, 1e-6));
 485 true;
 486
 487 (set_random_state(rnd_state),
 488  float_approx_equal_tolerance : save_tolerance,
 489  0);
 490 0;