R/RcppExports.R

In Rsubbotools: Fast Estimation of Subbottin and AEP Distributions (Generalized Error Distribution)

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Generates a Gamma-distributed sample
#'
#' Returns a sample from a gamma-distributed random variable.
#' The function uses the Marsaglia-Tsang fast gamma method used to generate
#' the random samples. See more details below. The name was chosen so it
#' doesn't clash with R's native method.
#'
#' The gamma distribution is given by the function:
#' \deqn{f(x) = \frac{1}{\Gamma(b)a^b}x^{b-1}e^{-x/a}, x > 0}
#' where \eqn{b} is a shape parameter and \eqn{a} is a scale parameter.
#' The RNG is given by Marsaglia and Tsang, "A Simple Method for
#' generating gamma variables", ACM Transactions on Mathematical
#' Software, Vol 26, No 3 (2000), p363-372.
#' Available at \doi{10.1145/358407.358414}.
#' The code is based on the original 'GSL' version, adapted to
#' use 'R' version of RNGs. All credits to the original authors.
#' Implemented by J.D.Lamb@btinternet.com, minor modifications for 'GSL'
#' by Brian Gough. Adapted to 'R' by Elias Haddad.
#'
#' @param n (int)
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a numeric vector containing a random sample with above parameters.
#'
#' @examples
#' sample_gamma <- rgamma_c(1000, 1, 1)
#' @export
#' @md
rgamma_c <- function(n, b = 2, a = 1/2) {
    .Call(`_Rsubbotools_rgamma_c`, n, b, a)
}

#' Generates a Laplace-distributed sample
#'
#' Returns a sample from a Laplace-distributed random variable.
#'
#' The Laplace distribution is given by the two-sided exponential distribution
#' given by the function:
#' \deqn{ f(x;a,m) = \frac{1}{2a} e^{- \left| \frac{x-m}{a} \right|} }
#' The random sampling is done by inverse transform sampling.
#' @param n (int) - the size of the sample.
#' @param m (numeric) - the location parameter.
#' @param a (numeric) - the scale parameter.
#' @return a numeric vector containing a random sample with above parameters.
#'
#' @examples
#' sample_gamma <- rlaplace(1000, 0, 1)
#' @export
#' @md
rlaplace <- function(n, m = 0, a = 1) {
    .Call(`_Rsubbotools_rlaplace`, n, m, a)
}

#' Generates an Asymmetric Laplace-distributed sample
#'
#' Returns a sample from an Asymetric Laplace distribution.
#'
#' The Asymmetric Laplace distribution is given by the two-sided exponential
#' distribution given by the function:
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_l}| }, x < m
#' }
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_r}| }, x > m
#' }
#' with:
#' \deqn{A = a_l + a_r}
#' The random sampling is done by inverse transform sampling.
#' @param n (int) - the size of the sample.
#' @param m (numeric) - the location parameter.
#' @param al,ar (numeric) - left and right scale parameters, respectively.
#' @return a numeric vector containing a random sample.
#'
#' @examples
#' sample_gamma <- ralaplace(1000)
#' @export
#' @md
ralaplace <- function(n, m = 0, al = 1, ar = 1) {
    .Call(`_Rsubbotools_ralaplace`, n, m, al, ar)
}

#' Generates a random sample from a Exponential Power distribution
#'
#' Returns a sample from a gamma-distributed random variable.
#'
#' The exponential power distribution (EP) is given by the function:
#' \deqn{ f(a,b) = \frac{1}{2a\Gamma(1+1/b)}e^{-|x/a|^b}, -\infty < x < \infty }.
#' where \eqn{b} is a shape parameter, \eqn{a} is a scale parameter and \eqn{\Gamma}
#' representes the gamma function. While not done here, the distribution can
#' be adapted to have non-zero location parameter.
#' The Exponential Power distribution is related to the gamma distribution by
#' the equation:
#' \deqn{E = a*G(1/b)^{1/b}}
#' where E and G are respectively EP and gamma random variables. This property
#' is used for cases where \eqn{b<1} and \eqn{b>4}. For \eqn{1 \leq b \leq 4}
#' rejection methods based on the Laplace and normal distributions are used,
#' which should be faster.
#' Technical details about this algorithm are available on:
#' P. R. Tadikamalla, "Random Sampling from the Exponential Power
#' Distribution", Journal of the American Statistical Association,
#' September 1980, Volume 75, Number 371, pages 683-686.
#' The code is based on the original 'GSL' version, adapted to
#' use 'R' version of RNGs by Elias Haddad. All credits to the original authors.
#' @param n (int) - size of the sample.
#' @param m (numeric) - the location parameter.
#' @param a (numeric) - scale parameter.
#' @param b (numeric) - shape parameter.
#' @return a numeric vector containing a random sample with above parameters.
#'
#' @examples
#' sample_gamma <- rpower(1000)
#' @export
#' @md
rpower <- function(n, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_rpower`, n, m, a, b)
}

#' Produces a random sample from a Subbotin distribution
#'
#' Generate pseudo random-number from a Subbotin distribution using the
#' Tadikamalla method.
#'
#' The Subbotin distribution is given by the function:
#' \deqn{ f(x;a,b,m) = \frac{1}{A} e^{- \frac{1}{b} \left|\frac{x-m}{a}\right|^b} }
#' where \eqn{m} is a location parameter, \eqn{b} is a shape parameter, \eqn{a}
#' is a scale parameter and \eqn{\Gamma} representes the gamma function.
#' Since the Subbotin distribution is basically the exponential distribution
#' with scale parameter \eqn{a = ab^{1/b}} and \eqn{m=0}, we use the same
#' method of the exponential power RNG and add the location parameter.
#' Details can be found on the documentation of the \code{rpower} function.
#' @param n (int) - the size of the sample.
#' @param m (numeric) - the location parameter.
#' @param a (numeric) - the scale parameter.
#' @param b (numeric) - the shape parameter.
#' @return a numeric vector containing a random sample.
#'
#' @examples
#' sample_gamma <- rsubbo(1000, 1, 1)
#' @export
#' @md
rsubbo <- function(n, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_rsubbo`, n, m, a, b)
}

