# R/VaRES.R In VaRES: Computes value at risk and expected shortfall for over 100 parametric distributions

```#exponential

dexponential=function (x, lambda=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=lambda*exp(-lambda*x)
pdf[log==TRUE]=log(lambda)-lambda*x
return(pdf)
}

pexponential=function (x, lambda=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-exp(-lambda*x)
cdf[log.p==FALSE&lower.tail==FALSE]=exp(-lambda*x)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-exp(-lambda*x))
cdf[log.p==TRUE&lower.tail==FALSE]=-lambda*x
return(cdf)
}

varexponential=function (p, lambda=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(-1/lambda)*log(1-p)
return(var)
}

esexponential=function (p, lambda=1)
{
f=function (x) {varexponential(x,lambda=lambda)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy exponential

dkumexp=function (x, lambda=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*b*lambda*exp(-lambda*x)*(1-exp(-lambda*x))**(a-1)*(1-(1-exp(-lambda*x))**a)**(b-1)
pdf[log==TRUE]=log(a*b*lambda)-lambda*x+(a-1)*log(1-exp(-lambda*x))+(b-1)*log(1-(1-exp(-lambda*x))**a)
return(pdf)
}

pkumexp=function (x, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-(1-exp(-lambda*x))**a)**b
cdf[log.p==FALSE&lower.tail==FALSE]=(1-(1-exp(-lambda*x))**a)**b
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-(1-exp(-lambda*x))**a)**b)
cdf[log.p==TRUE&lower.tail==FALSE]=b*log(1-(1-exp(-lambda*x))**a)
return(cdf)
}

varkumexp=function (p, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(-1/lambda)*log(1-(1-(1-p)**(1/b))**(1/a))
return(var)
}

eskumexp=function (p, lambda=1, a=1, b=1)
{
f=function (x) {varkumexp(x,lambda=lambda,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Exponentiated exponential

dexpexp=function (x, lambda=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*lambda*exp(-lambda*x)*(1-exp(-lambda*x))**(a-1)
pdf[log==TRUE]=log(a*lambda)-lambda*x+(a-1)*log(1-exp(-lambda*x))
return(pdf)
}

pexpexp=function (x, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=(1-exp(-lambda*x))**a
cdf[log.p==FALSE&lower.tail==FALSE]=1-(1-exp(-lambda*x))**a
cdf[log.p==TRUE&lower.tail==TRUE]=a*log(1-exp(-lambda*x))
cdf[log.p==TRUE&lower.tail==FALSE]=log(1-(1-exp(-lambda*x))**a)
return(cdf)
}

varexpexp=function (p, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(-1/lambda)*log(1-p**(1/a))
return(var)
}

esexpexp=function (p, lambda=1, a=1)
{
f=function (x) {varexpexp(x,lambda=lambda,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Inverse exponentiated exponential

dinvexpexp=function (x, lambda=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*lambda*(1/(x*x))*exp(-lambda/x)*(1-exp(-lambda/x))**(a-1)
pdf[log==TRUE]=log(a*lambda)-2*log(x)-lambda/x+(a-1)*log(1-exp(-lambda/x))
return(pdf)
}

pinvexpexp=function (x, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-exp(-lambda/x))**a
cdf[log.p==FALSE&lower.tail==FALSE]=(1-exp(-lambda/x))**a
cdf[log.p==TRUE&lower.tail==TRUE]=a*log(1-(1-exp(-lambda/x))**a)
cdf[log.p==TRUE&lower.tail==FALSE]=a*log(1-exp(-lambda/x))
return(cdf)
}

varinvexpexp=function (p, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=lambda/(-log(1-(1-p)**(1/a)))
return(var)
}

esinvexpexp=function (p, lambda=1, a=1)
{
f=function (x) {varinvexpexp(x,lambda=lambda,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Beta exponential

dbetaexp=function (x, lambda=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=lambda*exp(-b*lambda*x)*(1-exp(-lambda*x))**(a-1)/beta(a,b)
pdf[log==TRUE]=log(lambda)-b*lambda*x+(a-1)*log(1-exp(-lambda*x))-lbeta(a,b)
return(pdf)
}

pbetaexp=function (x, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(1-exp(-lambda*x),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varbetaexp=function (p, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(-1/lambda)*log(1-qbeta(p,shape1=a,shape2=b))
return(var)
}

esbetaexp=function (p, lambda=1, a=1, b=1)
{
f=function (x) {varbetaexp(x,lambda=lambda,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Logistic exponential

dlogisexp=function (x, lambda=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*lambda*exp(lambda*x)*(exp(lambda*x)-1)**(a-1)/(1+(exp(lambda*x)-1)**a)**2
pdf[log==TRUE]=log(a*lambda)+lambda*x+(a-1)*log(exp(lambda*x)-1)-2*log(1+(exp(lambda*x)-1)**a)
return(pdf)
}

plogisexp=function (x, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=(exp(lambda*x)-1)**a/(1+(exp(lambda*x)-1)**a)
cdf[log.p==FALSE&lower.tail==FALSE]=1/(1+(exp(lambda*x)-1)**a)
cdf[log.p==TRUE&lower.tail==TRUE]=a*log(exp(lambda*x)-1)-log(1+(exp(lambda*x)-1)**a)
cdf[log.p==TRUE&lower.tail==FALSE]=-log(1+(exp(lambda*x)-1)**a)
return(cdf)
}

varlogisexp=function (p, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/lambda)*log(1+(p/(1-p))**(1/a))
return(var)
}

eslogisexp=function (p, lambda=1, a=1)
{
f=function (x) {varlogisexp(x,lambda=lambda,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Exponential extension

dexpext=function (x, lambda=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*lambda*(1+lambda*x)**(a-1)*exp(1-(1+lambda*x)**a)
pdf[log==TRUE]=log(a*lambda)+(a-1)*log(1+lambda*x)+1-(1+lambda*x)**a
return(pdf)
}

pexpext=function (x, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-exp(1-(1+lambda*x)**a)
cdf[log.p==FALSE&lower.tail==FALSE]=exp(1-(1+lambda*x)**a)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-exp(1-(1+lambda*x)**a))
cdf[log.p==TRUE&lower.tail==FALSE]=1-(1+lambda*x)**a
return(cdf)
}

varexpext=function (p, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/lambda)*((1-log(1-p))**(1/a)-1)
return(var)
}

esexpext=function (p, lambda=1, a=1)
{
f=function (x) {varexpext(x,lambda=lambda,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Marshall Olkin exponential

dmoexp=function (x, lambda=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=lambda*exp(lambda*x)/(exp(lambda*x)-1+a)**2
pdf[log==TRUE]=log(lambda)+lambda*x-2*log(exp(lambda*x)-1+a)
return(pdf)
}

pmoexp=function (x, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=(exp(lambda*x)-2+a)/(exp(lambda*x)-1+a)
cdf[log.p==FALSE&lower.tail==FALSE]=1/(exp(lambda*x)-1+a)
cdf[log.p==TRUE&lower.tail==TRUE]=log(exp(lambda*x)-2+a)-log(exp(lambda*x)-1+a)
cdf[log.p==TRUE&lower.tail==FALSE]=-log(exp(lambda*x)-1+a)
return(cdf)
}

varmoexp=function (p, lambda=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/lambda)*(log(2-a-(1-a)*p)-log(1-p))
return(var)
}

esmoexp=function (p, lambda=1, a=1)
{
f=function (x) {varmoexp(x,lambda=lambda,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Perks

dperks=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*(a+1)*exp(b*x)/(1+a*exp(b*x))**2
pdf[log==TRUE]=log(a*(a+1))+b*x-2*log(1+a*exp(b*x))
return(pdf)
}

pperks=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=a*(exp(b*x)-1)/(1+a*exp(b*x))
cdf[log.p==FALSE&lower.tail==FALSE]=(1+a)/(1+a*exp(b*x))
cdf[log.p==TRUE&lower.tail==TRUE]=log(a)+log(exp(b*x)-1)-log(1+a*exp(b*x))
cdf[log.p==TRUE&lower.tail==FALSE]=log(1+a)-log(1+a*exp(b*x))
return(cdf)
}

varperks=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*(log(a+p)-log(a)-log(1-p))
return(var)
}

esperks=function (p, a=1, b=1)
{
f=function (x) {varperks(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Beard

dbeard=function (x, a=1, b=1, rho=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*exp(b*x)*(1+a*rho)**(rho**(-1/b))/(1+a*rho*exp(b*x))**(1+rho**(-1/b))
pdf[log==TRUE]=log(a)+b*x+(rho**(-1/b))*log(1+a*rho)-(1+rho**(-1/b))*log(1+a*rho*exp(b*x))
return(pdf)
}

pbeard=function (x, a=1, b=1, rho=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-((1+a*rho)/(1+a*rho*exp(b*x)))**(rho**(-1/b))
cdf[log.p==FALSE&lower.tail==FALSE]=((1+a*rho)/(1+a*rho*exp(b*x)))**(rho**(-1/b))
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-((1+a*rho)/(1+a*rho*exp(b*x)))**(rho**(-1/b)))
cdf[log.p==TRUE&lower.tail==FALSE]=(rho**(-1/b))*(log(1+a*rho)-log(1+a*rho*exp(b*x)))
return(cdf)
}

varbeard=function (p, a=1, b=1, rho=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*(log((1+a*rho)/(1-p)**(rho**(1/b))-1)-log(a*rho))
return(var)
}

esbeard=function (p, a=1, b=1, rho=1)
{
f=function (x) {varbeard(x,a=a,b=b,rho=rho)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Gompertz

dgompertz=function (x, b=1, eta=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=b*eta*exp(b*x)*exp(eta-eta*exp(b*x))
pdf[log==TRUE]=log(b*eta)+b*x+eta-eta*exp(b*x)
return(pdf)
}

pgompertz=function (x, b=1, eta=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-exp(eta-eta*exp(b*x))
cdf[log.p==FALSE&lower.tail==FALSE]=exp(eta-eta*exp(b*x))
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-exp(eta-eta*exp(b*x)))
cdf[log.p==TRUE&lower.tail==FALSE]=eta-eta*exp(b*x)
return(cdf)
}

vargompertz=function (p, b=1, eta=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*log(1-(1/eta)*log(1-p))
return(var)
}

esgompertz=function (p, b=1, eta=1)
{
f=function (x) {vargompertz(x,b=b,eta=eta)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Beta Gompertz

dbetagompertz=function (x, b=1, c=1, d=1, eta=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=b*eta*exp(d*eta-d*eta*exp(b*x)+b*x)*(1-exp(eta-eta*exp(b*x)))**(c-1)
pdf[log==TRUE]=log(b*eta)+d*eta-d*eta*exp(b*x)+b*x+(c-1)*log(1-exp(eta-eta*exp(b*x)))
return(pdf)
}

pbetagompertz=function (x, b=1, c=1, d=1, eta=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(1-exp(eta-eta*exp(b*x)),shape1=c,shape2=d,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varbetagompertz=function (p, b=1, c=1, d=1, eta=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*log(1-(1/eta)*log(1-qbeta(p,shape1=c,shape2=d)))
return(var)
}

esbetagompertz=function (p, b=1, c=1, d=1, eta=1)
{
f=function (x) {varbetagompertz(x,b=b,c=c,d=d,eta=eta)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Linear failure rate

dlfr=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(a+b*x)*exp(-a*x-b*x*x/2)
pdf[log==TRUE]=log(a+b*x)-a*x-b*x*x/2
return(pdf)
}

plfr=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-exp(-a*x-b*x*x/2)
cdf[log.p==FALSE&lower.tail==FALSE]=exp(-a*x-b*x*x/2)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-exp(-a*x-b*x*x/2))
cdf[log.p==TRUE&lower.tail==FALSE]=-a*x-b*x*x/2
return(cdf)
}

varlfr=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*(sqrt(a*a-2*b*log(1-p))-a)
return(var)
}

eslfr=function (p, a=1, b=1)
{
f=function (x) {varlfr(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Pareto

dpareto=function (x, K=1, c=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=c*K**c*x**(-c-1)
pdf[log==TRUE]=log(c)+c*log(K)-(c+1)*log(x)
return(pdf)
}

ppareto=function (x, K=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-K**c*x**(-c)
cdf[log.p==FALSE&lower.tail==FALSE]=K**c*x**(-c)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-K**c*x**(-c))
cdf[log.p==TRUE&lower.tail==FALSE]=c*(log(K)-log(x))
return(cdf)
}

varpareto=function (p, K=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=K*(1-p)**(-1/c)
return(var)
}

espareto=function (p, K=1, c=1)
{
f=function (x) {varpareto(x,K=K,c=c)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy Pareto

dkumpareto=function (x, K=1, a=1, b=1, c=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*b*c*K**c*x**(-c-1)*(1-K**c*x**(-c))**(a-1)*(1-(1-K**c*x**(-c))**a)**(b-1)
pdf[log==TRUE]=log(a*b*c)+c*log(K)-(c+1)*log(x)+(a-1)*log(1-K**c*x**(-c))+(b-1)*log(1-(1-K**c*x**(-c))**a)
return(pdf)
}

pkumpareto=function (x, K=1, a=1, b=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-(1-K**c*x**(-c))**a)**b
cdf[log.p==FALSE&lower.tail==FALSE]=(1-(1-K**c*x**(-c))**a)**b
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-(1-K**c*x**(-c))**a)**b)
cdf[log.p==TRUE&lower.tail==FALSE]=b*log(1-(1-K**c*x**(-c))**a)
return(cdf)
}

varkumpareto=function (p, K=1, a=1, b=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=K*(1-(1-(1-p)**(1/b))**(1/a))**(-1/c)
return(var)
}

eskumpareto=function (p, K=1, a=1, b=1, c=1)
{
f=function (x) {varkumpareto(x,K=K,a=a,b=b,c=c)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#F

dF=function (x, d1=1, d2=1, log=FALSE)
{
pdf=df(x,df1=d1,df2=d2,log=log)
return(pdf)
}

pF=function (x, d1=1, d2=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pf(x,df1=d1,df2=d2,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varF=function (p, d1=1, d2=1, log.p=FALSE, lower.tail=TRUE)
{
var=qf(p,df1=d1,df2=d2,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esF=function (p, d1=1, d2=1)
{
f=function (x) {varF(x,d1=d1,d2=d2)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized Pareto

dgenpareto=function (x, k=1, c=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/k)*(1-c*x/k)**(1/c-1)
pdf[log==TRUE]=-log(k)+(1/c-1)*log(1-c*x/k)
return(pdf)
}

pgenpareto=function (x, k=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-c*x/k)**(1/c)
cdf[log.p==FALSE&lower.tail==FALSE]=(1-c*x/k)**(1/c)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-c*x/k)**(1/c))
cdf[log.p==TRUE&lower.tail==FALSE]=(1/c)*log(1-c*x/k)
return(cdf)
}

vargenpareto=function (p, k=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(k/c)*(1-(1-p)**c)
return(var)
}

esgenpareto=function (p, k=1, c=1)
{
f=function (x) {vargenpareto(x,k=k,c=c)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Beta Pareto

dbetapareto=function (x, K=1, a=1, c=1, d=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*K**(a*d)*x**(-a*d-1)*(1-(K/x)**a)**(c-1)/beta(c,d)
pdf[log==TRUE]=log(a)+(a*d)*log(K)-(a*d+1)*log(x)+(c-1)*log(1-(K/x)**a)-lbeta(c,d)
return(pdf)
}

pbetapareto=function (x, K=1, a=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(1-(K/x)**a,shape1=c,shape2=d,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varbetapareto=function (p, K=1, a=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=K*(1-qbeta(p,shape1=c,shape2=d))**(-1/a)
return(var)
}

esbetapareto=function (p, K=1, a=1, c=1, d=1)
{
f=function (x) {varbetapareto(x,K=K,a=a,c=c,d=d)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Pareto positive stable

dparetostable=function (x, lambda=1, nu=1, sigma=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=nu*lambda*(1/x)*(log(x/sigma))**(nu-1)*exp(-lambda*(log(x/sigma))**nu)
pdf[log==TRUE]=log(nu*lambda)-log(x)+(nu-1)*log(log(x/sigma))-lambda*(log(x/sigma))**nu
return(pdf)
}

pparetostable=function (x, lambda=1, nu=1, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-exp(-lambda*(log(x/sigma))**nu)
cdf[log.p==FALSE&lower.tail==FALSE]=exp(-lambda*(log(x/sigma))**nu)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-exp(-lambda*(log(x/sigma))**nu))
cdf[log.p==TRUE&lower.tail==FALSE]=-lambda*(log(x/sigma))**nu
return(cdf)
}

varparetostable=function (p, lambda=1, nu=1, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=sigma*exp(((-1/lambda)*log(1-p))**(1/nu))
return(var)
}

esparetostable=function (p, lambda=1, nu=1, sigma=1)
{
f=function (x) {varparetostable(x,lambda=lambda,nu=nu,sigma=sigma)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Gamma

dGamma=function (x, a=1, b=1, log=FALSE)
{
pdf=dgamma(x,shape=a,rate=b,log=log)
return(pdf)
}

pGamma=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pgamma(x,shape=a,rate=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varGamma=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
var=qgamma(p,shape=a,rate=b,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esGamma=function (p, a=1, b=1)
{
f=function (x) {varGamma(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy gamma

dkumgamma=function (x, a=1, b=1, c=1, d=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=c*d*dgamma(x,shape=a,rate=b)*(pgamma(x,shape=a,rate=b))**(c-1)*(1-(pgamma(x,shape=a,rate=b))**c)**(d-1)
pdf[log==TRUE]=log(c*d)+dgamma(x,shape=a,rate=b,log=TRUE)+(c-1)*pgamma(x,shape=a,rate=b,log.p=TRUE)+(d-1)*log(1-(pgamma(x,shape=a,rate=b))**c)
return(pdf)
}

pkumgamma=function (x, a=1, b=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-(pgamma(x,shape=a,rate=b))**c)**d
cdf[log.p==FALSE&lower.tail==FALSE]=(1-(pgamma(x,shape=a,rate=b))**c)**d
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-(pgamma(x,shape=a,rate=b))**c)**d)
cdf[log.p==TRUE&lower.tail==FALSE]=d*log(1-(pgamma(x,shape=a,rate=b))**c)
return(cdf)
}

varkumgamma=function (p, a=1, b=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1/b)*qgamma((1-(1-p)**(1/d))**(1/c),shape=a)
return(var)
}

eskumgamma=function (p, a=1, b=1, c=1, d=1)
{
f=function (x) {varkumgamma(x,a=a,b=b,c=c,d=d)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Nakagami

dnakagami=function (x, m=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=2*(m/a)**m*x**(2*m-1)*exp(-m*x*x/a)/gamma(m)
pdf[log==TRUE]=log(2)+m*log(m/a)+(2*m-1)*log(x)-m*x*x/a-lgamma(m)
return(pdf)
}

pnakagami=function (x, m=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=pgamma(m*x*x/a,shape=m,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varnakagami=function (p, m=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=sqrt(a/m)*sqrt(qgamma(p,shape=m))
return(var)
}

esnakagami=function (p, m=1, a=1)
{
f=function (x) {varnakagami(x,m=m,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Reflected gamma

drgamma=function (x, a=1, theta=0, phi=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/(2*phi))*dgamma(abs(x-theta)/phi,shape=a)
pdf[log==TRUE]=dgamma(abs(x-theta)/phi,shape=a,log=TRUE)-log(2*phi)
return(pdf)
}

prgamma=function (x, a=1, theta=0, phi=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[x<=theta&log.p==FALSE&lower.tail==TRUE]=0.5-0.5*pgamma((theta-x)/phi,shape=a)
cdf[x>theta&log.p==FALSE&lower.tail==TRUE]=0.5+0.5*pgamma((x-theta)/phi,shape=a)
cdf[x<=theta&log.p==FALSE&lower.tail==FALSE]=0.5+0.5*pgamma((theta-x)/phi,shape=a)
cdf[x>theta&log.p==FALSE&lower.tail==FALSE]=0.5-0.5*pgamma((x-theta)/phi,shape=a)
cdf[x<=theta&log.p==TRUE&lower.tail==TRUE]=log(0.5-0.5*pgamma((theta-x)/phi,shape=a))
cdf[x>theta&log.p==TRUE&lower.tail==TRUE]=log(0.5+0.5*pgamma((x-theta)/phi,shape=a))
cdf[x<=theta&log.p==TRUE&lower.tail==FALSE]=log(0.5+0.5*pgamma((theta-x)/phi,shape=a))
cdf[x>theta&log.p==TRUE&lower.tail==FALSE]=log(0.5-0.5*pgamma((x-theta)/phi,shape=a))
return(cdf)
}

varrgamma=function (p, a=1, theta=0, phi=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<=0.5]=theta-phi*qgamma(1-2*p[p<=0.5],shape=a)
var[p>0.5]=theta+phi*qgamma(2*p[p>0.5]-1,shape=a)
return(var)
}

esrgamma=function (p, a=1, theta=0, phi=1)
{
f=function (x) {varrgamma(x,a=a,theta=theta,phi=phi)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Compound Laplace gamma

dclg=function (x, a=1, b=1, theta=0, log=FALSE)
{
pdf=x
pdf[log==FALSE]=0.5*a*b*(1+b*abs(x-theta))**(-a-1)
pdf[log==TRUE]=log(a*b)-log(2)-(a+1)*log(1+b*abs(x-theta))
return(pdf)
}

pclg=function (x, a=1, b=1, theta=0, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[x<=theta&log.p==FALSE&lower.tail==TRUE]=0.5*(1+b*abs(x-theta))**(-a)
cdf[x>theta&log.p==FALSE&lower.tail==TRUE]=1-0.5*(1+b*abs(x-theta))**(-a)
cdf[x<=theta&log.p==FALSE&lower.tail==FALSE]=1-0.5*(1+b*abs(x-theta))**(-a)
cdf[x>theta&log.p==FALSE&lower.tail==FALSE]=0.5*(1+b*abs(x-theta))**(-a)
cdf[x<=theta&log.p==TRUE&lower.tail==TRUE]=-log(2)-a*log(1+b*abs(x-theta))
cdf[x>theta&log.p==TRUE&lower.tail==TRUE]=log(1-0.5*(1+b*abs(x-theta))**(-a))
cdf[x<=theta&log.p==TRUE&lower.tail==FALSE]=log(1-0.5*(1+b*abs(x-theta))**(-a))
cdf[x>theta&log.p==TRUE&lower.tail==FALSE]=-log(2)-a*log(1+b*abs(x-theta))
return(cdf)
}

varclg=function (p, a=1, b=1, theta=0, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<=0.5]=theta-1/b-(1/b)*(2*p[p<=0.5])**(-1/a)
var[p>0.5]=theta-1/b+(1/b)*(2*(1-p[p>0.5]))**(-1/a)
return(var)
}

esclg=function (p, a=1, b=1, theta=0)
{
f=function (x) {varclg(x,a=a,b=b,theta=theta)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Log gamma

dloggamma=function (x, a=1, r=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/x)*dgamma(-log(x),shape=r,rate=a)
pdf[log==TRUE]=dgamma(-log(x),shape=r,rate=a,log=TRUE)-log(x)
return(pdf)
}

ploggamma=function (x, a=1, r=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-pgamma(-a*log(x),shape=r)
cdf[log.p==FALSE&lower.tail==FALSE]=pgamma(-a*log(x),shape=r)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-pgamma(-a*log(x),shape=r))
cdf[log.p==TRUE&lower.tail==FALSE]=pgamma(-a*log(x),shape=r,log.p=TRUE)
return(cdf)
}

varloggamma=function (p, a=1, r=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=exp(-(1/a)*qgamma(1-p,shape=r))
return(var)
}

esloggamma=function (p, a=1, r=1)
{
f=function (x) {varloggamma(x,a=a,r=r)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Inverse gamma

dinvgamma=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/x)*dgamma(1/x,shape=a,rate=b)
pdf[log==TRUE]=dgamma(1/x,shape=a,rate=b,log=TRUE)-log(x)
return(pdf)
}

pinvgamma=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-pgamma(b/x,shape=a)
cdf[log.p==FALSE&lower.tail==FALSE]=pgamma(b/x,shape=a)
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-pgamma(b/x,shape=a))
cdf[log.p==TRUE&lower.tail==FALSE]=pgamma(b/x,shape=a,log.p=TRUE)
return(cdf)
}

varinvgamma=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=b/qgamma(1-p,shape=a)
return(var)
}

esinvgamma=function (p, a=1, b=1)
{
f=function (x) {varinvgamma(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Stacy gamma

dstacygamma=function (x, gamma=1, c=1, theta=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(c*x**(c-1)/theta**c)*dgamma((x/theta)**c,shape=gamma)
pdf[log==TRUE]=log(c)+(c-1)*log(x)-c*log(theta)+dgamma((x/theta)**c,shape=gamma,log=TRUE)
return(pdf)
}

pstacygamma=function (x, gamma=1, c=1, theta=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pgamma((x/theta)**c,shape=gamma,log.p=log.p,lower.tail=TRUE)
return(cdf)
}

varstacygamma=function (p, gamma=1, c=1, theta=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=theta*(qgamma(p,shape=gamma))**(1/c)
return(var)
}

esstacygamma=function (p, gamma=1, c=1, theta=1)
{
f=function (x) {varstacygamma(x,gamma=gamma,c=c,theta=theta)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#beta

{
pdf=dbeta(x,shape1=a,shape2=b,log=log)
return(pdf)
}

pbetadist=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x,shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varbetadist=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
var=qbeta(p,shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(var)
}

{
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Uniform

duniform=function (x, a=0, b=1, log=FALSE)
{
pdf=dunif(x,min=a,max=b,log=log)
return(pdf)
}

puniform=function (x, a=0, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=punif(x,min=a,max=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varuniform=function (p, a=0, b=1, log.p=FALSE, lower.tail=TRUE)
{
var=qunif(p,min=a,max=b,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esuniform=function (p, a=0, b=1)
{
f=function (x) {varuniform(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized uniform

dgenunif=function (x, a=0, c=1, h=1, k=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=h*k*c*(x-a)**(c-1)*(1-k*(x-a)**c)**(h-1)
pdf[log==TRUE]=log(h*k*c)+(c-1)*log(x-a)+(h-1)*log(1-k*(x-a)**c)
return(pdf)
}

pgenunif=function (x, a=0, c=1, h=1, k=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-k*(x-a)**c)**h
cdf[log.p==FALSE&lower.tail==FALSE]=(1-k*(x-a)**c)**h
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-k*(x-a)**c)**h)
cdf[log.p==TRUE&lower.tail==FALSE]=h*log(1-k*(x-a)**c)
return(cdf)
}

vargenunif=function (p, a=0, c=1, h=1, k=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=a+k**(-1/c)*(1-(1-p)**(1/h))**(1/c)
return(var)
}

esgenunif=function (p, a=0, c=1, h=1, k=1)
{
f=function (x) {vargenunif(x,a=a,c=c,h=h,k=k)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#power function I

dpower1=function (x, a=1, log=FALSE)
{
pdf=dbeta(x,shape1=a,shape2=1,log=log)
return(pdf)
}

ppower1=function (x, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x,shape1=a,shape2=1,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varpower1=function (p, a=1, log.p=FALSE, lower.tail=TRUE)
{
var=qbeta(p,shape1=a,shape2=1,log.p=log.p,lower.tail=lower.tail)
return(var)
}

espower1=function (p, a=1)
{
f=function (x) {varpower1(x,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#power function II

dpower2=function (x, b=1, log=FALSE)
{
pdf=dbeta(x,shape1=1,shape2=b,log=log)
return(pdf)
}

ppower2=function (x, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x,shape1=1,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varpower2=function (p, b=1, log.p=FALSE, lower.tail=TRUE)
{
var=qbeta(p,shape1=1,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(var)
}

espower2=function (p, b=1)
{
f=function (x) {varpower2(x,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Log beta

dlogbeta=function (x, a=1, b=1, c=1, d=2, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(dbeta(log(x/c)/log(d/c),shape1=a,shape2=b))/(x*log(d/c))
pdf[log==TRUE]=dbeta(log(x/c)/log(d/c),shape1=a,shape2=b,log=TRUE)-log(x)-log(log(d/c))
return(pdf)
}

plogbeta=function (x, a=1, b=1, c=1, d=2, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(log(x/c)/log(d/c),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varlogbeta=function (p, a=1, b=1, c=1, d=2, log.p=FALSE, lower.tail=TRUE)
{
var=c*(d/c)**(qbeta(p,shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail))
return(var)
}

eslogbeta=function (p, a=1, b=1, c=1, d=2)
{
f=function (x) {varlogbeta(x,a=a,b=b,c=c,d=d)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Complementary beta

dcompbeta=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=beta(a,b)*(qbeta(x,shape1=a,shape2=b))**(1-a)*(1-qbeta(x,shape1=a,shape2=b))**(1-b)
pdf[log==TRUE]=lbeta(a,b)+(1-a)*log(qbeta(x,shape1=a,shape2=b))+(1-b)*log(1-qbeta(x,shape1=a,shape2=b))
return(pdf)
}

pcompbeta=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=qbeta(x,shape1=a,shape2=b)
cdf[log.p==FALSE&lower.tail==FALSE]=1-qbeta(x,shape1=a,shape2=b)
cdf[log.p==TRUE&lower.tail==TRUE]=log(qbeta(x,shape1=a,shape2=b))
cdf[log.p==TRUE&lower.tail==FALSE]=log(1-qbeta(x,shape1=a,shape2=b))
return(cdf)
}

varcompbeta=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=pbeta(p,shape1=a,shape2=b)
return(var)
}

escompbeta=function (p, a=1, b=1)
{
f=function (x) {varcompbeta(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Libby Novick beta

dLNbeta=function (x, lambda=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=lambda*(1+(lambda-1)*x)**(-2)*dbeta(lambda*x/(1+(lambda-1)*x),shape1=a,shape2=b)
pdf[log==TRUE]=log(lambda)-2*log(1+(lambda-1)*x)+dbeta(lambda*x/(1+(lambda-1)*x),shape1=a,shape2=b,log=TRUE)
return(pdf)
}

pLNbeta=function (x, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(lambda*x/(1+(lambda-1)*x),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varLNbeta=function (p, lambda=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=qbeta(p,shape1=a,shape2=b)/(lambda-(lambda-1)*qbeta(p,shape1=a,shape2=b))
return(var)
}

esLNbeta=function (p, lambda=1, a=1, b=1)
{
f=function (x) {varLNbeta(x,lambda=lambda,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#McDonald Richards beta

dMRbeta=function (x, a=1, b=1, r=1, q=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(r/b)*(x**(r-1)/q**r)*dbeta((1/b)*(x/q)**r,shape1=a,shape2=b)
pdf[log==TRUE]=log(r)-log(b)-r*log(q)+(r-1)*log(x)+dbeta((1/b)*(x/q)**r,shape1=a,shape2=b,log=TRUE)
return(pdf)
}

pMRbeta=function (x, a=1, b=1, r=1, q=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta((1/b)*(x/q)**r,shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varMRbeta=function (p, a=1, b=1, r=1, q=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=b**(1/r)*q*(qbeta(p,shape1=a,shape2=b))**(1/r)
return(var)
}

esMRbeta=function (p, a=1, b=1, r=1, q=1)
{
f=function (x) {varMRbeta(x,a=a,b=b,r=r,q=q)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized beta

dgenbeta=function (x, a=1, b=1, c=0, d=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/(d-c))*dbeta((x-c)/(d-c),shape1=a,shape2=b)
pdf[log==TRUE]=dbeta((x-c)/(d-c),shape1=a,shape2=b,log=TRUE)-log(d-c)
return(pdf)
}

pgenbeta=function (x, a=1, b=1, c=0, d=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta((x-c)/(d-c),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

vargenbeta=function (p, a=1, b=1, c=0, d=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=c+(d-c)*qbeta(p,shape1=a,shape2=b)
return(var)
}

esgenbeta=function (p, a=1, b=1, c=0, d=1)
{
f=function (x) {vargenbeta(x,a=a,b=b,c=c,d=d)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Arcsine

darcsine=function (x, a=0, b=1, log=FALSE)
{
pdf=dgenbeta(x,a=0.5,b=0.5,c=a,d=b,log=log)
return(pdf)
}

parcsine=function (x, a=0, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pgenbeta(x,a=0.5,b=0.5,c=a,d=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

vararcsine=function (p, a=0, b=1, log.p=FALSE, lower.tail=TRUE)
{
var=vargenbeta(p,a=0.5,b=0.5,c=a,d=b,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esarcsine=function (p, a=0, b=1)
{
f=function (x) {vararcsine(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Triangular

dtriangular=function (x, a=0, b=2, c=1, log=FALSE)
{
pdf=x
pdf[log==FALSE&x<a]=0
pdf[log==FALSE&x>=a&x<=c]=2*(x[x>=a&x<=c]-a)/((b-a)*(c-a))
pdf[log==FALSE&x>c&x<=b]=2*(b-x[x>c&x<=b])/((b-a)*(b-c))
pdf[log==FALSE&x>b]=1
pdf[log==TRUE&x<a]=-Inf
pdf[log==TRUE&x>=a&x<=c]=log(2*(x[x>=a&x<=c]-a)/((b-a)*(c-a)))
pdf[log==TRUE&x>c&x<=b]=log(2*(b-x[x>c&x<=b])/((b-a)*(b-c)))
pdf[log==TRUE&x>b]=0
return(pdf)
}

ptriangular=function (x, a=0, b=2, c=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE&x<a]=0
cdf[log.p==FALSE&lower.tail==TRUE&x>=a&x<=c]=(x[x>=a&x<=c]-a)**2/((b-a)*(c-a))
cdf[log.p==FALSE&lower.tail==TRUE&x>c&x<=b]=1-(b-x[x>c&x<=b])**2/((b-a)*(b-c))
cdf[log.p==FALSE&lower.tail==TRUE&x>b]=1
cdf[log.p==FALSE&lower.tail==FALSE&x<a]=1
cdf[log.p==FALSE&lower.tail==FALSE&x>=a&x<=c]=1-(x[x>=a&x<=c]-a)**2/((b-a)*(c-a))
cdf[log.p==FALSE&lower.tail==FALSE&x>c&x<=b]=(b-x[x>c&x<=b])**2/((b-a)*(b-c))
cdf[log.p==FALSE&lower.tail==FALSE&x>b]=0
cdf[log.p==TRUE&lower.tail==TRUE&x<a]=-Inf
cdf[log.p==TRUE&lower.tail==TRUE&x>=a&x<=c]=log((x[x>=a&x<=c]-a)**2/((b-a)*(c-a)))
cdf[log.p==TRUE&lower.tail==TRUE&x>c&x<=b]=log(1-(b-x[x>c&x<=b])**2/((b-a)*(b-c)))
cdf[log.p==TRUE&lower.tail==TRUE&x>b]=0
cdf[log.p==TRUE&lower.tail==FALSE&x<a]=0
cdf[log.p==TRUE&lower.tail==FALSE&x>=a&x<=c]=log(1-(x[x>=a&x<=c]-a)**2/((b-a)*(c-a)))
cdf[log.p==TRUE&lower.tail==FALSE&x>c&x<=b]=log((b-x[x>c&x<=b])**2/((b-a)*(b-c)))
cdf[log.p==TRUE&lower.tail==FALSE&x>b]=-Inf
return(cdf)
}

vartriangular=function (p, a=0, b=2, c=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<(c-a)/(b-a)]=a+sqrt(p[p<(c-a)/(b-a)]*(b-a)*(c-a))
var[p>=(c-a)/(b-a)]=b-sqrt((1-p[p>=(c-a)/(b-a)])*(b-a)*(b-c))
return(var)
}

estriangular=function (p, a=0, b=2, c=1)
{
f=function (x) {vartriangular(x,a=a,b=b,c=c)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized beta II

dgenbeta2=function (x, a=1, b=1, c=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=c*x**(c-1)*dbeta(x**c,shape1=a,shape2=b)
pdf[log==TRUE]=dbeta(x**c,shape1=a,shape2=b,log=TRUE)+log(c)+(c-1)*log(x)
return(pdf)
}

pgenbeta2=function (x, a=1, b=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x**c,shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

vargenbeta2=function (p, a=1, b=1, c=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(qbeta(p,shape1=a,shape2=b))**(1/c)
return(var)
}

esgenbeta2=function (p, a=1, b=1, c=1)
{
f=function (x) {vargenbeta2(x,a=a,b=b,c=c)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Inverse beta

dinvbeta=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1+x)**(-2)*dbeta(x/(1+x),shape1=a,shape2=b)
pdf[log==TRUE]=dbeta(x/(1+x),shape1=a,shape2=b,log=TRUE)-2*log(1+x)
return(pdf)
}

pinvbeta=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x/(1+x),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varinvbeta=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(qbeta(p,shape1=a,shape2=b))/(1-qbeta(p,shape1=a,shape2=b))
return(var)
}

esinvbeta=function (p, a=1, b=1)
{
f=function (x) {varinvbeta(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized inverse beta

dgeninvbeta=function (x, a=1, c=1, d=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*x**(a-1)*(1+x**a)**(-2)*dbeta(x**a/(1+x**a),shape1=c,shape2=d)
pdf[log==TRUE]=log(a)+(a-1)*log(x)-2*log(1+x**a)+dbeta(x**a/(1+x**a),shape1=c,shape2=d,log=TRUE)
return(pdf)
}

pgeninvbeta=function (x, a=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x**a/(1+x**a),shape1=c,shape2=d,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

vargeninvbeta=function (p, a=1, c=1, d=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=((qbeta(p,shape1=c,shape2=d))/(1-qbeta(p,shape1=c,shape2=d)))**(1/a)
return(var)
}

esgeninvbeta=function (p, a=1, c=1, d=1)
{
f=function (x) {vargeninvbeta(x,a=a,c=c,d=d)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Two sided power

dtsp=function (x, a=1, theta=0.5, log=FALSE)
{
pdf=x
pdf[log==FALSE&x<=theta]=a*(x[x<=theta]/theta)**(a-1)
pdf[log==FALSE&x>theta]=a*((1-x[x>theta])/(1-theta))**(a-1)
pdf[log==TRUE&x<=theta]=log(a)+(a-1)*log(x[x<=theta]/theta)
pdf[log==TRUE&x>theta]=log(a)+(a-1)*log((1-x[x>theta])/(1-theta))
return(pdf)
}

ptsp=function (x, a=1, theta=0.5, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE&x<=theta]=theta*(x[x<=theta]/theta)**a
cdf[log.p==FALSE&lower.tail==TRUE&x>theta]=1-(1-theta)*((1-x[x>theta])/(1-theta))**a
cdf[log.p==FALSE&lower.tail==FALSE&x<=theta]=1-theta*(x[x<=theta]/theta)**a
cdf[log.p==FALSE&lower.tail==FALSE&x>theta]=(1-theta)*((1-x[x>theta])/(1-theta))**a
cdf[log.p==TRUE&lower.tail==TRUE&x<=theta]=(1-a)*log(theta)+a*log(x[x<=theta])
cdf[log.p==TRUE&lower.tail==TRUE&x>theta]=log(1-(1-theta)*((1-x[x>theta])/(1-theta))**a)
cdf[log.p==TRUE&lower.tail==FALSE&x<=theta]=log(1-theta*(x[x<=theta]/theta)**a)
cdf[log.p==TRUE&lower.tail==FALSE&x>theta]=(1-a)*log(1-theta)+a*log(1-x[x>theta])
return(cdf)
}

vartsp=function (p, a=1, theta=0.5, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<=theta]=theta*(p[p<=theta]/theta)**(1/a)
var[p>theta]=1-(1-theta)*((1-p[p>theta])/(1-theta))**(1/a)
return(var)
}

estsp=function (p, a=1, theta=0.5)
{
f=function (x) {vartsp(x,a=a,theta=theta)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy

dkum=function (x, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*b*x**(a-1)*(1-x**a)**(b-1)
pdf[log==TRUE]=log(a*b)+(a-1)*log(x)+(b-1)*log(1-x**a)
return(pdf)
}

pkum=function (x, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-x**a)**b
cdf[log.p==FALSE&lower.tail==FALSE]=(1-x**a)**b
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-x**a)**b)
cdf[log.p==TRUE&lower.tail==FALSE]=b*log(1-x**a)
return(cdf)
}

varkum=function (p, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=(1-(1-p)**(1/b))**(1/a)
return(var)
}

eskum=function (p, a=1, b=1)
{
f=function (x) {varkum(x,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Normal

dnormal=function (x, mu=0, sigma=1, log=FALSE)
{
pdf=dnorm(x,mean=mu,sd=sigma,log=log)
return(pdf)
}

pnormal=function (x, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pnorm(x,mean=mu,sd=sigma,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varnormal=function (p, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
var=qnorm(p,mean=mu,sd=sigma,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esnormal=function (p, mu=0, sigma=1)
{
f=function (x) {varnormal(x,mu=mu,sigma=sigma)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy normal

dkumnormal=function (x, mu=0, sigma=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=a*b*dnorm(x,mean=mu,sd=sigma)*(pnorm(x,mean=mu,sd=sigma))**(a-1)*(1-(pnorm(x,mean=mu,sd=sigma))**a)**(b-1)
pdf[log==TRUE]=log(a*b)+dnorm(x,mean=mu,sd=sigma,log=TRUE)+(a-1)*pnorm(x,mean=mu,sd=sigma,log.p=TRUE)+(b-1)*log(1-(pnorm(x,mean=mu,sd=sigma))**a)
return(pdf)
}

pkumnormal=function (x, mu=0, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-(pnorm(x,mean=mu,sd=sigma))**a)**b
cdf[log.p==FALSE&lower.tail==FALSE]=(1-(pnorm(x,mean=mu,sd=sigma))**a)**b
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-(pnorm(x,mean=mu,sd=sigma))**a)**b)
cdf[log.p==TRUE&lower.tail==FALSE]=b*log(1-(pnorm(x,mean=mu,sd=sigma))**a)
return(cdf)
}

varkumnormal=function (p, mu=0, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=mu+sigma*qnorm((1-(1-p)**(1/b))**(1/a))
return(var)
}

eskumnormal=function (p, mu=0, sigma=1, a=1, b=1)
{
f=function (x) {varkumnormal(x,mu=mu,sigma=sigma,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Exponential power

dexppower=function (x, mu=0, sigma=1, a=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(exp(-(abs(x-mu))**a*sigma**(-a)/a))/(2*a**(1/a)*sigma*gamma(1+1/a))
pdf[log==TRUE]=-(abs(x-mu))**a*sigma**(-a)/a-log(2)-(1/a)*log(a)-log(sigma)-lgamma(1+1/a)
return(pdf)
}

pexppower=function (x, mu=0, sigma=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE&x<=mu]=0.5-0.5*pgamma((mu-x[x<=mu])**a/(a*sigma**a),shape=1/a)
cdf[log.p==FALSE&lower.tail==TRUE&x>mu]=0.5+0.5*pgamma((x[x>mu]-mu)**a/(a*sigma**a),shape=1/a)
cdf[log.p==FALSE&lower.tail==FALSE&x<=mu]=0.5+0.5*pgamma((mu-x[x<=mu])**a/(a*sigma**a),shape=1/a)
cdf[log.p==FALSE&lower.tail==FALSE&x>mu]=0.5-0.5*pgamma((x[x>mu]-mu)**a/(a*sigma**a),shape=1/a)
cdf[log.p==TRUE&lower.tail==TRUE&x<=mu]=log(0.5-0.5*pgamma((mu-x[x<=mu])**a/(a*sigma**a),shape=1/a))
cdf[log.p==TRUE&lower.tail==TRUE&x>mu]=log(0.5+0.5*pgamma((x[x>mu]-mu)**a/(a*sigma**a),shape=1/a))
cdf[log.p==TRUE&lower.tail==FALSE&x<=mu]=log(0.5+0.5*pgamma((mu-x[x<=mu])**a/(a*sigma**a),shape=1/a))
cdf[log.p==TRUE&lower.tail==FALSE&x>mu]=log(0.5-0.5*pgamma((x[x>mu]-mu)**a/(a*sigma**a),shape=1/a))
return(cdf)
}

varexppower=function (p, mu=0, sigma=1, a=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<=0.5]=mu-a**(1/a)*sigma*(qgamma(1-2*p[p<=0.5],shape=1/a))**(1/a)
var[p>0.5]=mu+a**(1/a)*sigma*(qgamma(2*p[p>0.5]-1,shape=1/a))**(1/a)
return(var)
}

esexppower=function (p, mu=0, sigma=1, a=1)
{
f=function (x) {varexppower(x,mu=mu,sigma=sigma,a=a)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Asymmetric exponential power

daep=function (x, q1=1, q2=1, alpha=0.5, log=FALSE)
{
K1= 1/(2*q1**(1/q1)*gamma(1+1/q1))
K2= 1/(2*q2**(1/q2)*gamma(1+1/q2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
pdf=x
pdf[log==FALSE&x<=0]=alpha*(exp(-(1/q1)*(abs(x[x<=0]/(2*alphastar)))**q1))/(2*alphastar*q1**(1/q1)*gamma(1+1/q1))
pdf[log==FALSE&x>0]=(1-alpha)*(exp(-(1/q2)*(abs(x[x>0]/(2-2*alphastar)))**q2))/(2*(1-alphastar)*q2**(1/q2)*gamma(1+1/q2))
pdf[log==TRUE&x<=0]=log(alpha)-(1/q1)*(abs(x[x<=0]/(2*alphastar)))**q1-log(2)-log(alphastar)-(1/q1)*log(q1)-lgamma(1+1/q1)
pdf[log==TRUE&x>0]=log(1-alpha)-(1/q2)*(abs(x[x>0]/(2-2*alphastar)))**q2-log(2)-log(1-alphastar)-(1/q2)*log(q2)-lgamma(1+1/q2)
return(pdf)
}

paep=function (x, q1=1, q2=1, alpha=0.5, log.p=FALSE, lower.tail=TRUE)
{
K1= 1/(2*q1**(1/q1)*gamma(1+1/q1))
K2= 1/(2*q2**(1/q2)*gamma(1+1/q2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE&x<=0]=alpha-alpha*pgamma((1/q1)*(-x[x<=0]/(2*alphastar))**q1,shape=1/q1)
cdf[log.p==FALSE&lower.tail==TRUE&x>0]=alpha+(1-alpha)*pgamma((1/q2)*(-x[x>0]/(2-2*alphastar))**q2,shape=1/q2)
cdf[log.p==FALSE&lower.tail==FALSE&x<=0]=1-alpha+alpha*pgamma((1/q1)*(-x[x<=0]/(2*alphastar))**q1,shape=1/q1)
cdf[log.p==FALSE&lower.tail==FALSE&x>0]=1-alpha-(1-alpha)*pgamma((1/q2)*(-x[x>0]/(2-2*alphastar))**q2,shape=1/q2)
cdf[log.p==TRUE&lower.tail==TRUE&x<=0]=log(alpha-alpha*pgamma((1/q1)*(-x[x<=0]/(2*alphastar))**q1,shape=1/q1))
cdf[log.p==TRUE&lower.tail==TRUE&x>0]=log(alpha+(1-alpha)*pgamma((1/q2)*(-x[x>0]/(2-2*alphastar))**q2,shape=1/q2))
cdf[log.p==TRUE&lower.tail==FALSE&x<=0]=log(1-alpha+alpha*pgamma((1/q1)*(-x[x<=0]/(2*alphastar))**q1,shape=1/q1))
cdf[log.p==TRUE&lower.tail==FALSE&x>0]=log(1-alpha-(1-alpha)*pgamma((1/q2)*(-x[x>0]/(2-2*alphastar))**q2,shape=1/q2))
return(cdf)
}

varaep=function (p, q1=1, q2=1, alpha=0.5, log.p=FALSE, lower.tail=TRUE)
{
K1= 1/(2*q1**(1/q1)*gamma(1+1/q1))
K2= 1/(2*q2**(1/q2)*gamma(1+1/q2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
var[p<=alpha]=-2*alphastar*(q1*qgamma(1-p[p<=alpha]/alpha,shape=1/q1))**(1/q1)
var[p>alpha]=-2*(1-alphastar)*(q2*qgamma(1-(1-p[p>alpha])/(1-alpha),shape=1/q2))**(1/q2)
return(var)
}

esaep=function (p, q1=1, q2=1, alpha=0.5)
{
f=function (x) {varaep(x,q1=q1,q2=q2,alpha=alpha)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Beta normal

dbetanorm=function (x, mu=0, sigma=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=dnorm(x,mean=mu,sd=sigma)*dbeta(pnorm(x,mean=mu,sd=sigma),shape1=a,shape2=b)
pdf[log==TRUE]=dnorm(x,mean=mu,sd=sigma,log=TRUE)+dbeta(pnorm(x,mean=mu,sd=sigma),shape1=a,shape2=b,log=TRUE)
return(pdf)
}

pbetanorm=function (x, mu=0, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(pnorm(x,mean=mu,sd=sigma),shape1=a,shape2=b,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varbetanorm=function (p, mu=0, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=mu+sigma*qnorm(qbeta(p,shape1=a,shape2=b))
return(var)
}

esbetanorm=function (p, mu=0, sigma=1, a=1, b=1)
{
f=function (x) {varbetanorm(x,mu=mu,sigma=sigma,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Half normal

dhalfnorm=function (x, sigma=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=2*dnorm(x,sd=sigma)
pdf[log==TRUE]=log(2)+dnorm(x,sd=sigma,log=TRUE)
return(pdf)
}

phalfnorm=function (x, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=2*pnorm(x,sd=sigma)-1
cdf[log.p==FALSE&lower.tail==FALSE]=2-2*pnorm(x,sd=sigma)
cdf[log.p==TRUE&lower.tail==TRUE]=log(2*pnorm(x,sd=sigma)-1)
cdf[log.p==TRUE&lower.tail==FALSE]=log(2-2*pnorm(x,sd=sigma))
return(cdf)
}

varhalfnorm=function (p, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=sigma*qnorm(0.5*(1+p))
return(var)
}

eshalfnorm=function (p, sigma=1)
{
f=function (x) {varhalfnorm(x,sigma=sigma)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Kumaraswamy half normal

dkumhalfnorm=function (x, sigma=1, a=1, b=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=2*a*b*dnorm(x,sd=sigma)*(2*pnorm(x,sd=sigma)-1)**(a-1)*(1-(2*pnorm(x,sd=sigma)-1)**a)**(b-1)
pdf[log==TRUE]=log(2*a*b)+dnorm(x,sd=sigma,log=TRUE)+(a-1)*log(2*pnorm(x,sd=sigma)-1)+(b-1)*log(1-(2*pnorm(x,sd=sigma)-1)**a)
return(pdf)
}

pkumhalfnorm=function (x, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=1-(1-(2*pnorm(x,sd=sigma)-1)**a)**b
cdf[log.p==FALSE&lower.tail==FALSE]=(1-(2*pnorm(x,sd=sigma)-1)**a)**b
cdf[log.p==TRUE&lower.tail==TRUE]=log(1-(1-(2*pnorm(x,sd=sigma)-1)**a)**b)
cdf[log.p==TRUE&lower.tail==FALSE]=b*log(1-(2*pnorm(x,sd=sigma)-1)**a)
return(cdf)
}

varkumhalfnorm=function (p, sigma=1, a=1, b=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=sigma*qnorm(0.5+0.5*(1-(1-p)**(1/b))**(1/a))
return(var)
}

eskumhalfnorm=function (p, sigma=1, a=1, b=1)
{
f=function (x) {varkumhalfnorm(x,sigma=sigma,a=a,b=b)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Student's t

dT=function (x, n=1, log=FALSE)
{
pdf=dt(x,df=n,log=log)
return(pdf)
}

pT=function (x, n=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pt(x,df=n,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varT=function (p, n=1, log.p=FALSE, lower.tail=TRUE)
{
var=qt(p,df=n,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esT=function (p, n=1)
{
f=function (x) {varT(x,n=n)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Generalized asymmetric Student's t

dast=function (x, nu1=1, nu2=1, alpha=0.5, log=FALSE)
{
K1= gamma((nu1+1)/2)/(sqrt(nu1/2)*gamma(nu1/2))
K2= gamma((nu2+1)/2)/(sqrt(nu2/2)*gamma(nu2/2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
pdf=x
pdf[log==FALSE&x<=0]=(alpha/alphastar)*dt(x[x<=0]/(2*alphastar),df=nu1)
pdf[log==FALSE&x>0]=((1-alpha)/(1-alphastar))*dt(x[x>0]/(2*(1-alphastar)),df=nu2)
pdf[log==TRUE&x<=0]=log(alpha/alphastar)+dt(x[x<=0]/(2*alphastar),df=nu1,log=TRUE)
pdf[log==TRUE&x>0]=log((1-alpha)/(1-alphastar))+dt(x[x>0]/(2*(1-alphastar)),df=nu2,log=TRUE)
return(pdf)
}

past=function (x, nu1=1, nu2=1, alpha=0.5, log.p=FALSE, lower.tail=TRUE)
{
K1= gamma((nu1+1)/2)/(sqrt(nu1/2)*gamma(nu1/2))
K2= gamma((nu2+1)/2)/(sqrt(nu2/2)*gamma(nu2/2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
cdf=x
mm=x
MM=x
for (i in 1:length(x))
{mm[i]=min(x[i],0)
MM[i]=max(x[i],0)}
cdf[log.p==FALSE&lower.tail==TRUE]=2*alpha*pt(mm/(2*alphastar),df=nu1)-1+alpha+2*(1-alpha)*pt(MM/(2-2*alphastar),df=nu2)
cdf[log.p==FALSE&lower.tail==FALSE]=1-2*alpha*pt(mm/(2*alphastar),df=nu1)+1-alpha-2*(1-alpha)*pt(MM/(2-2*alphastar),df=nu2)
cdf[log.p==TRUE&lower.tail==TRUE]=log(2*alpha*pt(mm/(2*alphastar),df=nu1)-1+alpha+2*(1-alpha)*pt(MM/(2-2*alphastar),df=nu2))
cdf[log.p==TRUE&lower.tail==FALSE]=log(1-2*alpha*pt(mm/(2*alphastar),df=nu1)+1-alpha-2*(1-alpha)*pt(MM/(2-2*alphastar),df=nu2))
return(cdf)
}

varast=function (p, nu1=1, nu2=1, alpha=0.5, log.p=FALSE, lower.tail=TRUE)
{
K1= gamma((nu1+1)/2)/(sqrt(nu1/2)*gamma(nu1/2))
K2= gamma((nu2+1)/2)/(sqrt(nu2/2)*gamma(nu2/2))
alphastar=alpha*K1/(alpha*K1+(1-alpha)*K2)
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=p
mm=p
MM=p
for (i in 1:length(p))
{mm[i]=min(p[i],alpha)
MM[i]=max(p[i],alpha)}
var=2*alphastar*qt(mm/(2*alpha),df=nu1)+2*(1-alphastar)*qt((MM+1-2*alpha)/(2-2*alpha),df=nu2)
return(var)
}

esast=function (p, nu1=1, nu2=1, alpha=0.5)
{
f=function (x) {varast(x,nu1=nu1,nu2=nu2,alpha=alpha)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Half Student's t

dhalfT=function (x, n=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=2*dt(x,df=n)
pdf[log==TRUE]=log(2)+dt(x,df=n,log=TRUE)
return(pdf)
}

phalfT=function (x, n=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pbeta(x*x/(n+x*x),shape1=0.5,shape2=n/2,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varhalfT=function (p, n=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=sqrt(n*qbeta(p,shape1=0.5,shape2=n/2)/(1-qbeta(p,shape1=0.5,shape2=n/2)))
return(var)
}

eshalfT=function (p, n=1)
{
f=function (x) {varhalfT(x,n=n)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Cauchy

dCauchy=function (x, mu=0, sigma=1, log=FALSE)
{
pdf=dcauchy(x,location=mu,scale=sigma,log=log)
return(pdf)
}

pCauchy=function (x, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=pcauchy(x,location=mu,scale=sigma,log.p=log.p,lower.tail=lower.tail)
return(cdf)
}

varCauchy=function (p, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
var=qcauchy(p,location=mu,scale=sigma,log.p=log.p,lower.tail=lower.tail)
return(var)
}

esCauchy=function (p, mu=0, sigma=1)
{
f=function (x) {varCauchy(x,mu=mu,sigma=sigma)}
es=p
for (i in 1:length(p)) {es[i]=(1/p[i])*integrate(f,lower=0,upper=p[i],stop.on.error=FALSE)\$value}
return(es)
}

#Log Cauchy

dlogcauchy=function (x, mu=0, sigma=1, log=FALSE)
{
pdf=x
pdf[log==FALSE]=(1/x)*dcauchy(log(x),location=mu,scale=sigma)
pdf[log==TRUE]=dcauchy(log(x),location=mu,scale=sigma,log=TRUE)-log(x)
return(pdf)
}

plogcauchy=function (x, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
cdf=x
cdf[log.p==FALSE&lower.tail==TRUE]=(1/pi)*atan((log(x)-mu)/sigma)
cdf[log.p==FALSE&lower.tail==FALSE]=1-(1/pi)*atan((log(x)-mu)/sigma)
cdf[log.p==TRUE&lower.tail==TRUE]=log((1/pi)*atan((log(x)-mu)/sigma))
cdf[log.p==TRUE&lower.tail==FALSE]=log(1-(1/pi)*atan((log(x)-mu)/sigma))
return(cdf)
}

varlogcauchy=function (p, mu=0, sigma=1, log.p=FALSE, lower.tail=TRUE)
{
if (log.p==TRUE) p=exp(p)
if (lower.tail==FALSE) p=1-p
var=exp(mu```