Help language development. Donate to The Perl Foundation

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

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

use Test;
use NativeCall;
use NativeHelpers::Blob;
use lib 'lib';
use Math::Libgsl::Raw::RandomDistribution :ALL;
use Math::Libgsl::Raw::Random;
use Math::Libgsl::Raw::Matrix :ALL;
use Math::Libgsl::Raw::BLAS :ALL;
use Math::Libgsl::Raw::LinearAlgebra :ALL;
use Math::Libgsl::Constants;

subtest 'gaussian distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gaussian($r, 3e0), 0.4017558243560277, :rel-tol(10⁻¹²), 'gaussian random variate';
  is-approx gsl_ran_gaussian_pdf(1e0, 3e0), 0.12579440923099774, :rel-tol(10⁻¹²), 'gaussian probability density';
  is-approx gsl_ran_gaussian_ziggurat($r, 3e0), 1.3570084962645521, :rel-tol(10⁻¹²), 'gaussian ziggurat';
  is-approx gsl_ran_gaussian_ratio_method($r, 3e0), 4.571692762640136, :rel-tol(10⁻¹²), 'gaussian ratio method';
  is-approx gsl_ran_ugaussian($r), 0.2712258328755056, :rel-tol(10⁻¹²), 'unit gaussian';
  is-approx gsl_ran_ugaussian_pdf(3e0), 0.0044318484119380075, :rel-tol(10⁻¹²), 'unit gaussian probability density';
  is-approx gsl_ran_ugaussian_ratio_method($r), 1.7149190276802055, :rel-tol(10⁻¹²), 'unit gaussian ratio method';
  is-approx gsl_cdf_gaussian_P(1e0, 3e0), 0.6305586598182363, :rel-tol(10⁻¹²), 'cumulative distribution P(x)';
  is-approx gsl_cdf_gaussian_Q(1e0, 3e0), 0.36944134018176367, :rel-tol(10⁻¹²), 'cumulative distribution Q(x)';
  is-approx gsl_cdf_gaussian_Pinv(0.6305586598182363e0, 3e0), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution P(x)';
  is-approx gsl_cdf_gaussian_Qinv(0.36944134018176367e0, 3e0), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution Q(x)';
  is-approx gsl_cdf_ugaussian_P(1e0), 0.8413447460685429, :rel-tol(10⁻¹²), 'unit cumulative distribution P(x)';
  is-approx gsl_cdf_ugaussian_Q(1e0), 0.15865525393145705, :rel-tol(10⁻¹²), 'unit cumulative distribution Q(x)';
  is-approx gsl_cdf_ugaussian_Pinv(0.8413447460685429e0), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution P(x)';
  is-approx gsl_cdf_ugaussian_Qinv(0.15865525393145705e0), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution Q(x)';
  gsl_rng_free($r);
}

subtest 'gaussian tail distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gaussian_tail($r, 1e0, 3e0), 5.023225218761322, :rel-tol(10⁻¹²), 'gaussian tail';
  is-approx gsl_ran_gaussian_tail_pdf(1e0, 1e0, 3e0), 0.3404990063350988, :rel-tol(10⁻¹²), 'gaussian tail probability density';
  is-approx gsl_ran_ugaussian_tail($r, 1e0), 1.198044133837379, :rel-tol(10⁻¹²), 'unit gaussian tail';
  is-approx gsl_ran_ugaussian_tail_pdf(1e0, 1e0), 1.525135276160981, :rel-tol(10⁻¹²), 'unit gaussian tail probability density';
  gsl_rng_free($r);
}

subtest 'bivariate gaussian distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my num64 ($x, $y) = (0e0, 0e0);
  gsl_ran_bivariate_gaussian($r, 3e0, 2e0, .3e0, $x, $y);
  ok (($x, $y) Z≅ (-0.1952914837346669e0, 0.21644212342026375e0)), 'bivariate gaussian';
  is-approx gsl_ran_bivariate_gaussian_pdf(0e0, 0e0, 3e0, 2e0, .3e0), 0.027806618922095617, :rel-tol(10⁻¹²), 'bivariate gaussian probability density';
  gsl_rng_free($r);
}