#' Produces a random sample from a Asymmetric Power Exponential distribution
#'
#'
#' Generate pseudo random-number from an asymmetric power exponential distribution
#' using the Tadikamalla method.
#' This codes is the original version of Bottazzi (2004)
#'
#' The AEP distribution is expressed by the function:
#' \deqn{f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{m} is a location parameter, \eqn{b*} are shape parameters, \eqn{a*}
#' are scale parameters and \eqn{\Gamma} representes the gamma function.
#' By a suitably transformation, it is possible to use the EP distribution with
#' the Tadikamalla method to sample from this distribution. We basically take
#' the absolute values of the numbers sampled from the \code{rpower} function,
#' which is equivalent from sampling from a half Exponential Power distribution.
#' These values are then weighted by a constant expressed in the parameters.
#' More details are available on the package vignette and on the
#' function \code{rpower}.
#' @param n (int) - size of the sample.
#' @param m (numeric) - location parameter.
#' @param bl,br (numeric) - shape parameters.
#' @param al,ar (numeric) - scale parameters.
#' @return a numeric vector containing a random sample.
#'
#' @examples
#' sample_gamma <- rasubbo(1000, 0, 0.5, 0.5,  1, 1)
#' @export
#' @md
rasubbo_orig <- function(n, m = 0, al = 1, ar = 1, bl = 2, br = 2) {
    .Call(`_Rsubbotools_rasubbo_orig`, n, m, al, ar, bl, br)
}

#' Produces a random sample from a Asymmetric Power Exponential distribution
#'
#'
#' Generate pseudo random-number from an asymmetric power exponential distribution
#' using the Tadikamalla method.
#' This version improves on Bottazzi (2004) by making the mass of each
#' distribution to depend on the ratio between the \eqn{al} and the \eqn{ar}
#' parameters.
#'
#' The AEP distribution is expressed by the function:
#' \deqn{f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{m} is a location parameter, \eqn{b*} are shape parameters, \eqn{a*}
#' are scale parameters and \eqn{\Gamma} representes the gamma function.
#' By a suitably transformation, it is possible to use the EP distribution with
#' the Tadikamalla method to sample from this distribution. We basically take
#' the absolute values of the numbers sampled from the \code{rpower} function,
#' which is equivalent from sampling from a half Exponential Power distribution.
#' These values are then weighted by a constant expressed in the parameters.
#' More details are available on the package vignette and on the
#' function \code{rpower}.
#' @param n (int) - size of the sample.
#' @param m (numeric) - location parameter.
#' @param bl,br (numeric) - shape parameters.
#' @param al,ar (numeric) - scale parameters.
#' @return a numeric vector containing a random sample.
#'
#' @examples
#' sample_gamma <- rasubbo(1000, 0, 0.5, 0.5,  1, 1)
#' @export
#' @md
rasubbo <- function(n, m = 0, al = 1, ar = 1, bl = 2, br = 2) {
    .Call(`_Rsubbotools_rasubbo`, n, m, al, ar, bl, br)
}

#' Fit an Asymmetric Laplace Distribution via maximum likelihood
#'
#' \code{alaplafit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the Asymmetric Laplace Distribution
#' for a sample. See details below.
#'
#' The Asymmetric Laplace distribution is a distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a_l,a_r,m) = \frac{1}{A} e^{-|\frac{x-m}{a_l}| }, x < m}
#' \deqn{f(x;a_l,a_r,m) = \frac{1}{A} e^{-|\frac{x-m}{a_r}| }, x > m}
#' with:
#' \deqn{A = a_l + a_r}
#' where \eqn{a*} are scale parameters, and \eqn{m} is a location parameter.
#' It is basically derived from the Asymmetric Exponential Power distribution
#' by setting \eqn{b_l = b_r = b}.
#' The estimations are produced by maximum likelihood, where
#' analytical formulas are available for the \eqn{a*} parameters.
#' The \eqn{m} parameter is found by an iterative method, using the median as
#' the initial guess. The method explore intervals around the last minimum
#' found, similar to the \code{subboafit} routine.
#' Details on the method can be found on the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  details of optim. routine
#' @param interv_step  int - the number of intervals to be explored after
#' the last minimum was found in the interval optimization. Default is 10.
#' @param provided_m_ NumericVector - if NULL, the m parameter is estimated
#' by the routine. If numeric, the estimation fixes m to the given value.
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#' * "matrix" - the covariance matrix for the parameters.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 1)
#' alaplafit(sample_subbo)
#' @export
#' @md
alaplafit <- function(data, verb = 0L, interv_step = 10L, provided_m_ = NULL) {
    .Call(`_Rsubbotools_alaplafit`, data, verb, interv_step, provided_m_)
}

