### ----------------------------------------------------
### -- Statistics
### -- Licenses: Artistic-2.0
### -- Authors: Suman Khanal <suman81765@gmail.com>
### -- File: t/04-prob_distributions.rakutest
### ----------------------------------------------------
### -- Chunk 1 of 2
#======================================================#
# Author: Dr Suman Khanal
# Date creation: Tuesday Oct 10 2023, 00:39:03
# ~: Tests for various prob distributions
# File: 04-prob_distributions.rakutest
#======================================================#

use Test;
use Statistics;
# use Probability::Distribution;

# Local function for testing two arrays with desired precision
sub element_wise_testing_with_precision(@functionoutput, @expectedoutput, $precision) {
    die "Provide equal elements!" if +@functionoutput != +@expectedoutput;
    loop (my $i = 0; $i < +@functionoutput; $i++)  {
        is-approx @functionoutput[$i], @expectedoutput[$i], $precision;
    }
}


####	Normal distribution    ####

diag("####    Testing raku_dnorm    ####");
is-approx raku_dnorm(5), 1.4867195147342977e-06, 1e-12;
is-approx raku_dnorm(3, mean => 2, sd => 1), 0.24197072451914337, 1e-12;
is raku_dnorm(3, mean => 2, sd => -1), NaN;
is raku_dnorm(3, mean => 2, sd => 0), 0;

diag("####    Testing raku_pnorm    ####");
is-approx raku_pnorm(5), 0.9999997133484281, 1e-12;
is-approx raku_pnorm(4, mean => 2, sd => 3), 0.7475074624530771, 1e-12;
is-approx raku_pnorm(4, mean => 2, sd => 3, log_p => True), 
                    -0.29101099055230867, 1e-12;

diag("####    Testing raku_qnorm    ####");
is raku_qnorm(5), NaN;
is-approx raku_qnorm(0.1), -1.2815515655446008, 1e-12;
is raku_qnorm(1), Inf;
is raku_qnorm(0.5), 0;
is-approx raku_qnorm(0.9, mean => 3, sd =>2), 5.563103131089202, 1e-12;
is-approx raku_qnorm(0.7),0.5244005127080407, 1e-12;


diag("####    Testing raku_rnorm    ####");
element_wise_testing_with_precision sub ($n) {set_seed(111, 222); raku_rnorm($n)}(12), 
                                    (0.2933064931910023, -0.43157564134676835, 0.1280404124560668, -1.001461811038476,
                                    -0.9219453227924342, 1.2988990178409578, -0.46867271131577315, -1.2678609176619775,
                                    0.3596981905325252, 0.35262105537769173, 0.5770092689090144, -1.0392300758070165),
                                    1e-12;

####	Gamma distribution    ####


diag("####    Testing raku_dgammma    ####");
is raku_dgamma(2,4, rate => 1), 0.1804470443154836;
is raku_dgamma(2,4, scale => 2), 0.030656620097620196;
is raku_dgamma(2,4, rate => 7), 0.0026620006060902226;
is raku_dgamma(2,4, scale => 1/7), 0.0026620006060902226;
is raku_dgamma(2,4, scale => 1/7), 0.0026620006060902226;
is raku_dgamma(2,4, scale => 2, log => True), -3.4849066497880004;

diag("####    Testing raku_pgamma    ####");
is raku_pgamma(2,4), 0.142876539501453;
is raku_pgamma(2,4, log_p => True), -1.9457743817619386;
is raku_pgamma(2,4, scale => 4, log_p => True), -6.34721274558465;
is raku_pgamma(2,4, rate => 4, log_p => True), -0.04330435635721676;
is raku_pgamma(2,4, scale => 1/4, log_p => True), -0.04330435635721676;

diag("####    Testing raku_qgamma    ####");
is raku_qgamma(2,4), NaN;
is raku_qgamma(1,0.5,log_p => True), NaN;
is raku_qgamma(1,2), Inf;


####	Beta distribution 	####

diag("####    Testing raku_dbeta    ####");
is raku_dbeta(2,4, 2), 0;
is raku_dbeta(2,4, 2, log => True), -Inf;
is raku_dbeta(1, 1, 0, ncp => 2,log => True), NaN;



