### ----------------------------------------------------
### -- Statistics::Distributions
### -- Licenses: Artistic-2.0
### -- Authors: Anton Antonov
### -- File: t/02-random-variate.rakutest
### ----------------------------------------------------
use Test;

#use lib <. lib>;

use Statistics::Distributions;

sub mean(@x) {
    [+](@x) / @x.elems
}
sub sd(@x) {
    sqrt([+]((@x X- mean(@x)) X** 2) / @x.elems)
}

plan *;

## 1
is random-variate(NormalDistribution.new(:mean(10), :sd(20)), 12).all ~~ Numeric,
        True,
        'Normal distribution';

## 2
is random-variate(NormalDistribution.new(:µ(10), :σ(20)), 12).all ~~ Numeric,
        True,
        'Normal distribution alternative spec';

## 3
is random-variate(NormalDistribution.new, 12).all ~~ Numeric,
        True,
        'Normal distribution default arguments';

## 4
my @res4 = random-variate(NormalDistribution.new, 1000);
is-approx abs(mean(@res4)), 0, 0.1,
        'Normal distribution mean with default arguments';

## 5
is-approx sd(@res4), 1.0, 0.3,
        'Normal distribution st. deviation with default arguments';

## 6
is random-variate(UniformDistribution.new(:min(10), :max(20)), 12).all ~~ Numeric,
        True,
        'Uniform distribution';

## 7
is random-variate(UniformDistribution.new, 12).all ~~ Numeric,
        True,
        'Uniform distribution default arguments';

## 8
my $res8 = random-variate(NormalDistribution.new, [4, 3]);

is ($res8 ~~ Array and $res8.shape eqv (4, 3)),
        True,
        'Normal distribution for 2D array';

## 9
isa-ok random-variate(UniformDistribution.new(:min(10), :max(20))), Numeric,
        'No size argument gives a number';

## 10
is random-variate(BernoulliDistribution.new, 12).all ~~ Numeric,
        True,
        'Bernoulli distribution default argument';

## 11
isa-ok random-variate(BernoulliDistribution.new(:p(0.1))), Numeric,
        'No size argument gives a number';

## 12
my @sample12 = random-variate(BernoulliDistribution.new(:p(0.4)), 10**4);
is-approx mean(@sample12), 0.4, 0.01;

## 13
is random-variate(UniformDistribution.new(3, 5), 100).all ~~ Numeric:D, True;

## 14
is random-variate(BetaDistribution.new(4, 4), 100).all ~~ Numeric:D, True;

## 15
is random-variate(MixtureDistribution.new([2, 5], [NormalDistribution.new(3, 4), NormalDistribution.new(6, 5)]), 100).all ~~ Numeric:D, True;

## 16
subtest {
    ok random-variate(ProductDistribution.new([NormalDistribution.new(3, 4), NormalDistribution.new(6, 5)]), 100);
    my @points = random-variate(ProductDistribution.new([NormalDistribution.new(3, 4), NormalDistribution.new(6, 5)]), 100);
    is @points.all ~~ Positional:D, True;

    is ([==] |@points».elems), True;
}

# 17
subtest {
    ok random-variate(DiscreteUniformDistribution.new(:min(10), :max(20)));
    isa-ok random-variate(DiscreteUniformDistribution.new(:min(10), :max(20))), Numeric:D,
            'No size argument gives a number for DiscreteUniformDistribution';
    is random-variate(DiscreteUniformDistribution.new(:min(-10), :max(20)), 10).all ~~ Int:D, True,
            'DiscreteUniformDistribution gives integers';
}

# 18
subtest {
    ok BinormalDistribution.new([10, 4], [4, 2], 0.5), 'three arguments';
    ok BinormalDistribution.new([1, 1], 0.5), 'two arguments';
    ok BinormalDistribution.new(0.5), 'one argument';

    is
            random-variate(BinormalDistribution.new([10, 4], [4, 2], 0.5), 10).all ~~ Positional:D,
            True,
            'list of positionals for three arugments';

    is
            random-variate(BinormalDistribution.new([4, 2], 0.5), 10).all ~~ Positional:D,
            True,
            'list of positionals for two arugments';
}

# 19
subtest {
    ok random-variate(ChiSquareDistribution.new(ν => 5));
    dies-ok { random-variate(ChiSquareDistribution.new(ν => -5)) };
    isa-ok random-variate(ChiSquareDistribution.new(3.3)), Numeric:D,
            'No size argument gives a number for ChiSquareDistribution';
    is random-variate(ChiSquareDistribution.new(20), 10).all ~~ Numeric:D, True,
            'ChiSquareDistribution gives numbers';
}

# 20
subtest {
    ok random-variate(StudentTDistribution.new(ν => 5));
    ok random-variate(StudentTDistribution.new(ν => 5, μ => 1, σ => 3));
    dies-ok { random-variate(StudentTDistribution.new(ν => -5)) };
    isa-ok random-variate(StudentTDistribution.new(3.3)), Numeric:D,
            'No size argument gives a number for StudentTDistribution';
    is random-variate(StudentTDistribution.new(20, 1, 3), 10).all ~~ Numeric:D, True,
            'StudentTDistribution gives numbers with three arguments';
}, 'StudentTDistribution';

# 21
subtest {
    ok random-variate(BenfordDistribution.new());
    ok random-variate(BenfordDistribution.new(b => 5));
    dies-ok { random-variate(BenfordDistribution.new(:2b)) };
    isa-ok random-variate(BenfordDistribution.new(:10b)),
            Int:D,
            'No size argument gives an integer for BenfordDistribution';

    is random-variate(BenfordDistribution.new(10), 12).all ~~ Int:D, True,
            'BenfordDistribution gives integers';

    # Expected distribution values check
    sub benford-law(UInt:D $d, UInt:D $b = 10) { log($d + 1, $b) - log($d, $b) };

    # In 10 differences runs at least 7 times we get expected distribution
    my @comp = do for ^10 {
        my %digitTally = random-variate(BenfordDistribution.new, 8_000).Bag;
        %digitTally = %digitTally.map({ $_.key => $_.value / %digitTally.values.sum });
        my @diffs = %digitTally.map({ abs($_.value - benford-law($_.key.Int)) });
        @diffs.map({ $_ < 0.01 }).sum == @diffs.elems
    }

    is @comp.sum ≥ 7, True, 'distribution check';

}, 'Benford disitribution';

done-testing;


