Help language development. Donate to The Perl Foundation

Math::Libgsl::RandomDistribution cpan:FRITH last updated on 2020-07-26

t/02-randist.rakutest
#!/usr/bin/env perl6

use Test;
use lib 'lib';
use Math::Libgsl::Random;
use Math::Libgsl::Matrix;
use Math::Libgsl::BLAS;
use Math::Libgsl::LinearAlgebra;
use Math::Libgsl::Constants;
use Math::Libgsl::RandomDistribution;

subtest 'gaussian distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gaussian($r, 3), 0.4017558243560277, :rel-tol(10⁻¹²), 'gaussian random variate';
  is-approx gaussian-pdf(1, 3), 0.12579440923099774, :rel-tol(10⁻¹²), 'gaussian probability density';
  is-approx gaussian-ziggurat($r, 3), 1.3570084962645521, :rel-tol(10⁻¹²), 'gaussian ziggurat';
  is-approx gaussian-ratio-method($r, 3), 4.571692762640136, :rel-tol(10⁻¹²), 'gaussian ratio method';
  is-approx ugaussian($r), 0.2712258328755056, :rel-tol(10⁻¹²), 'unit gaussian';
  is-approx ugaussian-pdf(3), 0.0044318484119380075, :rel-tol(10⁻¹²), 'unit gaussian probability density';
  is-approx ugaussian-ratio-method($r), 1.7149190276802055, :rel-tol(10⁻¹²), 'unit gaussian ratio method';
  is-approx gaussian-P(1, 3), 0.6305586598182363, :rel-tol(10⁻¹²), 'cumulative distribution P(x)';
  is-approx gaussian-Q(1, 3), 0.36944134018176367, :rel-tol(10⁻¹²), 'cumulative distribution Q(x)';
  is-approx gaussian-Pinv(0.6305586598182363, 3), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution P(x)';
  is-approx gaussian-Qinv(0.36944134018176367, 3), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution Q(x)';
  is-approx ugaussian-P(1), 0.8413447460685429, :rel-tol(10⁻¹²), 'unit cumulative distribution P(x)';
  is-approx ugaussian-Q(1), 0.15865525393145705, :rel-tol(10⁻¹²), 'unit cumulative distribution Q(x)';
  is-approx ugaussian-Pinv(0.8413447460685429), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution P(x)';
  is-approx ugaussian-Qinv(0.15865525393145705), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution Q(x)';
}

subtest 'gaussian tail distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gaussian-tail($r, 1, 3), 5.023225218761322, :rel-tol(10⁻¹²), 'gaussian tail';
  is-approx gaussian-tail-pdf(1, 1, 3), 0.3404990063350988, :rel-tol(10⁻¹²), 'gaussian tail probability density';
  is-approx ugaussian-tail($r, 1), 1.198044133837379, :rel-tol(10⁻¹²), 'unit gaussian tail';
  is-approx ugaussian-tail-pdf(1, 1), 1.525135276160981, :rel-tol(10⁻¹²), 'unit gaussian tail probability density';
}

subtest 'bivariate gaussian distribution' => {
  my $*TOLERANCE = 10⁻¹²;
  my Math::Libgsl::Random $r .= new;
  my ($x, $y) = bivariate-gaussian($r, 3, 2, .3);
  ok (($x, $y) Z≅ (-0.1952914837346669, 0.21644212342026375)), 'bivariate gaussian';
  is-approx bivariate-gaussian-pdf(0e0, 0e0, 3e0, 2e0, .3e0), 0.027806618922095617, :rel-tol(10⁻¹²), 'bivariate gaussian probability density';
}