diag("####    Testing raku_pbeta    ####");
is raku_pbeta(1,1,1), 1;
is raku_pbeta(0.5,1,1), 0.5;
is raku_pbeta(0.5, 2, 4, ncp => 2), 0.6590492649101354;
is raku_pbeta(0.5, 2, 4, ncp => 2, log_p => True), -0.41695699020227484;


diag("####    Testing raku_qbeta    ####");
is raku_qbeta(3, 4, 6), NaN;


####	Lognormal distribution 	####

diag("####    Testing raku_dlnorm    ####");
is-approx raku_dlnorm(3), 0.07272825613999472, 1e-12;
is-approx raku_dlnorm(3, meanlog => 3, log => True), -3.825188436274744, 1e-12;

diag("####    Testing raku_plnorm    ####");
is-approx raku_plnorm(3),0.8640313923585756, 1e-12;
is-approx raku_plnorm(5, meanlog => 2, sdlog => 2, log_p => True), -0.8613611948322474, 1e-12;

diag("####    Testing raku_qlnorm    ####");
is raku_qlnorm(2), NaN;
is raku_qlnorm(0.5), 1;
is raku_qlnorm(0), 0;
is-approx raku_qlnorm(0.1), 0.2776062418520097, 1e-12;


diag("####    Testing raku_rlnorm    ####");
element_wise_testing_with_precision sub ($n) {set_seed(1,3);raku_rlnorm($n)}(3), 
                                    (1.1751666576815853, 1.551169741442785, 0.35333870325756733),
                                    1e-12;
element_wise_testing_with_precision sub ($n) {set_seed(2,5);raku_rlnorm($n)}(3), 
                                    (0.32145646272210726, 0.6608107230044763, 0.588862748840785),
                                    1e-12;

####	Chi-Squared Distribution	####

diag("####    Testing raku_dchisq    ####");
is raku_dchisq(4,3), 0.10798193302637614;    
is raku_dchisq(4,3,ncp => 4), 0.09970211254391682;
is raku_dchisq(4,3,ncp => 4, log => True ), -2.3055684132326415;

diag("####    Testing raku_pchisq    ####");
is-approx raku_pchisq(3,4), 0.44217459962892536, 1e-12;
is raku_pchisq(0,0), 0;
is raku_pchisq(0,0, ncp => 4, lower_tail => False, log_p => True), -0.14541345786885906;

diag("####    Testing raku_qchisq    ####");
is raku_qchisq(5,4), NaN;
is raku_qchisq(4,4), NaN;
is raku_qchisq(1,3), Inf;
is raku_qchisq(0,4, ncp => 3), 0;
is raku_qchisq(0,0, ncp => 4, lower_tail => False, log_p => True), 0;
is raku_qchisq(75,25, ncp => 4, lower_tail => False, log_p => True), NaN;

####	F distribution 	####

diag("####    Testing raku_df    ####");
is-approx raku_df(1,3,2), 0.27885480092693404, 1e-12;
is-approx raku_df(1,3,2,4), 0.2255355699113039, 1e-12;

diag("####    Testing raku_pf    ####");
is-approx raku_pf(1,3,2), 0.46475800154489, 1e-12;
is-approx raku_pf(1,3,2, log_p => True), -0.7662384356489861, 1e-12;

diag("####    Testing raku_qf    ####");
is raku_qf(4, 2, 1), NaN;
is raku_qf(1, 2, 1), Inf;
is raku_qf(0, 2, 1), 0;
is raku_qf(0, 2, 1, log_p => True), Inf;

diag("####    Testing raku_rf    ####");
element_wise_testing_with_precision sub ($n) {set_seed(2,5); raku_rf($n, 1, 2)}(3), 
                                    (0.18435620618029944, 1.04131942262752, 14.825797754677964),
                                    1e-12;


####	Student t distribution	 ####

diag("####    Testing raku_dt    ####");
is-approx raku_dt(1,2), 0.19245008972987526, 1e-12;
is-approx raku_dt(1,2,4), 0.011764946752329405, 1e-12;

diag("####    Testing raku_pt    ####");
is-approx raku_pt(1,2), 0.7886751345948129, 1e-12;
is raku_pt(5,4), 0.9962547830593627;
is raku_pt(5,4, lower_tail => False, log_p => True), -5.587275735616685;

diag("####    Testing raku_qt    ####");
is raku_qt(1,2), Inf;
is-approx raku_qt(0.2,2), -1.060660171779821, 1e-12;

