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