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