#' Returns CDF from EP Distribution
#'
#' The \code{ppower} returns the Cumulative Distribution Function at point x for
#' the Exponential Power distribution with parameters \eqn{a}, \eqn{b} and \eqn{m}.
#'
#' The Exponential Power distribution (EP) is given by the function:
#' \deqn{ f(a,b) = \frac{1}{2a\Gamma(1+1/b)}e^{-|(x-m)/a|^b}, -\infty < x < \infty }.
#' where \eqn{b} is a shape parameter, \eqn{a} is a scale parameter, \eqn{m}
#' is a location parameter and \eqn{\Gamma} represents the gamma function.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter. Must be in the range
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' \eqn{(-\infty, \infty)}.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
ppower <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_ppower`, x, m, a, b)
}

#' Returns CDF from the Skewed Exponential Power distribution
#'
#' The \code{psep} returns the Cumulative Distribution Function at point x for
#' the Skewed Exponential Power distribution with parameters \eqn{a}, \eqn{b}.
#'
#' The  SEP is a exponential power distribution controlled
#' by four parameters, with formula:
#' \deqn{ f(x; m, b, a, \lambda) = 2 \Phi(w) e^{-|z|^b/b}/(c)}
#' where:
#' \deqn{z = (x-m)/a}
#' \deqn{w = sign(z) |z|^{(b/2)} \lambda \sqrt{2/b}}
#' \deqn{c = 2 ab^{(1/b)-1} \Gamma(1/b)}
#' with \eqn{\Phi} the cumulative normal distribution with mean zero and variance
#' one. The CDF is calculated through numerical integration using the 'GSL' suite.
#' @param x - vector with values to evaluate CDF.
#' @param m - the location parameter.
#' @param a - the scale parameter.
#' @param b - the shape parameter
#' @param lambda - the skewness parameter.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
psep <- function(x, m = 0, a = 2, b = 1, lambda = 0) {
    .Call(`_Rsubbotools_psep`, x, m, a, b, lambda)
}

#' Returns CDF from the Laplace Distribution
#'
#' The \code{plaplace} returns the Cumulative Distribution Function at point x
#' for the Laplace distribution with parameters \eqn{a} and \eqn{m}.
#'
#' The Laplace distribution is a distribution controlled
#' by two parameters, with formula:
#' \deqn{f(x;a,m) = \frac{1}{2a} e^{- \left| \frac{x-m}{a} \right| }}
#' where \eqn{a} is a scale parameter, and \eqn{m} is a location parameter.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param m (numeric) - location parameter.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
plaplace <- function(x, m = 0, a = 1) {
    .Call(`_Rsubbotools_plaplace`, x, m, a)
}

#' Returns CDF from Asymmetric Laplace Distribution
#'
#' The \code{palaplace} returns the Cumulative Distribution Function at point x
#' for the Asymmetric Laplace distribution with parameters \eqn{a*} and \eqn{m}.
#'
#' The Asymmetric Laplace distribution is a distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_l}| }, x < m
#' }
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_r}| }, x > m
#' }
#' with:
#' \deqn{A = a_l + a_r}
#' where \eqn{a*} are scale parameters, and \eqn{m} is a location parameter.
#' It is basically derived from the Asymmetric Exponential Power distribution
#' by setting \eqn{b_l = b_r = b}.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
palaplace <- function(x, m = 0, al = 1, ar = 1) {
    .Call(`_Rsubbotools_palaplace`, x, m, al, ar)
}

#' Returns CDF from Subbotin Distribution
#'
#' The \code{psubbo} returns the Cumulative Distribution Function (CDF) from the
#' the Subbotin evaluated at \eqn{a} and return \eqn{z}, such that
#' \eqn{P(X < a) = z}.
#'
#' The Subbotin cumulative distribution function is given by:
#' \deqn{F(x;a,b,m) = 0.5 + 0.5 \text{sign}(x -m)P(x, 1/b)}
#' where \eqn{P} is the normalized incomplete gamma function:
#' \deqn{P(x, 1/b) = 1 - \frac{1}{\Gamma(1/b)} \int_{0}^{x} t^{1/b -1}e^{-t} }
#' and \eqn{a} is a scale parameter, \eqn{b} controls the tails (lower values
#' represent fatter tails), and \eqn{m} is a location parameter.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
psubbo <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_psubbo`, x, m, a, b)
}

#' Returns CDF from the AEP Distribution
#'
#' The \code{pasubbo} returns the Cumulative Distribution Function at point x
#' for the AEP distribution with parameters \eqn{a*}, \eqn{b*}, \eqn{m}.
#'
#' The AEP is a exponential power distribution controlled
#' by five parameters, with formula:
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{l} and \eqn{r} represent left and right tails, \eqn{a*} are
#' scale parameters, \eqn{b*} control the tails (lower values represent
#' fatter tails), and \eqn{m} is a location parameter.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @param bl,br (numeric) - shape parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the probabilities.
#' @export
#' @md
pasubbo <- function(x, m = 0, al = 1, ar = 1, bl = 2, br = 2) {
    .Call(`_Rsubbotools_pasubbo`, x, m, al, ar, bl, br)
}

sortRcpp <- function(x) {
    invisible(.Call(`_Rsubbotools_sortRcpp`, x))
}

#' Returns density from EP Distribution
#'
#' The \code{dpower} returns the density at point x for the
#' Exponential Power distribution with parameters \eqn{a}, \eqn{b} and \eqn{m}.
#'
#' The Exponential Power distribution (EP) is given by the function:
#' \deqn{ f(a,b) = \frac{1}{2a\Gamma(1+1/b)}e^{-|(x-m)/a|^b}, -\infty < x < \infty }.
#' where \eqn{b} is a shape parameter, \eqn{a} is a scale parameter, \eqn{m}
#' is a location parameter and \eqn{\Gamma} represents the gamma function.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter. Must be in the range
#' \eqn{(-\infty, \infty)}.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dpower <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_dpower`, x, m, a, b)
}

#' Returns density from Skewed Exponential Power distribution
#'
#' The \code{dsep} returns the density at point x for the
#' Gamma distribution with parameters \eqn{a}, \eqn{b}.
#'
#' The  SEP is a exponential power distribution controlled
#' by four parameters, with formula:
#' \deqn{ f(x; m, b, a, \lambda) = 2 \Phi(w) e^{-|z|^b/b}/(c)}
#' where:
#' \deqn{z = (x-m)/a}
#' \deqn{w = sign(z) |z|^{(b/2)} \lambda \sqrt{2/b}}
#' \deqn{c = 2 ab^{(1/b)-1} \Gamma(1/b)}
#' with \eqn{\Phi} the cumulative normal distribution with mean zero and variance
#' one.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter. Must be in the range
#' \eqn{(-\infty, \infty)}.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @param lambda (numeric) - skewness parameter. Must be in the range \eqn{(-\infty, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dsep <- function(x, m = 0, a = 1, b = 1, lambda = 1) {
    .Call(`_Rsubbotools_dsep`, x, m, a, b, lambda)
}

