R/distributions.R
In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities
Defines functions
qcustom
qslash
qtnorm
ptnorm
dtnorm
phalft
rhalft
rsnorm
ral
qal
pal
dal
Documented in
dal
dtnorm
pal
phalft
ptnorm
qal
qslash
qtnorm
ral
rhalft
rsnorm
###########################################################################
######### Power Exponential (Generalized Normal) Distribution) ##########
###########################################################################
#' Power Exponential random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#'
#' @return vector
#' @export
#'
#' @examples
#' rpowexp()
rpowexp = function (n, mu = 0, sigma = 1, nu = 2)
{
if (any(sigma <= 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(n <= 0))
stop(paste("n must be a positive integer", "\n", ""))
n <- ceiling(n)
p <- runif(n)
r <- qpowexp(p, mu = mu, sigma = sigma, nu = nu)
r
}
#' Power Exponential quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qpowexp()
qpowexp = function (p, mu = 0, sigma = 1, nu = 2, lower.tail = TRUE, log = FALSE)
{
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
if (log == TRUE)
p <- exp(p)
else p <- p
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
if (any(p < 0) | any(p > 1))
stop(paste("p must be between 0 and 1", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
suppressWarnings(s <- qgamma((2 * p - 1) * sign(p - 0.5), shape = (1/nu), scale = 1))
z <- sign(p - 0.5) * ((2 * s)^(1/nu)) * c
ya <- mu + sigma * z
ya
}
#' Power Exponential probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dpowexp()
dpowexp = function (x, mu = 0, sigma = 1, nu = 2, log = FALSE) {
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
z <- (x - mu)/sigma
log.lik <- -log(sigma) + log(nu) - log.c - (0.5 * (abs(z/c)^nu)) -
(1 + (1/nu)) * log(2) - lgamma(1/nu)
if (log == FALSE)
fy <- exp(log.lik)
else fy <- log.lik
fy
}
#' Power Exponential cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' ppowexp()
ppowexp = function (q, mu = 0, sigma = 1, nu = 2, lower.tail = TRUE, log = FALSE) {
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
z <- (q - mu)/sigma
s <- 0.5 * ((abs(z/c))^nu)
cdf <- 0.5 * (1 + pgamma(s, shape = 1/nu, scale = 1) * sign(z))
if (length(nu) > 1)
cdf <- ifelse(nu > 10000, (q - (mu - sqrt(3) * sigma))/(sqrt(12) * sigma), cdf)
else cdf <- if (nu > 10000)
(q - (mu - sqrt(3) * sigma))/(sqrt(12) * sigma)
else cdf
if (lower.tail == TRUE)
cdf <- cdf
else cdf <- 1 - cdf
if (log == FALSE)
cdf <- cdf
else cdf <- log(cdf)
cdf
}
###########################################################################
######### Student's T Distribution ##########
###########################################################################
#' Student-T random number generator
#'
#' @param n number of observations
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#'
#' @return vector
#' @export
#'
#' @examples
#' rst()
rst = function (n, nu = 3, mu = 0, scale = 1) {
if (any(scale <= 0))
stop("the scale parameter must be positive.")
if (any(nu <= 0))
stop("the nu parameter must be positive.")
rGamma = function(n, shape = 1, rate = 1) {
return(qgamma(seq(0.0001, 0.9999, length.out = n),
shape, rate))
}
gamma.samps = rGamma(n, nu * 0.50, (scale^2) * (nu * 0.50))
sigmas = sqrt(1 / gamma.samps)
x = rnorm(n, mu, sigmas)
return(x)
}
#' Student-T probability density function
#'
#' @param x vector of quantiles
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dst()
dst = function (x, nu = 3, mu = 0, scale = 1, log = FALSE)
{
if (any(scale <= 0)) {
stop("the scale parameter must be positive.")
}
if (any(nu <= 0)) {
stop("the nu parameter must be positive.")
}
if (log) {
dt((x - mu)/scale, df = nu, log = TRUE) - log(scale)
}
else {
dt((x - mu)/scale, df = nu)/scale
}
}
#' Student-T quantile function
#'
#' @param p vector of probabilities.
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qst()
qst = function (p, nu = 3, mu = 0, scale = 1, lower.tail = TRUE, log = FALSE)
{
if (any(p < 0) || any(p > 1))
stop("p must be in [0,1].")
if (any(scale <= 0))
stop("the scale parameter must be positive.")
if (any(nu <= 0))
stop("the nu parameter must be positive.")
NN <- max(length(p), length(mu), length(scale), length(nu))
p <- rep(p, len = NN)
mu <- rep(mu, len = NN)
scale <- rep(scale, len = NN)
nu <- rep(nu, len = NN)
q <- mu + scale * qt(p, df = nu, lower.tail = lower.tail)
temp <- which(nu > 1e+06)
q[temp] <- qnorm(p[temp], mu[temp], scale[temp], lower.tail = lower.tail,
log.p = log)
return(q)
}
#' Student-T cumulative probability function
#'
#' @param q vector of quantiles
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' pst()
pst = function (q, nu = 3, mu = 0, scale = 1, lower.tail = TRUE, log = FALSE)
{
if (any(scale <= 0)) {
stop("the scale parameter must be positive.")
}
if (any(nu <= 0)) {
stop("the nu parameter must be positive.")
}
pt((q - mu)/scale, df = nu, lower.tail = lower.tail, log.p = log)
}
###########################################################################
######### Asymmetric Laplace Distribution ##########
###########################################################################
#' Asymmetric Laplace probability density function
#'
#' @param x vector of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
dal <- function(x, location=0, scale=1, tau=0.5, log=FALSE){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
x <- as.vector(x); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(x), length(location), length(scale), length(kappa))
x <- rep(x, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- rep(kappa, len=NN)
logconst <- 0.5 * log(2) - log(scale) + log(kappa) - log1p(kappa^2)
temp <- which(x < location); kappa[temp] <- 1/kappa[temp]
exponent <- -(sqrt(2) / scale) * abs(x - location) * kappa
dens <- logconst + exponent
if(log == FALSE) dens <- exp(logconst + exponent)
return(dens)
}
#' Asymmetric Laplace cumulative probability function
#'
#' @param q a vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
pal <- function(q, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
q <- as.vector(q); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if((kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(q), length(location), length(scale), length(kappa))
q <- rep(q, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- k2 <- rep(kappa, len=NN)
temp <- which(q < location); k2[temp] <- 1/kappa[temp]
exponent <- -(sqrt(2) / scale) * abs(q - location) * k2
temp <- exp(exponent) / (1 + kappa^2)
p <- 1 - temp
index1 <- (q < location)
p[index1] <- (kappa[index1])^2 * temp[index1]
return(p)
}
#' Asymmetric Laplace quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
qal <- function(p, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
p <- as.vector(p); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(p < 0) || any(p > 1)) stop("p must be in [0,1].")
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(p), length(location), length(scale), length(kappa))
p <- rep(p, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- rep(kappa, len=NN)
q <- p
temp <- kappa^2 / (1 + kappa^2)
index1 <- (p <= temp)
exponent <- p[index1] / temp[index1]
q[index1] <- location[index1] + (scale[index1] * kappa[index1]) *
log(exponent) / sqrt(2)
q[!index1] <- location[!index1] - (scale[!index1] / kappa[!index1]) *
(log1p((kappa[!index1])^2) + log1p(-p[!index1])) / sqrt(2)
q[p == 0] = -Inf
q[p == 1] = Inf
return(q)
}
#' Asymmetric Laplace random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#'
#' @return vector
#' @export
#'
ral <- function(n, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
location <- rep(location, len=n)
scale <- rep(scale, len=n)
kappa <- rep(kappa, len=n)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
x <- location + scale * log(runif(n)^kappa / runif(n)^(1/kappa)) / sqrt(2)
return(x)
}
###########################################################################
######### Sinh-Arcsinh (SHASH) Distribution ##########
###########################################################################
#' Sinh-Arcsinh (SHASH) random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#'
#' @return vector
#' @export
#'
#' @examples
#' rshash()
rshash = function (n, mu = 0, scale = 1, alpha = 0, nu = 1)
{
if (any(n <= 0))
stop(paste("n must be a positive integer", "\n", ""))
n <- ceiling(n)
p <- runif(n)
r <- qshash(p, mu = mu, scale = scale, alpha = alpha, nu = nu)
r
}
#' Sinh-Arcsinh (SHASH) probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dshash()
dshash = function (x, mu = 0, scale = 1, alpha = 0, nu = 1, log = FALSE){
if (any(scale < 0))
stop(paste("scale must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("nu must be positive", "\n", ""))
z <- (x - mu)/(scale * nu)
c <- cosh(nu * asinh(z) - alpha)
r <- sinh(nu * asinh(z) - alpha)
loglik <- -log(scale) - 0.5 * log(2 * pi) - 0.5 * log(1 + (z^2)) + log(c) - 0.5 * (r^2)
if (log == FALSE)
fy <- exp(loglik)
else fy <- loglik
fy
}
#' Sinh-Arcsinh (SHASH) cumulative probability function
#'
#' @param v vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' pshash()
pshash = function (q, mu = 0, scale = 1, alpha = 0, nu = 1, lower.tail = TRUE,
log.p = FALSE)
{
if (any(scale < 0))
stop(paste("scale must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("nu must be positive", "\n", ""))
z <- (q - mu)/(scale * nu)
c <- cosh(nu * asinh(z) - alpha)
r <- sinh(nu * asinh(z) - alpha)
p <- pNO(r)
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
if (log.p == FALSE)
p <- p
else p <- log(p)
p
}
#' Sinh-Arcsinh (SHASH) quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qshash()
qshash = function (p, mu = 0, scale = 1, alpha = 0, nu = 1, lower.tail = TRUE, log.p = FALSE)
{
if (log.p == TRUE)
p <- exp(p)
else p <- p
if (any(p <= 0) | any(p >= 1))
stop(paste("p must be between 0 and 1", "\n", ""))
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
y <- mu + (nu * scale) * sinh((1/nu) * asinh(qnorm(p)) +
(alpha/nu))
y
}
###########################################################################
############## Skew-Normal Distribution ##############
###########################################################################
#' Skew Normal random number generator
#'
#' @param n number of observations to generate
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#'
#' @return vector
#' @export
#'
#' @examples
#' rsnorm()
rsnorm <- function(n, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL){
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL)
{
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega)
with(args, {
z1 <- rnorm(n)
z2 <- rnorm(n)
id <- z2 > args$alpha * z1
z1[id] <- -z1[id]
xi + omega * z1
})
}
#' Skew Normal probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param log if TRUE, densities are returned as log(density)
#'
#' @return vector
#' @export
#'
#' @examples
#' dsnorm()
dsnorm <- function (x, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL,
log = FALSE) {
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL){
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, x = x)
out <- with(args, {
z <- (x - xi)/omega
if (length(alpha) == 1L) {
alpha <- rep(alpha, length(z))
}
logN <- -log(sqrt(2 * pi)) - log(omega) - z^2/2
logS <- ifelse(abs(alpha) < Inf, pnorm(alpha * z, log.p = TRUE),
log(as.numeric(sign(alpha) * z > 0)))
out <- logN + logS - pnorm(0, log.p = TRUE)
ifelse(abs(z) == Inf, -Inf, out)
})
if (!log) {
out <- exp(out)
}
out
}
#' Skew Normal cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' psnorm()
psnorm <- function (q, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL,
lower.tail = TRUE, log.p = FALSE)
{
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL) {
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, q = q)
out <- with(args, {
z <- (q - xi)/omega
nz <- length(z)
is_alpha_inf <- abs(alpha) == Inf
delta[is_alpha_inf] <- sign(alpha[is_alpha_inf])
out <- numeric(nz)
for (k in seq_len(nz)) {
if (is_alpha_inf[k]) {
if (alpha[k] > 0) {
out[k] <- 2 * (pnorm(pmax(z[k], 0)) - 0.5)
}
else {
out[k] <- 1 - 2 * (0.5 - pnorm(pmin(z[k], 0)))
}
}
else {
S <- matrix(c(1, -delta[k], -delta[k], 1), 2,
2)
out[k] <- 2 * mnormt::biv.nt.prob(0, lower = rep(-Inf, 2), upper = c(z[k], 0), mean = c(0, 0), S = S)
}
}
pmin(1, pmax(0, out))
})
if (!lower.tail) {
out <- 1 - out
}
if (log.p) {
out <- log(out)
}
out
}
#' Skew Normal quantile function
#'
#' @param p vector of probabilities
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#' @param tol tolerance for small values to be treated as equivalent to zero
#'
#' @return vector
#' @export
#'
#' @examples
#' qsnorm()
#'
qsnorm <- function (p, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, lower.tail = TRUE, log.p = FALSE, tol = 1e-08){
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
if (log.p) {
p <- exp(p)
}
if (!lower.tail) {
p <- 1 - p
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL){
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...) {
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, p = p)
snorm_cumulants <-function (xi = 0, omega = 1, alpha = 0, n = 4){
cumulants_half_norm <- function(n) {
n <- max(n, 2)
n <- as.integer(2 * ceiling(n/2))
half.n <- as.integer(n/2)
m <- 0:(half.n - 1)
a <- sqrt(2/pi)/(gamma(m + 1) * 2^m * (2 * m + 1))
signs <- rep(c(1, -1), half.n)[seq_len(half.n)]
a <- as.vector(rbind(signs * a, rep(0, half.n)))
coeff <- rep(a[1], n)
for (k in 2:n) {
ind <- seq_len(k - 1)
coeff[k] <- a[k] - sum(ind * coeff[ind] * a[rev(ind)]/k)
}
kappa <- coeff * gamma(seq_len(n) + 1)
kappa[2] <- 1 + kappa[2]
return(kappa)
}
expandfunc <- function (..., dots = list(), length = NULL)
{
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
args <- expandfunc(dots = nlist(xi, omega, alpha))
with(args, {
delta <- alpha/sqrt(1 + alpha^2)
kv <- cumulants_half_norm(n)
if (length(kv) > n) {
kv <- kv[-(n + 1)]
}
kv[2] <- kv[2] - 1
kappa <- outer(delta, 1:n, "^") * matrix(rep(kv, length(xi)),
ncol = n, byrow = TRUE)
kappa[, 2] <- kappa[, 2] + 1
kappa <- kappa * outer(omega, 1:n, "^")
kappa[, 1] <- kappa[, 1] + xi
kappa
})
}
out <- with(args, {
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
p <- replace(p, (na | zero | one), 0.5)
cum <- snorm_cumulants(0, 1, alpha, n = 4)
g1 <- cum[, 3]/cum[, 2]^(3/2)
g2 <- cum[, 4]/cum[, 2]^2
x <- qnorm(p)
x <- x + (x^2 - 1) * g1/6 + x * (x^2 - 3) * g2/24 - x *
(2 * x^2 - 5) * g1^2/36
x <- cum[, 1] + sqrt(cum[, 2]) * x
px <- psnorm(x, xi = 0, omega = 1, alpha = alpha)
max_err <- 1
while (max_err > tol) {
x1 <- x - (px - p)/dsnorm(x, xi = 0, omega = 1,
alpha = alpha)
x <- x1
px <- psnorm(x, xi = 0, omega = 1, alpha = alpha)
max_err <- max(abs(px - p))
if (is.na(max_err)) {
warning("Approximation in 'qsnorm' may have failed.")
}
}
x <- replace(x, na, NA)
x <- replace(x, zero, -Inf)
x <- replace(x, one, Inf)
as.numeric(xi + omega * x)
})
out
}
###########################################################################
############ Half Student's T Distribution ############
###########################################################################
#' Half Student's-Tdistribution random number generator
#'
#' @param n number of observations
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#'
#'
#' @return vector
#' @export
#'
#' @examples
#' rhalft()
#'
rhalft <- function(n, nu = 3, scale=1){
stopifnot(scale>0, nu>0)
return(abs(rst(n=n, nu=nu, mu = 0, scale = 1))*scale)
}
#' Half Student's-T distribution probability density function
#'
#' @param n number of observations
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dhalft()
#'
dhalft <- function (x, nu = 3, scale = 1, log = FALSE){
x <- as.vector(x)
scale <- as.vector(scale)
nu <- as.vector(nu)
if (any(scale <= 0))
stop("The scale parameter must be positive.")
NN <- max(length(x), length(scale), length(nu))
x <- rep(x, len = NN)
scale <- rep(scale, len = NN)
nu <- rep(nu, len = NN)
dens <- log(2) - log(scale) + lgamma((nu + 1)/2) - lgamma(nu/2) - 0.5 * log(pi * nu) - (nu + 1)/2 * log(1 + (1/nu) * (x/scale) * (x/scale))
if (log == FALSE)
dens <- exp(dens)
return(dens)
}
#' Half Student's-T distribution cumulative density function
#'
#' @param q vector of quantiles
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' phalft()
#'
phalft <- function(q, nu = 3, scale=1, lower.tail=TRUE, log.p=FALSE){
stopifnot(scale>0, nu>0)
result <- rep(0, length(q))
result[q>0] <- 2.0 * (pst(q[q>0]/scale, nu=nu) - 0.5)
if (!lower.tail) result <- 1 - result
if(log.p) return(log(result)) else return(result)
}
#' Half Student's-T distribution inverse cumulative density function
#'
#' @param p vector of probabilities.
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qhalft()
#'
qhalft <- function (p, nu = 3, scale = 1, lower.tail=TRUE,log.p=FALSE){
p <- as.vector(p)
scale <- as.vector(scale)
nu <- as.vector(nu)
if (any(p < 0) || any(p > 1))
stop("p must be in [0,1].")
if (any(scale <= 0))
stop("The scale parameter must be positive.")
NN <- max(length(p), length(scale), length(nu))
p <- rep(p, len = NN)
scale <- rep(scale, len = NN)
q <- qcustom(phalft, p=p, nu = nu, scale = scale, lower.tail=lower.tail,log.p=log.p)
return(q)
}
###########################################################################
############ Truncated Normal Distribution ############
###########################################################################
#' Truncated Gaussian probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper trunctation limit
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dtnorm()
#'
dtnorm <- function(x, mu=0, sigma=1, lower=-Inf, upper=Inf, log=FALSE)
{
ret <- numeric(length(x))
ret[x < lower | x > upper] <- if (log) -Inf else 0
ret[upper < lower] <- NaN
ind <- x >=lower & x <=upper
if (any(ind)) {
denom <- pnorm(upper, mu, sigma) - pnorm(lower, mu, sigma)
xtmp <- dnorm(x, mu, sigma, log)
if (log) xtmp <- xtmp - log(denom) else xtmp <- xtmp/denom
ret[x >=lower & x <=upper] <- xtmp[ind]
}
ret
}
#' Half Student's-T cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper truncation limit
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' ptnorm()
ptnorm <- function(q, mu=0, sigma=1, lower=-Inf, upper=Inf, lower.tail=TRUE, log.p=FALSE)
{
ret <- numeric(length(q))
if (lower.tail) {
ret[q < lower] <- 0
ret[q > upper] <- 1
}
else {
ret[q < lower] <- 1
ret[q > upper] <- 0
}
ret[upper < lower] <- NaN
ind <- q >=lower & q <=upper
if (any(ind)) {
denom <- pnorm(upper, mu, sigma) - pnorm(lower, mu, sigma)
if (lower.tail) qtmp <- pnorm(q, mu, sigma) - pnorm(lower, mu, sigma)
else qtmp <- pnorm(upper, mu, sigma) - pnorm(q, mu, sigma)
if (log.p) qtmp <- log(qtmp) - log(denom) else qtmp <- qtmp/denom
ret[q >=lower & q <=upper] <- qtmp[ind]
}
ret
}
#' Truncated Gaussian quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper truncation limit
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qtnorm()
qtnorm <- function(p, mu=0, sigma=1, lower=-Inf, upper=Inf, lower.tail=TRUE, log.p=FALSE)
{
qcustom(ptnorm, p=p, mu=mu, sigma=sigma, lower=lower, upper=upper, lbound=lower, ubound=upper, lower.tail=lower.tail, log.p=log.p)
}
#' Truncated Gaussian random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper trunctation limit
#'
#' @return vector
#' @export
#'
#' @examples
#' rtnorm()
#'
rtnorm <- function (n, mu = 0, sigma = 1, lower = -Inf, upper = Inf) {
if (length(n) > 1)
n <- length(n)
mu <- rep(mu, length=n)
sigma <- rep(sigma, length=n)
lower <- rep(lower, length=n)
upper <- rep(upper, length=n)
lower <- (lower - mu) / sigma ## Algorithm works on mu 0, sigma 1 scale
upper <- (upper - mu) / sigma
ind <- seq(length=n)
ret <- numeric(n)
nas <- is.na(mu) | is.na(sigma) | is.na(lower) | is.na(upper)
if (any(nas)) warning("NAs produced")
## Different algorithms depending on where upper/lower limits lie.
alg <- ifelse(
((lower > upper) | nas),
-1,# return NaN
ifelse(
((lower < 0 & upper == Inf) |
(lower == -Inf & upper > 0) |
(is.finite(lower) & is.finite(upper) & (lower < 0) & (upper > 0) & (upper-lower > sqrt(2*pi)))
),
0, # standard "simulate from normal and reject if outside limits" method. Use if bounds are wide.
ifelse(
(lower >= 0 & (upper > lower + 2*sqrt(exp(1)) /
(lower + sqrt(lower^2 + 4)) * exp((lower*2 - lower*sqrt(lower^2 + 4)) / 4))),
1, # rejection sampling with exponential proposal. Use if lower >> mu
ifelse(upper <= 0 & (-lower > -upper + 2*sqrt(exp(1)) /
(-upper + sqrt(upper^2 + 4)) * exp((upper*2 - -upper*sqrt(upper^2 + 4)) / 4)),
2, # rejection sampling with exponential proposal. Use if upper << mu.
3)))) # rejection sampling with uniform proposal. Use if bounds are narrow and central.
ind.nan <- ind[alg==-1]; ind.no <- ind[alg==0]; ind.expl <- ind[alg==1]; ind.expu <- ind[alg==2]; ind.u <- ind[alg==3]
ret[ind.nan] <- NaN
while (length(ind.no) > 0) {
y <- rnorm(length(ind.no))
done <- which(y >= lower[ind.no] & y <= upper[ind.no])
ret[ind.no[done]] <- y[done]
ind.no <- setdiff(ind.no, ind.no[done])
}
stopifnot(length(ind.no) == 0)
while (length(ind.expl) > 0) {
a <- (lower[ind.expl] + sqrt(lower[ind.expl]^2 + 4)) / 2
z <- rexp(length(ind.expl), a) + lower[ind.expl]
u <- runif(length(ind.expl))
done <- which((u <= exp(-(z - a)^2 / 2)) & (z <= upper[ind.expl]))
ret[ind.expl[done]] <- z[done]
ind.expl <- setdiff(ind.expl, ind.expl[done])
}
stopifnot(length(ind.expl) == 0)
while (length(ind.expu) > 0) {
a <- (-upper[ind.expu] + sqrt(upper[ind.expu]^2 +4)) / 2
z <- rexp(length(ind.expu), a) - upper[ind.expu]
u <- runif(length(ind.expu))
done <- which((u <= exp(-(z - a)^2 / 2)) & (z <= -lower[ind.expu]))
ret[ind.expu[done]] <- -z[done]
ind.expu <- setdiff(ind.expu, ind.expu[done])
}
stopifnot(length(ind.expu) == 0)
while (length(ind.u) > 0) {
z <- runif(length(ind.u), lower[ind.u], upper[ind.u])
rho <- ifelse(lower[ind.u] > 0,
exp((lower[ind.u]^2 - z^2) / 2), ifelse(upper[ind.u] < 0,
exp((upper[ind.u]^2 - z^2) / 2),
exp(-z^2/2)))
u <- runif(length(ind.u))
done <- which(u <= rho)
ret[ind.u[done]] <- z[done]
ind.u <- setdiff(ind.u, ind.u[done])
}
stopifnot(length(ind.u) == 0)
ret*sigma + mu
}
#' Slash Distribution inverse cumulative density function
#'
#' @name dslash
#' @title Slash distribution inverse cumulative density function
#' @description Inverse cumulative density function for the slash distribution.
#' @param x vector of quantiles
#' @param location vector of location parameters
#' @param scale vector of scale parameters
#' @param log if TRUE, the probabilites are given as log(p).
#' @returns a numeric vector
#' @export
#'
#'
#' @examples
#' qslash()
qslash <- function(p, location=0, scale=1, lower.tail=TRUE, log=FALSE)
{
qcustom(pslash, p=p, location=location, scale=scale, lower_tail=lower.tail, log=log)
}
qcustom <- function(pdist, p, ...){
args <- list(...)
if (!is.null(args$log)) args$log.p <- args$log
if (is.null(args$log.p)) args$log.p <- FALSE
if (is.null(args$lower.tail)) args$lower.tail <- TRUE
if (is.null(args$lbound)) args$lbound <- -Inf
if (is.null(args$ubound)) args$ubound <- Inf
if (args$log.p) p <- exp(p)
if (!args$lower.tail) p <- 1 - p
ret <- numeric(length(p))
ret[p == 0] <- args$lbound
ret[p == 1] <- args$ubound
args[c("lower.tail","log.p","lbound","ubound")] <- NULL
## Other args assumed to contain params of the distribution.
## Replicate all to their maximum length, along with p
maxlen <- max(sapply(c(args, p=list(p)), length))
for (i in seq(along=args))
args[[i]] <- rep(args[[i]], length.out=maxlen)
p <- rep(p, length.out=maxlen)
ret[p < 0 | p > 1] <- NaN
ind <- (p > 0 & p < 1)
if (any(ind)) {
hind <- seq(along=p)[ind]
h <- function(y) {
args <- lapply(args, function(x)x[hind[i]])
p <- p[hind[i]]
args$q <- y
(do.call(pdist, args) - p)
}
ptmp <- numeric(length(p[ind]))
for (i in 1:length(p[ind])) {
interval <- c(-1, 1)
while (h(interval[1])*h(interval[2]) >= 0) {
interval <- interval + c(-1,1)*0.5*(interval[2]-interval[1])
}
ptmp[i] <- uniroot(h, interval, tol=.Machine$double.eps)$root
}
ret[ind] <- ptmp
}
if (any(is.nan(ret))) warning("NaNs produced")
ret
}
abnormally-distributed/cvreg documentation built on May 3, 2020, 3:45 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions qcustom qslash qtnorm ptnorm dtnorm phalft rhalft rsnorm ral qal pal dal
Documented in dal dtnorm pal phalft ptnorm qal qslash qtnorm ral rhalft rsnorm
###########################################################################
######### Power Exponential (Generalized Normal) Distribution) ##########
###########################################################################
#' Power Exponential random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#'
#' @return vector
#' @export
#'
#' @examples
#' rpowexp()
rpowexp = function (n, mu = 0, sigma = 1, nu = 2)
{
if (any(sigma <= 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(n <= 0))
stop(paste("n must be a positive integer", "\n", ""))
n <- ceiling(n)
p <- runif(n)
r <- qpowexp(p, mu = mu, sigma = sigma, nu = nu)
r
}
#' Power Exponential quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qpowexp()
qpowexp = function (p, mu = 0, sigma = 1, nu = 2, lower.tail = TRUE, log = FALSE)
{
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
if (log == TRUE)
p <- exp(p)
else p <- p
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
if (any(p < 0) | any(p > 1))
stop(paste("p must be between 0 and 1", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
suppressWarnings(s <- qgamma((2 * p - 1) * sign(p - 0.5), shape = (1/nu), scale = 1))
z <- sign(p - 0.5) * ((2 * s)^(1/nu)) * c
ya <- mu + sigma * z
ya
}
#' Power Exponential probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dpowexp()
dpowexp = function (x, mu = 0, sigma = 1, nu = 2, log = FALSE) {
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
z <- (x - mu)/sigma
log.lik <- -log(sigma) + log(nu) - log.c - (0.5 * (abs(z/c)^nu)) -
(1 + (1/nu)) * log(2) - lgamma(1/nu)
if (log == FALSE)
fy <- exp(log.lik)
else fy <- log.lik
fy
}
#' Power Exponential cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param nu vector of shape parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' ppowexp()
ppowexp = function (q, mu = 0, sigma = 1, nu = 2, lower.tail = TRUE, log = FALSE) {
if (any(sigma < 0))
stop(paste("the sigma parameter must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("the nu parameter must be positive", "\n", ""))
log.c <- 0.5 * (-(2/nu) * log(2) + lgamma(1/nu) - lgamma(3/nu))
c <- exp(log.c)
z <- (q - mu)/sigma
s <- 0.5 * ((abs(z/c))^nu)
cdf <- 0.5 * (1 + pgamma(s, shape = 1/nu, scale = 1) * sign(z))
if (length(nu) > 1)
cdf <- ifelse(nu > 10000, (q - (mu - sqrt(3) * sigma))/(sqrt(12) * sigma), cdf)
else cdf <- if (nu > 10000)
(q - (mu - sqrt(3) * sigma))/(sqrt(12) * sigma)
else cdf
if (lower.tail == TRUE)
cdf <- cdf
else cdf <- 1 - cdf
if (log == FALSE)
cdf <- cdf
else cdf <- log(cdf)
cdf
}
###########################################################################
######### Student's T Distribution ##########
###########################################################################
#' Student-T random number generator
#'
#' @param n number of observations
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#'
#' @return vector
#' @export
#'
#' @examples
#' rst()
rst = function (n, nu = 3, mu = 0, scale = 1) {
if (any(scale <= 0))
stop("the scale parameter must be positive.")
if (any(nu <= 0))
stop("the nu parameter must be positive.")
rGamma = function(n, shape = 1, rate = 1) {
return(qgamma(seq(0.0001, 0.9999, length.out = n),
shape, rate))
}
gamma.samps = rGamma(n, nu * 0.50, (scale^2) * (nu * 0.50))
sigmas = sqrt(1 / gamma.samps)
x = rnorm(n, mu, sigmas)
return(x)
}
#' Student-T probability density function
#'
#' @param x vector of quantiles
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dst()
dst = function (x, nu = 3, mu = 0, scale = 1, log = FALSE)
{
if (any(scale <= 0)) {
stop("the scale parameter must be positive.")
}
if (any(nu <= 0)) {
stop("the nu parameter must be positive.")
}
if (log) {
dt((x - mu)/scale, df = nu, log = TRUE) - log(scale)
}
else {
dt((x - mu)/scale, df = nu)/scale
}
}
#' Student-T quantile function
#'
#' @param p vector of probabilities.
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qst()
qst = function (p, nu = 3, mu = 0, scale = 1, lower.tail = TRUE, log = FALSE)
{
if (any(p < 0) || any(p > 1))
stop("p must be in [0,1].")
if (any(scale <= 0))
stop("the scale parameter must be positive.")
if (any(nu <= 0))
stop("the nu parameter must be positive.")
NN <- max(length(p), length(mu), length(scale), length(nu))
p <- rep(p, len = NN)
mu <- rep(mu, len = NN)
scale <- rep(scale, len = NN)
nu <- rep(nu, len = NN)
q <- mu + scale * qt(p, df = nu, lower.tail = lower.tail)
temp <- which(nu > 1e+06)
q[temp] <- qnorm(p[temp], mu[temp], scale[temp], lower.tail = lower.tail,
log.p = log)
return(q)
}
#' Student-T cumulative probability function
#'
#' @param q vector of quantiles
#' @param nu vector of normality parameter values
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' pst()
pst = function (q, nu = 3, mu = 0, scale = 1, lower.tail = TRUE, log = FALSE)
{
if (any(scale <= 0)) {
stop("the scale parameter must be positive.")
}
if (any(nu <= 0)) {
stop("the nu parameter must be positive.")
}
pt((q - mu)/scale, df = nu, lower.tail = lower.tail, log.p = log)
}
###########################################################################
######### Asymmetric Laplace Distribution ##########
###########################################################################
#' Asymmetric Laplace probability density function
#'
#' @param x vector of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
dal <- function(x, location=0, scale=1, tau=0.5, log=FALSE){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
x <- as.vector(x); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(x), length(location), length(scale), length(kappa))
x <- rep(x, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- rep(kappa, len=NN)
logconst <- 0.5 * log(2) - log(scale) + log(kappa) - log1p(kappa^2)
temp <- which(x < location); kappa[temp] <- 1/kappa[temp]
exponent <- -(sqrt(2) / scale) * abs(x - location) * kappa
dens <- logconst + exponent
if(log == FALSE) dens <- exp(logconst + exponent)
return(dens)
}
#' Asymmetric Laplace cumulative probability function
#'
#' @param q a vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
pal <- function(q, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
q <- as.vector(q); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if((kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(q), length(location), length(scale), length(kappa))
q <- rep(q, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- k2 <- rep(kappa, len=NN)
temp <- which(q < location); k2[temp] <- 1/kappa[temp]
exponent <- -(sqrt(2) / scale) * abs(q - location) * k2
temp <- exp(exponent) / (1 + kappa^2)
p <- 1 - temp
index1 <- (q < location)
p[index1] <- (kappa[index1])^2 * temp[index1]
return(p)
}
#' Asymmetric Laplace quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
qal <- function(p, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
p <- as.vector(p); location <- as.vector(location)
scale <- as.vector(scale); kappa <- as.vector(kappa)
if(any(p < 0) || any(p > 1)) stop("p must be in [0,1].")
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
NN <- max(length(p), length(location), length(scale), length(kappa))
p <- rep(p, len=NN); location <- rep(location, len=NN)
scale <- rep(scale, len=NN); kappa <- rep(kappa, len=NN)
q <- p
temp <- kappa^2 / (1 + kappa^2)
index1 <- (p <= temp)
exponent <- p[index1] / temp[index1]
q[index1] <- location[index1] + (scale[index1] * kappa[index1]) *
log(exponent) / sqrt(2)
q[!index1] <- location[!index1] - (scale[!index1] / kappa[!index1]) *
(log1p((kappa[!index1])^2) + log1p(-p[!index1])) / sqrt(2)
q[p == 0] = -Inf
q[p == 1] = Inf
return(q)
}
#' Asymmetric Laplace random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param tau a vector of quantile cutpoints; value(s) between 0 and 1, where 0.5 corresponds to the median, for example.
#'
#' @return vector
#' @export
#'
ral <- function(n, location=0, scale=1, tau=0.5){
tau<-min(0.99, max(tau, 0.001));kappa<-(tau)/(1-tau)
location <- rep(location, len=n)
scale <- rep(scale, len=n)
kappa <- rep(kappa, len=n)
if(any(scale <= 0)) stop("The scale parameter must be positive.")
if(any(kappa<0)) stop("The kappa parameter must be positive.")
x <- location + scale * log(runif(n)^kappa / runif(n)^(1/kappa)) / sqrt(2)
return(x)
}
###########################################################################
######### Sinh-Arcsinh (SHASH) Distribution ##########
###########################################################################
#' Sinh-Arcsinh (SHASH) random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#'
#' @return vector
#' @export
#'
#' @examples
#' rshash()
rshash = function (n, mu = 0, scale = 1, alpha = 0, nu = 1)
{
if (any(n <= 0))
stop(paste("n must be a positive integer", "\n", ""))
n <- ceiling(n)
p <- runif(n)
r <- qshash(p, mu = mu, scale = scale, alpha = alpha, nu = nu)
r
}
#' Sinh-Arcsinh (SHASH) probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dshash()
dshash = function (x, mu = 0, scale = 1, alpha = 0, nu = 1, log = FALSE){
if (any(scale < 0))
stop(paste("scale must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("nu must be positive", "\n", ""))
z <- (x - mu)/(scale * nu)
c <- cosh(nu * asinh(z) - alpha)
r <- sinh(nu * asinh(z) - alpha)
loglik <- -log(scale) - 0.5 * log(2 * pi) - 0.5 * log(1 + (z^2)) + log(c) - 0.5 * (r^2)
if (log == FALSE)
fy <- exp(loglik)
else fy <- loglik
fy
}
#' Sinh-Arcsinh (SHASH) cumulative probability function
#'
#' @param v vector of quantiles
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' pshash()
pshash = function (q, mu = 0, scale = 1, alpha = 0, nu = 1, lower.tail = TRUE,
log.p = FALSE)
{
if (any(scale < 0))
stop(paste("scale must be positive", "\n", ""))
if (any(nu < 0))
stop(paste("nu must be positive", "\n", ""))
z <- (q - mu)/(scale * nu)
c <- cosh(nu * asinh(z) - alpha)
r <- sinh(nu * asinh(z) - alpha)
p <- pNO(r)
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
if (log.p == FALSE)
p <- p
else p <- log(p)
p
}
#' Sinh-Arcsinh (SHASH) quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param scale vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param nu vector of kurtosis parameter
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qshash()
qshash = function (p, mu = 0, scale = 1, alpha = 0, nu = 1, lower.tail = TRUE, log.p = FALSE)
{
if (log.p == TRUE)
p <- exp(p)
else p <- p
if (any(p <= 0) | any(p >= 1))
stop(paste("p must be between 0 and 1", "\n", ""))
if (lower.tail == TRUE)
p <- p
else p <- 1 - p
y <- mu + (nu * scale) * sinh((1/nu) * asinh(qnorm(p)) +
(alpha/nu))
y
}
###########################################################################
############## Skew-Normal Distribution ##############
###########################################################################
#' Skew Normal random number generator
#'
#' @param n number of observations to generate
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#'
#' @return vector
#' @export
#'
#' @examples
#' rsnorm()
rsnorm <- function(n, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL){
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL)
{
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega)
with(args, {
z1 <- rnorm(n)
z2 <- rnorm(n)
id <- z2 > args$alpha * z1
z1[id] <- -z1[id]
xi + omega * z1
})
}
#' Skew Normal probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param log if TRUE, densities are returned as log(density)
#'
#' @return vector
#' @export
#'
#' @examples
#' dsnorm()
dsnorm <- function (x, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL,
log = FALSE) {
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL){
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, x = x)
out <- with(args, {
z <- (x - xi)/omega
if (length(alpha) == 1L) {
alpha <- rep(alpha, length(z))
}
logN <- -log(sqrt(2 * pi)) - log(omega) - z^2/2
logS <- ifelse(abs(alpha) < Inf, pnorm(alpha * z, log.p = TRUE),
log(as.numeric(sign(alpha) * z > 0)))
out <- logN + logS - pnorm(0, log.p = TRUE)
ifelse(abs(z) == Inf, -Inf, out)
})
if (!log) {
out <- exp(out)
}
out
}
#' Skew Normal cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' psnorm()
psnorm <- function (q, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL,
lower.tail = TRUE, log.p = FALSE)
{
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL) {
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, q = q)
out <- with(args, {
z <- (q - xi)/omega
nz <- length(z)
is_alpha_inf <- abs(alpha) == Inf
delta[is_alpha_inf] <- sign(alpha[is_alpha_inf])
out <- numeric(nz)
for (k in seq_len(nz)) {
if (is_alpha_inf[k]) {
if (alpha[k] > 0) {
out[k] <- 2 * (pnorm(pmax(z[k], 0)) - 0.5)
}
else {
out[k] <- 1 - 2 * (0.5 - pnorm(pmin(z[k], 0)))
}
}
else {
S <- matrix(c(1, -delta[k], -delta[k], 1), 2,
2)
out[k] <- 2 * mnormt::biv.nt.prob(0, lower = rep(-Inf, 2), upper = c(z[k], 0), mean = c(0, 0), S = S)
}
}
pmin(1, pmax(0, out))
})
if (!lower.tail) {
out <- 1 - out
}
if (log.p) {
out <- log(out)
}
out
}
#' Skew Normal quantile function
#'
#' @param p vector of probabilities
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param alpha vector of skewness parameter
#' @param xi for internal use
#' @param omega for internal use
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#' @param tol tolerance for small values to be treated as equivalent to zero
#'
#' @return vector
#' @export
#'
#' @examples
#' qsnorm()
#'
qsnorm <- function (p, mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, lower.tail = TRUE, log.p = FALSE, tol = 1e-08){
if (any(sigma <= 0)) {
stop("sigma must be greater than 0.")
}
if (log.p) {
p <- exp(p)
}
if (!lower.tail) {
p <- 1 - p
}
cp2dp <- function(mu = 0, sigma = 1, alpha = 0, xi = NULL, omega = NULL, ...) {
delta <- alpha/sqrt(1 + alpha^2)
if (is.null(omega)) {
omega <- sigma/sqrt(1 - 2/pi * delta^2)
}
if (is.null(xi)) {
xi <- mu - omega * delta * sqrt(2/pi)
}
expandfunc <- function (..., dots = list(), length = NULL){
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...) {
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
expandfunc(dots = nlist(mu, sigma, alpha, xi, omega, delta, ...))
}
args <- cp2dp(mu, sigma, alpha, xi = xi, omega = omega, p = p)
snorm_cumulants <-function (xi = 0, omega = 1, alpha = 0, n = 4){
cumulants_half_norm <- function(n) {
n <- max(n, 2)
n <- as.integer(2 * ceiling(n/2))
half.n <- as.integer(n/2)
m <- 0:(half.n - 1)
a <- sqrt(2/pi)/(gamma(m + 1) * 2^m * (2 * m + 1))
signs <- rep(c(1, -1), half.n)[seq_len(half.n)]
a <- as.vector(rbind(signs * a, rep(0, half.n)))
coeff <- rep(a[1], n)
for (k in 2:n) {
ind <- seq_len(k - 1)
coeff[k] <- a[k] - sum(ind * coeff[ind] * a[rev(ind)]/k)
}
kappa <- coeff * gamma(seq_len(n) + 1)
kappa[2] <- 1 + kappa[2]
return(kappa)
}
expandfunc <- function (..., dots = list(), length = NULL)
{
dots <- c(dots, list(...))
max_dim <- NULL
if (is.null(length)) {
lengths <- lengths(dots)
length <- max(lengths)
max_dim <- dim(dots[[match(length, lengths)]])
}
out <- as.data.frame(lapply(dots, rep, length.out = length))
structure(out, max_dim = max_dim)
}
nlist <- function (...)
{
m <- match.call()
dots <- list(...)
no_names <- is.null(names(dots))
has_name <- if (no_names)
FALSE
else nzchar(names(dots))
if (all(has_name))
return(dots)
nms <- as.character(m)[-1]
if (no_names) {
names(dots) <- nms
}
else {
names(dots)[!has_name] <- nms[!has_name]
}
dots
}
args <- expandfunc(dots = nlist(xi, omega, alpha))
with(args, {
delta <- alpha/sqrt(1 + alpha^2)
kv <- cumulants_half_norm(n)
if (length(kv) > n) {
kv <- kv[-(n + 1)]
}
kv[2] <- kv[2] - 1
kappa <- outer(delta, 1:n, "^") * matrix(rep(kv, length(xi)),
ncol = n, byrow = TRUE)
kappa[, 2] <- kappa[, 2] + 1
kappa <- kappa * outer(omega, 1:n, "^")
kappa[, 1] <- kappa[, 1] + xi
kappa
})
}
out <- with(args, {
na <- is.na(p) | (p < 0) | (p > 1)
zero <- (p == 0)
one <- (p == 1)
p <- replace(p, (na | zero | one), 0.5)
cum <- snorm_cumulants(0, 1, alpha, n = 4)
g1 <- cum[, 3]/cum[, 2]^(3/2)
g2 <- cum[, 4]/cum[, 2]^2
x <- qnorm(p)
x <- x + (x^2 - 1) * g1/6 + x * (x^2 - 3) * g2/24 - x *
(2 * x^2 - 5) * g1^2/36
x <- cum[, 1] + sqrt(cum[, 2]) * x
px <- psnorm(x, xi = 0, omega = 1, alpha = alpha)
max_err <- 1
while (max_err > tol) {
x1 <- x - (px - p)/dsnorm(x, xi = 0, omega = 1,
alpha = alpha)
x <- x1
px <- psnorm(x, xi = 0, omega = 1, alpha = alpha)
max_err <- max(abs(px - p))
if (is.na(max_err)) {
warning("Approximation in 'qsnorm' may have failed.")
}
}
x <- replace(x, na, NA)
x <- replace(x, zero, -Inf)
x <- replace(x, one, Inf)
as.numeric(xi + omega * x)
})
out
}
###########################################################################
############ Half Student's T Distribution ############
###########################################################################
#' Half Student's-Tdistribution random number generator
#'
#' @param n number of observations
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#'
#'
#' @return vector
#' @export
#'
#' @examples
#' rhalft()
#'
rhalft <- function(n, nu = 3, scale=1){
stopifnot(scale>0, nu>0)
return(abs(rst(n=n, nu=nu, mu = 0, scale = 1))*scale)
}
#' Half Student's-T distribution probability density function
#'
#' @param n number of observations
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dhalft()
#'
dhalft <- function (x, nu = 3, scale = 1, log = FALSE){
x <- as.vector(x)
scale <- as.vector(scale)
nu <- as.vector(nu)
if (any(scale <= 0))
stop("The scale parameter must be positive.")
NN <- max(length(x), length(scale), length(nu))
x <- rep(x, len = NN)
scale <- rep(scale, len = NN)
nu <- rep(nu, len = NN)
dens <- log(2) - log(scale) + lgamma((nu + 1)/2) - lgamma(nu/2) - 0.5 * log(pi * nu) - (nu + 1)/2 * log(1 + (1/nu) * (x/scale) * (x/scale))
if (log == FALSE)
dens <- exp(dens)
return(dens)
}
#' Half Student's-T distribution cumulative density function
#'
#' @param q vector of quantiles
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' phalft()
#'
phalft <- function(q, nu = 3, scale=1, lower.tail=TRUE, log.p=FALSE){
stopifnot(scale>0, nu>0)
result <- rep(0, length(q))
result[q>0] <- 2.0 * (pst(q[q>0]/scale, nu=nu) - 0.5)
if (!lower.tail) result <- 1 - result
if(log.p) return(log(result)) else return(result)
}
#' Half Student's-T distribution inverse cumulative density function
#'
#' @param p vector of probabilities.
#' @param nu vector of kurtosis parameter values
#' @param scale vector of scale parameter values
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qhalft()
#'
qhalft <- function (p, nu = 3, scale = 1, lower.tail=TRUE,log.p=FALSE){
p <- as.vector(p)
scale <- as.vector(scale)
nu <- as.vector(nu)
if (any(p < 0) || any(p > 1))
stop("p must be in [0,1].")
if (any(scale <= 0))
stop("The scale parameter must be positive.")
NN <- max(length(p), length(scale), length(nu))
p <- rep(p, len = NN)
scale <- rep(scale, len = NN)
q <- qcustom(phalft, p=p, nu = nu, scale = scale, lower.tail=lower.tail,log.p=log.p)
return(q)
}
###########################################################################
############ Truncated Normal Distribution ############
###########################################################################
#' Truncated Gaussian probability density function
#'
#' @param x vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper trunctation limit
#' @param log if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' dtnorm()
#'
dtnorm <- function(x, mu=0, sigma=1, lower=-Inf, upper=Inf, log=FALSE)
{
ret <- numeric(length(x))
ret[x < lower | x > upper] <- if (log) -Inf else 0
ret[upper < lower] <- NaN
ind <- x >=lower & x <=upper
if (any(ind)) {
denom <- pnorm(upper, mu, sigma) - pnorm(lower, mu, sigma)
xtmp <- dnorm(x, mu, sigma, log)
if (log) xtmp <- xtmp - log(denom) else xtmp <- xtmp/denom
ret[x >=lower & x <=upper] <- xtmp[ind]
}
ret
}
#' Half Student's-T cumulative probability function
#'
#' @param q vector of quantiles
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper truncation limit
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' ptnorm()
ptnorm <- function(q, mu=0, sigma=1, lower=-Inf, upper=Inf, lower.tail=TRUE, log.p=FALSE)
{
ret <- numeric(length(q))
if (lower.tail) {
ret[q < lower] <- 0
ret[q > upper] <- 1
}
else {
ret[q < lower] <- 1
ret[q > upper] <- 0
}
ret[upper < lower] <- NaN
ind <- q >=lower & q <=upper
if (any(ind)) {
denom <- pnorm(upper, mu, sigma) - pnorm(lower, mu, sigma)
if (lower.tail) qtmp <- pnorm(q, mu, sigma) - pnorm(lower, mu, sigma)
else qtmp <- pnorm(upper, mu, sigma) - pnorm(q, mu, sigma)
if (log.p) qtmp <- log(qtmp) - log(denom) else qtmp <- qtmp/denom
ret[q >=lower & q <=upper] <- qtmp[ind]
}
ret
}
#' Truncated Gaussian quantile function
#'
#' @param p vector of probabilities.
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper truncation limit
#' @param lower.tail if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#' @param log.p if TRUE, probabilities p are given as log(p).
#'
#' @return vector
#' @export
#'
#' @examples
#' qtnorm()
qtnorm <- function(p, mu=0, sigma=1, lower=-Inf, upper=Inf, lower.tail=TRUE, log.p=FALSE)
{
qcustom(ptnorm, p=p, mu=mu, sigma=sigma, lower=lower, upper=upper, lbound=lower, ubound=upper, lower.tail=lower.tail, log.p=log.p)
}
#' Truncated Gaussian random number generator
#'
#' @param n number of observations
#' @param mu vector of location parameter values
#' @param sigma vector of scale parameter values
#' @param lower lower truncation limit
#' @param upper upper trunctation limit
#'
#' @return vector
#' @export
#'
#' @examples
#' rtnorm()
#'
rtnorm <- function (n, mu = 0, sigma = 1, lower = -Inf, upper = Inf) {
if (length(n) > 1)
n <- length(n)
mu <- rep(mu, length=n)
sigma <- rep(sigma, length=n)
lower <- rep(lower, length=n)
upper <- rep(upper, length=n)
lower <- (lower - mu) / sigma ## Algorithm works on mu 0, sigma 1 scale
upper <- (upper - mu) / sigma
ind <- seq(length=n)
ret <- numeric(n)
nas <- is.na(mu) | is.na(sigma) | is.na(lower) | is.na(upper)
if (any(nas)) warning("NAs produced")
## Different algorithms depending on where upper/lower limits lie.
alg <- ifelse(
((lower > upper) | nas),
-1,# return NaN
ifelse(
((lower < 0 & upper == Inf) |
(lower == -Inf & upper > 0) |
(is.finite(lower) & is.finite(upper) & (lower < 0) & (upper > 0) & (upper-lower > sqrt(2*pi)))
),
0, # standard "simulate from normal and reject if outside limits" method. Use if bounds are wide.
ifelse(
(lower >= 0 & (upper > lower + 2*sqrt(exp(1)) /
(lower + sqrt(lower^2 + 4)) * exp((lower*2 - lower*sqrt(lower^2 + 4)) / 4))),
1, # rejection sampling with exponential proposal. Use if lower >> mu
ifelse(upper <= 0 & (-lower > -upper + 2*sqrt(exp(1)) /
(-upper + sqrt(upper^2 + 4)) * exp((upper*2 - -upper*sqrt(upper^2 + 4)) / 4)),
2, # rejection sampling with exponential proposal. Use if upper << mu.
3)))) # rejection sampling with uniform proposal. Use if bounds are narrow and central.
ind.nan <- ind[alg==-1]; ind.no <- ind[alg==0]; ind.expl <- ind[alg==1]; ind.expu <- ind[alg==2]; ind.u <- ind[alg==3]
ret[ind.nan] <- NaN
while (length(ind.no) > 0) {
y <- rnorm(length(ind.no))
done <- which(y >= lower[ind.no] & y <= upper[ind.no])
ret[ind.no[done]] <- y[done]
ind.no <- setdiff(ind.no, ind.no[done])
}
stopifnot(length(ind.no) == 0)
while (length(ind.expl) > 0) {
a <- (lower[ind.expl] + sqrt(lower[ind.expl]^2 + 4)) / 2
z <- rexp(length(ind.expl), a) + lower[ind.expl]
u <- runif(length(ind.expl))
done <- which((u <= exp(-(z - a)^2 / 2)) & (z <= upper[ind.expl]))
ret[ind.expl[done]] <- z[done]
ind.expl <- setdiff(ind.expl, ind.expl[done])
}
stopifnot(length(ind.expl) == 0)
while (length(ind.expu) > 0) {
a <- (-upper[ind.expu] + sqrt(upper[ind.expu]^2 +4)) / 2
z <- rexp(length(ind.expu), a) - upper[ind.expu]
u <- runif(length(ind.expu))
done <- which((u <= exp(-(z - a)^2 / 2)) & (z <= -lower[ind.expu]))
ret[ind.expu[done]] <- -z[done]
ind.expu <- setdiff(ind.expu, ind.expu[done])
}
stopifnot(length(ind.expu) == 0)
while (length(ind.u) > 0) {
z <- runif(length(ind.u), lower[ind.u], upper[ind.u])
rho <- ifelse(lower[ind.u] > 0,
exp((lower[ind.u]^2 - z^2) / 2), ifelse(upper[ind.u] < 0,
exp((upper[ind.u]^2 - z^2) / 2),
exp(-z^2/2)))
u <- runif(length(ind.u))
done <- which(u <= rho)
ret[ind.u[done]] <- z[done]
ind.u <- setdiff(ind.u, ind.u[done])
}
stopifnot(length(ind.u) == 0)
ret*sigma + mu
}
#' Slash Distribution inverse cumulative density function
#'
#' @name dslash
#' @title Slash distribution inverse cumulative density function
#' @description Inverse cumulative density function for the slash distribution.
#' @param x vector of quantiles
#' @param location vector of location parameters
#' @param scale vector of scale parameters
#' @param log if TRUE, the probabilites are given as log(p).
#' @returns a numeric vector
#' @export
#'
#'
#' @examples
#' qslash()
qslash <- function(p, location=0, scale=1, lower.tail=TRUE, log=FALSE)
{
qcustom(pslash, p=p, location=location, scale=scale, lower_tail=lower.tail, log=log)
}
qcustom <- function(pdist, p, ...){
args <- list(...)
if (!is.null(args$log)) args$log.p <- args$log
if (is.null(args$log.p)) args$log.p <- FALSE
if (is.null(args$lower.tail)) args$lower.tail <- TRUE
if (is.null(args$lbound)) args$lbound <- -Inf
if (is.null(args$ubound)) args$ubound <- Inf
if (args$log.p) p <- exp(p)
if (!args$lower.tail) p <- 1 - p
ret <- numeric(length(p))
ret[p == 0] <- args$lbound
ret[p == 1] <- args$ubound
args[c("lower.tail","log.p","lbound","ubound")] <- NULL
## Other args assumed to contain params of the distribution.
## Replicate all to their maximum length, along with p
maxlen <- max(sapply(c(args, p=list(p)), length))
for (i in seq(along=args))
args[[i]] <- rep(args[[i]], length.out=maxlen)
p <- rep(p, length.out=maxlen)
ret[p < 0 | p > 1] <- NaN
ind <- (p > 0 & p < 1)
if (any(ind)) {
hind <- seq(along=p)[ind]
h <- function(y) {
args <- lapply(args, function(x)x[hind[i]])
p <- p[hind[i]]
args$q <- y
(do.call(pdist, args) - p)
}
ptmp <- numeric(length(p[ind]))
for (i in 1:length(p[ind])) {
interval <- c(-1, 1)
while (h(interval[1])*h(interval[2]) >= 0) {
interval <- interval + c(-1,1)*0.5*(interval[2]-interval[1])
}
ptmp[i] <- uniroot(h, interval, tol=.Machine$double.eps)$root
}
ret[ind] <- ptmp
}
if (any(is.nan(ret))) warning("NaNs produced")
ret
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.