# R/QBAD.R In QBAsyDist: Asymmetric Distributions and Quantile Estimation

#' @title Quantile-based asymmetric family of distributions
#'
#' @description Density, cumulative distribution function, quantile function and random sample generation
#' from the quantile-based asymmetric family of densities defined in Gijbels et al. (2019a).
#'
#' @param y,q These are each a vector of quantiles.
#' @param mu This is the location parameter \eqn{\mu}.
#' @param phi This is the scale parameter  \eqn{\phi}.
#' @param alpha This is the index parameter  \eqn{\alpha}.
#'
#' @param f This is the reference density function \eqn{f} which is a standard version of a unimodal and symmetric around 0 density.
#'
#'
#'
#' @references{
#'  Gijbels, I., Karim, R. and Verhasselt, A. (2019a). On quantile-based asymmetric family of distributions: properties and inference. \emph{International Statistical Review}, \url{https://doi.org/10.1111/insr.12324}.
#' }
#'
#'
NULL
#'
#'
#' @export

#' @param F This is the cumulative distribution function \eqn{F} of a unimodal and symmetric around 0 reference density function \eqn{f}.
#'
#' @export

#' @param beta This is a vector of probabilities.
#' @param QF This is the quantile function of the reference density \eqn{f}.
#' @import GoFKernel
#' @export
if (is.null(QF)){
QF<-GoFKernel::inverse(F, lower = -Inf, upper = Inf)
qf<-NA;
for (i in 1:length(beta)) {
qf[i]<-ifelse(beta[i]<alpha,(mu+(phi/(1-alpha))*QF(beta[i]/(2*alpha))),
mu+(phi/alpha)*QF((1+beta[i]-2*alpha)/(2*(1-alpha))))}
} else {

qf<-ifelse(beta<alpha,(mu+(phi/(1-alpha))*QF(beta/(2*alpha))),
mu+(phi/alpha)*QF((1+beta-2*alpha)/(2*(1-alpha))))
}
return(qf)
}

#' @param n This is the number of observations, which must be a positive integer that has length 1.
#' @import stats
#' The length of the result is determined by \eqn{n} for \code{\link{rQBAD}}, and is the maximum of the lengths of the numerical arguments for the other functions.
#'
#'
#' @examples
#' # Example 1: Let F be a standard normal cumulative distribution function then
#' f_N<-function(s){dnorm(s, mean = 0,sd = 1)} # density function of N(0,1)
#' F_N<-function(s){pnorm(s, mean = 0,sd = 1)} # distribution function of N(0,1)
#' QF_N<-function(beta){qnorm(beta, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)}
#' rnum<-rnorm(100)
#' beta=c(0.25,0.50,0.75)
#'
#' # Density
#'
#' # Distribution function
#'
#' # Quantile function
#'
#' # random sample generation
#'
#'
#'
#' # Example 2: Let F be a standard Laplace cumulative distribution function then
#' f_La<-function(s){0.5*exp(-abs(s))} # density function of Laplace(0,1)
#' F_La<-function(s){0.5+0.5*sign(s)*(1-exp(-abs(s)))} # distribution function of Laplace(0,1)
#' QF_La<-function(beta){-sign(beta-0.5)*log(1-2*abs(beta-0.5))}
#' rnum<-rnorm(100)
#' beta=c(0.25,0.50,0.75)
#'
#' # Density
#'
#' # Distribution function
#'
#' # Quantile function
#'
#' # random sample generation
#'
#' @export
u<-runif(n, min = 0, max = 1)
if (is.null(QF)){
QF<-GoFKernel::inverse(F, lower = -Inf, upper = Inf)
r<-NA;
for (i in 1:length(u)) {
r[i]<-ifelse(u[i]<alpha,(mu+(phi/(1-alpha))*QF(u[i]/(2*alpha))),
(mu+(phi/alpha)*QF((1+u[i]-2*alpha)/(2*(1-alpha)))))}
} else {
r<-ifelse(u<alpha,(mu+(phi/(1-alpha))*QF(u/(2*alpha))),
(mu+(phi/alpha)*QF((1+u-2*alpha)/(2*(1-alpha)))))
}
return(r)
}