#' Returns density from Laplace Distribution
#'
#' The \code{dlaplace} returns the density at point x for the
#' Laplace distribution with parameters \eqn{a} and \eqn{m}.
#'
#' The Laplace distribution is a distribution controlled
#' by two parameters, with formula:
#' \deqn{f(x;a,m) = \frac{1}{2a} e^{- \left| \frac{x-m}{a} \right| }}
#' where \eqn{a} is a scale parameter, and \eqn{m} is a location parameter.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dlaplace <- function(x, m = 0, a = 1) {
    .Call(`_Rsubbotools_dlaplace`, x, m, a)
}

#' Returns density from Asymmetric Laplace Distribution
#'
#' The \code{dalaplace} returns the density at point x for the
#' Asymmetric Laplace distribution with parameters \eqn{a*} and \eqn{m}.
#'
#' The Asymmetric Laplace distribution is a distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_l}| }, x < m
#' }
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_r}| }, x > m
#' }
#' with:
#' \deqn{A = a_l + a_r}
#' where \eqn{a*} are scale parameters, and \eqn{m} is a location parameter.
#' It is basically derived from the Asymmetric Exponential Power distribution
#' by setting \eqn{b_l = b_r = b}.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dalaplace <- function(x, m = 0, al = 1, ar = 1) {
    .Call(`_Rsubbotools_dalaplace`, x, m, al, ar)
}

#' Returns density from Subbotin Distribution
#'
#' The \code{dsubbo} returns the density at point x for the
#' Subbotin distribution with parameters \eqn{a}, \eqn{b}, \eqn{m}.
#'
#' The Subbotin distribution is a exponential power distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a,b,m) = \frac{1}{A} e^{-\frac{1}{b} |\frac{x-m}{a}|^b}}
#' with:
#' \deqn{A = 2ab^{1/b}\Gamma(1+1/b)}
#' where \eqn{a} is a scale parameter, \eqn{b} controls the tails (lower values
#' represent fatter tails), and \eqn{m} is a location parameter.
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dsubbo <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_dsubbo`, x, m, a, b)
}

#' Returns density from the AEP Distribution
#'
#' The \code{dasubbo} returns the density at point x for the
#' AEP distribution with parameters \eqn{a*}, \eqn{b*}, \eqn{m}. Notice
#' that the function can generate RNGs for both the \code{subboafit} and
#' \code{subbolafit} routines.
#'
#' The AEP is a exponential power distribution controlled
#' by five parameters, with formula:
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{l} and \eqn{r} represent left and right tails, \eqn{a*} are
#' scale parameters, \eqn{b*} control the tails (lower values represent
#' fatter tails), and \eqn{m} is a location parameter.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @param bl,br (numeric) - shape parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
dasubbo <- function(x, m = 0, al = 1, ar = 1, bl = 2, br = 2) {
    .Call(`_Rsubbotools_dasubbo`, x, m, al, ar, bl, br)
}

#' Fit a Laplace Distribution via maximum likelihood
#'
#' \code{laplafit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the Laplace Distribution for a
#' sample. See details below.
#'
#' The Laplace distribution is a distribution controlled
#' by two parameters, with formula:
#' \deqn{f(x;a,m) = \frac{1}{2a} e^{- \left| \frac{x-m}{a} \right| }}
#' where \eqn{a} is a scale parameter, and \eqn{m} is a location parameter.
#' The estimations are produced by maximum likelihood, where analytical
#' formulas are available. Details on the method can be found on
#' the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  details of optim. routine
#' @param interv_step  int - the number of intervals to be explored after
#' the last minimum was found in the interval optimization. Default is 10.
#' @param provided_m_ NumericVector - if NULL, the m parameter is estimated
#' by the routine. If numeric, the estimation fixes m to the given value.
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#' * "matrix" - the covariance matrix for the parameters.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 1)
#' laplafit(sample_subbo)
#' @export
#' @md
laplafit <- function(data, verb = 0L, interv_step = 10L, provided_m_ = NULL) {
    .Call(`_Rsubbotools_laplafit`, data, verb, interv_step, provided_m_)
}

#' Returns quantile from EP Distribution
#'
#' The \code{qpower} returns the density at point x for the
#' Exponential Power distribution with parameters \eqn{a}, \eqn{b} and \eqn{m}.
#'
#' The Exponential Power distribution (EP) is given by the function:
#' \deqn{ f(a,b) = \frac{1}{2a\Gamma(1+1/b)}e^{-|(x-m)/a|^b}, -\infty < x < \infty }.
#' where \eqn{b} is a shape parameter, \eqn{a} is a scale parameter, \eqn{m}
#' is a location parameter and \eqn{\Gamma} represents the gamma function.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter. Must be in the range
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' \eqn{(-\infty, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
qpower <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_qpower`, x, m, a, b)
}