diag("####    Testing raku_rt    ####");
element_wise_testing_with_precision sub ($n) {set_seed(3,10); raku_rt($n, 3, 4)}(5), 
                                    (3.2166139332581323, 2.9985780896164873, 3.4756218027009744, 5.833058113048503, 5.671377806342658),
                                    1e-12;

####	Binomial distribution 	####

diag("####    Testing raku_dbinom    ####");
is-approx raku_dbinom(1,2,0.5), 0.49999999999999994, 1e-12;
is raku_dbinom(1,2,1), 0;
is raku_dbinom(1,2,4), NaN;
is-approx raku_dbinom(1,2,0.2, log => True), -1.1394342831883648, 1e-12;

diag("####    Testing raku_pbinom    ####");
is raku_pbinom(1,2,0.2), 0.96;
is-approx raku_pbinom(1,2,0.2, log_p => True), 
                     -0.04082199452025514, 1e-12;
is-approx raku_pbinom(1,2,0.2, lower_tail => False, log_p => True), 
                     -3.2188758248682006, 1e-12;


diag("####    Testing raku_qbinom    ####");
is raku_qbinom(1,2,0.2), 2;
is raku_qbinom(1,2,3), NaN;
is raku_qbinom(1,2,3, log_p => True), NaN;

diag("####    Testing raku_rbinom    ####");
is sub ($n) {set_seed(1,9); raku_rbinom($n,2,0.2)}(3), (0, 0, 1);
is sub ($n) {set_seed(20,16); raku_rbinom($n,2,0.2)}(4), (0, 1, 0, 0);


####	Multinomial distribution 	####
# Test taken from here: 
#      https://stackoverflow.com/questions/63368639/how-to-use-multinomial-function-in-c-with-rmath-h

diag("####    Testing raku_rmultinom    ####");
is sub ($n) {set_seed(123, 456); my @prob = [0.1, 0.2, 0.4, 0.3]; raku_rmultinom(100, 4 , @prob)}(4), (9,14,43,34);
is -> $n {set_seed(123, 456); raku_rmultinom($n, 5, [0.1,0.2,0.3,0.2,0.2])}(0), (0, 0, 0, 0, 0);
is -> $n {set_seed(123, 456); raku_rmultinom($n, 5, [0.1,0.2,0.3,0.2,0.2])}(1), (0, 0, 0, 0, 1);
is -> $n {set_seed(123, 456); raku_rmultinom($n, 5, [0.1,0.2,0.3,0.2,0.2])}(2), (0, 0, 1, 0, 1);
is -> $n {set_seed(123, 456); raku_rmultinom($n, 5, [0.1,0.2,0.3,0.2,0.2])}(3), (0, 0, 1, 0, 2);
is -> $n {set_seed(123, 456); raku_rmultinom(100, $n, [0.1, 0.3, 0.5, 0.1])}(4), (9, 29, 49, 13);

####	Cauchy distribution 	####

diag("####    Testing raku_dcauchy    ####");
is-approx raku_dcauchy(1), 0.15915494309189535, 1e-12;
is-approx raku_dcauchy(1,2), 0.15915494309189535, 1e-12;
is-approx raku_dcauchy(1,2,3), 0.0954929658551372, 1e-12;

diag("####    Testing raku_pcauchy    ####");
is raku_pcauchy(1), 0.75;
is raku_pcauchy(1,2), 0.25;
is-approx raku_pcauchy(1,2,3), 0.39758361765043326, 1e-12;

diag("####    Testing raku_qcauchy    ####");
is raku_qcauchy(1), Inf;
is raku_qcauchy(1,2), Inf;
is raku_qcauchy(1,2,3), Inf;
is-approx raku_qcauchy(0.1), -3.077683537175254  , 1e-12;
is-approx raku_qcauchy(0.9),  3.0776835371752544 , 1e-12;
is raku_qcauchy(0.9, log_p => True),  NaN;
is raku_qcauchy(0), -Inf;

diag("####    Testing raku_rcauchy    ####");
element_wise_testing_with_precision sub ($n) {set_seed(2,5); raku_rcauchy($n)}(5), 
                                    (0.42607842113327926, -3.589933350285371, 1.8099808111776452, -0.5705944232265042, 1.3602101594780192),
                                    1e-12;

####	Exponential distribution 	####

