R/VaRES.R

#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

dbetadist=function (x, a=1, b=1, log=FALSE)
{
	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)
}

esbetadist=function (p, a=1, b=1)
{
	f=function (x) {varbetadist(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)
}


#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