#' Returns quantile from the Skewed Exponential Power distribution
#'
#' The \code{qsep} returns the Cumulative Distribution Function at point x for
#' the Skewed Exponential Power distribution with parameters \eqn{a}, \eqn{b}.
#'
#' The  SEP is a exponential power distribution controlled
#' by four parameters, with formula:
#' \deqn{ f(x; m, b, a, \lambda) = 2 \Phi(w) e^{-|z|^b/b}/(c)}
#' where:
#' \deqn{z = (x-m)/a}
#' \deqn{w = sign(z) |z|^{(b/2)} \lambda \sqrt{2/b}}
#' \deqn{c = 2 ab^{(1/b)-1} \Gamma(1/b)}
#' with \eqn{\Phi} the cumulative normal distribution with mean zero and variance
#' one. The CDF is calculated through numerical integration using the GSL suite
#' and the quantile is solved by inversion using a root-finding algorithm
#' (Newton-Raphson by default).
#' @param x (numeric) - vector with values to evaluate CDF.
#' @param m (numeric) - the location parameter.
#' @param a (numeric) - the scale parameter.
#' @param b (numeric) - the shape parameter
#' @param lambda (numeric) - the skewness parameter.
#' @param method (numeric) - If 0, uses the Newton-Raphson procedure for
#' optimization. If 1, uses Steffensen.
#' @param step_size (numeric) - the size of the step in the numerical
#' optimization (gradient descent). Default is 1e-4.
#' @param tol (numeric) - error tolerance (default is 1e-10).
#' @param max_iter (numeric) - maximum number of iterations for the
#' optimization procedure (default is 100).
#' @param verb (numeric) - verbosity level of the process (default 0).
#' @return a vector containing the values for the densities.
#' @export
#' @md
qsep <- function(x, m = 0, a = 2, b = 1, lambda = 0, method = 0L, step_size = 1e-4, tol = 1e-10, max_iter = 100L, verb = 0L) {
    .Call(`_Rsubbotools_qsep`, x, m, a, b, lambda, method, step_size, tol, max_iter, verb)
}

#' Returns quantile from Laplace Distribution
#'
#' The \code{qlaplace} returns the density at point x for the
#' Laplace distribution with parameters \eqn{a} and \eqn{m}.
#'
#' The Laplace distribution is a distribution controlled
#' by two parameters, with formula:
#' \deqn{f(x;a,m) = \frac{1}{2a} e^{- \left| \frac{x-m}{a} \right| }}
#' where \eqn{a} is a scale parameter, and \eqn{m} is a location parameter.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
qlaplace <- function(x, m = 0, a = 1) {
    .Call(`_Rsubbotools_qlaplace`, x, m, a)
}

#' Returns quantile from Asymmetric Laplace Distribution
#'
#' The \code{qalaplace} returns the density at point x for the
#' Asymmetric Laplace distribution with parameters \eqn{a*} and \eqn{m}.
#'
#' The Asymmetric Laplace distribution is a distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_l}| }, x < m
#' }
#' \deqn{f(x;a_l,a_r,m) =
#' \frac{1}{A} e^{-|\frac{x-m}{a_r}| }, x > m
#' }
#' with:
#' \deqn{A = a_l + a_r}
#' where \eqn{a*} are scale parameters, and \eqn{m} is a location parameter.
#' It is basically derived from the Asymmetric Exponential Power distribution
#' by setting \eqn{b_l = b_r = b}.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
qalaplace <- function(x, m = 0, al = 1, ar = 1) {
    .Call(`_Rsubbotools_qalaplace`, x, m, al, ar)
}

#' Returns CDF from Subbotin Distribution
#'
#' The \code{qsubbo} returns the Cumulative Distribution Function (CDF) from the
#' the Subbotin evaluated at \eqn{a} and return \eqn{z}, such that
#' \eqn{P(X < a) = z}.
#'
#' The Subbotin cumulative distribution function is given by:
#' \deqn{F(x;a,b,m) = 0.5 + 0.5 \text{sign}(x -m)P(x, 1/b)}
#' where \eqn{P} is the normalized incomplete gamma function:
#' \deqn{P(x, 1/b) = 1 - \frac{1}{\Gamma(1/b)} \int_{0}^{x} t^{1/b -1}e^{-t} }
#' and \eqn{a} is a scale parameter, \eqn{b} controls the tails (lower values
#' represent fatter tails), and \eqn{m} is a location parameter.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param a (numeric) - scale parameter. Must be in the range \eqn{(0, \infty)}.
#' @param b (numeric) - shape parameter. Must be in the range \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
qsubbo <- function(x, m = 0, a = 1, b = 2) {
    .Call(`_Rsubbotools_qsubbo`, x, m, a, b)
}

#' Returns CDF from the AEP Distribution
#'
#' The \code{qasubbo} returns the density at point x for the
#' AEP distribution with parameters \eqn{a*}, \eqn{b*}, \eqn{m}.
#'
#' The AEP is a exponential power distribution controlled
#' by five parameters, with formula:
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{l} and \eqn{r} represent left and right tails, \eqn{a*} are
#' scale parameters, \eqn{b*} control the tails (lower values represent
#' fatter tails), and \eqn{m} is a location parameter.
#'
#' @param x (numeric) - value in the range \eqn{(-\infty, \infty)} to evaluate
#' the density.
#' @param m (numeric) - location parameter.
#' @param al,ar (numeric) - scale parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @param bl,br (numeric) - shape parameters. Must be in the range
#' \eqn{(0, \infty)}.
#' @return a vector containing the values for the densities.
#' @export
#' @md
qasubbo <- function(x, m = 0, al = 1, ar = 1, bl = 2, br = 2) {
    .Call(`_Rsubbotools_qasubbo`, x, m, al, ar, bl, br)
}