diag("####    Testing raku_dexp    ####");
is-approx raku_dexp(1), 0.36787944117144233 , 1e-12;
is-approx raku_dexp(1,2), 0.2706705664732254 , 1e-12;
is-approx raku_dexp(1,2, log => True), -1.3068528194400546 , 1e-12;

diag("####    Testing raku_pexp    ####");
is-approx raku_pexp(1), 0.6321205588285577 , 1e-12;
is-approx raku_pexp(1,2), 0.8646647167633873 , 1e-12;
is-approx raku_pexp(1,2, log_p => True), -0.14541345786885906 , 1e-12;
is-approx raku_pexp(4,2), 0.9996645373720975, 1e-12;

diag("####    Testing raku_qexp    ####");
is raku_qexp(1), Inf;
is raku_qexp(1,2), Inf;
is raku_qexp(1,2, log_p => True), NaN;
is raku_qexp(4,2), NaN;
is-approx raku_qexp(0.5), 0.6931471805599453, 1e-12;
is raku_qexp(0.5, log_p => True), NaN;


diag("####    Testing raku_rexp    ####");
element_wise_testing_with_precision sub ($n) {set_seed(2,5); raku_rexp($n)}(2), 
                                    (1.411974694734631, 0.17295038773979732),
                                    1e-12;

element_wise_testing_with_precision sub ($n) {set_seed(2,5); raku_rexp($n, 4)}(2), 
                                    (0.35299367368365775, 0.04323759693494933),
                                    1e-12;

####	Geometric distribution    ####

diag("####    Testing raku_dgeom    ####");
is raku_dgeom(1,0.5), 0.25 ;
is raku_dgeom(1, 1), 0 ;
is raku_dgeom(1, 1, log => True), -Inf;

diag("####    Testing raku_pgeom    ####");
is raku_pgeom(1,0.5), 0.75 ;
is raku_pgeom(1, 1), 1 ;
is raku_pgeom(1, 1, log_p => True), 0;
is raku_pgeom(1,2), NaN;

diag("####    Testing raku_qgeom    ####");
is raku_qgeom(1,0.5), Inf ;
is raku_qgeom(1, 1), 0 ;
is raku_qgeom(1, 1, log_p => True), NaN;
is raku_qgeom(1,2), NaN;
is raku_qgeom(1,0), NaN;

diag("####    Testing raku_rgeom    ####");
is sub ($n) {set_seed(2,5); raku_rgeom($n,0.1)}(2), (13, 5);
is sub ($n) {set_seed(2,5); raku_rgeom($n,0.6)}(4), (1, 1, 1, 0);


####	Hypergeometric distribution 	####

diag("####    Testing raku_dhyper    ####");
is-approx raku_dhyper(1, 2, 3, 4), 0.4, 1e-12;  

diag("####    Testing raku_phyper    ####");
is-approx raku_phyper(1, 2, 3, 4), 0.4, 1e-12;  
is-approx raku_phyper(1, 2, 3, 4, lower_tail => False), 0.6, 1e-12;  
is raku_phyper(1, 2, 3, 6), NaN;  
is raku_phyper(1, 2, 3, 5, lower_tail => False), 1;  
is raku_phyper(11, 15, 0, 12, lower_tail => False), 1; 
is raku_phyper(11, 15, 0, 12, lower_tail => True), 0; 
is raku_phyper(11, 15, 0, 12, log_p => True), -Inf;  


diag("####    Testing raku_qhyper    ####");
is raku_qhyper(1, 2, 3, 4), 2;  
is raku_qhyper(1, 2, 3, 4, lower_tail => False), 1;  
is raku_qhyper(1, 2, 3, 6), NaN;  
is raku_qhyper(1, 2, 3, 5, lower_tail => False), 2;  
is raku_qhyper(0, 2, 3, 5, lower_tail => False), 2;  
is raku_qhyper(0.5, 2, 3, 5, lower_tail => False), 2;
is raku_qhyper(0, 20, 3, 6, lower_tail => False), 6; 
is raku_qhyper(0.7, 2, 3, 6, lower_tail => False), NaN;  

diag("####    Testing raku_rhyper    ####");
is sub ($n) {set_seed(1,9); raku_rhyper($n,2,3,4)}(1), (2);
is sub ($n) {set_seed(12,7); raku_rhyper($n,2,3,4)}(7), (1, 2, 2, 2, 1, 2, 2);



