### ----------------------------------------------------
### -- Math::Libgsl::RandomDistribution
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-randist.rakutest
### ----------------------------------------------------
### -- Chunk 2 of 2
-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;