#' Fit a Skewed Exponential Power density via maximum likelihood
#'
#' \code{sepfit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the skewed power exponential
#' for a sample. The process performs a global minimization over the negative
#' log-likelihood function. See details below.
#'
#' The  SEP is a exponential power distribution controlled
#' by four parameters, with formula:
#' \deqn{ f(x; m, b, a, \lambda) = 2 \Phi(w) e^{-|z|^b/b}/(c)}
#' where:
#' \deqn{z = (x-m)/a}
#' \deqn{w = sign(z) |z|^{(b/2)} \lambda \sqrt{2/b}}
#' \deqn{c = 2 ab^{(1/b)-1} \Gamma(1/b)}
#' with \eqn{\Phi} the cumulative normal distribution with mean zero and
#' variance one.
#' Details on the method are available on the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  headings and summary table
#' * 2  intermediate steps results
#' * 3  intermediate steps internals
#' * 4+  details of optim. routine
#' @param par NumericVector - vector containing the initial guess for parameters
#' m (location), a (scale), b (shape), lambda (skewness), respectively.
#' Default values of are c(0, 1, 2, 0), i.e. a normal distribution.
#' @param g_opt_par NumericVector - vector containing the global optimization
#' parameters.
#' The optimization parameters are:
#' * step  - (num) initial step size of the searching algorithm.
#' * tol   - (num) line search tolerance.
#' * iter  - (int) maximum number of iterations.
#' * eps   - (num) gradient tolerance. The stopping criteria is \eqn{||\text{gradient}||<\text{eps}}.
#' * msize - (num) simplex max size. stopping criteria given by \eqn{||\text{max edge}||<\text{msize}}
#' * algo  - (int) algorithm. the optimization method used:
#'   * 0 Fletcher-Reeves
#'   * 1 Polak-Ribiere
#'   * 2 Broyden-Fletcher-Goldfarb-Shanno
#'   * 3 Steepest descent
#'   * 4 Nelder-Mead simplex
#'   * 5 Broyden-Fletcher-Goldfarb-Shanno ver.2
#'
#' Details for each algorithm are available on the ['GSL' Manual](https://www.gnu.org/software/gsl/doc/html/multimin.html).
#' Default values are c(.1, 1e-2, 100, 1e-3, 1e-5, 2).
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#' * "matrix" - the covariance matrix for the parameters.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 2)
#' sepfit(sample_subbo)
#' @export
#' @md
sepfit <- function(data, verb = 0L, par = as.numeric( c(0., 1., 2., 0.)), g_opt_par = as.numeric( c(.1, 1e-2, 100, 1e-3, 1e-5, 2))) {
    .Call(`_Rsubbotools_sepfit`, data, verb, par, g_opt_par)
}

#' Returns the Fisher Information matrix and its inverse for the AEP
#'
#' Returns the Fisher Information matrix and its inverse for the Asymmetric
#' Power Exponential distribution for the given parameters.
#'
#' @param size (numeric) - number of observations (Default: 01)
#' @param bl (numeric) - set the left exponent (Default: 2.0)
#' @param br (numeric) - set the right exponent (Default: 2.0)
#' @param m  (numeric) - the location parameter (Default: 0.0)
#' @param al (numeric) - the left scale parameter (Default: 1.0)
#' @param ar (numeric) - the right scale parameter (Default: 1.0)
#' @param O_munknown (numeric) - if true assumes \eqn{m} is known
#' @return a list containing three elements:
#' * std_error - the standard error for the parameters
#' * infmatrix - the Fisher Information Matrix
#' * inv_infmatrix - the Inverse Fisher Information Matrix
#' @export
#' @md
subboafish <- function(size = 1L, bl = 2.0, br = 2.0, m = 0.0, al = 1.0, ar = 1.0, O_munknown = 0L) {
    .Call(`_Rsubbotools_subboafish`, size, bl, br, m, al, ar, O_munknown)
}

#' Fit an Asymmetric Power Exponential density via maximum likelihood
#'
#' \code{subboafit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the asymmetric power exponential for
#' a sample. The process can execute two steps, dependending on the level of
#' accuracy required. See details below.
#'
#' The AEP is a exponential power distribution controlled
#' by five parameters, with formula:
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a_l}|^{b_l} }, x < m
#' }
#' \deqn{ f(x;a_l,a_r,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a_r}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = a_lb_l^{1/b_l}\Gamma(1+1/b_l) + a_rb_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{l} and \eqn{r} represent left and right tails, \eqn{a*} are
#' scale parameters, \eqn{b*} control the tails (lower values represent
#' fatter tails), and \eqn{m} is a location parameter. Due to its lack of
#' simmetry, and differently from the Subbotin, there is no simple equations
#' available to use the method of moments, so we start directly by minimizing
#' the negative log-likelihood. This global optimization is executed without
#' restricting any parameters. If required (default), after the global
#' optimization is finished, the method proceeds to iterate over the intervals
#' between several two observations, iterating the same algorithm of the
#' global optimization. The last method happens because of the lack of
#' smoothness on the \eqn{m} parameter, and intervals must be used since the
#' likelihood function doesn't have a derivative whenever \eqn{m} equals a
#' sample observation. Due to the cost, these iterations are capped at most
#' \emph{interv_step} (default 10) from the last minimum observed.
#' Details on the method are available on the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  headings and summary table
#' * 2  intermediate steps results
#' * 3  intermediate steps internals
#' * 4+  details of optim. routine
#' @param method int - the steps that should be used to estimate the
#' parameters.
#' * 0 no optimization perform - just return the log-likelihood from initial guess.
#' * 1 global optimization not considering lack of smoothness in m
#' * 2 interval optimization taking non-smoothness in m into consideration
#' @param interv_step  int - the number of intervals to be explored after
#' the last minimum was found in the interval optimization. Default is 10.
#' @param provided_m_ NumericVector - if NULL, the m parameter is estimated
#' by the routine. If numeric, the estimation fixes m to the given value.
#' @param par NumericVector - vector containing the initial guess for
#' parameters bl, br, al, ar and m, respectively. Default values of are
#' c(2, 2, 1, 1, 0).
#' @param g_opt_par NumericVector - vector containing the global optimization
#' parameters.
#' The optimization parameters are:
#' * step  - (num) initial step size of the searching algorithm.
#' * tol   - (num) line search tolerance.
#' * iter  - (int) maximum number of iterations.
#' * eps   - (num) gradient tolerance. The stopping criteria is \eqn{||\text{gradient}||<\text{eps}}.
#' * msize - (num) simplex max size. stopping criteria given by \eqn{||\text{max edge}||<\text{msize}}
#' * algo  - (int) algorithm. the optimization method used:
#'   * 0 Fletcher-Reeves
#'   * 1 Polak-Ribiere
#'   * 2 Broyden-Fletcher-Goldfarb-Shanno
#'   * 3 Steepest descent
#'   * 4 Nelder-Mead simplex
#'   * 5 Broyden-Fletcher-Goldfarb-Shanno ver.2
#'
#' Details for each algorithm are available on the ['GSL' Manual](https://www.gnu.org/software/gsl/doc/html/multimin.html).
#' Default values are c(.1, 1e-2, 100, 1e-3, 1e-5, 2).
#' @param itv_opt_par NumericVector - interval optimization parameters. Fields
#' are the same as the ones for the global optimization. Default values
#' are c(.01, 1e-3, 200, 1e-3, 1e-5, 5).
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#' * "matrix" - the covariance matrix for the parameters.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 2)
#' subboafit(sample_subbo)
#' @export
#' @md
subboafit <- function(data, verb = 0L, method = 6L, interv_step = 10L, provided_m_ = NULL, par = as.numeric( c(2., 2., 1., 1., 0.)), g_opt_par = as.numeric( c(.1, 1e-2, 100, 1e-3, 1e-5, 2)), itv_opt_par = as.numeric( c(.01, 1e-3, 200, 1e-3, 1e-5, 5))) {
    .Call(`_Rsubbotools_subboafit`, data, verb, method, interv_step, provided_m_, par, g_opt_par, itv_opt_par)
}

