### ----------------------------------------------------
### -- Math::Libgsl::RandomDistribution
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: lib/Math/Libgsl/RandomDistribution.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 3
omial_Q($k, $p, $n);
}
# Pascal Distribution
sub pascal(Math::Libgsl::Random $r, Num() $p, UInt $n --> UInt) is export
{
  gsl_ran_pascal($r.r, $p, $n);
}
sub pascal-pdf(UInt $k, Num() $p, UInt $n --> Num) is export
{
  gsl_ran_pascal_pdf($k, $p, $n);
}
sub pascal-P(UInt $k, Num() $p, UInt $n --> Num) is export
{
  gsl_cdf_pascal_P($k, $p, $n);
}
sub pascal-Q(UInt $k, Num() $p, UInt $n --> Num) is export
{
  gsl_cdf_pascal_Q($k, $p, $n);
}
# Geometric Distribution
sub geometric(Math::Libgsl::Random $r, Num() $p --> UInt) is export
{
  gsl_ran_geometric($r.r, $p);
}
sub geometric-pdf(UInt $k where * ≥ 1, Num() $p --> Num) is export
{
  gsl_ran_geometric_pdf($k, $p);
}
sub geometric-P(UInt $k where * ≥ 1, Num() $p --> Num) is export
{
  gsl_cdf_geometric_P($k, $p);
}
sub geometric-Q(UInt $k where * ≥ 1, Num() $p --> Num) is export
{
  gsl_cdf_geometric_Q($k, $p);
}
# Hypergeometric Distribution
sub hypergeometric(Math::Libgsl::Random $r, UInt $n1, UInt $n2, UInt $t --> UInt) is export
{
  gsl_ran_hypergeometric($r.r, $n1, $n2, $t);
}
sub hypergeometric-pdf(UInt $k, UInt $n1, UInt $n2, UInt $t --> Num) is export
{
  gsl_ran_hypergeometric_pdf($k, $n1, $n2, $t);
}
sub hypergeometric-P(UInt $k, UInt $n1, UInt $n2, UInt $t --> Num) is export
{
  gsl_cdf_hypergeometric_P($k, $n1, $n2, $t);
}
sub hypergeometric-Q(UInt $k, UInt $n1, UInt $n2, UInt $t --> Num) is export
{
  gsl_cdf_hypergeometric_Q($k, $n1, $n2, $t);
}
# Logarithmic Distribution
sub logarithmic(Math::Libgsl::Random $r, Num() $p --> UInt) is export
{
  gsl_ran_logarithmic($r.r, $p);
}
sub logarithmic-pdf(UInt $k where * ≥ 1, Num() $p --> Num) is export
{
  gsl_ran_logarithmic_pdf($k, $p);
}
# Wishart Distribution
sub wishart(Math::Libgsl::Random $r,
            Num() $n,
            Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 && $n ≤ $L.matrix.size1 - 1 }
            --> Math::Libgsl::Matrix) is export
{
  my $N = $L.matrix.size1;
  my Math::Libgsl::Matrix $result .= new: $N, $N;
  my Math::Libgsl::Matrix $work   .= new: $N, $N;
  my $ret = gsl_ran_wishart($r.r, $n, $L.matrix, $result.matrix, $work.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in wishart" if $ret ≠ GSL_SUCCESS;
  return $result;
}
sub wishart-pdf(Math::Libgsl::Matrix $X where { $X.matrix.size1 == $X.matrix.size2 },
                Math::Libgsl::Matrix $L-X where { $L-X.matrix.size1 == $L-X.matrix.size2 },
                Num() $n where * > $X.matrix.size1 - 1,
                Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 &&
                                                $X.matrix.size1 == $L.matrix.size1 &&
                                                $L-X.matrix.size1 == $L.matrix.size1
                                              }
                --> Num) is export
{
  my $N = $L.matrix.size1;
  my Math::Libgsl::Matrix $work .= new: $N, $N;
  my num64 $result;
  my $ret = gsl_ran_wishart_pdf($X.matrix, $L-X.matrix, $n, $L.matrix, $result, $work.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in wishart-pdf" if $ret ≠ GSL_SUCCESS;
  return $result;
}
sub wishart-log-pdf(Math::Libgsl::Matrix $X where { $X.matrix.size1 == $X.matrix.size2 },
                Math::Libgsl::Matrix $L-X where { $L-X.matrix.size1 == $L-X.matrix.size2 },
                Num() $n where * > $X.matrix.size1 - 1,
                Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 &&
                                                $X.matrix.size1 == $L.matrix.size1 &&
                                                $L-X.matrix.size1 == $L.matrix.size1
                                              }
                --> Num) is export
{
  my $N = $L.matrix.size1;
  my Math::Libgsl::Matrix $work .= new: $N, $N;
  my num64 $result;
  my $ret = gsl_ran_wishart_log_pdf($X.matrix, $L-X.matrix, $n, $L.matrix, $result, $work.matrix);
  fail X::Libgsl.new: errno => $ret, error => "Error in wishart-log-pdf" if $ret ≠ GSL_SUCCESS;
  return $result;
}
# Shuffling and Sampling
sub shuffle(Math::Libgsl::Random $r, @base --> List) is export
{
  my $type = @base[0].WHAT;

  given $type {
    when Int|UInt {
      my CArray[uint64] $base .= new: @base;
      gsl_ran_shuffle($r.r, pointer-to($base), @base.elems, nativesizeof(uint64));
      return $base.list;
    }
    when Num|Rat {
      my CArray[num64] $base .= new: @base».Num;
      gsl_ran_shuffle($r.r, pointer-to($base), @base.elems, nativesizeof(num64));
      return $base.list;
    }
  }
}
sub choose(Math::Libgsl::Random $r, Int $k, @src where *.elems ≥ $k --> List) is export
{
  my $type = @src[0].WHAT;

  given $type {
    when Int|UInt {
      my CArray[uint64] $src .= new: @src;
      my $out = CArray[uint64].allocate: $k;
      gsl_ran_choose($r.r, pointer-to($out), $k, pointer-to($src), @src.elems, nativesizeof(uint64));
      return $out.list;
    }
    when Num|Rat {
      my CArray[num64] $src .= new: @src».Num;
      my $out = CArray[num64].allocate: $k;
      gsl_ran_choose($r.r, pointer-to($out), $k, pointer-to($src), @src.elems, nativesizeof(num64));
      return $out.list;
    }
  }
}
sub sample(Math::Libgsl::Random $r, Int $k, @src --> List) is export
{
  my $type = @src[0].WHAT;

  given $type {
    when Int|UInt {
      my CArray[uint64] $src .= new: @src;
      my $out = CArray[uint64].allocate: $k;
      gsl_ran_sample($r.r, pointer-to($out), $k, pointer-to($src), @src.elems, nativesizeof(uint64));
      return $out.list;
    }
    when Num|Rat {
      my CArray[num64] $src .= new: @src».Num;
      my $out = CArray[num64].allocate: $k;
      gsl_ran_sample($r.r, pointer-to($out), $k, pointer-to($src), @src.elems, nativesizeof(num64));
      return $out.list;
    }
  }
}

