#!/usr/bin/perl

use constant pi => 4 * atan2(1, 1);
use constant pi_inverse => 0.25 / atan2(1, 1);
use constant two_pi_sqrt_inverse => 1 / sqrt(8 * atan2(1, 1));
use constant sqrt_twopi => sqrt(8 * atan2(1, 1));


# bernoulli($x, $p) returns 1-$p if $x is 0, $p if $x is 1, 0 otherwise.
#
sub bernoulli {
    my ($x, $p) = @_;
    return unless $p > 0 && $p < 1;
    return $x ? ( ($x == 1) ? $p : 0 ) : (1 - $p);
}
sub bernoulli_expected { $_[0] }
sub bernoulli_variance { $_[0] * (1 - $_[0]) }


# beta( $x, $a, $b ) returns the Beta distribution for $x given the
# Beta parameters $a and $b.
#
sub beta {
    my ($x, $a, $b) = @_;
    return unless $a > 0 and $b > 0;
    return factorial ($a + $b - 1) / factorial ($a - 1) /
        factorial ($b - 1) * ($x ** ($a - 1)) * ((1 - $x) ** ($b - 1));
}

sub beta_expected { $_[0] / ($_[0] + $_[1]) }
sub beta_variance { ($_[0] * $_[1]) / (($_[0] + $_[1]) ** 2) /
                        ($_[0] + $_[1] + 1) }

# binomial($x, $n, $p);
# binomial_expected($n, $p);
#
sub binomial {
    my ($x, $n, $p) = @_;
    return unless $x >= 0 && $x == int $x && $n > 0 &&
        $n == int $n && $p > 0 && $p < 1;
    return factorial($n) / factorial($x) / factorial($n - $x) *
        ($p ** $x) * ((1 - $p) ** ($n - $x));
}

sub binomial_expected { $_[0] * $_[1] }
sub binomial_variance { $_[0] * $_[1] * (1 - $_[1]) }

sub cauchy {
    my ($x, $a, $b) = @_;
    return unless $a > 0;
    return pi_inverse * $a / (($a ** 2) + (($x - $b) ** 2));
}
sub cauchy_expected { $_[1] }

sub chi_square {
    my ($x, $n) = @_;
    return 0 unless $x > 0;
    return 1 / factorial($n/2 - 1) * (2 ** (-$n / 2)) *
        ($x ** (($n / 2) - 1)) * exp(-$x / 2);
}
sub chi_square_expected { $_[0] }
sub chi_square_variance { 2 * $_[0] }


sub erlang {
    my ($x, $a, $n) = @_;
    return unless $a > 0 && $n > 0 && $n == int($n);
    return 0 unless $x > 0;
    return ($a ** $n) * ($x ** ($n-1)) * exp(-$a * $x) / factorial($n-1);
}

sub erlang_expected { $_[1] / $_[0] }
sub erlang_variance { $_[1] / ($_[0] ** 2) }

sub exponential {
    my ($x, $a) = @_;
    return unless $a > 0;
    return 0 unless $x > 0;
    return $a * exp(-$a * $x);
}
sub exponential_expected { 1 / $_[0] }
sub exponential_variance { 1 / ($_[0] ** 2) }

sub gamma {
    my ($x, $a, $b) = @_;
    return unless $a > -1 && $b > 0;
    return 0 unless $x > 0;
    return ($x ** $a) * exp(-$x / $b) / factorial($a) / ($b ** ($a + 1));
}

sub gamma_expected { ($_[0] + 1) * $_[1] }
sub gamma_variance { ($_[0] + 1) * ($_[1] ** 2) }


sub gaussian {
    my ($x, $mean, $variance) = @_;
    return two_pi_sqrt_inverse *

        exp( -( ($x - $mean) ** 2 ) / (2 * $variance) ) /
            sqrt $variance;
}

sub geometric {
    my ($x, $p) = @_;
    return unless $p > 0 && $p < 1;
    return 0 unless $x == int($x);
    return $p * ((1 - $p) ** ($x - 1));
}

sub geometric_expected { 1 / $_[0] }
sub geometric_variance { (1 - $_[0]) / ($_[0] ** 2) }

sub hypergeometric {
    my ($x, $k, $m, $n) = @_;
    return unless $m > 0 && $m == int($m) && $n > 0 && $n == int($n) &&
        $k > 0 && $k <= $m + $n;
    return 0 unless $x <= $k && $x == int($x);
    return choose($m, $x) * choose($n, $k - $x) / choose($m + $n, $k);
}

sub hypergeometric_expected { $_[0] * $_[1] / ($_[1] + $_[2]) }
sub hypergeometric_variance {
    my ($k, $m, $n) = @_;
    return $m * $n * $k * ($m + $n - $k) / (($m + $n) ** 2) /
        ($m + $n - 1);
}


# laplace($x, $a, $b)
sub laplace {
    return unless $_[1] > 0;
    return $_[1] / 2 * exp( -$_[1] * abs($_[0] - $_[2]) );
}

sub laplace_expected { $_[1] }
sub laplace_variance { 2 / ($_[0] ** 2) }

sub lognormal {
    my ($x, $a, $b, $std) = @_;
    return unless $std > 0;
    return 0 unless $x > $a;
    return (exp -(((log($x - $a) - $b) ** 2) / (2 * ($std ** 2)))) /
        (sqrt_twopi * $std * ($x - $a));
}

sub lognormal_expected { $_[0] + exp($_[1] + 0.5 * ($_[2] ** 2)) }
sub lognormal_variance { exp(2 * $_[1] + ($_[2] ** 2)) * (exp($_[2] ** 2)
                         - 1) }

sub maxwell {
    my ($x, $a) = @_;
    return unless $a > 0;
    return 0 unless $x > 0;
    return sqrt(2 / pi) * ($a ** 3) * ($x ** 2) *
        exp($a * $a * $x * $x / -2);
}

sub maxwell_expected { sqrt( 8/pi ) / $_[0] }
sub maxwell_variance { (3 - 8/pi) / ($_[0] ** 2) }

sub pascal {
    my ($x, $n, $p) = @_;
    return unless $p > 0 && $p < 1 && $n > 0 && $n == int($n);
    return 0 unless $x >= $n && $x == int($x);
    return choose($x - 1, $n - 1) * ($p ** $n) * ((1 - $p) ** ($x - $n));
}

sub pascal_expected { $_[0] / $_[1] }
sub pascal_variance { $_[0] * (1 - $_[1]) / ($_[1] ** 2) }

sub poisson {
    my ($x, $a) = @_;
    return unless $a >= 0 && $x >= 0 && $x == int($x);
    return ($a ** $x) * exp(-$a) / factorial($x);
}

sub poisson_expected { $_[0] }
sub poisson_variance { $_[0] }

sub rayleigh {
    my ($x, $a) = @_;
    return unless $a > 0;
    return 0 unless $x > 0;
    return ($a ** 2) * $x * exp( -($a ** 2) * ($x ** 2) / 2 );
}

sub rayleigh_expected { sqrt(pi / 2) / $_[0] }
sub rayleigh_variance { (2 - pi / 2) / ($_[0] ** 2) }

sub uniform {
    my ($x, $a, $b) = @_;
    return unless $b > $a;
    return 0 unless $x > $a && $x < $b;
    return 1 / ($b - $a);
}

sub uniform_expected { ($_[0] + $_[1]) / 2 }
sub uniform_variance { (($_[1] - $_[0]) ** 2) / 12 }