### ----------------------------------------------------
### -- Statistics
### -- Licenses: Artistic-2.0
### -- Authors: Suman Khanal <suman81765@gmail.com>
### -- File: lib/Probability/Distribution.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 2
, :$log_p = False) 
    is export {
        with $ncp {
            return qnf($p.Num, $df1.Num, $df2.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else {
            return qf($p.Num, $df1.Num, $df2.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }


sub raku_rf($n, $df1!, $df2!, $ncp?) 
    is export {
        with $ncp {
            return (1..$n).map: {(rnchisq($df1.Num, $ncp.Num)/$df1)/(rchisq($df2.Num)/$df2)}
        } else {
            return (1..$n).map: {rf($df1.Num, $df2.Num)}
        }
    }



####	Student t distribution	####

# dt(x, df, ncp, log = FALSE)
# pt(q, df, ncp, lower.tail = TRUE, log.p = FALSE)
# qt(p, df, ncp, lower.tail = TRUE, log.p = FALSE)
# rt(n, df, ncp)

sub raku_dt($x!, $df!, $ncp?, :$log = False) 
    is export {
        with $ncp {
            dnt($x.Num, $df.Num, $ncp.Num, $log ?? 1 !! 0)
        } else {
            dt($x.Num, $df.Num, $log ?? 1 !! 0)
        }
    }


sub raku_pt($q!, $df!, $ncp?, :$lower_tail = True, :$log_p = False) 
    is export {
        with $ncp {
            pnt($q.Num, $df.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else{
            pt($q.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }



sub raku_qt($p!, $df!, $ncp?, :$lower_tail = True, :$log_p = False) 
    is export {
        with $ncp {
            qnt($p.Num, $df.Num, $ncp.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else{
            qt($p.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }



sub raku_rt($n, $df, $ncp?) 
    is export {
        with $ncp {
            (1..$n).map: {rnorm($ncp.Num, 1.Num)/sqrt(rchisq($df.Num)/$df)}
        } else {
            (1..$n).map: {rt($df.Num)}
        }
    }


####	Binomial distribution	####

# dbinom(x, size, prob, log = FALSE)
# pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
# qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)
# rbinom(n, size, prob)


sub raku_dbinom($x!, $size!, $prob!, :$log = False)
    is export {
        dbinom($x.Num, $size.Num, $prob.Num, $log ?? 1 !! 0)
    }


sub raku_pbinom($q!, $size!, $prob!, :$lower_tail = True, :$log_p = False)
    is export {
        pbinom($q.Num, $size.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_qbinom($p, $size, $prob, :$lower_tail = True, :$log_p = False)
    is export {
        qbinom($p.Num, $size.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_rbinom($n!, $size!, $prob!)
    is export {
        (1..$n).map: {rbinom($size.Num, $prob.Num)}
    }

#===========================================================
#===========================================================
####	Multinomial distribution	####


# https://stackoverflow.com/questions/63368639/how-to-use-multinomial-function-in-c-with-rmath-h

sub raku_rmultinom($n , $size, @prob is copy) is export {
    @prob.append(0 xx $size - @prob);
    my $prob = CArray[num64].new(@prob».Num);
    my $ints = CArray[int32].allocate($size);
    rmultinom($n, $prob, $size, $ints);
    return $ints.list
}

# sub raku_dmultinom($x, $size = NULL, $prob, $log = False)
#     is export {

#     }


####	Cauchy distribution	####

# dcauchy(x, location = 0, scale = 1, log = FALSE)
# pcauchy(q, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
# qcauchy(p, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
# rcauchy(n, location = 0, scale = 1)

sub raku_dcauchy($x, $location = 0, $scale = 1, :$log = False)
    is export {
        dcauchy($x.Num, $location.Num, $scale.Num, $log ?? 1 !! 0);
    }


sub raku_pcauchy($q, $location = 0, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        pcauchy($q.Num, $location.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }



sub raku_qcauchy($p, $location = 0, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qcauchy($p.Num, $location.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }


sub raku_rcauchy($n, $location = 0, $scale = 1)
    is export {
        (1..$n).map: {rcauchy($location.Num, $scale.Num)}
    }

####	Exponential distribution	####

# dexp(x, rate = 1, log = FALSE)
# pexp(q, rate = 1, lower.tail = TRUE, log.p = FALSE)
# qexp(p, rate = 1, lower.tail = TRUE, log.p = FALSE)
# rexp(n, rate = 1)

sub raku_dexp($x!, $rate = 1, :$log = False)
    is export {
        dexp($x.Num, (1/$rate).Num, $log ?? 1 !! 0 )
    }


sub raku_pexp($q, $rate = 1, :$lower_tail = True, :$log_p = False)
    is export {
        pexp($q.Num, (1/$rate).Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }

sub raku_qexp($p, $rate = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qexp($p.Num, (1/$rate).Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_rexp($n, $rate = 1)
    is export {
        (1..$n).map: {rexp((1/$rate).Num)}
    }

####	Geometric distribution	####

# dgeom(x, prob, log = FALSE)
# pgeom(q, prob, lower.tail = TRUE, log.p = FALSE)
# qgeom(p, prob, lower.tail = TRUE, log.p = FALSE)
# rgeom(n, prob)

sub raku_dgeom($x!, $prob!, :$log = False)
    is export {
        dgeom($x.Num, $prob.Num, $log ?? 1 !! 0);
    }

sub raku_pgeom($q!, $prob!, :$lower_tail = True, :$log_p = False)
    is export {
        pgeom($q.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }

sub raku_qgeom($p!, $prob!, :$lower_tail = True, :$log_p = False)
    is export {
        qgeom($p.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }

sub raku_rgeom($n!, $prob!)
    is export {
        (1..$n).map: {rgeom($prob.Num)}
    }

####	Hypergeometric distribution	####

# dhyper(x, m, n, k, log = FALSE)
# phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
# qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
# rhyper(nn, m, n, k)

sub raku_dhyper($x!, $m!, $n!, $k!, :$log = False)
    is export {
        dhyper($x.Num, $m.Num, $n.Num, $k.Num, $log ?? 1 !! 0);
    }

 
sub raku_phyper($q!, $m!, $n!, $k!, :$lower_tail = True, :$log_p = False)
    is export {
        phyper($q.Num, $m.Num, $n.Num, $k.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }


sub raku_qhyper($p!, $m!, $n!, $k!, :$lower_tail = True, :$log_p = False)
    is export {
        qhyper($p.Num, $m.Num, $n.Num, $k.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }

 
sub raku_rhyper($nn, $m, $n, $k)
    is export {
        (1..$nn).map: {rhyper($m.Num, $n.Num, $k.Num)}
    }

####	Negative binomial distribution	####

# dnbinom(x, size, prob, mu, log = FALSE)
# pnbinom(q, size, prob, mu, lower.tail = TRUE, log.p = FALSE)
# qnbinom(p, size, prob, mu, lower.tail = TRUE, log.p = FALSE)
# rnbinom(n, size, prob, mu)

sub raku_dnbinom($x, $size, :$prob, :$mu, :$log = False)
    is export {
        if $mu.defined {
            if $prob.defined {
                die("'prob' and 'mu' both specified")
            }
            dnbinom_mu($x.Num, $size.Num, $mu.Num, $log ?? 1 !! 0)
        } else {
            dnbinom($x.Num, $size.Num, $prob.Num, $log ?? 1 !! 0)
        }
    }




sub raku_pnbinom($q, $size, :$prob, :$mu, :$lower_tail = True, :$log_p = False)
    is export {
        if $mu.defined {
            if $prob.defined {
                die("'prob' and 'mu' both specified")
            }
            pnbinom_mu($q.Num, $size.Num, $mu.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else {
            pnbinom($q.Num, $size.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }





sub raku_qnbinom($p, $size, :$prob, :$mu, :$lower_tail = True, :$log_p = False)
    is export {
        if $mu.defined {
            if $prob.defined {
                die("'prob' and 'mu' both specified")
            }
            qnbinom_mu($p.Num, $size.Num, $mu.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        } else {
            qnbinom($p.Num, $size.Num, $prob.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
        }
    }





sub raku_rnbinom($n, $size, :$prob, :$mu)
    is export {
        if !$prob.defined && !$mu.defined {
            die "Supply atleast \$prob or \$mu"
        }
        if $mu.defined {
            if $prob.defined {
                die("'prob' and 'mu' both specified")
            }
            (1..$n).map: {rnbinom_mu($size.Num, $mu.Num)}
        } else {
            (1..$n).map: {rnbinom($size.Num, $prob.Num)}
        }
    }

####	Poisson distribution	####

# dpois(x, lambda, log = FALSE)
# ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)
# qpois(p, lambda, lower.tail = TRUE, log.p = FALSE)
# rpois(n, lambda)

sub raku_dpois($x, $lambda, :$log = False)
    is export {
        dpois($x.Num, $lambda.Num, $log ?? 1 !! 0)
    }



sub raku_ppois($q, $lambda, :$lower_tail = True, :$log_p = False)
    is export {
        ppois($q.Num, $lambda.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }




sub raku_qpois($p, $lambda, :$lower_tail = True, :$log_p = False)
    is export {
        qpois($p.Num, $lambda.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_rpois($n, $lambda)
    is export {
        (1..$n).map: {rpois($lambda.Num)}
    }

####	Weibull Distribution	####

# dweibull(x, shape, scale = 1, log = FALSE)
# pweibull(q, shape, scale = 1, lower.tail = TRUE, log.p = FALSE)
# qweibull(p, shape, scale = 1, lower.tail = TRUE, log.p = FALSE)
# rweibull(n, shape, scale = 1)

sub raku_dweibull($x, $shape, $scale = 1, :$log = False)
    is export {
        dweibull($x.Num, $shape.Num, $scale.Num, $log ?? 1 !! 0)
    }


sub raku_pweibull($q, $shape, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        pweibull($q.Num, $shape.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_qweibull($p, $shape, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qweibull($p.Num, $shape.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_rweibull($n, $shape, $scale = 1)
    is export {
        (1..$n).map: {rweibull($shape.Num, $scale.Num)}
    }

####	Logistic distribution	####

# dlogis(x, location = 0, scale = 1, log = FALSE)
# plogis(q, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
# qlogis(p, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
# rlogis(n, location = 0, scale = 1)



sub raku_dlogis($x, $location = 0, $scale = 1, :$log = False)
    is export {
        dlogis($x.Num, $location.Num, $scale.Num, $log ?? 1 !! 0)
    }


sub raku_plogis($q, $location = 0, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        plogis($q.Num, $location.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_qlogis($p, $location = 0, $scale = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qlogis($p.Num, $location.Num, $scale.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)

    }


sub raku_rlogis($n, $location = 0, $scale = 1)
    is export {
        (1..$n).map: {rlogis($location.Num, $scale.Num)}
    }

####	Studentized Range Distribution	####

# ptukey(q, nmeans, df, nranges = 1, lower.tail = TRUE, log.p = FALSE)
# qtukey(p, nmeans, df, nranges = 1, lower.tail = TRUE, log.p = FALSE)

sub raku_ptukey($q, $nmeans, $df, $nranges = 1, :$lower_tail = True, :$log_p = False)
    is export {
        ptukey($q.Num, $nranges.Num, $nmeans.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }


sub raku_qtukey($p, $nmeans, $df, $nranges = 1, :$lower_tail = True, :$log_p = False)
    is export {
        qtukey($p.Num, $nranges.Num, $nmeans.Num, $df.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0)
    }

####	Wilcoxon Rank Sum Distribution	####

# dwilcox(x, m, n, log = FALSE)
# pwilcox(q, m, n, lower.tail = TRUE, log.p = FALSE)
# qwilcox(p, m, n, lower.tail = TRUE, log.p = FALSE)
# rwilcox(nn, m, n)

sub raku_dwilcox($x, $m, $n, :$log = False)
    is export {
        LEAVE { wilcox_free() };
        dwilcox($x.Num, $m.Num, $n.Num, $log ?? 1 !! 0);
    }



sub raku_pwilcox($q, $m, $n, :$lower_tail = True, :$log_p = False)
    is export {
        LEAVE { wilcox_free() };
        pwilcox($q.Num, $m.Num, $n.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }


sub raku_qwilcox($p, $m, $n, :$lower_tail = True, :$log_p = False)
    is export {
        LEAVE { wilcox_free() };
        qwilcox($p.Num, $m.Num, $n.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }


sub raku_rwilcox($nn, $m, $n)
    is export {
        (1..$nn).map: {rwilcox($m.Num, $n.Num)}
    }


####	Wilcoxon Signed Rank Distribution	####

# dsignrank(x, n, log = FALSE)
# psignrank(q, n, lower.tail = TRUE, log.p = FALSE)
# qsignrank(p, n, lower.tail = TRUE, log.p = FALSE)
# rsignrank(nn, n)

sub raku_dsignrank($x, $n, :$log = False)
    is export {
        LEAVE { signrank_free() };
        dsignrank($x.Num, $n.Num, $log ?? 1 !! 0);
    }


sub raku_psignrank($q, $n, :$lower_tail = True, :$log_p = False)
    is export {
        LEAVE { signrank_free() };
        psignrank($q.Num, $n.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }


sub raku_qsignrank($p, $n, :$lower_tail = True, :$log_p = False)
    is export {
        LEAVE { signrank_free() };
        qsignrank($p.Num, $n.Num, $lower_tail ?? 1 !! 0, $log_p ?? 1 !! 0);
    }
 
sub raku_rsignrank($nn, $n)
    is export {
        (1..$nn).map: {rsignrank($n.Num)}
    }


# Other functions

multi beta(Num() $a, Num() $b) is export 
{
    return beta(my num64 $ = $a, my num64 $ = $b)
}

multi lbeta(Num() $a, Num() $b) is export 
{
    return lbeta(my num64 $ = $a, my num64 $ = $b)
}

sub gamma($x) is export 
{
    return gammafn($x.Num)
}

sub lgamma($x) is export 
{
    return lgammafn($x.Num)
}

multi psigamma(Num() $x, Num() $deriv = 0) is export 
{
    return psigamma(my num64 $ = $x, my num64 $ = $deriv)
}

multi digamma(Num() $x) is export 
{
    return digamma(my num64 $ = $x)
}

multi trigamma(Num() $x) is export 
{
    return trigamma(my num64 $ = $x)
}

multi choose(Num() $n, Num() $k) is export
{
    return choose(my num64 $ = $n, my num64 $ = $k)
}

multi lchoose(Num() $n, Num() $k) is export 
{
    return lchoose(my num64 $ = $n, my num64 $ = $k)
}

sub factorial($x) is export
{
    gammafn(($x + 1).Num)
}

sub lfactorial($x) is export
{
    lgammafn(($x + 1).Num)
}


