#### Functions for continuous truncated Pareto distributions

# by Lai jiangshan

# dtparetoProbability density

# ptparetoProbability distribution (CDF)

# qtparetoQuantile function

# rtparetoRandom variable generation

# Probability density of truncated Pareto distributions

# Gives NA on values below the lower

# Input: Data vector, lower threshold,upper threshold, scaling exponent, "log" flag

# Output: Vector of (log) probability densities

dtpareto <- function(x, lower,upper, exponent, log=FALSE) {

 # Avoid doing limited-precision arithmetic followed by logs if we want

 # the log!

 if (!log) {

   prefactor <- (1-exponent)/(upper^(1-exponent)-lower^(1-exponent))

   f <- function(x) {prefactor*(x)^(-exponent)}

 } else {

   prefactor.log <- log((1-exponent)/(upper^(1-exponent)-lower^(1-exponent)))

   f <- function(x) {prefactor.log -exponent*log(x)}

 }

 d <- ifelse(x<lower,NA,f(x))

 return(d)

}

# Cumulative distribution function of the Pareto distributions

# Gives NA on values < threshold

# Input: Data vector, lower threshold,upper threshold scaling exponent, usual flags

# Output: Vector of (log) probabilities

ptpareto <- function(x, lower, upper,exponent, lower.tail=TRUE, log.p=FALSE)

{

 if ((!lower.tail) && (!log.p)) {

   f <- function(x) {1-(x^(1-exponent)-lower^(1-exponent))/(upper^(1-exponent)-lower^(1-exponent))}

}

 if ((lower.tail) && (!log.p)) {

   f <- function(x) { (x^(1-exponent)-lower^(1-exponent))/(upper^(1-exponent)-lower^(1-exponent))}

}

 if ((!lower.tail) && (log.p)) {

   f <- function(x) {log(1-(x^(1-exponent)-lower^(1-exponent))/(upper^(1-exponent)-lower^(1-exponent)))}

}

 if ((lower.tail) && (log.p)) {

   f <- function(x) {log((x^(1-exponent)-lower^(1-exponent))/(upper^(1-exponent)-lower^(1-exponent)))}

}

 p <- ifelse(x < lower, NA, f(x))

 return(p)

}

# Quantiles of Pareto distributions

# Input: vector of probabilities, lower threshold,upper threshold,scaling exponent, usual flags

# Output: Vector of quantile values

qtpareto <- function(p, lower,upper, exponent, lower.tail=TRUE, log.p=FALSE) {

 if (log.p) {

   p <- exp(p)

 }

 if (lower.tail) {

   p <- 1-p

 }

 # This works, via the recycling rule

 # q<-(p^(1/(1-exponent)))*lower

 q <- (p*(upper^(1-exponent)-lower^(1-exponent))+lower^(1-exponent))^(1/(1-exponent))

 return(q)

}

# Generate Pareto-distributed random variates

# Input: Integer size, lower threshold,upper threshold scaling exponent

# Output: Vector of real-valued random variates

rtpareto <- function(n, lower,upper, exponent) {

 # Using the transformation method, because we know the quantile function

 # analytically

 # Consider replacing with a non-R implementation of transformation method

 ru <- runif(n)

 r<-qtpareto(ru,lower,upper,exponent)

 return(r)

}