subtest 'bivariate gaussian distribution' => {
  my $*TOLERANCE = 10⁻¹²;
  my $N = 100;
  my $d = 2;
  my Math::Libgsl::Random $r .= new;
  my Math::Libgsl::Matrix $samples .= new: $N, $d;
  my Math::Libgsl::Matrix $L       .= new: $d, $d;
  my Math::Libgsl::Vector $mu      .= new: $d;
  my Math::Libgsl::Vector $x       .= new: $d;
  $mu[0] = 1; $mu[1] = 2;
  $L[0;0] = 4; $L[0;1] = 2; $L[1;0] = 2; $L[1;1] = 3;
  cholesky-decomp1($L);

  for ^$N -> $i {
    my $sample = multivariate-gaussian($r, $mu, $L);
    $samples.set-row($i, $sample);
  }
  ok ($samples.get-row(0) Z≅ (1.2678372162373517, 2.0093249906121344)), 'multivariate gaussian';
  my $mu-hat = multivariate-gaussian-mean($samples);
  my $L-hat  = multivariate-gaussian-vcov($samples);
  ok ($mu-hat[^2] Z≅ (0.9849998488450891, 1.987771398882359)), 'multivariate gaussian mean';
  ok ((gather for ^2 X ^2 -> ($i, $j) { take $L-hat[$i;$j] })
      Z≅
      (3.981679494074944, 1.981330609263434, 1.981330609263434, 2.9937183347090186)), 'multivariate gaussian vcov';
  is-approx multivariate-gaussian-pdf($x, $mu, $L), 0.028294217120391e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';
  is-approx multivariate-gaussian-log-pdf($x, $mu, $L), -3.565097837249263e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';
}

subtest 'exponential distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx exponential($r, 2), 16.523156432740787, :abs-tol(10⁻¹⁰), 'exponential';
  is-approx exponential-pdf(2, 2), 0.18393972058572117, :abs-tol(10⁻¹⁰), 'exponential with pdf';
  is-approx exponential-P(2, 2), 0.6321205588285577, :abs-tol(10⁻¹⁰), 'P cumulative exponential';
  is-approx exponential-Q(2, 2), 0.36787944117144233, :abs-tol(10⁻¹⁰), 'Q cumulative exponential';
  is-approx exponential-Pinv(0.6321205588285577, 2), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative exponential';
  is-approx exponential-Qinv(0.36787944117144233, 2), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative exponential';
}

subtest 'Laplace distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx laplace($r, 2.75), 0.001420747954609721, :abs-tol(10⁻¹⁰), 'laplace';
  is-approx laplace-pdf(2, 2.75), 0.08785910567087736, :abs-tol(10⁻¹⁰), 'laplace with pdf';
  is-approx laplace-P(2, 2.75), 0.7583874594050872, :abs-tol(10⁻¹⁰), 'P cumulative laplace';
  is-approx laplace-Q(2, 2.75), 0.24161254059491272, :abs-tol(10⁻¹⁰), 'Q cumulative laplace';
  is-approx laplace-Pinv(0.7583874594050872, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative laplace';
  is-approx laplace-Qinv(0.24161254059491272, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative laplace';
}

subtest 'exponential power distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx exppow($r, 3.7, 0.3), 96.18418180115754, :abs-tol(10⁻¹⁰), 'exppow';
  is-approx exppow-pdf(2, 3.7, 0.3), 0.006353720526045815, :abs-tol(10⁻¹⁰), 'exppow with pdf';
  is-approx exppow-P(2, 3.7, 0.3), 0.5155820414218526, :abs-tol(10⁻¹⁰), 'P cumulative expow';
  is-approx exppow-Q(2, 3.7, 0.3), 0.4844179585781474, :abs-tol(10⁻¹⁰), 'Q cumulative expow';
}

subtest 'Cauchy distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx cauchy($r, 2), -0.001622639831190868, :abs-tol(10⁻¹⁰), 'cauchy';
  is-approx cauchy-pdf(2, 2), 0.07957747154594767, :abs-tol(10⁻¹⁰), 'cauchy with pdf';
  is-approx cauchy-P(2, 2), 0.75, :abs-tol(10⁻¹⁰), 'P cumulative cauchy';
  is-approx cauchy-Q(2, 2), 0.25, :abs-tol(10⁻¹⁰), 'Q cumulative cauchy';
  is-approx cauchy-Pinv(.75, 2), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative cauchy';
  is-approx cauchy-Qinv(.25, 2), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative cauchy';
}