#' Return the Fisher Information Matrix for the Subbotin Distribution
#'
#' Calculate the standard errors, the correlation, the
#' Fisher Information matrix and its inverse for a power exponential density
#' with given parameters
#'
#' @param size numeric - number of observations (Default: 01)
#' @param b numeric - the exponent b (Default: 2.0)
#' @param m numeric - the location parameter (Default: 0.0)
#' @param a numeric - the scale parameter (Default: 1.0)
#' @param O_munknown numeric - if true assumes m known
#' @return a list containing four elements:
#' * std_error - the standard error for the parameters
#' * cor_ab    - the correlation between parameters a and b
#' * infmatrix - the Fisher Information Matrix
#' * inv_infmatrix - the Inverse Fisher Information Matrix
#' @export
#' @md
subbofish <- function(size = 1L, b = 2.0, m = 0.0, a = 1.0, O_munknown = 0L) {
    .Call(`_Rsubbotools_subbofish`, size, b, m, a, O_munknown)
}

mm_subbotin <- function(std_over_aad, verb = 0L) {
    .Call(`_Rsubbotools_mm_subbotin`, std_over_aad, verb)
}

optim_method_moments <- function(data, fmin, provided_m_ = NULL, verb = 0L) {
    .Call(`_Rsubbotools_optim_method_moments`, data, fmin, provided_m_, verb)
}

#' Fit a power exponential density via maximum likelihood
#'
#' \code{subbofit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the Subbotin Distribution for a
#' sample. The process can execute three steps, dependending on the level of
#' accuracy required. See details below.
#'
#'
#' The Subbotin distribution is a exponential power distribution controlled
#' by three parameters, with formula:
#' \deqn{f(x;a,b,m) = \frac{1}{A} e^{-\frac{1}{b} |\frac{x-m}{a}|^b}}
#' with:
#' \deqn{A = 2ab^{1/b}\Gamma(1+1/b)}
#' where \eqn{a} is a scale parameter, \eqn{b} controls the tails (lower values
#' represent fatter tails), and \eqn{m} is a location parameter. Due to its
#' simmetry, the equations are simple enough to be estimated by the method of
#' moments, which produce rough estimations that should be used only for first
#' explorations. The maximum likelihood global estimation improves on this
#' initial guess by using a optimization routine, defaulting to the
#' Broyden-Fletcher-Goldfarb-Shanno method. However, due to the lack of
#' smoothness of this function on the \eqn{m} parameter (derivatives are zero
#' whenever \eqn{m} equals a sample observation), an exhaustive search must be
#' done by redoing the previous step in all intervals between two observations.
#' For a sample of \eqn{n} observations, this would lead to \eqn{n-1}
#' optimization problems. Given the computational cost of such procedure,
#' an interval search is used, where the optimization is repeated in the
#' intervals at most the value of the \emph{interv_step} from the last
#' minimum found. Details on the method are available on the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  headings and summary table
#' * 2  intermediate steps results
#' * 3  intermediate steps internals
#' * 4+  details of optim. routine
#' @param method int - the steps that should be used to estimate the
#' parameters.
#' * 0 no optimization perform - just return the log-likelihood from initial guess.
#' * 1 initial estimation based on method of moments
#' * 2 global optimization not considering lack of smoothness in m
#' * 3 interval optimization taking non-smoothness in m into consideration
#' @param interv_step  int - the number of intervals to be explored after
#' the last minimum was found in the interval optimization. Default is 10.
#' @param provided_m_ NumericVector - if NULL, the m parameter is estimated
#' by the routine. If numeric, the estimation fixes m to the given value.
#' @param par NumericVector - vector containing the initial guess for
#' parameters b, a and m, respectively. Default values of are c(2, 1, 0).
#' @param g_opt_par NumericVector - vector containing the global optimization
#' parameters.
#' The optimization parameters are:
#' * step  - (num) initial step size of the searching algorithm.
#' * tol   - (num) line search tolerance.
#' * iter  - (int) maximum number of iterations.
#' * eps   - (num) gradient tolerance. The stopping criteria is \eqn{||\text{gradient}||<\text{eps}}.
#' * msize - (num) simplex max size. stopping criteria given by \eqn{||\text{max edge}||<\text{msize}}
#' * algo  - (int) algorithm. the optimization method used:
#'   * 0 Fletcher-Reeves
#'   * 1 Polak-Ribiere
#'   * 2 Broyden-Fletcher-Goldfarb-Shanno
#'   * 3 Steepest descent
#'   * 4 Nelder-Mead simplex
#'   * 5 Broyden-Fletcher-Goldfarb-Shanno ver.2
#'
#' Details for each algorithm are available on the ['GSL' Manual](https://www.gnu.org/software/gsl/doc/html/multimin.html).
#' Default values are c(.1, 1e-2, 100, 1e-3, 1e-5, 3,0).
#' @param itv_opt_par NumericVector - interval optimization parameters. Fields
#' are the same as the ones for the global optimization. Default values
#' are c(.01, 1e-3, 200, 1e-3, 1e-5, 5, 0).
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#' * "matrix" - the covariance matrix for the parameters.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 2)
#' subbofit(sample_subbo)
#' @export
#' @md
subbofit <- function(data, verb = 0L, method = 3L, interv_step = 10L, provided_m_ = NULL, par = as.numeric( c(2.,1.,0.)), g_opt_par = as.numeric( c(.1, 1e-2, 100, 1e-3, 1e-5, 3)), itv_opt_par = as.numeric( c(.01, 1e-3, 200, 1e-3, 1e-5, 5))) {
    .Call(`_Rsubbotools_subbofit`, data, verb, method, interv_step, provided_m_, par, g_opt_par, itv_opt_par)
}