subtest 'multivariate gaussian distribution' => {
  my $N = 100_000;
  my $d = 2;
  my num64 $t2        = 0e0;
  my num64 $threshold = 0e0;
  my num64 $alpha     = 0.05e0;
  my num64 $pvalue    = 0e0;
  my gsl_vector $mu            = gsl_vector_calloc($d);
  my gsl_matrix $sigma         = gsl_matrix_calloc($d, $d);
  my gsl_matrix $L             = gsl_matrix_calloc($d, $d);
  my gsl_vector $sample        = gsl_vector_calloc($d);
  my gsl_matrix $samples       = gsl_matrix_calloc($N, $d);
  my gsl_vector $mu_hat        = gsl_vector_calloc($d);
  my gsl_matrix $sigma_hat     = gsl_matrix_calloc($d, $d);
  my gsl_vector $mu_hat_ctr    = gsl_vector_calloc($d);
  my gsl_matrix $sigma_hat_inv = gsl_matrix_calloc($d, $d);
  my gsl_vector $tmp           = gsl_vector_calloc($d);
  my gsl_rng    $r             = mgsl_rng_setup(DEFAULT);

  gsl_vector_set($mu, 0, 1e0);
  gsl_vector_set($mu, 1, 2e0);
  gsl_matrix_set($sigma, 0, 0, 4e0);
  gsl_matrix_set($sigma, 0, 1, 2e0);
  gsl_matrix_set($sigma, 1, 0, 2e0);
  gsl_matrix_set($sigma, 1, 1, 3e0);

  gsl_matrix_memcpy($L, $sigma);
  gsl_linalg_cholesky_decomp1($L);

  for ^$N -> $i {
    gsl_ran_multivariate_gaussian($r, $mu, $L, $sample);
    gsl_matrix_set_row($samples, $i, $sample);
  }

  gsl_ran_multivariate_gaussian_mean($samples, $mu_hat);
  gsl_ran_multivariate_gaussian_vcov($samples, $sigma_hat);

  gsl_vector_memcpy($mu_hat_ctr, $mu_hat);
  gsl_vector_sub($mu_hat_ctr, $mu);
  gsl_matrix_memcpy($sigma_hat_inv, $sigma_hat);
  gsl_linalg_cholesky_decomp1($sigma_hat_inv);
  gsl_linalg_cholesky_invert($sigma_hat_inv);
  gsl_blas_dgemv(CblasNoTrans, 1e0, $sigma_hat_inv, $mu_hat_ctr, 0e0, $tmp);
  gsl_blas_ddot($mu_hat_ctr, $tmp, $t2);
  $t2 *= $N;

  my num64 $diff = ($N - $d).Num;
  my num64 $dd   = $d.Num;
  $threshold = ($N - 1e0) * $d / $diff * gsl_cdf_fdist_Pinv(1e0 - $alpha, $dd, $diff);

  ok $t2 > $threshold && $pvalue < $alpha, 'multivariate gaussian, mean, vcov';

  my num64 $obs_res = 0e0;
  my gsl_vector $x  = gsl_vector_calloc($d);

  gsl_ran_multivariate_gaussian_pdf($x, $mu, $L, $obs_res, $tmp);

  is-approx $obs_res, 0.028294217120391e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';

  gsl_ran_multivariate_gaussian_log_pdf($x, $mu, $L, $obs_res, $tmp);

  is-approx $obs_res, -3.565097837249263e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian log with pdf';

  gsl_vector_free($mu);
  gsl_matrix_free($sigma);
  gsl_matrix_free($L);
  gsl_vector_free($sample);
  gsl_matrix_free($samples);
  gsl_vector_free($mu_hat);
  gsl_matrix_free($sigma_hat);
  gsl_vector_free($mu_hat_ctr);
  gsl_matrix_free($sigma_hat_inv);
  gsl_vector_free($tmp);
  gsl_rng_free($r);
}