=begin pod

=head1 NAME

Math::Libgsl::RandomDistribution - An interface to libgsl, the Gnu Scientific Library - Random Number Distributions

=head1 SYNOPSIS

=begin code :lang<raku>

use Math::Libgsl::Random;
use Math::Libgsl::RandomDistribution;

my Math::Libgsl::Random $r .= new;
say poisson($r, 3) for ^10;

=end code

=head1 DESCRIPTION

Math::Libgsl::RandomDistribution is an interface to the Random Number Distributions section of libgsl, the Gnu Scientific Library.

=head3 gaussian(Math::Libgsl::Random $r, Num() $sigma --> Num)

This function returns a Gaussian random variate, with mean zero and standard deviation $sigma.

=head3 gaussian-pdf(Num() $x, Num() $sigma --> Num)

This function computes the probability density at $x for a Gaussian distribution with standard deviation $sigma.

=head3 gaussian-ziggurat(Math::Libgsl::Random $r, Num() $sigma --> Num)
=head3 gaussian-ratio-method(Math::Libgsl::Random $r, Num() $sigma --> Num)

This function computes a Gaussian random variate using the alternative Marsaglia-Tsang ziggurat and Kinderman-Monahan-Leva ratio methods.

=head3 ugaussian(Math::Libgsl::Random $r --> Num)
=head3 ugaussian-pdf(Num() $x --> Num)
=head3 ugaussian-ratio-method(Math::Libgsl::Random $r --> Num)