#' @title Log-likelihood function for the quantile-based asymmetric family of distributions.
#' @description Log-Likelihood function \eqn{\ell_n(\mu,\phi,\alpha)=\ln[L_n(\mu,\phi,\alpha)]}
#' in the three parameter quantile-based asymmetric family of densities defined in Section 3.2 of Gijbels et al. (2019a).
#' @param y This is a vector of quantiles.
#' @param mu This is the location parameter \eqn{\mu}.
#' @param phi This is the scale parameter  \eqn{\phi}.
#' @param alpha This is the index parameter  \eqn{\alpha}.
#'
#' @param f This is the reference density function \eqn{f} which is a standard version of a unimodal and symmetric around 0 density.
#' @return \code{\link{LogLikQBAD}} provides the realized value of the Log-likelihood function of quantile-based asymmetric family of distributions.
#'
#'
#' @references{
#'  Gijbels, I., Karim, R. and Verhasselt, A. (2019a). On quantile-based asymmetric family of distributions: properties and inference. \emph{International Statistical Review}, \url{https://doi.org/10.1111/insr.12324}.
#' }
#'
#'
#' @examples
#' # Example 1: Let F be a standard normal cumulative distribution function then
#' f_N<-function(s){dnorm(s, mean = 0,sd = 1)} # density function of N(0,1)
#' y<-rnorm(100)
#'
#' # Example 2: Let F be a standard Laplace cumulative distribution function then
#' f_La<-function(s){0.5*exp(-abs(s))} # density function of Laplace(0,1)
#' @export
LL<-(log(ifelse(y>mu,((2*alpha*(1-alpha)/phi)*f(alpha*(y-mu)/phi)), ((2*alpha*(1-alpha)/phi)*f((1-alpha)*(mu-y)/phi)))))
return(sum(LL[!is.infinite(LL)]))
}

#' @title Moment estimation for the quantile-based asymmetric family of distributions.
#' @description Mean, variance, skewness, kurtosis and moments about the location parameter (i.e., \eqn{\alpha}th quantile) of the quantile-based asymmetric family of densities defined in Gijbels et al. (2019a) useful for quantile regression with location parameter equal to \eqn{\mu}, scale parameter \eqn{\phi} and index parameter \eqn{\alpha}.
#' @param f This is the reference density function \eqn{f} which is a standard version of a unimodal and symmetric around 0 density.
#' @param r This is a value which is used to calculate the \eqn{r}th moment about \eqn{\mu}.
#' @param k This is an integer value (\eqn{k=1,2,3,\ldots}) for calculating \eqn{\mu_k=\int_{0}^{\infty} 2s^k f(s) ds} and \eqn{\gamma_k=\int_{0}^{\infty}s^{k-1}\frac{[f'(s)]^2}{f(s)}ds}.
#' @import stats
#'
#'
#'
#' @references{
#'  Gijbels, I., Karim, R. and Verhasselt, A. (2019a). On quantile-based asymmetric family of distributions: properties and inference. \emph{International Statistical Review}, \url{https://doi.org/10.1111/insr.12324}.
#' }
#'
#'
#'
NULL
#' @export
mu_k<-function(f,k){
integrand <- function(x) {x^k*f(x)}
return(2*stats::integrate(integrand, lower = 0, upper = Inf)$value)} #' @import Deriv #' @rdname momentQBAD #' @export gamma_k<-function(f,k){ fn <- function(x) {1/f(x)} df<- Deriv::Deriv(f) gamma_fn<-function(x){x^(k-1)*(df(x))^2*(fn(x))} return(stats::integrate(gamma_fn, lower = 0, upper = 30, stop.on.error=FALSE)$ value)
}

#' @param mu This is the location parameter \eqn{\mu}.
#' @param phi This is the scale parameter  \eqn{\phi}.
#' @param alpha This is the index parameter  \eqn{\alpha}.
#'
#' @param mu_1 This is the quantity \eqn{\int_{0}^{\infty} 2 s f(s) ds}.
#' @export

#' @param mu_2 This is the quantity \eqn{\int_{0}^{\infty} 2 s^2 f(s) ds}.
#' @export