subtest 'Rayleigh distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx rayleigh($r, 1.9), 0.04318348873449323, :abs-tol(10⁻¹⁰), 'rayleigh';
  is-approx rayleigh-pdf(2, 1.9), 0.3183584873991547, :abs-tol(10⁻¹⁰), 'rayleigh with pdf';
  is-approx rayleigh-P(2, 1.9), 0.4253629302445258, :abs-tol(10⁻¹⁰), 'P cumulative rayleigh';
  is-approx rayleigh-Q(2, 1.9), 0.5746370697554741, :abs-tol(10⁻¹⁰), 'Q cumulative rayleigh';
  is-approx rayleigh-Pinv(0.4253629302445258, 1.9), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative rayleigh';
  is-approx rayleigh-Qinv(0.5746370697554741, 1.9), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative rayleigh';
}

subtest 'Rayleigh tail distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx rayleigh-tail($r, 2.7, 1.9), 2.7003453137884574, :abs-tol(10⁻¹⁰), 'rayleigh tail';
  is-approx rayleigh-tail-pdf(2, 2.7, 1.9), 0, :abs-tol(10⁻¹⁰), 'rayleigh tail with pdf';
}

subtest 'Landau distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx landau($r), 3880.037426254597, :abs-tol(10⁻¹⁰), 'landau';
  is-approx landau-pdf(2), 0.10491298949080555, :abs-tol(10⁻¹⁰), 'landau with pdf';
}

subtest 'Levy alpha-Stable distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx levy($r, 5, 1), 6162.79707164862, :abs-tol(10⁻¹⁰), 'levy';
}

subtest 'Levy skew alpha-Stable distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx levy-skew($r, 5, 1, 0), 6162.79707164862, :abs-tol(10⁻¹⁰), 'levy skew';
}

subtest 'gamma distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gamma($r, 2.5, 2.17), 7.592395030767021, :abs-tol(10⁻¹⁰), 'gamma';
  is-approx gamma-knuth($r, 2.5, 2.17), 7.554777177884447, :abs-tol(10⁻¹⁰), 'gamma knuth';
  is-approx gamma-pdf(2, 2.5, 2.17), 0.12203602332361416, :abs-tol(10⁻¹⁰), 'gamma with pdf';
  is-approx gamma-P(2, 2.5, 2.17), 0.12962769871792748, :abs-tol(10⁻¹⁰), 'P cumulative gamma';
  is-approx gamma-Q(2, 2.5, 2.17), 0.8703723012820725, :abs-tol(10⁻¹⁰), 'Q cumulative gamma';
  is-approx gamma-Pinv(0.12962769871792748, 2.5, 2.17), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative gamma';
  is-approx gamma-Qinv(0.8703723012820725, 2.5, 2.17), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gamma';
}

subtest 'flat distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx flat($r, 3, 4), 3.999741748906672, :abs-tol(10⁻¹⁰), 'flat';
  is-approx flat-pdf(2, 3, 4), 0, :abs-tol(10⁻¹⁰), 'flat with pdf';
  is-approx flat-P(2, 3, 4), 0, :abs-tol(10⁻¹⁰), 'P cumulative flat';
  is-approx flat-Q(2, 3, 4), 1, :abs-tol(10⁻¹⁰), 'Q cumulative flat';
  is-approx flat-Pinv(0, 3, 4), 3, :abs-tol(10⁻¹⁰), 'inverse P cumulative flat';
  is-approx flat-Qinv(1, 3, 4), 3, :abs-tol(10⁻¹⁰), 'inverse Q cumulative flat';
}