These functions compute results for the unit Gaussian distribution. They are equivalent to the functions above with $sigma = 1.

=head3 gaussian-P(Num() $x, Num() $sigma --> Num)
=head3 gaussian-Q(Num() $x, Num() $sigma --> Num)
=head3 gaussian-Pinv(Num() $P, Num() $sigma --> Num)
=head3 gaussian-Qinv(Num() $Q, Num() $sigma --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses.

=head3 ugaussian-P(Num() $x --> Num)
=head3 ugaussian-Q(Num() $x --> Num)
=head3 ugaussian-Pinv(Num() $P --> Num)
=head3 ugaussian-Qinv(Num() $Q --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the unit Gaussian distribution.

=head3 gaussian-tail(Math::Libgsl::Random $r, Num() $a where * > 0, Num() $sigma --> Num)

This function provides random variates from the upper tail of a Gaussian distribution with standard deviation $sigma.

=head3 gaussian-tail-pdf(Num() $x, Num() $a where * > 0, Num() $sigma --> Num)

This function computes the probability density p(x) at $x for a Gaussian tail distribution with standard deviation $sigma and lower limit $a.

=head3 ugaussian-tail(Math::Libgsl::Random $r, Num() $a where * > 0 --> Num)
=head3 ugaussian-tail-pdf(Num() $x, Num() $a where * > 0 --> Num)

These functions compute results for the tail of a unit Gaussian distribution.

=head3 bivariate-gaussian(Math::Libgsl::Random $r, Num() $sigma-x, Num() $sigma-y, Num() $rho where -1 < * < 1 --> List)

This function generates a pair of correlated Gaussian variates, with mean zero, correlation coefficient $rho and
standard deviations $sigma_x and $sigma_y.
The function returns two Num values: the Gauusian correlates along the x and y directions.

=head3 bivariate-gaussian-pdf(Num() $x, Num() $y, Num() $sigma-x, Num() $sigma-y, Num() $rho where -1 < * < 1 --> Num)

This function computes the probability density p(x, y) at ($x, $y) for a bivariate Gaussian distribution with standard deviations $sigma_x, $sigma_y and correlation coefficient $rho.

=head3 multivariate-gaussian(Math::Libgsl::Random $r, Math::Libgsl::Vector $mu, Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 && $L.matrix.size1 == $mu.vector.size } --> Math::Libgsl::Vector)

This function generates a random vector satisfying the k-dimensional multivariate Gaussian distribution with mean μ and variance-covariance matrix Σ. On input, the k-vector μ is given in $mu, and the Cholesky factor of the k-by-k matrix Σ = LT(L) is given in the lower triangle of $L, as output from Math::Libgsl::LinearAlgebra::cholesky-decomp().
The function returns the random vector as a Math::Libgsl::Vector object.

=head3 multivariate-gaussian-pdf(Math::Libgsl::Vector $x, Math::Libgsl::Vector $mu where *.vector.size == $x.vector.size, Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 && $L.matrix.size1 == $mu.vector.size } --> Num)
=head3 multivariate-gaussian-log-pdf(Math::Libgsl::Vector $x, Math::Libgsl::Vector $mu where *.vector.size == $x.vector.size, Math::Libgsl::Matrix $L where { $L.matrix.size1 == $L.matrix.size2 && $L.matrix.size1 == $mu.vector.size } --> Num)

These functions compute p(x) or log p(x) at the point $x, using mean vector $mu and variance-covariance matrix specified by its Cholesky factor $L.