####	Negative binomial distribution  	####

diag("####    Testing raku_dnbinom    ####");
is-approx raku_dnbinom(2, 3, prob => 0.2), 0.030719999999999997, 1e-12;
is-approx raku_dnbinom(2, 3, mu => 1), 0.1582031249999999, 1e-12;
is-approx raku_dnbinom(5, 10, mu => 1), 0.004792628828370662, 1e-12;
is raku_dnbinom(5, 10, prob => 1), 0;


diag("####    Testing raku_pnbinom    ####");
is-approx raku_pnbinom(2, 3, prob => 0.2), 0.05792, 1e-12;
is-approx raku_pnbinom(2, 3, mu => 1), 0.8964843749999999, 1e-12;
is-approx raku_pnbinom(5, 10, mu => 1), 0.9986311157049472, 1e-12;
is raku_pnbinom(5, 10, prob => 1), 1;

diag("####    Testing raku_qnbinom    ####");
is raku_qnbinom(2, 3, prob => 0.2), NaN;
is raku_qnbinom(2, 3, mu => 1), NaN;
is raku_qnbinom(5, 10, mu => 1), NaN;
is raku_qnbinom(5, 10, prob => 1), 0;
is raku_qnbinom(1, 10, prob => 0.5), Inf;

diag("####    Testing raku_rnbinom    ####");

# Run these commands you will get same random values in C and R
# RNGkind("Marsaglia-Multicarry")
# .Random.seed <- c(401L, 111L, 222L)
is -> $n {set_seed(111,222); raku_rnbinom($n, 2, mu => 3)}(5), (2,8,1,4,3);
is -> $n {set_seed(111,222); raku_rnbinom($n, 2, prob => 0.2)}(10), (6,17,8,12,0,9,18,20,9,6);


####	Poisson distribution 	####

diag("####    Testing raku_dpois    ####");
is-approx raku_dpois(1, 2), 0.2706705664732254, 1e-12;
is-approx raku_dpois(1, 0.2), 0.16374615061559641, 1e-12;
is-approx raku_dpois(1, 20, log => True), -17.00426772644601, 1e-12;
is-approx raku_dpois(5, 2), 0.03608940886309672, 1e-12;
is-approx raku_dpois(10, 2), 3.8189850648779595e-05, 1e-12;

diag("####    Testing raku_ppois    ####");
is-approx raku_ppois(1, 2), 0.4060058497098381, 1e-12;
is-approx raku_ppois(1, 0.2), 0.9824769036935782, 1e-12;
is-approx raku_ppois(1, 20, log_p => True), -16.955477562276577, 1e-12;
is-approx raku_ppois(5, 2), 0.9834363915193856, 1e-12;
is-approx raku_ppois(10, 2), 0.9999916917756315, 1e-12;


diag("####    Testing raku_qpois    ####");
is raku_qpois(1, 2), Inf;
is raku_qpois(1, 0.2), Inf;
is raku_qpois(0.5, 20, log_p => True), NaN;
is raku_qpois(0.1, 2), 0;
is raku_qpois(10, 2), NaN;
is raku_qpois(0.2, 2), 1;

diag("####    Testing raku_rpois    ####");
is -> $n {set_seed(111,222); raku_rpois($n, 2)}(5), (2, 2, 1, 3, 2);
is -> $n {set_seed(111,222); raku_rpois($n, 0.2)}(10), (0, 0, 0, 1, 0, 1, 0, 0, 0, 0);

####	Weibull distribution 	####

diag("####    Testing raku_dweibull    ####");
is-approx raku_dweibull(3,1,2), 0.11156508007421491, 1e-12;
is-approx raku_dweibull(3,1,2,log => True), -2.1931471805599454, 1e-12;
is raku_dweibull(3,7), 0;
is-approx raku_dweibull(30,1,0.2), 3.5875479865822053e-65, 1e-12;

diag("####    Testing raku_pweibull    ####");
is-approx raku_pweibull(3, 1, 2), 0.7768698398515702, 1e-12;
is-approx raku_pweibull(3, 1, 2,log_p => True), -0.25248245892545396, 1e-12;
is raku_pweibull(3, 7), 1;
is raku_pweibull(30, 1, 0.2), 1;



diag("####    Testing raku_qweibull    ####");
is raku_qweibull(1, 1, 2), Inf;
is raku_qweibull(