subtest 'lognormal distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx lognormal($r, 2.7, 3), 12.239990687580857, :abs-tol(10⁻¹⁰), 'lognormal';
  is-approx lognormal-pdf(2, 2.7, 3), 0.053160178764929525, :abs-tol(10⁻¹⁰), 'lognormal with pdf';
  is-approx lognormal-P(2, 2.7, 3), 0.25176338701440004, :abs-tol(10⁻¹⁰), 'P cumulative lognormal';
  is-approx lognormal-Q(2, 2.7, 3), 0.7482366129856, :abs-tol(10⁻¹⁰), 'Q cumulative lognormal';
  is-approx lognormal-Pinv(0.25176338701440004, 2.7, 3), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative lognormal';
  is-approx lognormal-Qinv(0.7482366129856, 2.7, 3), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative lognormal';
}

subtest 'chi-squared distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx chisq($r, 13), 16.53548783112717, :abs-tol(10⁻¹⁰), 'chisq';
  is-approx chisq-pdf(2, 13), 0.0006389340989638759, :abs-tol(10⁻¹⁰), 'chisq with pdf';
  is-approx chisq-P(2, 13), 0.00022625008465604695, :abs-tol(10⁻¹⁰), 'P cumulative chisq';
  is-approx chisq-Q(2, 13), 0.999773749915344, :abs-tol(10⁻¹⁰), 'Q cumulative chisq';
  is-approx chisq-Pinv(0.00022625008465604695, 13), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative chisq';
  is-approx chisq-Qinv(0.999773749915344, 13), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative chisq';
}

subtest 'f-distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx fdist($r, 3, 4), 1.265383001140288, :abs-tol(10⁻¹⁰), 'fdist';
  is-approx fdist-pdf(2, 3, 4), 0.13942740046346694, :abs-tol(10⁻¹⁰), 'fdist with pdf';
  is-approx fdist-P(2, 3, 4), 0.743612802471824, :abs-tol(10⁻¹⁰), 'P cumulative fdist';
  is-approx fdist-Q(2, 3, 4), 0.2563871975281759, :abs-tol(10⁻¹⁰), 'Q cumulative fdist';
  is-approx fdist-Pinv(0.743612802471824, 3, 4), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative fdist';
  is-approx fdist-Qinv(0.2563871975281759, 3, 4), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative fdist';
}

subtest 't-distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx tdist($r, 1.75), 0.2825086397013947, :abs-tol(10⁻¹⁰), 'tdist';
  is-approx tdist-pdf(2, 1.75), 0.06778157292044809, :abs-tol(10⁻¹⁰), 'tdist with pdf';
  is-approx tdist-P(2, 1.75), 0.8993101012459747, :abs-tol(10⁻¹⁰), 'P cumulative tdist';
  is-approx tdist-Q(2, 1.75), 0.1006898987540253, :abs-tol(10⁻¹⁰), 'Q cumulative tdist';
  is-approx tdist-Pinv(0.8993101012459747, 1.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative tdist';
  is-approx tdist-Qinv(0.1006898987540253, 1.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative tdist';
}

subtest 'beta distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx beta($r, 2, 3), 0.45158505658193976, :abs-tol(10⁻¹⁰), 'beta';
  is-approx beta-pdf(2, 2, 3), 0, :abs-tol(10⁻¹⁰), 'beta with pdf';
  is-approx beta-P(2, 2, 3), 1, :abs-tol(10⁻¹⁰), 'P cumulative beta';
  is-approx beta-Q(2, 2, 3), 0, :abs-tol(10⁻¹⁰), 'Q cumulative beta';
  is-approx beta-Pinv(1, 2, 3), 1, :abs-tol(10⁻¹⁰), 'inverse P cumulative beta';
  is-approx beta-Qinv(0, 2, 3), 1, :abs-tol(10⁻¹⁰), 'inverse Q cumulative beta';
}