subtest 'exponential distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_exponential($r, 2e0), 16.523156432740787e0, :abs-tol(10⁻¹⁰), 'exponential';
  is-approx gsl_ran_exponential_pdf(2e0, 2e0), 0.18393972058572117e0, :abs-tol(10⁻¹⁰), 'exponential with pdf';
  is-approx gsl_cdf_exponential_P(2e0, 2e0), 0.6321205588285577e0, :abs-tol(10⁻¹⁰), 'P cumulative exponential';
  is-approx gsl_cdf_exponential_Q(2e0, 2e0), 0.36787944117144233e0, :abs-tol(10⁻¹⁰), 'Q cumulative exponential';
  is-approx gsl_cdf_exponential_Pinv(0.6321205588285577e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative exponential';
  is-approx gsl_cdf_exponential_Qinv(0.36787944117144233e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative exponential';
  gsl_rng_free($r);
}

subtest 'Laplace distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_laplace($r, 2.75e0), 0.001420747954609721e0, :abs-tol(10⁻¹⁰), 'laplace';
  is-approx gsl_ran_laplace_pdf(2e0, 2.75e0), 0.08785910567087736e0, :abs-tol(10⁻¹⁰), 'laplace with pdf';
  is-approx gsl_cdf_laplace_P(2e0, 2.75e0), 0.7583874594050872e0, :abs-tol(10⁻¹⁰), 'P cumulative laplace';
  is-approx gsl_cdf_laplace_Q(2e0, 2.75e0), 0.24161254059491272e0, :abs-tol(10⁻¹⁰), 'Q cumulative laplace';
  is-approx gsl_cdf_laplace_Pinv(0.7583874594050872e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative laplace';
  is-approx gsl_cdf_laplace_Qinv(0.24161254059491272e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative laplace';
  gsl_rng_free($r);
}

subtest 'exponential power distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_exppow($r, 3.7e0, 0.3e0), 96.18418180115754e0, :abs-tol(10⁻¹⁰), 'exppow';
  is-approx gsl_ran_exppow_pdf(2e0, 3.7e0, 0.3e0), 0.006353720526045815e0, :abs-tol(10⁻¹⁰), 'exppow with pdf';
  is-approx gsl_cdf_exppow_P(2e0, 3.7e0, 0.3e0), 0.5155820414218526e0, :abs-tol(10⁻¹⁰), 'P cumulative expow';
  is-approx gsl_cdf_exppow_Q(2e0, 3.7e0, 0.3e0), 0.4844179585781474e0, :abs-tol(10⁻¹⁰), 'Q cumulative expow';
  gsl_rng_free($r);
}

subtest 'Cauchy distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_cauchy($r, 2e0), -0.001622639831190868e0, :abs-tol(10⁻¹⁰), 'cauchy';
  is-approx gsl_ran_cauchy_pdf(2e0, 2e0), 0.07957747154594767e0, :abs-tol(10⁻¹⁰), 'cauchy with pdf';
  is-approx gsl_cdf_cauchy_P(2e0, 2e0), 0.75e0, :abs-tol(10⁻¹⁰), 'P cumulative cauchy';
  is-approx gsl_cdf_cauchy_Q(2e0, 2e0), 0.25e0, :abs-tol(10⁻¹⁰), 'Q cumulative cauchy';
  is-approx gsl_cdf_cauchy_Pinv(.75e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative cauchy';
  is-approx gsl_cdf_cauchy_Qinv(.25e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative cauchy';
  gsl_rng_free($r);
}

subtest 'Rayleigh distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_rayleigh($r, 1.9e0), 0.04318348873449323e0, :abs-tol(10⁻¹⁰), 'rayleigh';
  is-approx gsl_ran_rayleigh_pdf(2e0, 1.9e0), 0.3183584873991547e0, :abs-tol(10⁻¹⁰), 'rayleigh with pdf';
  is-approx gsl_cdf_rayleigh_P(2e0, 1.9e0), 0.4253629302445258e0, :abs-tol(10⁻¹⁰), 'P cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Q(2e0, 1.9e0), 0.5746370697554741e0, :abs-tol(10⁻¹⁰), 'Q cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Pinv(0.4253629302445258e0, 1.9e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Qinv(0.5746370697554741e0, 1.9e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative rayleigh';
  gsl_rng_free($r);
}

subtest 'Rayleigh tail distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_rayleigh_tail($r, 2.7e0, 1.9e0), 2.7003453137884574e0, :abs-tol(10⁻¹⁰), 'rayleigh tail';
  is-approx gsl_ran_rayleigh_tail_pdf(2e0, 2.7e0, 1.9e0), 0e0, :abs-tol(10⁻¹⁰), 'rayleigh tail with pdf';
  gsl_rng_free($r);
}

subtest 'Landau distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_landau($r), 3880.037426254597e0, :abs-tol(10⁻¹⁰), 'landau';
  is-approx gsl_ran_landau_pdf(2e0), 0.10491298949080555e0, :abs-tol(10⁻¹⁰), 'landau with pdf';
  gsl_rng_free($r);
}

subtest 'Levy alpha-Stable distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_levy($r, 5e0, 1e0), 6162.79707164862e0, :abs-tol(10⁻¹⁰), 'levy';
  gsl_rng_free($r);
}

