### ----------------------------------------------------
### -- Statistics
### -- Licenses: Artistic-2.0
### -- Authors: Suman Khanal <suman81765@gmail.com>
### -- File: lib/Probability/Distribution.rakumod
### ----------------------------------------------------
### -- Chunk 1 of 2
#======================================================#
# Author: Dr Suman Khanal
# Date creation: Monday Oct 09 2023, 20:39:32
# ~: Probability Distribution module
# File: Distribution.rakumod
#======================================================#


unit module Probability::Distribution;

use NativeCall;
constant RMATH = %?RESOURCES<libraries/Rmath>;


##========================================================##
# NOTE:
# cmake naively generates lib<libraryname>.dll in Windows
# and lib<libraryname>.so in Linux
# and lib<libraryname>.dylib in MacOS

##========================================================##
# Despite libRmath.so in Linux and libRmath.dylib in macos
# constant RMATH = %?RESOURCES<libraries/Rmath>
# libraries/Rmath in META6.json works
# in both OSes

##==========================================##
# BUT it does not in Windows, instead requires
# constant RMATH = %?RESOURCES<libraries/libRmath>
# libraries/libRmath in META6.json 
#===========================================##

##=====================================================##
# To avoid this and to make the module cross-platform
# there is 
# if (WIN32)
#   set_target_properties(Rmath PROPERTIES PREFIX "")
# endif()
# in CMakeLists.txt, so that in Windows it builds only 
# Rmath.dll that is compatible with
##======================================================##
# constant RMATH = %?RESOURCES<libraries/Rmath>
# libraries/Rmath in META6.json in Windows
##======================================================##

##======================================================##
# The module needs gcc, cmake and ninja for install. I find it better
# to avoid using **make** because it has so many different names 
# (make, mingw32-make, gmake, nmake), and might run into issues.
# Best would be to use **ninja** 
##======================================================##


####    Function signatures     ####


#	/* Normal Distribution */