subtest 'logistic distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx logistic($r, 3.1), 25.61009178896598, :abs-tol(10⁻¹⁰), 'logistic';
  is-approx logistic-pdf(2, 3.1), 0.07280296434630028, :abs-tol(10⁻¹⁰), 'logistic with pdf';
  is-approx logistic-P(2, 3.1), 0.6559192436053651, :abs-tol(10⁻¹⁰), 'P cumulative logistic';
  is-approx logistic-Q(2, 3.1), 0.344080756394635, :abs-tol(10⁻¹⁰), 'Q cumulative logistic';
  is-approx logistic-Pinv(0.6559192436053651, 3.1), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative logistic';
  is-approx logistic-Qinv(0.344080756394635, 3.1), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative logistic';
}

subtest 'Pareto distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx pareto($r, 1.9, 2.75), 2.7503738581610317, :abs-tol(10⁻¹⁰), 'pareto';
  is-approx pareto-pdf(3, 1.9, 2.75), 0.5368266659600216, :abs-tol(10⁻¹⁰), 'pareto with pdf';
  is-approx pareto-P(3, 1.9, 2.75), 0.15237894848417666, :abs-tol(10⁻¹⁰), 'P cumulative pareto';
  is-approx pareto-Q(3, 1.9, 2.75), 0.8476210515158233, :abs-tol(10⁻¹⁰), 'Q cumulative pareto';
  is-approx pareto-Pinv(0.15237894848417666, 1.9, 2.75), 3, :abs-tol(10⁻¹⁰), 'inverse P cumulative pareto';
  is-approx pareto-Qinv(0.8476210515158233, 1.9, 2.75), 3, :abs-tol(10⁻¹⁰), 'inverse Q cumulative pareto';
}

subtest 'spherical vector distribution' => {
  my Math::Libgsl::Random $r .= new;
  my $*TOLERANCE = 10⁻¹⁰;
  my ($x, $y) = d2dir($r);
  ok $x ≅ -0.617745613497854 && $y ≅ -0.7863779988047479, '2D random vector';
  ($x, $y) = d2dir-trig-method($r);
  ok $x ≅ 0.11500033916619544 && $y ≅ 0.9933654523848008, '2D random vector trig method';
  my $z;
  ($x, $y, $z) = d3dir($r);
  ok $x ≅ -0.024188741233946875 && $y ≅ 0.7364240456943486 && $z ≅ 0.6760876642275653, '3D random vector';
  my @ret = dndir($r, 5);
  ok (@ret[^5]
    Z≅
    (0.1717392856031707, 0.547345960645299, -0.814904936760339, 0.05741225862147077, -0.059596927348619114)),
    'nD random vector';
}

subtest 'Weibull distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx weibull($r, 3.14, 2.75), 0.15568189304354582, :abs-tol(10⁻¹⁰), 'weibull';
  is-approx weibull-pdf(2, 3.14, 2.75), 0.2978229948509798, :abs-tol(10⁻¹⁰), 'weibull with pdf';
  is-approx weibull-P(2, 3.14, 2.75), 0.25117631509661126, :abs-tol(10⁻¹⁰), 'P cumulative weibull';
  is-approx weibull-Q(2, 3.14, 2.75), 0.7488236849033888, :abs-tol(10⁻¹⁰), 'Q cumulative weibull';
  is-approx weibull-Pinv(0.25117631509661126, 3.14, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative weibull';
  is-approx weibull-Qinv(0.7488236849033888, 3.14, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative weibull';
}

subtest 'type-1 Gumbel distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gumbel1($r, 3.12, 4.56), 3.134221698863258, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gumbel1-pdf(2, 3.12, 4.56), 0.02749542323890301, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gumbel1-P(2, 3.12, 4.56), 0.9911480698975662, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gumbel1-Q(2, 3.12, 4.56), 0.008851930102433793, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gumbel1-Pinv(0.9911480698975662, 3.12, 4.56), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gumbel1-Qinv(0.008851930102433793, 3.12, 4.56), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel1';
}