#' @param mu_3 This is the quantity \eqn{\int_{0}^{\infty} 2 s^3 f(s) ds}.
#' @export
nu<-(1-2*alpha)*((1-2*alpha)^2*(mu_3-3*mu_1*mu_2+2*mu_1^3)+alpha*(1-alpha)*(2*mu_3-3*mu_1*mu_2))
deno<-((1-2*alpha)^2*(mu_2-mu_1^2)+alpha*(1-alpha)*mu_2)^(3/2)
skew<-(nu/deno)
return(skew)
}

#' @param mu_4 This is the quantity \eqn{\int_{0}^{\infty} 2 s^4 f(s) ds}.
#'
#' @export
nu<-((1-alpha)^5+alpha^5)*mu_4-(1-2*alpha)^2*(4*(1-2*alpha+2*alpha^2)*mu_1*mu_3-6*(1-3*alpha+3*alpha^2)*mu_1^2*mu_2+3*(1-2*alpha)^2*mu_1^4)
deno<-((1-2*alpha)^2*(mu_2-mu_1^2)+alpha*(1-alpha)*mu_2)^(2)
kurto<-(nu/deno)
return(kurto)
}

#' @import stats
#' @examples
#' # Example 1: Let F be a standard normal cumulative distribution function then
#' f_N<-function(s){dnorm(s, mean = 0,sd = 1)} # density function of N(0,1)
#' mu_k(f=f_N,k=1)
#' gamma_k(f=f_N,k=1)
#' mu.1_N=sqrt(2/pi)
#' mu.2_N=1
#' mu.3_N=2*sqrt(2/pi)
#' mu.4_N=4
#'
#'
#' # Example 2: Let F be a standard Laplace cumulative distribution function then
#' f_La<-function(s){0.5*exp(-abs(s))} # density function of Laplace(0,1)
#' mu_k(f=f_La,k=1)
#' gamma_k(f=f_La,k=1)
#' mu.1_La=1
#' mu.2_La=2
#' mu.3_La=6
#' mu.4_La=24
#' @export
integrand <- function(x) {x^r*f(x)}
mu.r<-2*stats::integrate(integrand, lower = 0, upper = Inf)$value return(phi^r*((1-alpha)^(r+1)+(-1)^r*alpha^(r+1))*mu.r/(alpha^r*(1-alpha)^r)) } #' @title Method of moments (MoM) estimation for the quantile-based asymmetric family of distributions. #' @description Parameter estimation in the quantile-based asymmetric family of densities by using method of moments are discussed in Section 3.1 of Gijbels et al. (2019a). #' @param y This is a vector of quantiles. #' @param alpha This is the index parameter \eqn{\alpha}. If \eqn{\alpha} is unknown, indicate NULL which is default option. In this case, the sample skewness will be used to estimate \eqn{\alpha}. If \eqn{\alpha} is known, then the value of \eqn{\alpha} has to be specified in the function. #' #' #' @param f This is the reference density function \eqn{f} which is a standard version of a unimodal and symmetric around 0 density. #' #' @import stats #' @return \code{\link{momQBAD}} provides the method of moments estimates of the unknown parameters of the distribution. #' #' #' #' #' @references{ #' Gijbels, I., Karim, R. and Verhasselt, A. (2019a). On quantile-based asymmetric family of distributions: properties and inference. \emph{International Statistical Review}, \url{https://doi.org/10.1111/insr.12324}. #' } #' #' #' #' @name momQBAD NULL #' @rdname momQBAD #' @examples #' # Example 1: Let F be a standard normal cumulative distribution function then #' f_N<-function(s){dnorm(s, mean = 0,sd = 1)} # density function of N(0,1) #' y=rnorm(100) #' momQBAD(y=y,f=f_N,alpha=0.5) # If alpha is known with alpha=0.5 #' momQBAD(y=y,f=f_N) # If alpha is unknown #' #' #' # Example 2: Let F be a standard Laplace cumulative distribution function then #' f_La<-function(s){0.5*exp(-abs(s))} # density function of Laplace(0,1) #' momQBAD(y=y,f=f_La,alpha=0.5) # If alpha is known with alpha=0.5 #' momQBAD(y=y,f=f_La) # If alpha is unknown #' @export momQBAD<-function(y,f,alpha=NULL){ muMoM<-function(y,f,alpha=NULL){ if (is.null(alpha)){ m1<-mean(y) m2<-sum(y^2)/length(y) cm2<-sum((y-m1)^2)/length(y) cm3<-sum((y-m1)^3)/length(y) cm4<-sum((y-m1)^4)/length(y) integrand1 <- function(x) {x*f(x)} integrand2 <- function(x) {x^2*f(x)} integrand3 <- function(x) {x^3*f(x)} mu1<-2*stats::integrate(integrand1, lower = 0, upper = Inf)$ value
mu2<-2*stats::integrate(integrand2, lower = 0, upper = Inf)$value mu3<-2*stats::integrate(integrand3, lower = 0, upper = Inf)$ value
diffskew<-matrix(NA,99,2)
for (al.i in 1:99){
al=al.i/100
nu<-(1-2*al)*((1-2*al)^2*(mu3-3*mu1*mu2+2*mu1^3)+al*(1-al)*(2*mu3-3*mu1*mu2))
deno<-((1-2*al)^2*(mu2-mu1^2)+al*(1-al)*mu2)^(3/2)
skewPop<-(nu/deno)
diffskew[al.i,1]<-abs(skewPop- cm3/cm2^(3/2))
diffskew[al.i,2]<-al}
alpha<-diffskew[diffskew[,1]== min( diffskew[,1]), ][2]
k1<-(1-2*alpha)*mu1/(alpha*(1-alpha))
k2<-((1-alpha)^3+alpha^3)*mu2/(alpha^2*(1-alpha)^2)
mu<-m1-(k1/sqrt(k2-k1^2))*sqrt(m2-m1^2)
} else {
m1<-mean(y)
m2<-sum(y^2)/length(y)
integrand1 <- function(x) {x*f(x)}
integrand2 <- function(x) {x^2*f(x)}
mu1<-2*stats::integrate(integrand1, lower = 0, upper = Inf)$value mu2<-2*stats::integrate(integrand2, lower = 0, upper = Inf)$ value

k1<-(1-2*alpha)*mu1/(alpha*(1-alpha))
k2<-((1-alpha)^3+alpha^3)*mu2/(alpha^2*(1-alpha)^2)
mu<-m1-(k1/sqrt(k2-k1^2))*sqrt(m2-m1^2)
}
return(mu)
}
phiMoM<-function(y,f,alpha=NULL){
if (is.null(alpha)){
m1<-mean(y)
m2<-sum(y^2)/length(y)
cm2<-sum((y-m1)^2)/length(y)
cm3<-sum((y-m1)^3)/length(y)
integrand1 <- function(x) {x*f(x)}
integrand2 <- function(x) {x^2*f(x)}
integrand3 <- function(x) {x^3*f(x)}
mu1<-2*stats::integrate(integrand1, lower = 0, upper = Inf)$value mu2<-2*stats::integrate(integrand2, lower = 0, upper = Inf)$ value
mu3<-2*stats::integrate(integrand3, lower = 0, upper = Inf)$value diffskew<-matrix(NA,99,2) for (al.i in 1:99){ al=al.i/100 nu<-(1-2*al)*((1-2*al)^2*(mu3-3*mu1*mu2+2*mu1^3)+al*(1-al)*(2*mu3-3*mu1*mu2)) deno<-((1-2*al)^2*(mu2-mu1^2)+al*(1-al)*mu2)^(3/2) skewPop<-(nu/deno) diffskew[al.i,1]<-abs(skewPop- cm3/cm2^(3/2)) diffskew[al.i,2]<-al} alpha<-diffskew[diffskew[,1]== min( diffskew[,1]), ][2] k1<-(1-2*alpha)*mu1/(alpha*(1-alpha)) k2<-((1-alpha)^3+alpha^3)*mu2/(alpha^2*(1-alpha)^2) phi<-(1/sqrt(k2-k1^2))*sqrt(m2-m1^2) } else { m1<-mean(y) m2<-sum(y^2)/length(y) integrand1 <- function(x) {x*f(x)} integrand2 <- function(x) {x^2*f(x)} mu1<-2*stats::integrate(integrand1, lower = 0, upper = Inf)$ value
mu2<-2*stats::integrate(integrand2, lower = 0, upper = Inf)$value k1<-(1-2*alpha)*mu1/(alpha*(1-alpha)) k2<-((1-alpha)^3+alpha^3)*mu2/(alpha^2*(1-alpha)^2) phi<-(1/sqrt(k2-k1^2))*sqrt(m2-m1^2) } return(phi) } alphaMoM<-function(y,f){ m1<-mean(y) cm2<-sum((y-m1)^2)/length(y) cm3<-sum((y-m1)^3)/length(y) integrand1 <- function(x) {x*f(x)} integrand2 <- function(x) {x^2*f(x)} integrand3 <- function(x) {x^3*f(x)} mu1<-2*stats::integrate(integrand1, lower = 0, upper = Inf)$ value
mu2<-2*stats::integrate(integrand2, lower = 0, upper = Inf)$value mu3<-2*stats::integrate(integrand3, lower = 0, upper = Inf)$ value
diffskew<-matrix(NA,99,2)
for (al.i in 1:99){
al=al.i/100
nu<-(1-2*al)*((1-2*al)^2*(mu3-3*mu1*mu2+2*mu1^3)+al*(1-al)*(2*mu3-3*mu1*mu2))
deno<-((1-2*al)^2*(mu2-mu1^2)+al*(1-al)*mu2)^(3/2)
skewPop<-(nu/deno)
diffskew[al.i,1]<-abs(skewPop- cm3/cm2^(3/2))
diffskew[al.i,2]<-al
}
return(diffskew[diffskew[,1]== min( diffskew[,1]), ][2])
}