sub dnorm4(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * };
sub pnorm5(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub qnorm5(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub rnorm(num64, num64) returns num64
    is native( RMATH ) { * };   


#	/* Uniform Distribution */

sub dunif(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * };    
sub punif(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };    
sub qunif(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };    
sub runif(num64, num64) returns num64
    is native( RMATH ) { * };    

#	/* Gamma Distribution */

sub dgamma(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * };    
sub pgamma(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };    
sub qgamma(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };    
sub rgamma(num64, num64) returns num64
    is native( RMATH ) { * };    


#	/* Beta Distribution */


sub dbeta(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pbeta(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qbeta(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rbeta(num64, num64) returns num64
    is native( RMATH ) { * }; 



#	/* Lognormal Distribution */

sub dlnorm(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub plnorm(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qlnorm(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rlnorm(num64, num64) returns num64
    is native( RMATH ) { * }; 


#	/* Chi-squared Distribution */

sub dchisq(num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pchisq(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qchisq(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rchisq(num64) returns num64
    is native( RMATH ) { * }; 


#	/* Non-central Chi-squared Distribution */

sub dnchisq(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pnchisq(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qnchisq(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rnchisq(num64, num64) returns num64
    is native( RMATH ) { * }; 

#	/* F Distibution */


sub df(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pf(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qf(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rf(num64, num64) returns num64
    is native( RMATH ) { * }; 


#	/* Student t Distibution */


sub dt(num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pt(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qt(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rt(num64) returns num64
    is native( RMATH ) { * }; 


#	/* Binomial Distribution */

sub dbinom(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pbinom(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qbinom(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rbinom(num64, num64) returns num64
    is native( RMATH ) { * }; 

#	/* Multnomial Distribution */

sub rmultinom(int32, CArray[num64], int32, CArray[int32])
    is native( RMATH ) { * }; 




#	/* Cauchy Distribution */

sub dcauchy(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pcauchy(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qcauchy(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rcauchy(num64, num64) returns num64
    is native( RMATH ) { * }; 




#	/* Exponential Distribution */

sub dexp(num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub pexp(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qexp(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rexp(num64) returns num64
    is native( RMATH ) { * }; 




#	/* Geometric Distribution */


sub dgeom(num64, num64, int32) returns num64
    is native( RMATH ) { * };
sub pgeom(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub qgeom(num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub rgeom(num64) returns num64
    is native( RMATH ) { * };


#	/* Hypergeometric Distibution */


sub dhyper(num64, num64, num64, num64, int32) returns num64
    is native( RMATH ) { * }; 
sub phyper(num64, num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub qhyper(num64, num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * }; 
sub rhyper(num64, num64, num64) returns num64
    is native( RMATH ) { * }; 




#	/* Negative Binomial Distribution */

sub dnbinom(num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pnbinom(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qnbinom(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rnbinom(num64, num64) returns num64 
    is native( RMATH ) { * };


sub dnbinom_mu(num64, num64, num64, int32) returns num64
    is native( RMATH ) { * };
sub pnbinom_mu(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub qnbinom_mu(num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub rnbinom_mu(num64, num64) returns num64
    is native( RMATH ) { * };




#	/* Poisson Distribution */

sub dpois(num64, num64, int32) returns num64
    is native( RMATH ) { * };
sub ppois(num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qpois(num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rpois(num64) returns num64 
    is native( RMATH ) { * };

#	/* Weibull Distribution */

sub dweibull(num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pweibull(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qweibull(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rweibull(num64, num64) returns num64 
    is native( RMATH ) { * };


#	/* Logistic Distribution */

sub dlogis(num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub plogis(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qlogis(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rlogis(num64, num64) returns num64 
    is native( RMATH ) { * };



#	/* Non-central Beta Distribution */

sub dnbeta(num64, num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pnbeta(num64, num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qnbeta(num64, num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rnbeta(num64, num64, num64) returns num64 
    is native( RMATH ) { * };

#	/* Non-central F Distribution */

sub dnf(num64, num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pnf(num64, num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qnf(num64, num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };




#	/* Non-central Student t Distribution */

sub dnt(num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pnt(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qnt(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };

#	/* Studentized Range Distribution */

sub ptukey(num64, num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };
sub qtukey(num64, num64, num64, num64, int32, int32) returns num64
    is native( RMATH ) { * };




#	/* Wilcoxon Rank Sum Distribution */

sub dwilcox(num64, num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub pwilcox(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qwilcox(num64, num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rwilcox(num64, num64) returns num64 
    is native( RMATH ) { * };
sub wilcox_free()
    is native( RMATH ) { * };

#	/* Wilcoxon Signed Rank Distribution */

sub dsignrank(num64, num64, int32) returns num64 
    is native( RMATH ) { * };
sub psignrank(num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub qsignrank(num64, num64, int32, int32) returns num64 
    is native( RMATH ) { * };
sub rsignrank(num64) returns num64 
    is native( RMATH ) { * };
sub signrank_free()
    is native( RMATH ) { * };

#	/* Gamma and Related Functions */

sub gammafn(num64) returns num64
    is native( RMATH ) { * }; # gamma in R
sub lgammafn(num64) returns num64 
    is native( RMATH ) { * }; # lgamma in R
multi psigamma(num64, num64) returns num64 
    is native( RMATH ) { * };
multi digamma(num64) returns num64 
    is native( RMATH ) { * };
multi trigamma(num64) returns num64 
    is native( RMATH ) { * };



multi beta(num64, num64) returns num64 
    is native( RMATH ) { * };
multi lbeta(num64, num64) returns num64 
    is native( RMATH ) { * };
multi choose(num64, num64) returns num64
    is native( RMATH ) { * };
multi lchoose(num64, num64) returns num64 
    is native( RMATH ) { * };

#	/* Bessel Functions */

sub bessel_i(num64, num64, num64) returns num64 
    is native( RMATH ) { * };
sub bessel_j(num64, num64) returns num64 
    is native( RMATH ) { * };
sub bessel_k(num64, num64, num64) returns num64 
    is native( RMATH ) { * };
sub bessel_y(num64, num64) returns num64 
    is native( RMATH ) { * };




####	Random number generators	####

sub norm_rand() returns num64
    is export
    is native( RMATH ) { * };

sub unif_rand() returns num64
    is export
    is native( RMATH ) { * };

sub exp_rand() returns num64
    is export
    is native( RMATH ) { * };

multi set_seed(uint32, uint32)
    is native( RMATH ) { * };

multi set_seed(UInt $a, UInt $b) is export {
    # Don't give negative parameters, only unsigned Ints
    set_seed(my uint32 $ = $a, my uint32 $ = $b)
}

sub get_seed(uint32 is rw, uint32 is rw)
    is native( RMATH ) { * };

multi set_N01_kind(int32) is native( RMATH ) { * };

multi set_N01_kind(UInt $value) is export { # this one is exported 
    # create error
    my $error = qq:to/END/;
    The N01_type can only be one of: \n 
    BUGGY_KINDERMAN_RAMAGE => 0,
    AHRENS_DIETER => 1,
    BOX_MULLER => 2,
    INVERSION => 4,
    KINDERMAN_RAMAGE => 5

    USER_NORM not supported in standalone.
    END
    # check for acceptable value
    if $value !~~ (0|1|2|4|5) {die $error};
    set_N01_kind(my int32 $ = $value)        
}



####	Normal distribution	####

# dnorm(x, mean = 0, sd = 1, log = FALSE)
# pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
# qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
# rnorm(n, mean = 0, sd = 1)


sub raku_dnorm($x, :$mean = 0 , :$sd = 1, :$log = False ) is export {
    return dnorm4($x.Num, $mean.Num, $sd.Num, $log ?? 1 !! 0);
}

sub raku_pnorm($q, :$mean = 0 , :$sd = 1, :$lower_tail = True, :$log_p = False) 
    is export {
    return pnorm5($q.Num, $mean.Num, $sd.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
}

sub raku_qnorm($p, :$mean = 0 , :$sd = 1, :$lower_tail = True, :$log_p = False) 
    is export {
    return qnorm5($p.Num, $mean.Num, $sd.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
}

sub raku_rnorm($n, :$mean = 0, :$sd = 1) is export {
    (1..$n).map: {rnorm($mean.Num, $sd.Num)}
}

####	Uniform distribution	####

# dunif(x, min = 0, max = 1, log = FALSE)
# punif(q, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
# qunif(p, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
# runif(n, min = 0, max = 1)

sub raku_dunif($x, :$min = 0 , :$max = 1, :$log = False) 
    is export {
    dunif($x.Num, $min.Num, $max.Num, $log ?? 1 !! 0);
}

sub raku_punif($q, :$min = 0 , :$max = 1, :$lower_tail = True, :$log_p = False) 
    is export {
    punif($q.Num, $min.Num, $max.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
}

sub raku_qunif($p, :$min = 0 , :$max = 1, :$lower_tail = True, :$log_p = False) 
    is export {
    qunif($p.Num, $min.Num, $max.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);

}

sub raku_runif($n, :$min = 0 , :$max = 1) 
    is export {
    (1..$n).map: {runif($min.Num, $max.Num)};
}


####	Gamma distribution 	####

# dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
# pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,
#        log.p = FALSE)
# qgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,
#        log.p = FALSE)
# rgamma(n, shape, rate = 1, scale = 1/rate)



sub raku_dgamma($x, $shape , :$rate is copy, :$scale is copy, :$log = False) 
    is export {
    if $rate.defined && $scale.defined {
        if abs($rate*$scale - 1) < 1e-15 {
            warn("Warning:\n\tspecify 'rate' or 'scale' but not both")
        }
        else {
            die "specify 'rate' or 'scale' but not both"
        }
    } 
    $rate //= 1;
    $scale //= 1/$rate;
    return dgamma($x.Num, $shape.Num, $scale.Num, $log ?? 1 !! 0);    
}


sub raku_pgamma($q, $shape , :$rate is copy, :$scale is copy, :$lower_tail = True, :$log_p = False) 
    is export {
    if $rate.defined && $scale.defined {
        if abs($rate*$scale - 1) < 1e-15 {
            warn("Warning:\n\tspecify 'rate' or 'scale' but not both")
        }
        else {
            die "specify 'rate' or 'scale' but not both"
        }
    } 
    $rate //= 1;
    $scale //= 1/$rate;
    return pgamma($q.Num, $shape.Num, $scale.Num, $lower_tail ?? 1 !! 0,  $log_p ?? 1 !! 0);    
}


sub raku_qgamma($p, $shape , :$rate is copy, :$scale is copy, :$lower_tail = True, :$log_p = False) 
    is export {
    if $rate.defined && $scale.defined {
        if abs($rate*$scale - 1) < 1e-15 {
            warn("Warning:\n\tspecify 'rate' or 'scale' but not both")
        }
        else {
            die "specify 'rate' or 'scale' but not both"
        }
    } 
    $rate //= 1;
    $scale //= 1/$rate;
    return qgamma($p.Num, $shape.Num, $scale.Num, $lower_tail ?? 1 !! 0,  $log_p ?? 1 !! 0);    
}


sub raku_rgamma($n, $shape , :$rate is copy, :$scale is copy) 
    is export {
    if $rate.defined && $scale.defined {
        if abs($rate*$scale - 1) < 1e-15 {
            warn("Warning:\n\tspecify 'rate' or 'scale' but not both")
        }
        else {
            die "specify 'rate' or 'scale' but not both"
        }
    } 
    $rate //= 1;
    $scale //= 1/$rate;
    (1..$n).map: {rgamma($shape.Num, $scale.Num)}
}



####	Beta distribution	####

# dbeta(x, shape1, shape2, ncp = 0, log = FALSE)
# pbeta(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
# qbeta(p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
# rbeta(n, shape1, shape2, ncp = 0)


sub raku_dbeta($x, $shape1, $shape2, :$ncp is copy, :$log = False)
    is export {
        if !$ncp.defined {
            return dbeta($x.Num, $shape1.Num, $shape2.Num, $log ?? 1 !! 0)
        } else {
            $ncp //= 0;
            return dnbeta($x.Num, $shape1.Num, $shape2.Num, $ncp.Num, $log ?? 1 !! 0)
        }
    }

sub raku_pbeta($q, $shape1, $shape2, :$ncp is copy, :$lower_tail = True, :$log_p = False)
    is export {
        if !$ncp.defined {
            return pbeta($q.Num, $shape1.Num, $shape2.Num, $lower_tail ?? 1 !! 0 ,$log_p ?? 1 !! 0)
        } else {
            $ncp //= 0;
            return pnbeta($q.Num, $shape1.Num, $shape2.Num, $ncp.Num, $lower_tail ?? 1 !! 0 ,$log_p ?? 1 !! 0)
        }
    }

sub raku_qbeta($p, $shape1, $shape2, :$ncp is copy, :$lower_tail = True, :$log_p = False)
    is export {
        if !$ncp.defined {
            return qbeta($p.Num, $shape1.Num, $shape2.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else {
            $ncp //= 0;
            return qnbeta($p.Num, $shape1.Num, $shape2.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }


# needs raku_chisq function
sub raku_rbeta($n, $shape1, $shape2, :$ncp is copy)
    is export {
        if !$ncp.defined {
            return (1..$n).map: {rbeta($shape1.Num, $shape2.Num)}
        } else {
            $ncp //= 0;
            my $x = raku_rchisq($n, 2 * $shape1, ncp => $ncp);
            $x/($x + raku_rchisq($n, 2 * $shape2))
        }
    }

####	Lognormal distribution	####

# dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)
# plnorm(q, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# qlnorm(p, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# rlnorm(n, meanlog = 0, sdlog = 1)

sub raku_dlnorm($x, :$meanlog = 0, :$sdlog = 1, :$log = False)
    is export {
        dlnorm($x.Num, $meanlog.Num, $sdlog.Num, $log ?? 1 !! 0)
    }

sub raku_plnorm($q, :$meanlog = 0, :$sdlog = 1, :$lower_tail = True, :$log_p = False)
    is export {
        plnorm($q.Num, $meanlog.Num, $sdlog.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }

sub raku_qlnorm($p, :$meanlog = 0, :$sdlog = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qlnorm($p.Num, $meanlog.Num, $sdlog.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }

sub raku_rlnorm($n, :$meanlog = 0, :$sdlog = 1)
    is export {
        (1..$n).map: {rlnorm($meanlog.Num, $sdlog.Num)}
    }


####	Chi-squared distribution	####

# dchisq(x, df, ncp = 0, log = FALSE)
# pchisq(q, df, ncp = 0, lower.tail = TRUE, log.p = FALSE)
# qchisq(p, df, ncp = 0, lower.tail = TRUE, log.p = FALSE)
# rchisq(n, df, ncp = 0)

sub raku_dchisq($x, $df, :$ncp is copy, :$log = False) is export { 
    if !$ncp.defined {
        dchisq($x.Num, $df.Num, $log ?? 1 !! 0 )
    } else {
        $ncp //= 0;
        return dnchisq($x.Num, $df.Num, $ncp.Num, $log ?? 1 !! 0 )
    }
}

sub raku_pchisq($q, $df, :$ncp is copy, :$lower_tail = True, :$log_p = False) is export { 
    if !$ncp.defined {
        pchisq($q.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0 )
    } else {
        $ncp //= 0;
        return pnchisq($q.Num, $df.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0 )
    }
}

sub raku_qchisq($p, $df, :$ncp is copy, :$lower_tail = True, :$log_p = False) is export { 
    if !$ncp.defined {
        qchisq($p.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0 )
    } else {
        $ncp //= 0;
        return qnchisq($p.Num, $df.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0 )
    }
 }

sub raku_rchisq($n, $df, :$ncp is copy) is export { 
    if !$ncp.defined {
        return (1..$n).map: {rchisq($df.Num)}
    } else {
        $ncp //= 0;
        return (1..$n).map: {rnchisq($df.Num, $ncp.Num)}
    }
}


####	F Distribution	####

# df(x, df1, df2, ncp, log = FALSE)
# pf(q, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)
# qf(p, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)
# rf(n, df1, df2, ncp)


sub raku_df($x!, $df1!, $df2!, $ncp?, :$log = False) 
    is export {
        with $ncp {
            return dnf($x.Num, $df1.Num, $df2.Num, $ncp.Num, $log ?? 1 !! 0)
        } else {
            return df($x.Num, $df1.Num, $df2.Num, $log ?? 1 !! 0)
        }
    }


sub raku_pf($q!, $df1!, $df2!, $ncp?, :$lower_tail = True, :$log_p = False) 
    is export {
        with $ncp {
            return pnf($q.Num, $df1.Num, $df2.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else {
            return pf($q.Num, $df1.Num, $df2.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }


sub raku_qf($p!, $df1!, $df2!, $ncp?, :$lower_tail = True