subtest 'type-2 Gumbel distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gumbel2($r, 3.12, 4.56), 22.970750721534436, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gumbel2-pdf(2, 3.12, 4.56), 0.4842675579059279, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gumbel2-P(2, 3.12, 4.56), 0.5918470962288894, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gumbel2-Q(2, 3.12, 4.56), 0.4081529037711105, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gumbel2-Pinv(0.5918470962288894, 3.12, 4.56), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gumbel2-Qinv(0.4081529037711105, 3.12, 4.56), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel2';
}

subtest 'Dirichlet distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my Math::Libgsl::Random $r .= new;
  my $K = 2;
  my @alpha = 2.5, 5;

  my @theta = dirichlet($r, $K, @alpha);
  ok (@theta[^$K] Z≅ (0.37979184993320686, 0.6202081500667932)), 'dirichlet';

  @theta = .2, .8;
  ok dirichlet-pdf($K, @alpha, @theta) ≅ 2.148771883658201, 'dirichlet with pdf';
  ok dirichlet-lnpdf($K, @alpha, @theta) ≅ 0.7648964620298799, 'ln dirichlet with pdf';
}

subtest 'general discrete distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my Math::Libgsl::Random $r .= new;
  my Int $size = 3;
  my @probability = .59, .4, .01;

  my Math::Libgsl::RandomDistribution::Discrete $d .= new: :$size, :@probability;
  isa-ok $d, Math::Libgsl::RandomDistribution::Discrete, 'crete a discrete distribution object';

  ok $d.discrete($r) ~~ 0|1|2, 'get discrete random number';
  ok $d.discrete-pdf(2) ≅ .01e0, 'probability of a value';
  ok $d.size == $size, 'read size from object';
  is-deeply $d.probability, [0.59, 0.4, 0.01], 'read probability from object';
}

subtest 'Poisson distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok poisson($r, 30) == 34, 'poisson';
  is-approx poisson-pdf(2, 5), 0.08422433748856832, :abs-tol(10⁻¹⁰), 'poisson with pdf';
  is-approx poisson-P(2, 5), 0.12465201948308077, :abs-tol(10⁻¹⁰), 'P cumulative poisson';
  is-approx poisson-Q(2, 5), 0.8753479805169192, :abs-tol(10⁻¹⁰), 'Q cumulative poisson';
}

subtest 'Bernoulli distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok bernoulli($r, .3) == 0, 'bernoulli';
  is-approx bernoulli-pdf(1, .3), 0.3, :abs-tol(10⁻¹⁰), 'bernoulli with pdf';
}

subtest 'binomial distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok binomial($r, .3, 5) == 5, 'binomial';
  is-approx binomial-pdf(2, .3, 5), 0.3086999999999999, :abs-tol(10⁻¹⁰), 'binomial with pdf';
  is-approx binomial-P(2, .3, 5), 0.83692, :abs-tol(10⁻¹⁰), 'P cumulative binomial';
  is-approx binomial-Q(2, .3, 5), 0.16307999999999997, :abs-tol(10⁻¹⁰), 'Q cumulative binomial';
}

subtest 'multinomial distribution' => {
  my Math::Libgsl::Random $r .= new;
  my $K = 3;
  my $sum_n = 100;
  my @p = 2, 7, 1;

  my @n = multinomial($r, $K, $sum_n, @p);
  is-deeply @n, [20, 72, 8], 'multinomial';
  is-approx multinomial-pdf($K, @p, @n), 0.011455694618856253, :abs-tol(10⁻¹⁰), 'multinomial with pdf';
  is-approx multinomial-lnpdf($K, @p, @n), -4.469268325992729, :abs-tol(10⁻¹⁰), 'ln multinomial with pdf';
}