subtest 'Levy skew alpha-Stable distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_levy_skew($r, 5e0, 1e0, 0e0), 6162.79707164862e0, :abs-tol(10⁻¹⁰), 'levy skew';
  gsl_rng_free($r);
}

subtest 'gamma distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gamma($r, 2.5e0, 2.17e0), 7.592395030767021e0, :abs-tol(10⁻¹⁰), 'gamma';
  is-approx gsl_ran_gamma_knuth($r, 2.5e0, 2.17e0), 7.554777177884447e0, :abs-tol(10⁻¹⁰), 'gamma knuth';
  is-approx gsl_ran_gamma_pdf(2e0, 2.5e0, 2.17e0), 0.12203602332361416e0, :abs-tol(10⁻¹⁰), 'gamma with pdf';
  is-approx gsl_cdf_gamma_P(2e0, 2.5e0, 2.17e0), 0.12962769871792748e0, :abs-tol(10⁻¹⁰), 'P cumulative gamma';
  is-approx gsl_cdf_gamma_Q(2e0, 2.5e0, 2.17e0), 0.8703723012820725e0, :abs-tol(10⁻¹⁰), 'Q cumulative gamma';
  is-approx gsl_cdf_gamma_Pinv(0.12962769871792748e0, 2.5e0, 2.17e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gamma';
  is-approx gsl_cdf_gamma_Qinv(0.8703723012820725e0, 2.5e0, 2.17e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gamma';
  gsl_rng_free($r);
}

subtest 'flat distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_flat($r, 3e0, 4e0), 3.999741748906672e0, :abs-tol(10⁻¹⁰), 'flat';
  is-approx gsl_ran_flat_pdf(2e0, 3e0, 4e0), 0e0, :abs-tol(10⁻¹⁰), 'flat with pdf';
  is-approx gsl_cdf_flat_P(2e0, 3e0, 4e0), 0e0, :abs-tol(10⁻¹⁰), 'P cumulative flat';
  is-approx gsl_cdf_flat_Q(2e0, 3e0, 4e0), 1e0, :abs-tol(10⁻¹⁰), 'Q cumulative flat';
  is-approx gsl_cdf_flat_Pinv(0e0, 3e0, 4e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative flat';
  is-approx gsl_cdf_flat_Qinv(1e0, 3e0, 4e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative flat';
  gsl_rng_free($r);
}

subtest 'lognormal distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_lognormal($r, 2.7e0, 3e0), 12.239990687580857e0, :abs-tol(10⁻¹⁰), 'lognormal';
  is-approx gsl_ran_lognormal_pdf(2e0, 2.7e0, 3e0), 0.053160178764929525e0, :abs-tol(10⁻¹⁰), 'lognormal with pdf';
  is-approx gsl_cdf_lognormal_P(2e0, 2.7e0, 3e0), 0.25176338701440004e0, :abs-tol(10⁻¹⁰), 'P cumulative lognormal';
  is-approx gsl_cdf_lognormal_Q(2e0, 2.7e0, 3e0), 0.7482366129856e0, :abs-tol(10⁻¹⁰), 'Q cumulative lognormal';
  is-approx gsl_cdf_lognormal_Pinv(0.25176338701440004e0, 2.7e0, 3e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative lognormal';
  is-approx gsl_cdf_lognormal_Qinv(0.7482366129856e0, 2.7e0, 3e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative lognormal';
  gsl_rng_free($r);
}

subtest 'chi-squared distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_chisq($r, 13e0), 16.53548783112717e0, :abs-tol(10⁻¹⁰), 'chisq';
  is-approx gsl_ran_chisq_pdf(2e0, 13e0), 0.0006389340989638759e0, :abs-tol(10⁻¹⁰), 'chisq with pdf';
  is-approx gsl_cdf_chisq_P(2e0, 13e0), 0.00022625008465604695e0, :abs-tol(10⁻¹⁰), 'P cumulative chisq';
  is-approx gsl_cdf_chisq_Q(2e0, 13e0), 0.999773749915344e0, :abs-tol(10⁻¹⁰), 'Q cumulative chisq';
  is-approx gsl_cdf_chisq_Pinv(0.00022625008465604695e0, 13e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative chisq';
  is-approx gsl_cdf_chisq_Qinv(0.999773749915344e0, 13e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative chisq';
  gsl_rng_free($r);
}

subtest 'f-distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_fdist($r, 3e0, 4e0), 1.265383001140288e0, :abs-tol(10⁻¹⁰), 'fdist';
  is-approx gsl_ran_fdist_pdf(2e0, 3e0, 4e0), 0.13942740046346694e0, :abs-tol(10⁻¹⁰), 'fdist with pdf';
  is-approx gsl_cdf_fdist_P(2e0, 3e0, 4e0), 0.743612802471824e0, :abs-tol(10⁻¹⁰), 'P cumulative fdist';
  is-approx gsl_cdf_fdist_Q(2e0, 3e0, 4e0), 0.2563871975281759e0, :abs-tol(10⁻¹⁰), 'Q cumulative fdist';
  is-approx gsl_cdf_fdist_Pinv(0.743612802471824e0, 3e0, 4e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative fdist';
  is-approx gsl_cdf_fdist_Qinv(0.2563871975281759e0, 3e0, 4e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative fdist';
  gsl_rng_free($r);
}

subtest 't-distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_tdist($r, 1.75e0), 0.2825086397013947e0, :abs-tol(10⁻¹⁰), 'tdist';
  is-approx gsl_ran_tdist_pdf(2e0, 1.75e0), 0.06778157292044809e0, :abs-tol(10⁻¹⁰), 'tdist with pdf';
  is-approx gsl_cdf_tdist_P(2e0, 1.75e0), 0.8993101012459747e0, :abs-tol(10⁻¹⁰), 'P cumulative tdist';
  is-approx gsl_cdf_tdist_Q(2e0, 1.75e0), 0.1006898987540253e0, :abs-tol(10⁻¹⁰), 'Q cumulative tdist';
  is-approx gsl_cdf_tdist_Pinv(0.8993101012459747e0, 1.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative tdist';
  is-approx gsl_cdf_tdist_Qinv(0.1006898987540253e0, 1.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative tdist';
  gsl_rng_free($r);
}

subtest 'beta distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_beta($r, 2e0, 3e0), 0.45158505658193976e0, :abs-tol(10⁻¹⁰), 'beta';
  is-approx gsl_ran_beta_pdf(2e0, 2e0, 3e0), 0e0, :abs-tol(10⁻¹⁰), 'beta with pdf';
  is-approx gsl_cdf_beta_P(2e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'P cumulative beta';
  is-approx gsl_cdf_beta_Q(2e0, 2e0, 3e0), 0e0, :abs-tol(10⁻¹⁰), 'Q cumulative beta';
  is-approx gsl_cdf_beta_Pinv(1e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative beta';
  is-approx gsl_cdf_beta_Qinv(0e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative beta';
  gsl_rng_free($r);
}

subtest 'logistic distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_logistic($r, 3.1e0), 25.61009178896598e0, :abs-tol(10⁻¹⁰), 'logistic';
  is-approx gsl_ran_logistic_pdf(2e0, 3.1e0), 0.07280296434630028e0, :abs-tol(10⁻¹⁰), 'logistic with pdf';
  is-approx gsl_cdf_logistic_P(2e0, 3.1e0), 0.6559192436053651e0, :abs-tol(10⁻¹⁰), 'P cumulative logistic';
  is-approx gsl_cdf_logistic_Q(2e0, 3.1e0), 0.344080756394635e0, :abs-tol(10⁻¹⁰), 'Q cumulative logistic';
  is-approx gsl_cdf_logistic_Pinv(0.6559192436053651e0, 3.1e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative logistic';
  is-approx gsl_cdf_logistic_Qinv(0.344080756394635e0, 3.1e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative logistic';
  gsl_rng_free($r);
}

subtest 'Pareto distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_pareto($r, 1.9e0, 2.75e0), 2.7503738581610317e0, :abs-tol(10⁻¹⁰), 'pareto';
  is-approx gsl_ran_pareto_pdf(3e0, 1.9e0, 2.75e0), 0.5368266659600216e0, :abs-tol(10⁻¹⁰), 'pareto with pdf';
  is-approx gsl_cdf_pareto_P(3e0, 1.9e0, 2.75e0), 0.15237894848417666e0, :abs-tol(10⁻¹⁰), 'P cumulative pareto';
  is-approx gsl_cdf_pareto_Q(3e0, 1.9e0, 2.75e0), 0.8476210515158233e0, :abs-tol(10⁻¹⁰), 'Q cumulative pareto';
  is-approx gsl_cdf_pareto_Pinv(0.15237894848417666e0, 1.9e0, 2.75e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative pareto';
  is-approx gsl_cdf_pareto_Qinv(0.8476210515158233e0, 1.9e0, 2.75e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative pareto';
  gsl_rng_free($r);
}

subtest 'spherical vector distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  my num64 ($x, $y) = (0e0, 0e0);
  gsl_ran_dir_2d($r, $x, $y);
  ok $x ≅ -0.617745613497854e0 && $y ≅ -0.7863779988047479e0, '2D random vector';

  ($x, $y) = (0e0, 0e0);
  gsl_ran_dir_2d_trig_method($r, $x, $y);
  ok $x ≅ 0.11500033916619544e0 && $y ≅ 0.9933654523848008e0, '2D random vector trig method';

  my num64 $z;
  ($x, $y, $z) = (0e0, 0e0, 0e0);
  gsl_ran_dir_3d($r, $x, $y, $z);
  ok $x ≅ -0.024188741233946875e0 && $y ≅ 0.7364240456943486e0 && $z ≅ 0.6760876642275653e0, '3D random vector';

  my $n = 5;
  my $xarr = CArray[num64].new: 0e0 xx $n;
  gsl_ran_dir_nd($r, $n, $xarr);
  ok ($xarr[^$n].list
    Z≅
    (0.1717392856031707e0, 0.547345960645299e0, -0.814904936760339e0, 0.05741225862147077e0, -0.059596927348619114e0)),
    'nD random vector';

  gsl_rng_free($r);
}

subtest 'Weibull distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_weibull($r, 3.14e0, 2.75e0), 0.15568189304354582e0, :abs-tol(10⁻¹⁰), 'weibull';
  is-approx gsl_ran_weibull_pdf(2e0, 3.14e0, 2.75e0), 0.2978229948509798e0, :abs-tol(10⁻¹⁰), 'weibull with pdf';
  is-approx gsl_cdf_weibull_P(2e0, 3.14e0, 2.75e0), 0.25117631509661126e0, :abs-tol(10⁻¹⁰), 'P cumulative weibull';
  is-approx gsl_cdf_weibull_Q(2e0, 3.14e0, 2.75e0), 0.7488236849033888e0, :abs-tol(10⁻¹⁰), 'Q cumulative weibull';
  is-approx gsl_cdf_weibull_Pinv(0.25117631509661126e0, 3.14e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative weibull';
  is-approx gsl_cdf_weibull_Qinv(0.7488236849033888e0, 3.14e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative weibull';
  gsl_rng_free($r);
}

subtest 'type-1 Gumbel distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gumbel1($r, 3.12e0, 4.56e0), 3.134221698863258e0, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gsl_ran_gumbel1_pdf(2e0, 3.12e0, 4.56e0), 0.02749542323890301e0, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gsl_cdf_gumbel1_P(2e0, 3.12e0, 4.56e0), 0.9911480698975662e0, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Q(2e0, 3.12e0, 4.56e0), 0.008851930102433793e0, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Pinv(0.9911480698975662e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Qinv(0.008851930102433793e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel1';
  gsl_rng_free($r);
}

subtest 'type-2 Gumbel distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gumbel2($r, 3.12e0, 4.56e0), 22.970750721534436e0, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gsl_ran_gumbel2_pdf(2e0, 3.12e0, 4.56e0), 0.4842675579059279e0, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gsl_cdf_gumbel2_P(2e0, 3.12e0, 4.56e0), 0.5918470962288894e0, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Q(2e0, 3.12e0, 4.56e0), 0.4081529037711105e0, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Pinv(0.5918470962288894e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Qinv(0.4081529037711105e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel2';
  gsl_rng_free($r);
}

subtest 'Dirichlet distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $K = 2;
  my $alpha = CArray[num64].new: 2.5e0, 5e0;
  my $theta = CArray[num64].new: 0e0, 0e0;

  gsl_ran_dirichlet($r, $K, $alpha, $theta);
  ok ($theta[^$K].list Z≅ (0.37979184993320686e0, 0.6202081500667932e0)), 'dirichlet';

  $theta[0] = .2e0;
  $theta[1] = .8e0;
  ok gsl_ran_dirichlet_pdf($K, $alpha, $theta) ≅ 2.148771883658201e0, 'dirichlet with pdf';

  ok gsl_ran_dirichlet_lnpdf($K, $alpha, $theta) ≅ 0.7648964620298799e0, 'ln dirichlet with pdf';

  gsl_rng_free($r);
}

subtest 'general discrete distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $P = CArray[num64].new: .59e0, .4e0, .01e0;

  my $g = gsl_ran_discrete_preproc(3, $P);
  isa-ok $g, gsl_ran_discrete_t, 'create lookup table';
  ok gsl_ran_discrete($r, $g) ~~ 0|1|2, 'get discrete random number';
  ok gsl_ran_discrete_pdf(2, $g) ≅ .01e0, 'probability of a value';
  lives-ok { gsl_ran_discrete_free($g) }, 'destroy lookup table ';

  gsl_rng_free($r);
}

subtest 'Poisson distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_poisson($r, 30e0) == 34, 'poisson';
  is-approx gsl_ran_poisson_pdf(2, 5e0), 0.08422433748856832e0, :abs-tol(10⁻¹⁰), 'poisson with pdf';
  is-approx gsl_cdf_poisson_P(2, 5e0), 0.12465201948308077e0, :abs-tol(10⁻¹⁰), 'P cumulative poisson';
  is-approx gsl_cdf_poisson_Q(2, 5e0), 0.8753479805169192e0, :abs-tol(10⁻¹⁰), 'Q cumulative poisson';

  gsl_rng_free($r);
}

subtest 'Bernoulli distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_bernoulli($r, .3e0) == 0, 'bernoulli';
  is-approx gsl_ran_bernoulli_pdf(1, .3e0), 0.3e0, :abs-tol(10⁻¹⁰), 'bernoulli with pdf';

  gsl_rng_free($r);
}

subtest 'binomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_binomial($r, .3e0, 5) == 5, 'binomial';
  is-approx gsl_ran_binomial_pdf(2, .3e0, 5), 0.3086999999999999e0, :abs-tol(10⁻¹⁰), 'binomial with pdf';
  is-approx gsl_cdf_binomial_P(2, .3e0, 5), 0.83692e0, :abs-tol(10⁻¹⁰), 'P cumulative binomial';
  is-approx gsl_cdf_binomial_Q(2, .3e0, 5), 0.16307999999999997e0, :abs-tol(10⁻¹⁰), 'Q cumulative binomial';

  gsl_rng_free($r);
}

subtest 'multinomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $K = 3;
  my $sum_n = 100;
  my $p = CArray[num64].new: 2e0, 7e0, 1e0;
  my $n = CArray[uint32].new: 0 xx 3;

  gsl_ran_multinomial($r, $K, $sum_n, $p, $n);
  is-deeply $n.list, (20, 72, 8), 'multinomial';
  is-approx gsl_ran_multinomial_pdf($K, $p, $n), 0.011455694618856253, :abs-tol(10⁻¹⁰), 'multinomial with pdf';
  is-approx gsl_ran_multinomial_lnpdf($K, $p, $n), -4.469268325992729, :abs-tol(10⁻¹⁰), 'ln multinomial with pdf';

  gsl_rng_free($r);
}

subtest 'negative binomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_negative_binomial($r, .3e0, 20e0) == 45, 'negative binomial';
  is-approx gsl_ran_negative_binomial_pdf(20, .3e0, 20e0), 0.0019175722376507298e0, :abs-tol(10⁻¹⁰), 'negative binomial with pdf';
  is-approx gsl_cdf_negative_binomial_P(20, .3e0, 20e0), 0.006254504372434949e0, :abs-tol(10⁻¹⁰), 'P cumulative negative binomial';
  is-approx gsl_cdf_negative_binomial_Q(20, .3e0, 20e0), 0.9937454956275651e0, :abs-tol(10⁻¹⁰), 'Q cumulative negative binomial';

  gsl_rng_free($r);
}

subtest 'Pascal distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_pascal($r, .8e0, 3) == 0, 'pascal';
  is-approx gsl_ran_pascal_pdf(2, .8e0, 3), 0.12287999999999971e0, :abs-tol(10⁻¹⁰), 'pascal with pdf';
  is-approx gsl_cdf_pascal_P(2, .8e0, 3), 0.94208e0, :abs-tol(10⁻¹⁰), 'P cumulative pascal';
  is-approx gsl_cdf_pascal_Q(2, .8e0, 3), 0.057920000000000006e0, :abs-tol(10⁻¹⁰), 'Q cumulative pascal';

  gsl_rng_free($r);
}