#' Fit a (Less) Asymmetric Power Exponential density via maximum likelihood
#'
#' \code{subbolafit} returns the parameters, standard errors. negative
#' log-likelihood and covariance matrix of the (less) asymmetric power exponential
#' for a sample. The main difference from \code{subboafit} is that
#' \eqn{a_l = a_r = a}. The process can execute two steps, dependending on the
#' level of accuracy required. See details below.
#'
#' The  LAPE is a exponential power distribution controlled
#' by four parameters, with formula:
#' \deqn{ f(x;a,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_l} |\frac{x-m}{a}|^{b_l} }, x < m
#' }
#' \deqn{ f(x;a,b_l,b_r,m) =
#' \frac{1}{A} e^{- \frac{1}{b_r} |\frac{x-m}{a}|^{b_r} }, x > m
#' }
#' with:
#' \deqn{A = ab_l^{1/b_l}\Gamma(1+1/b_l) + ab_r^{1/b_r}\Gamma(1+1/b_r)}
#' where \eqn{l} and \eqn{r} represent left and right tails, \eqn{a} is a
#' scale parameter, \eqn{b*} control the tails (lower values represent
#' fatter tails), and \eqn{m} is a location parameter. Due to its lack of
#' simmetry, and differently from the Subbotin, there is no simple equations
#' available to use the method of moments, so we start directly by minimizing
#' the negative log-likelihood. This global optimization is executed without
#' restricting any parameters. If required (default), after the global
#' optimization is finished, the method proceeds to iterate over the intervals
#' between several two observations, iterating the same algorithm of the
#' global optimization. The last method happens because of the lack of
#' smoothness on the \eqn{m} parameter, and intervals must be used since the
#' likelihood function doesn't have a derivative whenever \eqn{m} equals a
#' sample observation. Due to the cost, these iterations are capped at most
#' \emph{interv_step} (default 10) from the last minimum observed.
#' Details on the method are available on the package vignette.
#'
#' @param data (NumericVector) - the sample used to fit the distribution.
#' @param verb (int) - the level of verbosity. Select one of:
#' * 0  just the final result
#' * 1  headings and summary table
#' * 2  intermediate steps results
#' * 3  intermediate steps internals
#' * 4+  details of optim. routine
#' @param method int - the steps that should be used to estimate the
#' parameters.
#' * 0 no optimization perform - just return the log-likelihood from initial guess.
#' * 1 global optimization not considering lack of smoothness in m
#' * 2 interval optimization taking non-smoothness in m into consideration
#' @param interv_step  int - the number of intervals to be explored after
#' the last minimum was found in the interval optimization. Default is 10.
#' @param provided_m_ NumericVector - if NULL, the m parameter is estimated
#' by the routine. If numeric, the estimation fixes m to the given value.
#' @param par NumericVector - vector containing the initial guess for
#' parameters bl, br, a and m, respectively. Default values of are
#' c(2, 2, 1, 0).
#' @param g_opt_par NumericVector - vector containing the global optimization
#' parameters.
#' The optimization parameters are:
#' * step  - (num) initial step size of the searching algorithm.
#' * tol   - (num) line search tolerance.
#' * iter  - (int) maximum number of iterations.
#' * eps   - (num) gradient tolerance. The stopping criteria is \eqn{||\text{gradient}||<\text{eps}}.
#' * msize - (num) simplex max size. stopping criteria given by \eqn{||\text{max edge}||<\text{msize}}
#' * algo  - (int) algorithm. the optimization method used:
#'   * 0 Fletcher-Reeves
#'   * 1 Polak-Ribiere
#'   * 2 Broyden-Fletcher-Goldfarb-Shanno
#'   * 3 Steepest descent
#'   * 4 Nelder-Mead simplex
#'   * 5 Broyden-Fletcher-Goldfarb-Shanno ver.2
#'
#' Details for each algorithm are available on the ['GSL' Manual](https://www.gnu.org/software/gsl/doc/html/multimin.html).
#' Default values are c(.1, 1e-2, 100, 1e-3, 1e-5, 2).
#' @param itv_opt_par NumericVector - interval optimization parameters. Fields
#' are the same as the ones for the global optimization. Default values
#' are c(.01, 1e-3, 200, 1e-3, 1e-5, 2).
#' @return a list containing the following items:
#' * "dt" - dataset containing parameters estimations and standard deviations.
#' * "log-likelihood" - negative log-likelihood value.
#'
#' @examples
#' sample_subbo <- rpower(1000, 1, 2)
#' subbolafit(sample_subbo)
#' @export
#' @md
subbolafit <- function(data, verb = 0L, method = 2L, interv_step = 10L, provided_m_ = NULL, par = as.numeric( c(2., 2., 1., 0.)), g_opt_par = as.numeric( c(.1, 1e-2, 100, 1e-3, 1e-5, 2)), itv_opt_par = as.numeric( c(.01, 1e-3, 200, 1e-3, 1e-5, 2))) {
    .Call(`_Rsubbotools_subbolafit`, data, verb, method, interv_step, provided_m_, par, g_opt_par, itv_opt_par)
}