subtest 'negative binomial distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok negative-binomial($r, .3, 20) == 45, 'negative binomial';
  is-approx negative-binomial-pdf(20, .3, 20), 0.0019175722376507298, :abs-tol(10⁻¹⁰), 'negative binomial with pdf';
  is-approx negative-binomial-P(20, .3, 20), 0.006254504372434949, :abs-tol(10⁻¹⁰), 'P cumulative negative binomial';
  is-approx negative-binomial-Q(20, .3, 20), 0.9937454956275651, :abs-tol(10⁻¹⁰), 'Q cumulative negative binomial';
}

subtest 'Pascal distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok pascal($r, .8, 3) == 0, 'pascal';
  is-approx pascal-pdf(2, .8, 3), 0.12287999999999971, :abs-tol(10⁻¹⁰), 'pascal with pdf';
  is-approx pascal-P(2, .8, 3), 0.94208, :abs-tol(10⁻¹⁰), 'P cumulative pascal';
  is-approx pascal-Q(2, .8, 3), 0.057920000000000006, :abs-tol(10⁻¹⁰), 'Q cumulative pascal';
}

subtest 'geometric distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok geometric($r, .5) == 1, 'geometric';
  is-approx geometric-pdf(2, .5), .25, :abs-tol(10⁻¹⁰), 'geometric with pdf';
  is-approx geometric-P(2, .5), .75, :abs-tol(10⁻¹⁰), 'P cumulative geometric';
  is-approx geometric-Q(2, .5), .25, :abs-tol(10⁻¹⁰), 'Q cumulative geometric';
}

subtest 'hypergeometric distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok hypergeometric($r, 5, 7, 4) == 2, 'hypergeometric';
  is-approx hypergeometric-pdf(2, 5, 7, 4), 0.4242424242424253, :abs-tol(10⁻¹⁰), 'hypergeometric with pdf';
  is-approx hypergeometric-P(2, 5, 7, 4), 0.8484848484848483, :abs-tol(10⁻¹⁰), 'P cumulative hypergeometric';
  is-approx hypergeometric-Q(2, 5, 7, 4), 0.1515151515151517, :abs-tol(10⁻¹⁰), 'Q cumulative hypergeometric';
}

subtest 'logarithmic distribution' => {
  my Math::Libgsl::Random $r .= new;
  ok logarithmic($r, .4) == 1, 'logarithmic';
  is-approx logarithmic-pdf(2, .4), 0.15660921511769743, :abs-tol(10⁻¹⁰), 'logarithmic with pdf';
}

subtest 'Wishart distribution' => {
  my Math::Libgsl::Random $r   .= new;
  my Math::Libgsl::Matrix $L   .= new: 2, 2;
  my Math::Libgsl::Matrix $X   .= new: 2, 2;
  my Math::Libgsl::Matrix $L-X .= new: 2, 2;
  my $df = 3;
  my $*TOLERANCE = 10⁻¹⁰;

  $X[0;0] = 2.213322; $X[1;1] = 3.285779; $X[0;1] = 1.453357; $X[1;0] = 1.453357;
  $L[0;0] = 1;        $L[1;1] = 1;        $L[0;1] = .3;       $L[1;0] = .3;
  $L-X.copy($X);
  cholesky-decomp1($L);
  cholesky-decomp1($L-X);

  my $res = wishart-log-pdf($X, $L-X, $df, $L);
  ok $res ≅ -4.931913612377813e0, 'log wishart with pdf';

  $res = wishart-pdf($X, $L-X, $df, $L);
  ok $res ≅ 0.007212687778224e0, 'wishart with pdf';
}

subtest 'shuffling and sampling' => {
  my Math::Libgsl::Random $r .= new;
  my @x = ^10;

  my @y = shuffle($r, @x) for ^10;
  is-deeply @y, [5, 7, 4, 0, 9, 3, 2, 8, 6, 1], 'shuffle';

  my @elems = ^3;
  @y = choose($r, 3, @x);
  is-deeply @y, [2, 3, 5], 'choose';

  @y = sample($r, 3, @x);
  is-deeply @y, [2, 3, 4], 'sample';
}

done-testing;