subtest 'geometric distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_geometric($r, .5e0) == 1, 'geometric';
  is-approx gsl_ran_geometric_pdf(2, .5e0), .25e0, :abs-tol(10⁻¹⁰), 'geometric with pdf';
  is-approx gsl_cdf_geometric_P(2, .5e0), .75e0, :abs-tol(10⁻¹⁰), 'P cumulative geometric';
  is-approx gsl_cdf_geometric_Q(2, .5e0), .25e0, :abs-tol(10⁻¹⁰), 'Q cumulative geometric';

  gsl_rng_free($r);
}

subtest 'hypergeometric distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_hypergeometric($r, 5, 7, 4) == 2, 'hypergeometric';
  is-approx gsl_ran_hypergeometric_pdf(2, 5, 7, 4), 0.4242424242424253e0, :abs-tol(10⁻¹⁰), 'hypergeometric with pdf';
  is-approx gsl_cdf_hypergeometric_P(2, 5, 7, 4), 0.8484848484848483e0, :abs-tol(10⁻¹⁰), 'P cumulative hypergeometric';
  is-approx gsl_cdf_hypergeometric_Q(2, 5, 7, 4), 0.1515151515151517e0, :abs-tol(10⁻¹⁰), 'Q cumulative hypergeometric';

  gsl_rng_free($r);
}