=head3 multivariate-gaussian-mean(Math::Libgsl::Matrix $X --> Math::Libgsl::Vector)

Given a set of n samples Xⱼ from a k-dimensional multivariate Gaussian distribution, this function computes the maximum likelihood estimate of the mean of the distribution.
The function returns the maximum likelihood estimate as a Math::Libgsl::Vector object.

=head3 multivariate-gaussian-vcov(Math::Libgsl::Matrix $X --> Math::Libgsl::Matrix)

Given a set of n samples Xⱼ from a k-dimensional multivariate Gaussian distribution, this function computes the maximum likelihood estimate of the variance-covariance matrix of the distribution.
The function returns the maximum likelihood estimate of the variance-covariance matrix as a Math::Libgsl::Matrix object.

=head3 exponential(Math::Libgsl::Random $r, Num() $mu --> Num)

This function returns a random variate from the exponential distribution with mean $mu.

=head3 exponential-pdf(Num() $x, Num() $mu --> Num)

This function computes the probability density p(x) at $x for an exponential distribution with mean $mu.

=head3 exponential-P(Num() $x, Num() $mu --> Num)
=head3 exponential-Q(Num() $x, Num() $mu --> Num)
=head3 exponential-Pinv(Num() $P, Num() $mu --> Num)
=head3 exponential-Qinv(Num() $Q, Num() $mu --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the exponential distribution with mean $mu.

=head3 laplace(Math::Libgsl::Random $r, Num() $a --> Num)

This function returns a random variate from the Laplace distribution with width $a.

=head3 laplace-pdf(Num() $x, Num() $a --> Num)

This function computes the probability density p(x) at $x for a Laplace distribution with width $a.

=head3 laplace-P(Num() $x, Num() $a --> Num)
=head3 laplace-Q(Num() $x, Num() $a --> Num)
=head3 laplace-Pinv(Num() $P, Num() $a --> Num)
=head3 laplace-Qinv(Num() $Q, Num() $a --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Laplace distribution with width $a.

=head3 exppow(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the exponential power distribution with scale parameter $a and exponent $b.

=head3 exppow-pdf(Num() $x, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for an exponential power distribution with scale parameter $a and exponent $b.

=head3 exppow-P(Num() $x, Num() $a, Num() $b --> Num)
=head3 exppow-Q(Num() $x, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) for the exponential power distribution with parameters $a and $b.

=head3 cauchy(Math::Libgsl::Random $r, Num() $a --> Num)

This function returns a random variate from the Cauchy distribution with scale parameter $a.

=head3 cauchy-pdf(Num() $x, Num() $a --> Num)

This function computes the probability density p(x) at $x for a Cauchy distribution with scale parameter $a.

=head3 cauchy-P(Num() $x, Num() $a --> Num)
=head3 cauchy-Q(Num() $x, Num() $a --> Num)
=head3 cauchy-Pinv(Num() $P, Num() $a --> Num)
=head3 cauchy-Qinv(Num() $Q, Num() $a --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Cauchy distribution with scale parameter $a.

=head3 rayleigh(Math::Libgsl::Random $r, Num() $sigma --> Num)

This function returns a random variate from the Rayleigh distribution with scale parameter $sigma.

=head3 rayleigh-pdf(Num() $x, Num() $sigma --> Num)

This function computes the probability density p(x) at $x for a Rayleigh distribution with scale parameter $sigma.

=head3 rayleigh-P(Num() $x, Num() $sigma --> Num)
=head3 rayleigh-Q(Num() $x, Num() $sigma --> Num)
=head3 rayleigh-Pinv(Num() $P, Num() $sigma --> Num)
=head3 rayleigh-Qinv(Num() $Q, Num() $sigma --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Rayleigh distribution with scale parameter $sigma.

=head3 rayleigh-tail(Math::Libgsl::Random $r, Num() $a, Num() $sigma --> Num)

This function returns a random variate from the tail of the Rayleigh distribution with scale parameter $sigma and a lower limit of $a.

=head3 rayleigh-tail-pdf(Num() $x, Num() $a, Num() $sigma --> Num)

This function computes the probability density p(x) at $x for a Rayleigh tail distribution with scale parameter $sigma and lower limit $a.

=head3 landau(Math::Libgsl::Random $r --> Num)

This function returns a random variate from the Landau distribution.

=head3 landau-pdf(Num() $x --> Num)

This function computes the probability density p(x) at $x for the Landau distribution.

=head3 levy(Math::Libgsl::Random $r, Num() $c, Num() $alpha where 0 < * ≤ 2 --> Num)

This function returns a random variate from the Levy symmetric stable distribution with scale $c and exponent $alpha.

=head3 levy-skew(Math::Libgsl::Random $r, Num() $c, Num() $alpha where 0 < * ≤ 2, Num() $beta --> Num)

This function returns a random variate from the Levy skew stable distribution with scale $c, exponent $alpha and skewness parameter $beta.

=head3 gamma(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the gamma distribution.

=head3 gamma-knuth(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a gamma variate using the algorithms from Knuth.

=head3 gamma-pdf(Num() $x, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a gamma distribution with parameters $a and $b.

=head3 gamma-P(Num() $x, Num() $a, Num() $b --> Num)
=head3 gamma-Q(Num() $x, Num() $a, Num() $b --> Num)
=head3 gamma-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 gamma-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the gamma distribution with parameters $a and $b.

=head3 flat(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the flat (uniform) distribution from $a to $b.

=head3 flat-pdf(Num() $x, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a uniform distribution from $a to $b.

=head3 flat-P(Num() $x, Num() $a, Num() $b --> Num)
=head3 flat-Q(Num() $x, Num() $a, Num() $b --> Num)
=head3 flat-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 flat-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for a uniform distribution from $a to $b.

=head3 lognormal(Math::Libgsl::Random $r, Num() $zeta, Num() $sigma --> Num)

This function returns a random variate from the lognormal distribution.

=head3 lognormal-pdf(Num() $x, Num() $zeta, Num() $sigma --> Num)

This function computes the probability density p(x) at $x for a lognormal distribution with parameters $zeta and $sigma.

=head3 lognormal-P(Num() $x, Num() $zeta, Num() $sigma --> Num)
=head3 lognormal-Q(Num() $x, Num() $zeta, Num() $sigma --> Num)
=head3 lognormal-Pinv(Num() $P, Num() $zeta, Num() $sigma --> Num)
=head3 lognormal-Qinv(Num() $Q, Num() $zeta, Num() $sigma --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the lognormal distribution with parameters $zeta and $sigma.

=head3 chisq(Math::Libgsl::Random $r, Num() $nu --> Num)

This function returns a random variate from the chi-squared distribution with $nu degrees of freedom.

=head3 chisq-pdf(Num() $x where * ≥ 0, Num() $nu --> Num)

This function computes the probability density p(x) at $x for a chi-squared distribution with $nu degrees of freedom.

=head3 chisq-P(Num() $x where * ≥ 0, Num() $nu --> Num)
=head3 chisq-Q(Num() $x where * ≥ 0, Num() $nu --> Num)
=head3 chisq-Pinv(Num() $P, Num() $nu --> Num)
=head3 chisq-Qinv(Num() $Q, Num() $nu --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the chi-squared distribution with $nu degrees of freedom.

=head3 fdist(Math::Libgsl::Random $r, Num() $nu1, Num() $nu2 --> Num)

This function returns a random variate from the F-distribution with degrees of freedom $nu1 and $nu2.

=head3 fdist-pdf(Num() $x where * ≥ 0, Num() $nu1, Num() $nu2 --> Num)

This function computes the probability density p(x) at $x for an F-distribution with $nu1 and $nu2 degrees of freedom.

=head3 fdist-P(Num() $x where * ≥ 0, Num() $nu1, Num() $nu2 --> Num)
=head3 fdist-Q(Num() $x where * ≥ 0, Num() $nu1, Num() $nu2 --> Num)
=head3 fdist-Pinv(Num() $P, Num() $nu1, Num() $nu2 --> Num)
=head3 fdist-Qinv(Num() $Q, Num() $nu1, Num() $nu2 --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the F-distribution with $nu1 and $nu2 degrees of freedom.

=head3 tdist(Math::Libgsl::Random $r, Num() $nu --> Num)

This function returns a random variate from the t-distribution.

=head3 tdist-pdf(Num() $x, Num() $nu --> Num)

This function computes the probability density p(x) at $x for a t-distribution with $nu degrees of freedom.

=head3 tdist-P(Num() $x, Num() $nu --> Num)
=head3 tdist-Q(Num() $x, Num() $nu --> Num)
=head3 tdist-Pinv(Num() $P, Num() $nu --> Num)
=head3 tdist-Qinv(Num() $Q, Num() $nu --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the t-distribution with $nu degrees of freedom.

=head3 beta(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the beta distribution.

=head3 beta-pdf(Num() $x, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a beta distribution with parameters $a and $b.

=head3 beta-P(Num() $x, Num() $a, Num() $b --> Num)
=head3 beta-Q(Num() $x, Num() $a, Num() $b --> Num)
=head3 beta-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 beta-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the beta distribution with parameters $a and $b.

=head3 logistic(Math::Libgsl::Random $r, Num() $a --> Num)

This function returns a random variate from the logistic distribution.

=head3 logistic-pdf(Num() $x, Num() $a --> Num)

This function computes the probability density p(x) at $x for a logistic distribution with scale parameter $a.

=head3 logistic-P(Num() $x, Num() $a --> Num)
=head3 logistic-Q(Num() $x, Num() $a --> Num)
=head3 logistic-Pinv(Num() $P, Num() $a --> Num)
=head3 logistic-Qinv(Num() $Q, Num() $a --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the logistic distribution with scale parameter $a.

=head3 pareto(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the Pareto distribution of order $a and scale $b.

=head3 pareto-pdf(Num() $x, Num() $a, Num() $b where { $x ≥ $b } --> Num)

This function computes the probability density p(x) at $x for a Pareto distribution with exponent $a and scale $b.

=head3 pareto-P(Num() $x, Num() $a, Num() $b where { $x ≥ $b } --> Num)
=head3 pareto-Q(Num() $x, Num() $a, Num() $b where { $x ≥ $b } --> Num)
=head3 pareto-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 pareto-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Pareto distribution with exponent $a and scale $b.

=head3 d2dir(Math::Libgsl::Random $r --> List)
=head3 d2dir-trig-method(Math::Libgsl::Random $r --> List)

This function returns a random direction vector v = (x, y) in two dimensions.
The return value is a List of two Num(s): the x and y components of the 2D vector.

=head3 d3dir(Math::Libgsl::Random $r --> List)

This function returns a random direction vector v = (x, y, z) in three dimensions.
The return value is a List of three Num(s): the x, y, and z components of the 3D vector.

=head3 dndir(Math::Libgsl::Random $r, Int $n --> List)

This function returns a random direction vector v = (x 1 , x 2 , . . . , x n ) in $n dimensions.
The return value is a List of $n Num(s): the components of the n-D vector.

=head3 weibull(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the Weibull distribution.

=head3 weibull-pdf(Num() $x where $x ≥ 0, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a Weibull distribution with scale $a and exponent $b.

=head3 weibull-P(Num() $x where $x ≥ 0, Num() $a, Num() $b --> Num)
=head3 weibull-Q(Num() $x where $x ≥ 0, Num() $a, Num() $b --> Num)
=head3 weibull-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 weibull-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Weibull distribution with scale $a and exponent $b.

=head3 gumbel1(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the Type-1 Gumbel distribution.

=head3 gumbel1-pdf(Num() $x, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a Type-1 Gumbel distribution with parameters $a and $b.

=head3 gumbel1-P(Num() $x, Num() $a, Num() $b --> Num)
=head3 gumbel1-Q(Num() $x, Num() $a, Num() $b --> Num)
=head3 gumbel1-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 gumbel1-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-1 Gumbel distribution with parameters $a and $b.

=head3 gumbel2(Math::Libgsl::Random $r, Num() $a, Num() $b --> Num)

This function returns a random variate from the Type-2 Gumbel distribution.

=head3 gumbel2-pdf(Num() $x where * > 0, Num() $a, Num() $b --> Num)

This function computes the probability density p(x) at $x for a Type-2 Gumbel distribution with parameters $a and $b.

=head3 gumbel2-P(Num() $x where * > 0, Num() $a, Num() $b --> Num)
=head3 gumbel2-Q(Num() $x where * > 0, Num() $a, Num() $b --> Num)
=head3 gumbel2-Pinv(Num() $P, Num() $a, Num() $b --> Num)
=head3 gumbel2-Qinv(Num() $Q, Num() $a, Num() $b --> Num)

These functions compute the cumulative distribution functions P(x), Q(x) and their inverses for the Type-2 Gumbel distribution with parameters $a and $b.

=head3 dirichlet(Math::Libgsl::Random $r, Int $K, @alpha where *.all > 0 --> List)

This function returns an array of K random variates as a List of Num(s) from a Dirichlet distribution of order K-1.

=head3 dirichlet-pdf(Int $K, @alpha where *.all > 0, @theta where *.all ≥ 0 --> Num)

This function computes the probability density p(θ₁, …, θₖ) at theta[K] for a Dirichlet distribution with parameters alpha[K].

=head3 dirichlet-lnpdf(Int $K, @alpha where *.all > 0, @theta where *.all ≥ 0 --> Num)

This function computes the logarithm of the probability density p(θ₁, …, θₖ) at theta[K] for a Dirichlet distribution with parameters alpha[K].

=head3 poisson(Math::Libgsl::Random $r, Num() $mu --> UInt)

This function returns a random integer from the Poisson distribution with mean $mu.

=head3 poisson-pdf(UInt $k, Num() $mu --> Num)

This function computes the probability p(k) of obtaining $k from a Poisson distribution with mean $mu.

=head3 poisson-P(UInt $k, Num() $mu --> Num)
=head3 poisson-Q(UInt $k, Num() $mu --> Num)

These functions compute the cumulative distribution functions P(k), Q(k) for the Poisson distribution with parameter $mu.

=head3 bernoulli(Math::Libgsl::Random $r, Num() $p --> Int)

This function returns either 0 or 1, the result of a Bernoulli trial with probability $p.

=head3 bernoulli-pdf(UInt $k, Num() $p --> Num)

This function computes the probability p(k) of obtaining $k from a Bernoulli distribution with probability parameter $p.

=head3 binomial(Math::Libgsl::Random $r, Num() $p, UInt $n --> Int)

This function returns a random integer from the binomial distribution, the number of successes in $n independent trials with probability $p.

=head3 binomial-pdf(UInt $k, Num() $p, UInt $n where * ≥ $k --> Num)

This function computes the probability p(k) of obtaining $k from a binomial distribution with parameters $p and $n.

=head3 binomial-P(UInt $k, Num() $p, UInt $n where * ≥ $k --> Num)
=head3 binomial-Q(UInt $k, Num() $p, UInt $n where * ≥ $k --> Num)

These functions compute the cumulative distribution functions P (k), Q(k) for the binomial distribution with parameters $p and $n.

=head3 multinomial(Math::Libgsl::Random $r, Int $K, UInt $N, *@p where { @p.all ~~ Numeric } --> List)

This function returns a random sample from the multinomial distribution formed by $N trials from an underlying distribution p[K] as a List of UInt(s).

=head3 multinomial-pdf(Int $K, @p where { @p.all ~~ Numeric }, *@n where { @n.all ~~ UInt } --> Num)

This function computes the probability P(n₁, n₂, …, nₖ) of sampling n[K] from a multinomial distribution with parameters p[K].

=head3 multinomial-lnpdf(Int $K, @p where { @p.all ~~ Numeric }, *@n where { @n.all ~~ UInt } --> Num)

This function computes the logarithm of the probability P(n₁, n₂, …, n