if (is.null(alpha)){
list(mu.MoM=muMoM(y,f),phi.MoM=phiMoM(y,f),alpha.MoM=alphaMoM(y,f))
} else {
list(mu.MoM=muMoM(y,f),phi.MoM=phiMoM(y,f))}
}

#' @title Maximum likelihood estimation (MLE) for the quantile-based asymmetric family of distributions.
#' @description The log-likelihood function \eqn{\ell_n(\mu,\phi,\alpha)=\ln[L_n(\mu,\phi,\alpha)]}
#' and parameter estimation of \eqn{ \theta=(\mu,\phi,\alpha)} in the three parameter quantile-based asymmetric family of densities
#' by using the maximum likelihood estimation are discussed in Section 3.2 of Gijbels et al. (2019a).
#' @param y This is a vector of quantiles.
#' @param f This is the reference density function \eqn{f} which is a standard version of a unimodal and symmetric around 0 density.
#' @param alpha This is the index parameter  \eqn{\alpha}.
#'
#' @return The maximum likehood estimate of paramter \eqn{\theta=(\mu,\phi,\alpha)} of the quantile-based asymmetric family of densities
#'
#' @import ald
#' @import nloptr
#'
#'
#' @references{
#'  Gijbels, I., Karim, R. and Verhasselt, A. (2019a). On quantile-based asymmetric family of distributions: properties and inference. \emph{International Statistical Review}, \url{https://doi.org/10.1111/insr.12324}.
#' }
#'
#'
#'
NULL
#' @examples
#' # Example 1: Let F be a standard normal cumulative distribution function then
#' f_N<-function(s){dnorm(s, mean = 0,sd = 1)} # density function of N(0,1)
#' rnum=rnorm(100)
#'
#' # Example 2: Let F be a standard Laplace cumulative distribution function then
#' f_La<-function(s){0.5*exp(-abs(s))} # density function of Laplace(0,1)
#' @export
if (is.null(alpha)){
theta0<-ald::mleALD(y, initial = NA)$par fn <- function(theta){ mu=theta[1] phi=theta[2] alpha=theta[3] LL<-log(dQBAD(y,mu,phi,alpha,f)) return(-sum(LL[!is.infinite(LL)]))} Est.par=nloptr::bobyqa(x0=theta0,fn = fn ,lower = c(min(y),00.001,0.001), upper = c(max(y),Inf,1))$par
list(mu.MLE=Est.par[1], phi.MLE=Est.par[2],alpha.MLE=Est.par[3],LogLikelihood=estimated.LL)
} else {
theta0<-c(median(y),sd(y),alpha)
fn <- function(theta){
mu=theta[1]
phi=theta[2]
alpha=alpha