subtest 'logarithmic distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_logarithmic($r, .4e0) == 1, 'logarithmic';
  is-approx gsl_ran_logarithmic_pdf(2, .4e0), 0.15660921511769743e0, :abs-tol(10⁻¹⁰), 'logarithmic with pdf';

  gsl_rng_free($r);
}

subtest 'Wishart distribution' => {
  my gsl_matrix $V    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $L    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $X    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $L_X  = gsl_matrix_calloc(2, 2);
  my gsl_matrix $work = gsl_matrix_calloc(2, 2);
  my num64 $obs_res;
  my num64 $df = 3e0;
  my $*TOLERANCE = 10⁻¹⁰;

  gsl_matrix_set($V, 0, 0, 1e0);
  gsl_matrix_set($V, 1, 1, 1e0);
  gsl_matrix_set($V, 0, 1, .3e0);
  gsl_matrix_set($V, 1, 0, .3e0);

  gsl_matrix_memcpy($L, $V);
  gsl_linalg_cholesky_decomp1($L);

  gsl_matrix_set($X, 0, 0, 2.213322e0);
  gsl_matrix_set($X, 1, 1, 3.285779e0);
  gsl_matrix_set($X, 0, 1, 1.453357e0);
  gsl_matrix_set($X, 1, 0, 1.453357e0);

  gsl_matrix_memcpy($L_X, $X);
  gsl_linalg_cholesky_decomp1($L_X);

  gsl_ran_wishart_log_pdf($X, $L_X, $df, $L, $obs_res, $work);
  ok $obs_res ≅ -4.931913612377813e0, 'log wishart with pdf';

  gsl_ran_wishart_pdf($X, $L_X, $df, $L, $obs_res, $work);
  ok $obs_res ≅ 0.007212687778224e0, 'wishart with pdf';

  gsl_matrix_free($V);
  gsl_matrix_free($L);
  gsl_matrix_free($X);
  gsl_matrix_free($L_X);
  gsl_matrix_free($work);
}

subtest 'shuffling and sampling' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my CArray[int32] $x .= new: ^10;
  my CArray[int32] $y .= new: ^3;

  gsl_ran_shuffle($r, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $x[^10], (5, 7, 2, 3, 4, 1, 0, 6, 9, 8), 'shuffle';

  gsl_ran_choose($r, pointer-to($y), 3, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $y[^3], (2, 3, 9), 'choose';

  gsl_ran_sample($r, pointer-to($y), 3, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $y[^3], (3, 3, 8), 'sample';

  gsl_rng_free($r);
}

done-testing;