Nothing
R/egp.R
In mev: Modelling of Extreme Values
Defines functions
extgp.fit.ml.boot
extgp.fit.ml
extgp.nll
extgp.fit.pwm.boot
extgp.fit.pwm
extgp.pwm
rextgp
qextgp
dextgp
pextgp
rextgp.G
qextgp.G
dextgp.G
pextgp.G
egp2.fit
fit.extgp
Documented in
dextgp
dextgp.G
egp2.fit
fit.extgp
pextgp
pextgp.G
qextgp
qextgp.G
rextgp
rextgp.G
### Code by Raphael Huser with contributions from P. Naveau for the methods described in @references Naveau, P., R. Huser, P.
### Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection,
### \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
#' Fit an extended generalized Pareto distribution of Naveau et al.
#'
#' This is a wrapper function to obtain PWM or MLE estimates for
#' the extended GP models of Naveau et al. (2016) for rainfall intensities. The function calculates confidence intervals
#' by means of nonparametric percentile bootstrap and returns histograms and QQ plots of
#' the fitted distributions. The function handles both censoring and rounding.
#'
#' @details
#' The different models include the following transformations:
#' \itemize{
#' \item \code{model} 0 corresponds to uniform carrier, \eqn{G(u)=u}.
#' \item \code{model} 1 corresponds to a three parameters family, with carrier \eqn{G(u)=u^\kappa}.
#' \item \code{model} 2 corresponds to a three parameters family, with carrier \eqn{G(u)=1-V_\delta((1-u)^\delta)}.
#' \item \code{model} 3 corresponds to a four parameters family, with carrier \deqn{G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}}.
#' \item \code{model} 4 corresponds to a five parameter model (a mixture of \code{type} 2, with \eqn{G(u)=pu^\kappa + (1-p)*u^\delta}
#' }
#'
#' @author Raphael Huser and Philippe Naveau
#' @param data data vector.
#' @param model integer ranging from 0 to 4 indicating the model to select (see \code{\link{extgp}}).
#' @param method string; either \code{'mle'} for maximum likelihood, or \code{'pwm'} for probability weighted moments, or both.
#' @param init vector of initial values, comprising of \eqn{p}, \eqn{\kappa}, \eqn{\delta},\eqn{\sigma},\eqn{\xi} (in that order) for the optimization. All parameters may not appear depending on \code{model}.
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring; \code{censoring=c(0,Inf)} is equivalent to no censoring.
#' @param rounded numeric giving the instrumental precision (and rounding of the data), with default of 0.
#' @param confint logical; should confidence interval be returned (percentile bootstrap).
#' @param R integer; number of bootstrap replications.
#' @param ncpus integer; number of CPUs for parallel calculations (default: 1).
#' @param plots logical; whether to produce histogram and density plots.
#' @export
#' @importFrom grDevices rgb
#' @importFrom graphics hist layout
#' @seealso \code{\link{egp.fit}}, \code{\link{egp}}, \code{\link{extgp}}
#' @references Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
#' @examples
#' \dontrun{
#' data(rain, package = "ismev")
#' fit.extgp(rain[rain>0], model=1, method = 'mle', init = c(0.9, gp.fit(rain, 0)$est),
#' rounded = 0.1, confint = TRUE, R = 20)
#' }
fit.extgp <- function(data,
model = 1,
method = c("mle", "pwm"),
init,
censoring = c(0, Inf),
rounded = 0,
confint = FALSE,
R = 1000,
ncpus = 1,
plots = TRUE) {
method <- match.arg(method, c("mle", "pwm"), several.ok = TRUE)
if("pwm" %in% method){
if (!requireNamespace("gmm", quietly = TRUE)) {
stop("Package \"gmm\" needed for this function to work. Please install it.",
call. = FALSE)
}
}
if(confint){
if (!requireNamespace("boot", quietly = TRUE)) {
stop("Package \"boot\" needed for this function to work. Please install it.",
call. = FALSE)
}
}
# Sanity checks
initsize <- switch(model, 3, 3, 4, 5)
if (length(init) != initsize) {
stop("Invalid starting values in \"init\"; incorrect length.")
}
data = data[data > 0]
if (model == 3 && "pwm" %in% method) {
stop("No explicit formula for PWMs for this model...")
}
x <- seq(0, max(data, na.rm = TRUE), by = 0.05)
dextgps <- c()
qextgps <- c()
if (any(method == "pwm")) {
if (model == 1) {
fit.pwm <- extgp.fit.pwm(x = data, type = 1, kappa0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 1, kappa0 = init[1], sigma0 = init[2],
ncpus = ncpus, xi0 = init[3], censoring = censoring, parallel = "multicore")
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.kappa, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 1, kappa = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 1, kappa = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(1),
kappa = as.list(fit.pwm.boot$t[,
1]), sigma = as.list(fit.pwm.boot$t[, 2]), xi = as.list(fit.pwm.boot$t[, 3]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 2) {
fit.pwm <- extgp.fit.pwm(x = data, type = 2, delta0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 2, delta0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 2, delta = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 2, delta = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)),
type = list(2), delta = as.list(fit.pwm.boot$t[,
1]), sigma = as.list(fit.pwm.boot$t[, 2]), xi = as.list(fit.pwm.boot$t[, 3]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 3) {
fit.pwm <- extgp.fit.pwm(x = data, type = 3, kappa0 = init[1], delta0 = init[2], sigma0 = init[3], xi0 = init[4], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 3, kappa0 = init[1], delta0 = init[2],
sigma0 = init[3], xi0 = init[4], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 4, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.kappa, CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 3, kappa = fit.pwm[1], delta = fit.pwm[2], sigma = fit.pwm[3], xi = fit.pwm[4])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 3, kappa = fit.pwm[1], delta = fit.pwm[2], sigma = fit.pwm[3],
xi = fit.pwm[4])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(3), kappa = as.list(fit.pwm.boot$t[,
1]), delta = as.list(fit.pwm.boot$t[, 2]), sigma = as.list(fit.pwm.boot$t[, 3]), xi = as.list(fit.pwm.boot$t[,
4]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 4) {
fit.pwm <- extgp.fit.pwm(x = data, type = 4, prob0 = init[1], kappa0 = init[2], delta0 = init[3], sigma0 = init[4], xi0 = init[5],
censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 4, prob0 = init[1], kappa0 = init[2],
delta0 = init[3], sigma0 = init[4], xi0 = init[5], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.prob <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 4, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 5, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.prob, CI.pwm.kappa, CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 4, prob = fit.pwm[1], kappa = fit.pwm[2], delta = fit.pwm[3], sigma = fit.pwm[4], xi = fit.pwm[5])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 4, prob = fit.pwm[1], kappa = fit.pwm[2], delta = fit.pwm[3],
sigma = fit.pwm[4], xi = fit.pwm[5])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(4), prob = as.list(fit.pwm.boot$t[,
1]), kappa = as.list(fit.pwm.boot$t[, 2]), delta = as.list(fit.pwm.boot$t[, 3]), sigma = as.list(fit.pwm.boot$t[,
4]), xi = as.list(fit.pwm.boot$t[, 5]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
}
}
if (any(method == "mle")) {
if (model == 1) {
fit.mle <- extgp.fit.ml(x = data, type = 1, kappa0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 1, kappa0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.kappa, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 1, kappa = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 1, kappa = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(1), kappa = as.list(fit.mle.boot$t[,
1]), sigma = as.list(fit.mle.boot$t[, 2]), xi = as.list(fit.mle.boot$t[, 3]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 2) {
fit.mle <- extgp.fit.ml(x = data, type = 2, delta0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 2, delta0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 2, delta = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 2, delta = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(2), delta = as.list(fit.mle.boot$t[,
1]), sigma = as.list(fit.mle.boot$t[, 2]), xi = as.list(fit.mle.boot$t[, 3]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 3) {
fit.mle <- extgp.fit.ml(x = data, type = 3, kappa0 = init[1], delta0 = init[2], sigma0 = init[3], xi0 = init[4], censoring = censoring,
rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 3, kappa0 = init[1], delta0 = init[2],
sigma0 = init[3], xi0 = init[4], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 4, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.kappa, CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 3, kappa = fit.mle[1], delta = fit.mle[2], sigma = fit.mle[3], xi = fit.mle[4])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 3, kappa = fit.mle[1], delta = fit.mle[2], sigma = fit.mle[3],
xi = fit.mle[4])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(3), kappa = as.list(fit.mle.boot$t[,
1]), delta = as.list(fit.mle.boot$t[, 2]), sigma = as.list(fit.mle.boot$t[, 3]), xi = as.list(fit.mle.boot$t[,
4]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 4) {
fit.mle <- extgp.fit.ml(x = data, type = 4, prob0 = init[1], kappa0 = init[2], delta0 = init[3], sigma0 = init[4], xi0 = init[5],
censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 4, prob0 = init[1], kappa0 = init[2],
delta0 = init[3], sigma0 = init[4], xi0 = init[5], censoring = censoring, rounded = rounded, parallel = "multicore",
ncpus = ncpus)
CI.mle.prob <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 4, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 5, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.prob, CI.mle.kappa, CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 4, prob = fit.mle[1], kappa = fit.mle[2], delta = fit.mle[3], sigma = fit.mle[4], xi = fit.mle[5])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 4, prob = fit.mle[1], kappa = fit.mle[2], delta = fit.mle[3],
sigma = fit.mle[4], xi = fit.mle[5])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(4), prob = as.list(fit.mle.boot$t[,
1]), kappa = as.list(fit.mle.boot$t[, 2]), delta = as.list(fit.mle.boot$t[, 3]), sigma = as.list(fit.mle.boot$t[,
4]), xi = as.list(fit.mle.boot$t[, 5]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
}
}
if (any(method == "pwm")) {
if (any(method == "mle")) {
fits <- list(pwm = fit.pwm, mle = fit.mle)
if (confint) {
CIs <- list(pwm = CIs.pwm, mle = CIs.mle)
} else {
CIs <- NULL
}
} else {
fits <- list(pwm = fit.pwm)
if (confint & plots) {
CIs <- list(pwm = CIs.pwm)
} else {
CIs <- NULL
}
}
} else {
if (any(method == "mle")) {
fits <- list(mle = fit.mle)
if (confint) {
CIs <- list(mle = CIs.mle)
} else {
CIs <- NULL
}
} else {
fits <- NULL
CIs <- NULL
}
}
if (plots) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
mat <- matrix(c(1:(1 + plots)), nrow = 1, ncol = 1 + plots, byrow = TRUE)
layout(mat)
par(mar = c(4, 4, 1, 1))
hist(data, breaks = c(0:30) * max(data + 10^-1)/30, freq = FALSE, xlim = range(x),
xlab = "Rainfall amounts [mm]", ylab = "Density",
main = "", col = "lightgrey")
dens <- density(data, from = 0)
lines(dens$x, dens$y, col = "black")
abline(v = censoring, col = "lightgrey")
if (any(method == "pwm")) {
lines(x, dextgp.pwm, col = "red", lty = 2)
}
if (any(method == "mle")) {
lines(x, dextgp.mle, col = "blue")
}
if (any(method == "pwm")) {
plot(data, qextgp.pwm, asp = 1, xlab = "Empirical quantiles", ylab = "Fitted quantiles", ylim = range(qextgps, na.rm = TRUE),
type = "n")
} else {
if (any(method == "mle")) {
plot(data, qextgp.mle, asp = 1, xlab = "Empirical quantiles", ylab = "Fitted quantiles", ylim = range(qextgps, na.rm = TRUE),
type = "n")
}
}
if (confint) {
if (any(method == "pwm")) {
polygon(x = c(sort(data), sort(data, decreasing = TRUE)), y = c(q.pwm.L, q.pwm.U[length(q.pwm.U):1]), col = rgb(1,
0, 0, alpha = 0.1), lty = 2, border = rgb(1, 0, 0, alpha = 0.5))
}
if (any(method == "mle")) {
polygon(x = c(sort(data), sort(data, decreasing = TRUE)), y = c(q.mle.L, q.mle.U[length(q.mle.U):1]), col = rgb(0,
0, 1, alpha = 0.1), lty = 1, border = rgb(0, 0, 1, alpha = 0.5))
}
}
if (any(method == "pwm")) {
lines(sort(data), qextgp.pwm, lty = 2, type = "b", pch = 20, col = "red")
}
if (any(method == "mle")) {
lines(sort(data), qextgp.mle, lty = 1, type = "b", pch = 20, col = "blue")
}
abline(0, 1, col = "lightgrey")
}
return(list(fit = fits, confint = CIs))
}
#' Fit an extended generalized Pareto distribution of Naveau et al.
#'
#' Deprecated function name to fit an extended generalized Pareto family. The user should call \code{\link[mev]{fit.extgp}} instead.
#'
#' @seealso \code{\link[mev]{fit.extgp}}
#' @export
#' @keywords internal
#' @inheritParams fit.extgp
egp2.fit <- function(data, model = 1, method = c("mle", "pwm"), init, censoring = c(0, Inf), rounded = 0, CI = FALSE, R = 1000, ncpus = 1, plots = TRUE) {
fit.extgp(data = data, model = model, method = method, init = init, censoring = censoring, rounded = rounded,
confint = CI, R = R, ncpus = ncpus, plots = plots)
}
################################################################ ### Distribution/Density/Quantile/Random generator functions ### ###
### different types of distribution G(u):
### type=0: G(u)=u type=1: G(u)=u^kappa
### type=2: G(u)=1-V_delta((1-u)^delta)
### type=3: G(u)=(1-V_delta((1-u)^delta))^(kappa/2)
### type=4: G(u)=p*u^kappa + (1-p)*u^delta
#' @title Carrier distribution for the extended GP distributions of Naveau et al.
#' @description Density, distribution function, quantile function and random number
#' generation for the carrier distributions of the extended Generalized Pareto distributions.
#' @name extgp.G
#' @seealso \code{\link{extgp}}
#' @param u vector of observations (\code{dextgp.G}), probabilities (\code{qextgp.G}) or quantiles (\code{pextgp.G}), in \eqn{[0,1]}
#' @param prob mixture probability for model \code{type} \code{4}
#' @param kappa shape parameter for \code{type} \code{1}, \code{3} and \code{4}
#' @param delta additional parameter for \code{type} \code{2}, \code{3} and \code{4}
#' @param type integer between 0 to 5 giving the model choice
#' @param log logical; should the log-density be returned (default to \code{FALSE})?
#' @param n sample size
#' @param unifsamp sample of uniform; if provided, the data will be used in place of new uniform random variates
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring
#' @param direct logical; which method to use for sampling in model of \code{type} \code{4}?
#' @author Raphael Huser and Philippe Naveau
#'
#' @section Usage: \code{pextgp.G(u, type=1, prob, kappa, delta)}
#' @section Usage: \code{dextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA, log=FALSE)}
#' @section Usage: \code{qextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA)}
#' @section Usage: \code{rextgp.G(n, prob=NA, kappa=NA, delta=NA,
#' type=1, unifsamp=NULL, direct=FALSE, censoring=c(0,1))}
NULL
#' @rdname extgp-functions
#' @export
#' @keywords internal
pextgp.G <- function(u, type = 1, prob, kappa, delta) {
type <- as.integer(type[1])
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
if (type %in% c(1, 3, 4)) {
if(missing(kappa)){
stop("Argument \"kappa\" missing.")
} else{
if(!isTRUE(all(is.numeric(kappa),
length(kappa) == 1L,
kappa > 0))){
stop("Invalid \"kappa\" parameter.")
}
}
}
if (type %in% c(2, 3, 4)){
if(missing(delta)) {
stop("Argument \"delta\" missing.")
} else{
if(!isTRUE(all(is.numeric(delta),
length(delta) == 1L,
delta > 0))){
stop("Invalid \"delta\" parameter.")
}
}
}
if (type == 4){
if(missing(prob)) {
stop("Argument \"prob\" missing.")
} else{
if(!isTRUE(all(is.numeric(prob),
length(prob) == 1L,
prob >= 0,
prob <= 1))){
stop("Invalid \"prob\" parameter.")
}
}
}
u <- pmax(pmin(u, 1), 0)
if (type == 0) {
return(u)
} else if (type == 1) {
return(u^kappa)
} else if (type == 2) {
return(1 - pbeta((1 - u)^delta, 1/delta, 2))
} else if (type == 3) {
return((1 - pbeta((1 - u)^delta, 1/delta, 2))^(kappa/2))
} else if (type == 4) {
return(prob * u^kappa + (1 - prob) * u^delta)
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
dextgp.G <- function(u, type = 1, prob = NA, kappa = NA, delta = NA, log = FALSE) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob\" missing.")
}
if (log == FALSE) {
if (type == 0) {
return(1)
} else if (type == 1) {
return(kappa * u^(kappa - 1))
} else if (type == 2) {
return(dbeta((1 - u)^delta, 1/delta, 2) * delta * (1 - u)^(delta - 1))
} else if (type == 3) {
return((kappa/2) * (1 - pbeta((1 - u)^delta, 1/delta, 2))^(kappa/2 - 1) * dbeta((1 - u)^delta, 1/delta, 2) * delta * (1 -
u)^(delta - 1))
} else if (type == 4) {
return(prob * kappa * u^(kappa - 1) + (1 - prob) * delta * u^(delta - 1))
}
} else {
if (type == 0) {
return(0)
} else if (type == 1) {
return(log(kappa) + (kappa - 1) * log(u))
} else if (type == 2) {
return(dbeta((1 - u)^delta, 1/delta, 2, log = TRUE) + log(delta) + (delta - 1) * log(1 - u))
} else if (type == 3) {
return(log(kappa/2) + (kappa/2 - 1) * log(1 - pbeta((1 - u)^delta, 1/delta, 2)) + dbeta((1 - u)^delta, 1/delta, 2, log = TRUE) +
log(delta) + (delta - 1) * log(1 - u))
} else if (type == 4) {
return(log(prob * kappa * u^(kappa - 1) + (1 - prob) * delta * u^(delta - 1)))
}
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
qextgp.G <- function(u, type = 1, prob = NA, kappa = NA, delta = NA) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob\" missing.")
}
if (type == 0) {
return(u)
} else if (type == 1) {
return(u^(1/kappa))
} else if (type == 2) {
return(1 - qbeta(1 - u, 1/delta, 2)^(1/delta))
} else if (type == 3) {
return(1 - qbeta(1 - u^(2/kappa), 1/delta, 2)^(1/delta))
} else if (type == 4) {
dummy.func <- function(u, p, prob = NA, kappa = NA, delta = NA) {
return(pextgp.G(u = u, prob = prob, kappa = kappa, delta = delta, type = 4) - p)
}
find.root <- function(u, prob = NA, kappa = NA, delta = NA) {
return(uniroot(dummy.func, interval = c(0, 1), p = u, prob = prob, kappa = kappa, delta = delta)$root)
}
return(sapply(u, FUN = find.root, prob = prob, kappa = kappa, delta = delta))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
rextgp.G <- function(n, prob = NA, kappa = NA, delta = NA, type = 1, unifsamp = NULL, direct = FALSE, censoring = c(0, 1)) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob missing.")
}
if (is.null(unifsamp)) {
unifsamp <- runif(n)
}
if (type != 4 | (type == 4 & direct)) {
if (censoring[1] == 0 & censoring[2] == 1) {
return(qextgp.G(unifsamp, prob = prob, kappa = kappa, delta = delta, type = type))
} else {
x.L <- pextgp.G(censoring[1], prob = prob, kappa = kappa, delta = delta, type = type)
x.U <- pextgp.G(censoring[2], prob = prob, kappa = kappa, delta = delta, type = type)
return(qextgp.G(x.L + (x.U - x.L) * unifsamp, prob = prob, kappa = kappa, delta = delta, type = type))
}
} else if (type == 4 & !direct) {
if (censoring[1] == 0 & censoring[2] == 1) {
components <- sample(x = c(1, 2), size = n, replace = TRUE, prob = c(prob, 1 - prob))
res <- c()
res[components == 1] <- qextgp.G(unifsamp[components == 1], prob = NA, kappa = kappa, delta = NA, type = 1)
res[components == 2] <- qextgp.G(unifsamp[components == 2], prob = NA, kappa = delta, delta = NA, type = 1)
return(res)
} else {
x.L <- pextgp.G(censoring[1], prob = prob, kappa = kappa, delta = delta, type = type)
x.U <- pextgp.G(censoring[2], prob = prob, kappa = kappa, delta = delta, type = type)
prop1 <- prob * (pextgp.G(censoring[2], prob = NA, kappa = kappa, delta = NA, type = 1) - pextgp.G(censoring[1], prob = NA,
kappa = kappa, delta = NA, type = 1))
prop2 <- (1 - prob) * (pextgp.G(censoring[2], prob = NA, kappa = delta, delta = NA, type = 1) - pextgp.G(censoring[1], prob = NA,
kappa = delta, delta = NA, type = 1))
new.prob <- prop1/(prop1 + prop2)
components <- sample(x = c(1, 2), size = n, replace = TRUE, prob = c(new.prob, 1 - new.prob))
res <- c()
res[components == 1] <- qextgp.G(x.L + (x.U - x.L) * unifsamp[components == 1], prob = NA, kappa = kappa, delta = NA, type = 1)
res[components == 2] <- qextgp.G(x.L + (x.U - x.L) * unifsamp[components == 2], prob = NA, kappa = delta, delta = NA, type = 1)
return(res)
}
}
}
### different types of extended GP type=0:
### exact GP ---> 2 parameters (sigma,xi)
### type=1: extgp with G(u)=u^kappa ---> 3 parameters(kappa,sigma,xi)
### type=2: extgp with G(u)=1-V_delta((1-u)^delta) ---> 3 parameters (delta,sigma,xi)
### type=3: extgp with G(u)=(1-V_delta((1-u)^delta))^(kappa/2) ---> 4 parameters (kappa,delta,sigma,xi)
### type=4: extgp with mixture G(u)=p*u^kappa + (1-p)*u^delta ---> 5 parameters (p,kappa,delta,sigma,xi)
#' @title Extended generalised Pareto families of Naveau et al. (2016)
#'
#' @description Density function, distribution function, quantile function and random generation for the extended generalized Pareto distribution (GPD) with scale and shape parameters.
#'
#' @details The extended generalized Pareto families proposed in Naveau \emph{et al.} (2016)
#'retain the tail index of the distribution while being compliant with the theoretical behavior of extreme
#'low rainfall. There are five proposals, the first one being equivalent to the GP distribution.
#'\itemize{
#' \item \code{type} 0 corresponds to uniform carrier, \eqn{G(u)=u}.
#' \item \code{type} 1 corresponds to a three parameters family, with carrier \eqn{G(u)=u^\kappa}.
#' \item \code{type} 2 corresponds to a three parameters family, with carrier \eqn{G(u)=1-V_\delta((1-u)^\delta)}.
#' \item \code{type} 3 corresponds to a four parameters family, with carrier \deqn{G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}}.
#' \item \code{type} 4 corresponds to a five parameter model (a mixture of \code{type} 2, with \eqn{G(u)=pu^\kappa + (1-p)*u^\delta}
#' }
#' @name extgp
#' @param q vector of quantiles
#' @param x vector of observations
#' @param p vector of probabilities
#' @param n sample size
#' @param prob mixture probability for model \code{type} \code{4}
#' @param kappa shape parameter for \code{type} \code{1}, \code{3} and \code{4}
#' @param delta additional parameter for \code{type} \code{2}, \code{3} and \code{4}
#' @param sigma scale parameter
#' @param xi shape parameter
#' @param type integer between 0 to 5 giving the model choice
#' @param step function of step size for discretization with default \code{0}, corresponding to continuous quantiles
#' @param log logical; should the log-density be returned (default to \code{FALSE})?
#' @param unifsamp sample of uniform; if provided, the data will be used in place of new uniform random variates
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring
#'
#' @section Usage: \code{pextgp(q, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)}
#' @section Usage: \code{dextgp(x, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, log=FALSE)}
#' @section Usage: \code{qextgp(p, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)}
#' @section Usage: \code{rextgp(n, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, unifsamp=NULL, censoring=c(0,Inf))}
#' @author Raphael Huser and Philippe Naveau
#' @references Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
NULL
#' Extended GP functions
#'
#' These functions are documented in \code{\link{extgp}} and in \code{\link{extgp.G}} for the carrier distributions supported in the unit interval.
#' @name extgp-functions
#' @export
#' @seealso \code{\link{extgp}}, \code{\link{extgp.G}}
#' @keywords internal
pextgp <- function(q, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1) {
return(pextgp.G(mev::pgp(q, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type))
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
dextgp <- function(x, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, log = FALSE) {
if (log == FALSE) {
return(dextgp.G(mev::pgp(x, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type) * mev::dgp(x, scale = sigma,
shape = xi))
} else {
return(dextgp.G(mev::pgp(x, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type, log = TRUE) + mev::dgp(x,
scale = sigma, shape = xi, log = TRUE))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
qextgp <- function(p, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, step = 0) {
stopifnot(length(step) == 1L, step >= 0, is.finite(step))
if(isTRUE(all.equal(step, 0, check.attributes = FALSE))){
return(mev::qgp(qextgp.G(p, prob = prob, kappa = kappa, delta = delta, type = type), scale = sigma, shape = xi))
} else{
qcont <- mev::qgp(qextgp.G(p, prob = prob, kappa = kappa, delta = delta, type = type), scale = sigma, shape = xi)
return(step*floor((qcont-step)/step))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
rextgp <- function(n, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, unifsamp = NULL, censoring = c(0, Inf)) {
if (censoring[1] == 0 & censoring[2] == Inf) {
return(mev::qgp(rextgp.G(n, prob = prob, kappa = kappa, delta = delta, type = type, unifsamp, censoring = c(0, 1)), scale = sigma,
shape = xi))
} else {
H.L <- mev::pgp(censoring[1], loc = 0, scale = sigma, shape = xi)
H.U <- mev::pgp(censoring[2], loc = 0, scale = sigma, shape = xi)
return(mev::qgp(rextgp.G(n, prob = prob, kappa = kappa, delta = delta, type = type, unifsamp, censoring = c(H.L, H.U)), scale = sigma,
shape = xi))
}
}
######################### ### Inference via PWM ### ###
extgp.pwm <- function(orders, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, censoring = c(0, Inf), empiric = FALSE,
unifsamp = NULL, NbSamples = 10^4, N = 200) {
H.L <- ifelse(censoring[1] > 0, mev::pgp(censoring[1], loc = 0, scale = sigma, shape = xi), 0)
H.U <- ifelse(censoring[2] < Inf, mev::pgp(censoring[2], loc = 0, scale = sigma, shape = xi), 1)
F.L <- ifelse(censoring[1] > 0, pextgp(censoring[1], prob = prob, kappa = kappa, delta = delta, sigma = sigma, xi = xi, type = type),
0)
F.U <- ifelse(censoring[2] < Inf, pextgp(censoring[2], prob = prob, kappa = kappa, delta = delta, sigma = sigma, xi = xi, type = type),
1)
prob.LU <- F.U - F.L
if (!empiric) {
E2 <- ((1 - F.L)^(orders + 1) - (1 - F.U)^(orders + 1))/(prob.LU * (orders + 1))
if (type == 0) {
E1 <- (H.U - H.L)^(orders - xi)/(orders - xi + 1)
return((sigma/xi) * (E1 - E2))
} else if (type == 1) {
beta.coefs <- beta(c(1:(max(orders) + 1)) * kappa, 1 - xi)
probs.beta <- pbeta(H.U, shape1 = c(1:(max(orders) + 1)) * kappa, shape2 = 1 - xi) - pbeta(H.L, shape1 = c(1:(max(orders) +
1)) * kappa, shape2 = 1 - xi)
E1 <- c()
for (i in 1:length(orders)) {
E1[i] <- (kappa/prob.LU) * sum((-1)^c(0:orders[i]) * choose(orders[i], c(0:orders[i])) * beta.coefs[1:(orders[i] +
1)] * probs.beta[1:(orders[i] + 1)])
}
return((sigma/xi) * (E1 - E2))
} else if (type == 2) {
E1 <- c()
for (i in 1:length(orders)) {
expos1 <- orders[i] - xi + delta * c(0:orders[i]) + 1
expos2 <- orders[i] - xi + delta * (c(0:orders[i]) + 1) + 1
numer.j <- expos1 * ((1 - H.U)^expos2 - (1 - H.L)^expos2) + expos2 * ((1 - H.L)^expos1 - (1 - H.U)^expos1)
denom.j <- expos1 * expos2
E1[i] <- (1/(delta^(orders[i] + 1) * prob.LU)) * sum((-1)^c(0:orders[i]) * choose(orders[i], c(0:orders[i])) * (1 +
delta)^(orders[i] - c(0:orders[i]) + 1) * numer.j/denom.j)
}
return((sigma/xi) * (E1 - E2))
} else if (type == 3) {
print("No explicit formula for PWMs for this model...")
return(NA)
} else if (type == 4) {
E1 <- c()
for (i in 1:length(orders)) {
l.vals <- matrix(c(0:orders[i]), nrow = orders[i] + 1, ncol = orders[i] + 1, byrow = TRUE)
m.vals <- matrix(c(0:orders[i]), nrow = orders[i] + 1, ncol = orders[i] + 1)
inds <- l.vals < m.vals
l.vals[inds] <- m.vals[inds] <- NA
const <- choose(orders[i], l.vals) * choose(l.vals, m.vals) * (-1)^l.vals * prob^m.vals * (1 - prob)^(l.vals - m.vals)
beta.coefs1 <- beta(kappa * (m.vals + 1) + delta * (l.vals - m.vals), 1 - xi)
beta.coefs2 <- beta(kappa * m.vals + delta * (l.vals - m.vals + 1), 1 - xi)
probs.beta1 <- pbeta(H.U, shape1 = kappa * (m.vals + 1) + delta * (l.vals - m.vals), shape2 = 1 - xi) - pbeta(H.L,
shape1 = kappa * (m.vals + 1) + delta * (l.vals - m.vals), shape2 = 1 - xi)
probs.beta2 <- pbeta(H.U, shape1 = kappa * m.vals + delta * (l.vals - m.vals + 1), shape2 = 1 - xi) - pbeta(H.L, shape1 = kappa *
m.vals + delta * (l.vals - m.vals + 1), shape2 = 1 - xi)
E1[i] <- (1/prob.LU) * sum(const * (prob * kappa * beta.coefs1 * probs.beta1 + (1 - prob) * delta * beta.coefs2 *
probs.beta2), na.rm = TRUE)
}
return((sigma/xi) * (E1 - E2))
}
} else {
if (is.null(unifsamp)) {
unifsamp <- runif(NbSamples)
}
if (censoring[1] == 0 & censoring[2] == Inf) {
V <- rextgp.G(length(unifsamp), prob, kappa, delta, type, unifsamp, direct = TRUE, censoring = c(0, 1))
} else {
V <- rextgp.G(length(unifsamp), prob, kappa, delta, type, unifsamp, direct = TRUE, censoring = c(H.L, H.U))
}
X <- mev::qgp(V, scale = sigma, shape = xi)
res <- c()
for (i in 1:length(orders)) {
res[i] <- mean(X * (1 - (F.L + (F.U - F.L) * unifsamp))^orders[i], na.rm = TRUE)
}
return(res)
}
}
extgp.fit.pwm <- function(x, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf), empiric = FALSE,
unifsamp = NULL, NbSamples = 10^4) {
Fn = ecdf(x)
inds <- x > censoring[1] & x < censoring[2]
mu0hat = mean(x[inds])
mu1hat = mean(x[inds] * (1 - Fn(x[inds])))
mu2hat = mean(x[inds] * (1 - Fn(x[inds]))^2)
mu3hat = mean(x[inds] * (1 - Fn(x[inds]))^3)
mu4hat = mean(x[inds] * (1 - Fn(x[inds]))^4)
if (type == 0) {
if (censoring[1] == 0 & censoring[2] == Inf) {
xihat <- (mu0hat - 4 * mu1hat)/(mu0hat - 2 * mu1hat)
sigmahat <- mu0hat * (1 - xihat)
thetahat <- c(sigmahat, xihat)
} else {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:1), sigma = theta[1], xi = theta[2], type = 0, censoring = censoring, empiric = empiric,
unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat)
return(matrix(pwm.theor - pwm.empir, ncol = 2))
}
theta0 <- c(sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 10^(-6)), upper = c(Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
}
names(thetahat) <- c("sigma", "xi")
return(thetahat)
} else if (type == 1) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:2), kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat)
return(matrix(pwm.theor - pwm.empir, ncol = 3))
}
theta0 <- c(kappa0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 10^(-6)), upper = c(Inf, Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("kappa", "sigma", "xi")
return(thetahat)
} else if (type == 2) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:2), delta = theta[1], sigma = theta[2], xi = theta[3], type = 2, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat)
return(matrix(pwm.theor - pwm.empir, ncol = 3))
}
theta0 <- c(delta0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 10^(-6)), upper = c(100, Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("delta", "sigma", "xi")
return(thetahat)
} else if (type == 3) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:3), kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4], type = 3,
censoring = censoring, empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat, mu3hat)
return(matrix(pwm.theor - pwm.empir, ncol = 4))
}
theta0 <- c(kappa0, delta0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 1e-04, 10^(-6)), upper = c(Inf, 100, Inf, 0.99),
vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("kappa", "delta", "sigma", "xi")
return(thetahat)
} else if (type == 4) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:4), prob = theta[1], kappa = theta[2], delta = theta[2] + theta[3], sigma = theta[4],
xi = theta[5], type = 4, censoring = censoring, empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat, mu3hat, mu4hat)
return(matrix(pwm.theor - pwm.empir, ncol = 5))
}
res0 <- extgp.fit.pwm(x[x > quantile(x, 0.9)] - quantile(x, 0.9), type = 0)
if (is.na(prob0)) {
prob0 <- 0.5
}
if (is.na(kappa0)) {
kappa0 <- mean(x[x < quantile(x, 0.1)])/(quantile(x, 0.1) - mean(x[x < quantile(x, 0.1)]))
}
if (is.na(delta0)) {
delta0 <- kappa0 + 0.01
}
Ddelta0 <- delta0 - kappa0
if (is.na(sigma0)) {
sigma0 <- res0[1]
}
if (is.na(xi0)) {
xi0 <- res0[2]
}
theta0 <- c(prob0, kappa0, Ddelta0, sigma0, xi0)
names(theta0) <- c("prob", "kappa", "Ddelta", "sigma", "xi")
# print(theta0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(0, 1e-04, 0, 1e-04, 10^(-6)), upper = c(1, Inf, Inf, Inf, 0.99),
vcov = "iid")
thetahat <- res$coefficients
thetahat[3] <- thetahat[2] + thetahat[3]
names(thetahat) <- c("prob", "kappa", "delta", "sigma", "xi")
return(thetahat)
}
}
extgp.fit.pwm.boot <- function(data, i, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf),
empiric = FALSE, unifsamp = NULL, NbSamples = 10^4) {
return(extgp.fit.pwm(data[i], type = type, prob0 = prob0, kappa0 = kappa0, delta0 = delta0, sigma0 = sigma0, xi0 = xi0, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples))
}
######################## ### Inference via ML ### ###
extgp.nll <- function(theta, x, censoring, rounded, type) {
x.cens1 <- x[x < censoring[1]]
x.cens2 <- x[x > censoring[2]]
x.not.cens <- x[x >= censoring[1] & x <= censoring[2]]
if (type == 0) {
if (theta[1] <= 0 | theta[2] <= 10^(-6) | theta[2] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], sigma = theta[1], xi = theta[2], type = 0), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], sigma = theta[1], xi = theta[2], type = 0), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, sigma = theta[1], xi = theta[2], type = 0, log = TRUE),
na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, sigma = theta[1], xi = theta[2], type = 0) - pextgp(x.not.cens,
sigma = theta[1], xi = theta[2], type = 0))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 1) {
if (theta[1] <= 0 | theta[2] <= 0 | theta[3] <= 10^(-6) | theta[3] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1),
0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1),
1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1,
log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, kappa = theta[1], sigma = theta[2], xi = theta[3],
type = 1) - pextgp(x.not.cens, kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 2) {
if (theta[1] <= 0 | theta[1] > 100 | theta[2] <= 0 | theta[3] <= 10^(-6) | theta[3] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], delta = theta[1], sigma = theta[2], xi = theta[3], type = 2),
0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], delta = theta[1], sigma = theta[2], xi = theta[3], type = 2),
1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, delta = theta[1], sigma = theta[2], xi = theta[3], type = 2,
log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, delta = theta[1], sigma = theta[2], xi = theta[3],
type = 2) - pextgp(x.not.cens, delta = theta[1], sigma = theta[2], xi = theta[3], type = 2))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 3) {
if (theta[1] <= 0 | theta[2] <= 0 | theta[2] > 100 | theta[3] <= 0 | theta[4] <= 10^(-6) | theta[4] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3, log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, kappa = theta[1], delta = theta[2], sigma = theta[3],
xi = theta[4], type = 3) - pextgp(x.not.cens, kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 4) {
if (theta[1] < 0 | theta[1] > 1 | theta[2] <= 0 | theta[3] <= 0 | theta[4] <= 0 | theta[5] <= 10^(-6) | theta[5] > 0.99 |
theta[3] < theta[2]) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4, log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, prob = theta[1], kappa = theta[2],
delta = theta[3], sigma = theta[4], xi = theta[5], type = 4) - pextgp(x.not.cens, prob = theta[1], kappa = theta[2],
delta = theta[3], sigma = theta[4], xi = theta[5], type = 4))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
}
}
extgp.fit.ml <- function(x, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf), rounded = 0.1,
print = FALSE) {
if (type == 0) {
theta0 <- c(sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("sigma", "xi")
return(thetahat)
} else if (type == 1) {
theta0 <- c(kappa0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("kappa", "sigma", "xi")
return(thetahat)
} else if (type == 2) {
theta0 <- c(delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("delta", "sigma", "xi")
return(thetahat)
} else if (type == 3) {
theta0 <- c(kappa0, delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("kappa", "delta", "sigma", "xi")
return(thetahat)
} else if (type == 4) {
theta0 <- c(prob0, kappa0, delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("prob", "kappa", "delta", "sigma", "xi")
return(thetahat)
}
}
extgp.fit.ml.boot <- function(data, i, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf),
rounded = 0.1, print = FALSE) {
return(extgp.fit.ml(data[i], type = type, prob0 = prob0, kappa0 = kappa0, delta0 = delta0, sigma0 = sigma0, xi0 = xi0, censoring = censoring,
rounded = rounded, print = print))
}
Try the mev package in your browser
Any scripts or data that you put into this service are public.
mev documentation built on April 20, 2023, 5:10 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 extgp.fit.ml.boot extgp.fit.ml extgp.nll extgp.fit.pwm.boot extgp.fit.pwm extgp.pwm rextgp qextgp dextgp pextgp rextgp.G qextgp.G dextgp.G pextgp.G egp2.fit fit.extgp
Documented in dextgp dextgp.G egp2.fit fit.extgp pextgp pextgp.G qextgp qextgp.G rextgp rextgp.G
### Code by Raphael Huser with contributions from P. Naveau for the methods described in @references Naveau, P., R. Huser, P.
### Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection,
### \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
#' Fit an extended generalized Pareto distribution of Naveau et al.
#'
#' This is a wrapper function to obtain PWM or MLE estimates for
#' the extended GP models of Naveau et al. (2016) for rainfall intensities. The function calculates confidence intervals
#' by means of nonparametric percentile bootstrap and returns histograms and QQ plots of
#' the fitted distributions. The function handles both censoring and rounding.
#'
#' @details
#' The different models include the following transformations:
#' \itemize{
#' \item \code{model} 0 corresponds to uniform carrier, \eqn{G(u)=u}.
#' \item \code{model} 1 corresponds to a three parameters family, with carrier \eqn{G(u)=u^\kappa}.
#' \item \code{model} 2 corresponds to a three parameters family, with carrier \eqn{G(u)=1-V_\delta((1-u)^\delta)}.
#' \item \code{model} 3 corresponds to a four parameters family, with carrier \deqn{G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}}.
#' \item \code{model} 4 corresponds to a five parameter model (a mixture of \code{type} 2, with \eqn{G(u)=pu^\kappa + (1-p)*u^\delta}
#' }
#'
#' @author Raphael Huser and Philippe Naveau
#' @param data data vector.
#' @param model integer ranging from 0 to 4 indicating the model to select (see \code{\link{extgp}}).
#' @param method string; either \code{'mle'} for maximum likelihood, or \code{'pwm'} for probability weighted moments, or both.
#' @param init vector of initial values, comprising of \eqn{p}, \eqn{\kappa}, \eqn{\delta},\eqn{\sigma},\eqn{\xi} (in that order) for the optimization. All parameters may not appear depending on \code{model}.
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring; \code{censoring=c(0,Inf)} is equivalent to no censoring.
#' @param rounded numeric giving the instrumental precision (and rounding of the data), with default of 0.
#' @param confint logical; should confidence interval be returned (percentile bootstrap).
#' @param R integer; number of bootstrap replications.
#' @param ncpus integer; number of CPUs for parallel calculations (default: 1).
#' @param plots logical; whether to produce histogram and density plots.
#' @export
#' @importFrom grDevices rgb
#' @importFrom graphics hist layout
#' @seealso \code{\link{egp.fit}}, \code{\link{egp}}, \code{\link{extgp}}
#' @references Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
#' @examples
#' \dontrun{
#' data(rain, package = "ismev")
#' fit.extgp(rain[rain>0], model=1, method = 'mle', init = c(0.9, gp.fit(rain, 0)$est),
#' rounded = 0.1, confint = TRUE, R = 20)
#' }
fit.extgp <- function(data,
model = 1,
method = c("mle", "pwm"),
init,
censoring = c(0, Inf),
rounded = 0,
confint = FALSE,
R = 1000,
ncpus = 1,
plots = TRUE) {
method <- match.arg(method, c("mle", "pwm"), several.ok = TRUE)
if("pwm" %in% method){
if (!requireNamespace("gmm", quietly = TRUE)) {
stop("Package \"gmm\" needed for this function to work. Please install it.",
call. = FALSE)
}
}
if(confint){
if (!requireNamespace("boot", quietly = TRUE)) {
stop("Package \"boot\" needed for this function to work. Please install it.",
call. = FALSE)
}
}
# Sanity checks
initsize <- switch(model, 3, 3, 4, 5)
if (length(init) != initsize) {
stop("Invalid starting values in \"init\"; incorrect length.")
}
data = data[data > 0]
if (model == 3 && "pwm" %in% method) {
stop("No explicit formula for PWMs for this model...")
}
x <- seq(0, max(data, na.rm = TRUE), by = 0.05)
dextgps <- c()
qextgps <- c()
if (any(method == "pwm")) {
if (model == 1) {
fit.pwm <- extgp.fit.pwm(x = data, type = 1, kappa0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 1, kappa0 = init[1], sigma0 = init[2],
ncpus = ncpus, xi0 = init[3], censoring = censoring, parallel = "multicore")
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.kappa, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 1, kappa = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 1, kappa = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(1),
kappa = as.list(fit.pwm.boot$t[,
1]), sigma = as.list(fit.pwm.boot$t[, 2]), xi = as.list(fit.pwm.boot$t[, 3]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 2) {
fit.pwm <- extgp.fit.pwm(x = data, type = 2, delta0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 2, delta0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 2, delta = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 2, delta = fit.pwm[1], sigma = fit.pwm[2], xi = fit.pwm[3])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)),
type = list(2), delta = as.list(fit.pwm.boot$t[,
1]), sigma = as.list(fit.pwm.boot$t[, 2]), xi = as.list(fit.pwm.boot$t[, 3]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 3) {
fit.pwm <- extgp.fit.pwm(x = data, type = 3, kappa0 = init[1], delta0 = init[2], sigma0 = init[3], xi0 = init[4], censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 3, kappa0 = init[1], delta0 = init[2],
sigma0 = init[3], xi0 = init[4], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 4, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.kappa, CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 3, kappa = fit.pwm[1], delta = fit.pwm[2], sigma = fit.pwm[3], xi = fit.pwm[4])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 3, kappa = fit.pwm[1], delta = fit.pwm[2], sigma = fit.pwm[3],
xi = fit.pwm[4])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(3), kappa = as.list(fit.pwm.boot$t[,
1]), delta = as.list(fit.pwm.boot$t[, 2]), sigma = as.list(fit.pwm.boot$t[, 3]), xi = as.list(fit.pwm.boot$t[,
4]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 4) {
fit.pwm <- extgp.fit.pwm(x = data, type = 4, prob0 = init[1], kappa0 = init[2], delta0 = init[3], sigma0 = init[4], xi0 = init[5],
censoring = censoring)
if (confint) {
fit.pwm.boot <- boot::boot(data = data, statistic = extgp.fit.pwm.boot, R = R, type = 4, prob0 = init[1], kappa0 = init[2],
delta0 = init[3], sigma0 = init[4], xi0 = init[5], censoring = censoring, parallel = "multicore", ncpus = ncpus)
CI.pwm.prob <- boot::boot.ci(boot.out = fit.pwm.boot, index = 1, type = "perc")$perc[4:5]
CI.pwm.kappa <- boot::boot.ci(boot.out = fit.pwm.boot, index = 2, type = "perc")$perc[4:5]
CI.pwm.delta <- boot::boot.ci(boot.out = fit.pwm.boot, index = 3, type = "perc")$perc[4:5]
CI.pwm.sigma <- boot::boot.ci(boot.out = fit.pwm.boot, index = 4, type = "perc")$perc[4:5]
CI.pwm.xi <- boot::boot.ci(boot.out = fit.pwm.boot, index = 5, type = "perc")$perc[4:5]
CIs.pwm <- cbind(CI.pwm.prob, CI.pwm.kappa, CI.pwm.delta, CI.pwm.sigma, CI.pwm.xi)
}
dextgp.pwm <- dextgp(x = x, type = 4, prob = fit.pwm[1], kappa = fit.pwm[2], delta = fit.pwm[3], sigma = fit.pwm[4], xi = fit.pwm[5])
dextgps <- c(dextgps, dextgp.pwm)
if (plots) {
qextgp.pwm <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 4, prob = fit.pwm[1], kappa = fit.pwm[2], delta = fit.pwm[3],
sigma = fit.pwm[4], xi = fit.pwm[5])
qextgps <- c(qextgps, qextgp.pwm)
if (confint) {
q.pwm.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(4), prob = as.list(fit.pwm.boot$t[,
1]), kappa = as.list(fit.pwm.boot$t[, 2]), delta = as.list(fit.pwm.boot$t[, 3]), sigma = as.list(fit.pwm.boot$t[,
4]), xi = as.list(fit.pwm.boot$t[, 5]))
q.pwm.L <- apply(q.pwm.boot, 1, quantile, 0.025, na.rm = TRUE)
q.pwm.U <- apply(q.pwm.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
}
}
if (any(method == "mle")) {
if (model == 1) {
fit.mle <- extgp.fit.ml(x = data, type = 1, kappa0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 1, kappa0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.kappa, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 1, kappa = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 1, kappa = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(1), kappa = as.list(fit.mle.boot$t[,
1]), sigma = as.list(fit.mle.boot$t[, 2]), xi = as.list(fit.mle.boot$t[, 3]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 2) {
fit.mle <- extgp.fit.ml(x = data, type = 2, delta0 = init[1], sigma0 = init[2], xi0 = init[3], censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 2, delta0 = init[1], sigma0 = init[2],
xi0 = init[3], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 2, delta = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 2, delta = fit.mle[1], sigma = fit.mle[2], xi = fit.mle[3])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(2), delta = as.list(fit.mle.boot$t[,
1]), sigma = as.list(fit.mle.boot$t[, 2]), xi = as.list(fit.mle.boot$t[, 3]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 3) {
fit.mle <- extgp.fit.ml(x = data, type = 3, kappa0 = init[1], delta0 = init[2], sigma0 = init[3], xi0 = init[4], censoring = censoring,
rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 3, kappa0 = init[1], delta0 = init[2],
sigma0 = init[3], xi0 = init[4], censoring = censoring, rounded = rounded, parallel = "multicore", ncpus = ncpus)
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 4, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.kappa, CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 3, kappa = fit.mle[1], delta = fit.mle[2], sigma = fit.mle[3], xi = fit.mle[4])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 3, kappa = fit.mle[1], delta = fit.mle[2], sigma = fit.mle[3],
xi = fit.mle[4])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(3), kappa = as.list(fit.mle.boot$t[,
1]), delta = as.list(fit.mle.boot$t[, 2]), sigma = as.list(fit.mle.boot$t[, 3]), xi = as.list(fit.mle.boot$t[,
4]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
} else if (model == 4) {
fit.mle <- extgp.fit.ml(x = data, type = 4, prob0 = init[1], kappa0 = init[2], delta0 = init[3], sigma0 = init[4], xi0 = init[5],
censoring = censoring, rounded = rounded)
if (confint) {
fit.mle.boot <- boot::boot(data = data, statistic = extgp.fit.ml.boot, R = R, type = 4, prob0 = init[1], kappa0 = init[2],
delta0 = init[3], sigma0 = init[4], xi0 = init[5], censoring = censoring, rounded = rounded, parallel = "multicore",
ncpus = ncpus)
CI.mle.prob <- boot::boot.ci(boot.out = fit.mle.boot, index = 1, type = "perc")$perc[4:5]
CI.mle.kappa <- boot::boot.ci(boot.out = fit.mle.boot, index = 2, type = "perc")$perc[4:5]
CI.mle.delta <- boot::boot.ci(boot.out = fit.mle.boot, index = 3, type = "perc")$perc[4:5]
CI.mle.sigma <- boot::boot.ci(boot.out = fit.mle.boot, index = 4, type = "perc")$perc[4:5]
CI.mle.xi <- boot::boot.ci(boot.out = fit.mle.boot, index = 5, type = "perc")$perc[4:5]
CIs.mle <- cbind(CI.mle.prob, CI.mle.kappa, CI.mle.delta, CI.mle.sigma, CI.mle.xi)
}
dextgp.mle <- dextgp(x = x, type = 4, prob = fit.mle[1], kappa = fit.mle[2], delta = fit.mle[3], sigma = fit.mle[4], xi = fit.mle[5])
dextgps <- c(dextgps, dextgp.mle)
if (plots) {
qextgp.mle <- qextgp(p = c(1:length(data))/(length(data) + 1), type = 4, prob = fit.mle[1], kappa = fit.mle[2], delta = fit.mle[3],
sigma = fit.mle[4], xi = fit.mle[5])
qextgps <- c(qextgps, qextgp.mle)
if (confint) {
q.mle.boot <- mapply(FUN = qextgp, p = list(c(1:length(data))/(length(data) + 1)), type = list(4), prob = as.list(fit.mle.boot$t[,
1]), kappa = as.list(fit.mle.boot$t[, 2]), delta = as.list(fit.mle.boot$t[, 3]), sigma = as.list(fit.mle.boot$t[,
4]), xi = as.list(fit.mle.boot$t[, 5]))
q.mle.L <- apply(q.mle.boot, 1, quantile, 0.025, na.rm = TRUE)
q.mle.U <- apply(q.mle.boot, 1, quantile, 0.975, na.rm = TRUE)
}
}
}
}
if (any(method == "pwm")) {
if (any(method == "mle")) {
fits <- list(pwm = fit.pwm, mle = fit.mle)
if (confint) {
CIs <- list(pwm = CIs.pwm, mle = CIs.mle)
} else {
CIs <- NULL
}
} else {
fits <- list(pwm = fit.pwm)
if (confint & plots) {
CIs <- list(pwm = CIs.pwm)
} else {
CIs <- NULL
}
}
} else {
if (any(method == "mle")) {
fits <- list(mle = fit.mle)
if (confint) {
CIs <- list(mle = CIs.mle)
} else {
CIs <- NULL
}
} else {
fits <- NULL
CIs <- NULL
}
}
if (plots) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
mat <- matrix(c(1:(1 + plots)), nrow = 1, ncol = 1 + plots, byrow = TRUE)
layout(mat)
par(mar = c(4, 4, 1, 1))
hist(data, breaks = c(0:30) * max(data + 10^-1)/30, freq = FALSE, xlim = range(x),
xlab = "Rainfall amounts [mm]", ylab = "Density",
main = "", col = "lightgrey")
dens <- density(data, from = 0)
lines(dens$x, dens$y, col = "black")
abline(v = censoring, col = "lightgrey")
if (any(method == "pwm")) {
lines(x, dextgp.pwm, col = "red", lty = 2)
}
if (any(method == "mle")) {
lines(x, dextgp.mle, col = "blue")
}
if (any(method == "pwm")) {
plot(data, qextgp.pwm, asp = 1, xlab = "Empirical quantiles", ylab = "Fitted quantiles", ylim = range(qextgps, na.rm = TRUE),
type = "n")
} else {
if (any(method == "mle")) {
plot(data, qextgp.mle, asp = 1, xlab = "Empirical quantiles", ylab = "Fitted quantiles", ylim = range(qextgps, na.rm = TRUE),
type = "n")
}
}
if (confint) {
if (any(method == "pwm")) {
polygon(x = c(sort(data), sort(data, decreasing = TRUE)), y = c(q.pwm.L, q.pwm.U[length(q.pwm.U):1]), col = rgb(1,
0, 0, alpha = 0.1), lty = 2, border = rgb(1, 0, 0, alpha = 0.5))
}
if (any(method == "mle")) {
polygon(x = c(sort(data), sort(data, decreasing = TRUE)), y = c(q.mle.L, q.mle.U[length(q.mle.U):1]), col = rgb(0,
0, 1, alpha = 0.1), lty = 1, border = rgb(0, 0, 1, alpha = 0.5))
}
}
if (any(method == "pwm")) {
lines(sort(data), qextgp.pwm, lty = 2, type = "b", pch = 20, col = "red")
}
if (any(method == "mle")) {
lines(sort(data), qextgp.mle, lty = 1, type = "b", pch = 20, col = "blue")
}
abline(0, 1, col = "lightgrey")
}
return(list(fit = fits, confint = CIs))
}
#' Fit an extended generalized Pareto distribution of Naveau et al.
#'
#' Deprecated function name to fit an extended generalized Pareto family. The user should call \code{\link[mev]{fit.extgp}} instead.
#'
#' @seealso \code{\link[mev]{fit.extgp}}
#' @export
#' @keywords internal
#' @inheritParams fit.extgp
egp2.fit <- function(data, model = 1, method = c("mle", "pwm"), init, censoring = c(0, Inf), rounded = 0, CI = FALSE, R = 1000, ncpus = 1, plots = TRUE) {
fit.extgp(data = data, model = model, method = method, init = init, censoring = censoring, rounded = rounded,
confint = CI, R = R, ncpus = ncpus, plots = plots)
}
################################################################ ### Distribution/Density/Quantile/Random generator functions ### ###
### different types of distribution G(u):
### type=0: G(u)=u type=1: G(u)=u^kappa
### type=2: G(u)=1-V_delta((1-u)^delta)
### type=3: G(u)=(1-V_delta((1-u)^delta))^(kappa/2)
### type=4: G(u)=p*u^kappa + (1-p)*u^delta
#' @title Carrier distribution for the extended GP distributions of Naveau et al.
#' @description Density, distribution function, quantile function and random number
#' generation for the carrier distributions of the extended Generalized Pareto distributions.
#' @name extgp.G
#' @seealso \code{\link{extgp}}
#' @param u vector of observations (\code{dextgp.G}), probabilities (\code{qextgp.G}) or quantiles (\code{pextgp.G}), in \eqn{[0,1]}
#' @param prob mixture probability for model \code{type} \code{4}
#' @param kappa shape parameter for \code{type} \code{1}, \code{3} and \code{4}
#' @param delta additional parameter for \code{type} \code{2}, \code{3} and \code{4}
#' @param type integer between 0 to 5 giving the model choice
#' @param log logical; should the log-density be returned (default to \code{FALSE})?
#' @param n sample size
#' @param unifsamp sample of uniform; if provided, the data will be used in place of new uniform random variates
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring
#' @param direct logical; which method to use for sampling in model of \code{type} \code{4}?
#' @author Raphael Huser and Philippe Naveau
#'
#' @section Usage: \code{pextgp.G(u, type=1, prob, kappa, delta)}
#' @section Usage: \code{dextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA, log=FALSE)}
#' @section Usage: \code{qextgp.G(u, type=1, prob=NA, kappa=NA, delta=NA)}
#' @section Usage: \code{rextgp.G(n, prob=NA, kappa=NA, delta=NA,
#' type=1, unifsamp=NULL, direct=FALSE, censoring=c(0,1))}
NULL
#' @rdname extgp-functions
#' @export
#' @keywords internal
pextgp.G <- function(u, type = 1, prob, kappa, delta) {
type <- as.integer(type[1])
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
if (type %in% c(1, 3, 4)) {
if(missing(kappa)){
stop("Argument \"kappa\" missing.")
} else{
if(!isTRUE(all(is.numeric(kappa),
length(kappa) == 1L,
kappa > 0))){
stop("Invalid \"kappa\" parameter.")
}
}
}
if (type %in% c(2, 3, 4)){
if(missing(delta)) {
stop("Argument \"delta\" missing.")
} else{
if(!isTRUE(all(is.numeric(delta),
length(delta) == 1L,
delta > 0))){
stop("Invalid \"delta\" parameter.")
}
}
}
if (type == 4){
if(missing(prob)) {
stop("Argument \"prob\" missing.")
} else{
if(!isTRUE(all(is.numeric(prob),
length(prob) == 1L,
prob >= 0,
prob <= 1))){
stop("Invalid \"prob\" parameter.")
}
}
}
u <- pmax(pmin(u, 1), 0)
if (type == 0) {
return(u)
} else if (type == 1) {
return(u^kappa)
} else if (type == 2) {
return(1 - pbeta((1 - u)^delta, 1/delta, 2))
} else if (type == 3) {
return((1 - pbeta((1 - u)^delta, 1/delta, 2))^(kappa/2))
} else if (type == 4) {
return(prob * u^kappa + (1 - prob) * u^delta)
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
dextgp.G <- function(u, type = 1, prob = NA, kappa = NA, delta = NA, log = FALSE) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob\" missing.")
}
if (log == FALSE) {
if (type == 0) {
return(1)
} else if (type == 1) {
return(kappa * u^(kappa - 1))
} else if (type == 2) {
return(dbeta((1 - u)^delta, 1/delta, 2) * delta * (1 - u)^(delta - 1))
} else if (type == 3) {
return((kappa/2) * (1 - pbeta((1 - u)^delta, 1/delta, 2))^(kappa/2 - 1) * dbeta((1 - u)^delta, 1/delta, 2) * delta * (1 -
u)^(delta - 1))
} else if (type == 4) {
return(prob * kappa * u^(kappa - 1) + (1 - prob) * delta * u^(delta - 1))
}
} else {
if (type == 0) {
return(0)
} else if (type == 1) {
return(log(kappa) + (kappa - 1) * log(u))
} else if (type == 2) {
return(dbeta((1 - u)^delta, 1/delta, 2, log = TRUE) + log(delta) + (delta - 1) * log(1 - u))
} else if (type == 3) {
return(log(kappa/2) + (kappa/2 - 1) * log(1 - pbeta((1 - u)^delta, 1/delta, 2)) + dbeta((1 - u)^delta, 1/delta, 2, log = TRUE) +
log(delta) + (delta - 1) * log(1 - u))
} else if (type == 4) {
return(log(prob * kappa * u^(kappa - 1) + (1 - prob) * delta * u^(delta - 1)))
}
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
qextgp.G <- function(u, type = 1, prob = NA, kappa = NA, delta = NA) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob\" missing.")
}
if (type == 0) {
return(u)
} else if (type == 1) {
return(u^(1/kappa))
} else if (type == 2) {
return(1 - qbeta(1 - u, 1/delta, 2)^(1/delta))
} else if (type == 3) {
return(1 - qbeta(1 - u^(2/kappa), 1/delta, 2)^(1/delta))
} else if (type == 4) {
dummy.func <- function(u, p, prob = NA, kappa = NA, delta = NA) {
return(pextgp.G(u = u, prob = prob, kappa = kappa, delta = delta, type = 4) - p)
}
find.root <- function(u, prob = NA, kappa = NA, delta = NA) {
return(uniroot(dummy.func, interval = c(0, 1), p = u, prob = prob, kappa = kappa, delta = delta)$root)
}
return(sapply(u, FUN = find.root, prob = prob, kappa = kappa, delta = delta))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
rextgp.G <- function(n, prob = NA, kappa = NA, delta = NA, type = 1, unifsamp = NULL, direct = FALSE, censoring = c(0, 1)) {
if (!type %in% 1:5) {
stop("Invalid \"type\" argument")
}
type <- type[1]
if (type %in% c(1, 3, 4) && missing(kappa)) {
stop("Argument \"kappa\" missing.")
}
if (type %in% c(2, 3, 4) && missing(delta)) {
stop("Argument \"delta\" missing.")
}
if (type == 4 && missing(prob)) {
stop("Argument \"prob missing.")
}
if (is.null(unifsamp)) {
unifsamp <- runif(n)
}
if (type != 4 | (type == 4 & direct)) {
if (censoring[1] == 0 & censoring[2] == 1) {
return(qextgp.G(unifsamp, prob = prob, kappa = kappa, delta = delta, type = type))
} else {
x.L <- pextgp.G(censoring[1], prob = prob, kappa = kappa, delta = delta, type = type)
x.U <- pextgp.G(censoring[2], prob = prob, kappa = kappa, delta = delta, type = type)
return(qextgp.G(x.L + (x.U - x.L) * unifsamp, prob = prob, kappa = kappa, delta = delta, type = type))
}
} else if (type == 4 & !direct) {
if (censoring[1] == 0 & censoring[2] == 1) {
components <- sample(x = c(1, 2), size = n, replace = TRUE, prob = c(prob, 1 - prob))
res <- c()
res[components == 1] <- qextgp.G(unifsamp[components == 1], prob = NA, kappa = kappa, delta = NA, type = 1)
res[components == 2] <- qextgp.G(unifsamp[components == 2], prob = NA, kappa = delta, delta = NA, type = 1)
return(res)
} else {
x.L <- pextgp.G(censoring[1], prob = prob, kappa = kappa, delta = delta, type = type)
x.U <- pextgp.G(censoring[2], prob = prob, kappa = kappa, delta = delta, type = type)
prop1 <- prob * (pextgp.G(censoring[2], prob = NA, kappa = kappa, delta = NA, type = 1) - pextgp.G(censoring[1], prob = NA,
kappa = kappa, delta = NA, type = 1))
prop2 <- (1 - prob) * (pextgp.G(censoring[2], prob = NA, kappa = delta, delta = NA, type = 1) - pextgp.G(censoring[1], prob = NA,
kappa = delta, delta = NA, type = 1))
new.prob <- prop1/(prop1 + prop2)
components <- sample(x = c(1, 2), size = n, replace = TRUE, prob = c(new.prob, 1 - new.prob))
res <- c()
res[components == 1] <- qextgp.G(x.L + (x.U - x.L) * unifsamp[components == 1], prob = NA, kappa = kappa, delta = NA, type = 1)
res[components == 2] <- qextgp.G(x.L + (x.U - x.L) * unifsamp[components == 2], prob = NA, kappa = delta, delta = NA, type = 1)
return(res)
}
}
}
### different types of extended GP type=0:
### exact GP ---> 2 parameters (sigma,xi)
### type=1: extgp with G(u)=u^kappa ---> 3 parameters(kappa,sigma,xi)
### type=2: extgp with G(u)=1-V_delta((1-u)^delta) ---> 3 parameters (delta,sigma,xi)
### type=3: extgp with G(u)=(1-V_delta((1-u)^delta))^(kappa/2) ---> 4 parameters (kappa,delta,sigma,xi)
### type=4: extgp with mixture G(u)=p*u^kappa + (1-p)*u^delta ---> 5 parameters (p,kappa,delta,sigma,xi)
#' @title Extended generalised Pareto families of Naveau et al. (2016)
#'
#' @description Density function, distribution function, quantile function and random generation for the extended generalized Pareto distribution (GPD) with scale and shape parameters.
#'
#' @details The extended generalized Pareto families proposed in Naveau \emph{et al.} (2016)
#'retain the tail index of the distribution while being compliant with the theoretical behavior of extreme
#'low rainfall. There are five proposals, the first one being equivalent to the GP distribution.
#'\itemize{
#' \item \code{type} 0 corresponds to uniform carrier, \eqn{G(u)=u}.
#' \item \code{type} 1 corresponds to a three parameters family, with carrier \eqn{G(u)=u^\kappa}.
#' \item \code{type} 2 corresponds to a three parameters family, with carrier \eqn{G(u)=1-V_\delta((1-u)^\delta)}.
#' \item \code{type} 3 corresponds to a four parameters family, with carrier \deqn{G(u)=1-V_\delta((1-u)^\delta))^{\kappa/2}}.
#' \item \code{type} 4 corresponds to a five parameter model (a mixture of \code{type} 2, with \eqn{G(u)=pu^\kappa + (1-p)*u^\delta}
#' }
#' @name extgp
#' @param q vector of quantiles
#' @param x vector of observations
#' @param p vector of probabilities
#' @param n sample size
#' @param prob mixture probability for model \code{type} \code{4}
#' @param kappa shape parameter for \code{type} \code{1}, \code{3} and \code{4}
#' @param delta additional parameter for \code{type} \code{2}, \code{3} and \code{4}
#' @param sigma scale parameter
#' @param xi shape parameter
#' @param type integer between 0 to 5 giving the model choice
#' @param step function of step size for discretization with default \code{0}, corresponding to continuous quantiles
#' @param log logical; should the log-density be returned (default to \code{FALSE})?
#' @param unifsamp sample of uniform; if provided, the data will be used in place of new uniform random variates
#' @param censoring numeric vector of length 2 containing the lower and upper bound for censoring
#'
#' @section Usage: \code{pextgp(q, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)}
#' @section Usage: \code{dextgp(x, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, log=FALSE)}
#' @section Usage: \code{qextgp(p, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1)}
#' @section Usage: \code{rextgp(n, prob=NA, kappa=NA, delta=NA, sigma=NA, xi=NA, type=1, unifsamp=NULL, censoring=c(0,Inf))}
#' @author Raphael Huser and Philippe Naveau
#' @references Naveau, P., R. Huser, P. Ribereau, and A. Hannart (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, \emph{Water Resour. Res.}, 52, 2753-2769, \code{doi:10.1002/2015WR018552}.
NULL
#' Extended GP functions
#'
#' These functions are documented in \code{\link{extgp}} and in \code{\link{extgp.G}} for the carrier distributions supported in the unit interval.
#' @name extgp-functions
#' @export
#' @seealso \code{\link{extgp}}, \code{\link{extgp.G}}
#' @keywords internal
pextgp <- function(q, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1) {
return(pextgp.G(mev::pgp(q, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type))
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
dextgp <- function(x, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, log = FALSE) {
if (log == FALSE) {
return(dextgp.G(mev::pgp(x, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type) * mev::dgp(x, scale = sigma,
shape = xi))
} else {
return(dextgp.G(mev::pgp(x, scale = sigma, shape = xi), prob = prob, kappa = kappa, delta = delta, type = type, log = TRUE) + mev::dgp(x,
scale = sigma, shape = xi, log = TRUE))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
qextgp <- function(p, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, step = 0) {
stopifnot(length(step) == 1L, step >= 0, is.finite(step))
if(isTRUE(all.equal(step, 0, check.attributes = FALSE))){
return(mev::qgp(qextgp.G(p, prob = prob, kappa = kappa, delta = delta, type = type), scale = sigma, shape = xi))
} else{
qcont <- mev::qgp(qextgp.G(p, prob = prob, kappa = kappa, delta = delta, type = type), scale = sigma, shape = xi)
return(step*floor((qcont-step)/step))
}
}
#' @rdname extgp-functions
#' @export
#' @keywords internal
rextgp <- function(n, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, unifsamp = NULL, censoring = c(0, Inf)) {
if (censoring[1] == 0 & censoring[2] == Inf) {
return(mev::qgp(rextgp.G(n, prob = prob, kappa = kappa, delta = delta, type = type, unifsamp, censoring = c(0, 1)), scale = sigma,
shape = xi))
} else {
H.L <- mev::pgp(censoring[1], loc = 0, scale = sigma, shape = xi)
H.U <- mev::pgp(censoring[2], loc = 0, scale = sigma, shape = xi)
return(mev::qgp(rextgp.G(n, prob = prob, kappa = kappa, delta = delta, type = type, unifsamp, censoring = c(H.L, H.U)), scale = sigma,
shape = xi))
}
}
######################### ### Inference via PWM ### ###
extgp.pwm <- function(orders, prob = NA, kappa = NA, delta = NA, sigma = NA, xi = NA, type = 1, censoring = c(0, Inf), empiric = FALSE,
unifsamp = NULL, NbSamples = 10^4, N = 200) {
H.L <- ifelse(censoring[1] > 0, mev::pgp(censoring[1], loc = 0, scale = sigma, shape = xi), 0)
H.U <- ifelse(censoring[2] < Inf, mev::pgp(censoring[2], loc = 0, scale = sigma, shape = xi), 1)
F.L <- ifelse(censoring[1] > 0, pextgp(censoring[1], prob = prob, kappa = kappa, delta = delta, sigma = sigma, xi = xi, type = type),
0)
F.U <- ifelse(censoring[2] < Inf, pextgp(censoring[2], prob = prob, kappa = kappa, delta = delta, sigma = sigma, xi = xi, type = type),
1)
prob.LU <- F.U - F.L
if (!empiric) {
E2 <- ((1 - F.L)^(orders + 1) - (1 - F.U)^(orders + 1))/(prob.LU * (orders + 1))
if (type == 0) {
E1 <- (H.U - H.L)^(orders - xi)/(orders - xi + 1)
return((sigma/xi) * (E1 - E2))
} else if (type == 1) {
beta.coefs <- beta(c(1:(max(orders) + 1)) * kappa, 1 - xi)
probs.beta <- pbeta(H.U, shape1 = c(1:(max(orders) + 1)) * kappa, shape2 = 1 - xi) - pbeta(H.L, shape1 = c(1:(max(orders) +
1)) * kappa, shape2 = 1 - xi)
E1 <- c()
for (i in 1:length(orders)) {
E1[i] <- (kappa/prob.LU) * sum((-1)^c(0:orders[i]) * choose(orders[i], c(0:orders[i])) * beta.coefs[1:(orders[i] +
1)] * probs.beta[1:(orders[i] + 1)])
}
return((sigma/xi) * (E1 - E2))
} else if (type == 2) {
E1 <- c()
for (i in 1:length(orders)) {
expos1 <- orders[i] - xi + delta * c(0:orders[i]) + 1
expos2 <- orders[i] - xi + delta * (c(0:orders[i]) + 1) + 1
numer.j <- expos1 * ((1 - H.U)^expos2 - (1 - H.L)^expos2) + expos2 * ((1 - H.L)^expos1 - (1 - H.U)^expos1)
denom.j <- expos1 * expos2
E1[i] <- (1/(delta^(orders[i] + 1) * prob.LU)) * sum((-1)^c(0:orders[i]) * choose(orders[i], c(0:orders[i])) * (1 +
delta)^(orders[i] - c(0:orders[i]) + 1) * numer.j/denom.j)
}
return((sigma/xi) * (E1 - E2))
} else if (type == 3) {
print("No explicit formula for PWMs for this model...")
return(NA)
} else if (type == 4) {
E1 <- c()
for (i in 1:length(orders)) {
l.vals <- matrix(c(0:orders[i]), nrow = orders[i] + 1, ncol = orders[i] + 1, byrow = TRUE)
m.vals <- matrix(c(0:orders[i]), nrow = orders[i] + 1, ncol = orders[i] + 1)
inds <- l.vals < m.vals
l.vals[inds] <- m.vals[inds] <- NA
const <- choose(orders[i], l.vals) * choose(l.vals, m.vals) * (-1)^l.vals * prob^m.vals * (1 - prob)^(l.vals - m.vals)
beta.coefs1 <- beta(kappa * (m.vals + 1) + delta * (l.vals - m.vals), 1 - xi)
beta.coefs2 <- beta(kappa * m.vals + delta * (l.vals - m.vals + 1), 1 - xi)
probs.beta1 <- pbeta(H.U, shape1 = kappa * (m.vals + 1) + delta * (l.vals - m.vals), shape2 = 1 - xi) - pbeta(H.L,
shape1 = kappa * (m.vals + 1) + delta * (l.vals - m.vals), shape2 = 1 - xi)
probs.beta2 <- pbeta(H.U, shape1 = kappa * m.vals + delta * (l.vals - m.vals + 1), shape2 = 1 - xi) - pbeta(H.L, shape1 = kappa *
m.vals + delta * (l.vals - m.vals + 1), shape2 = 1 - xi)
E1[i] <- (1/prob.LU) * sum(const * (prob * kappa * beta.coefs1 * probs.beta1 + (1 - prob) * delta * beta.coefs2 *
probs.beta2), na.rm = TRUE)
}
return((sigma/xi) * (E1 - E2))
}
} else {
if (is.null(unifsamp)) {
unifsamp <- runif(NbSamples)
}
if (censoring[1] == 0 & censoring[2] == Inf) {
V <- rextgp.G(length(unifsamp), prob, kappa, delta, type, unifsamp, direct = TRUE, censoring = c(0, 1))
} else {
V <- rextgp.G(length(unifsamp), prob, kappa, delta, type, unifsamp, direct = TRUE, censoring = c(H.L, H.U))
}
X <- mev::qgp(V, scale = sigma, shape = xi)
res <- c()
for (i in 1:length(orders)) {
res[i] <- mean(X * (1 - (F.L + (F.U - F.L) * unifsamp))^orders[i], na.rm = TRUE)
}
return(res)
}
}
extgp.fit.pwm <- function(x, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf), empiric = FALSE,
unifsamp = NULL, NbSamples = 10^4) {
Fn = ecdf(x)
inds <- x > censoring[1] & x < censoring[2]
mu0hat = mean(x[inds])
mu1hat = mean(x[inds] * (1 - Fn(x[inds])))
mu2hat = mean(x[inds] * (1 - Fn(x[inds]))^2)
mu3hat = mean(x[inds] * (1 - Fn(x[inds]))^3)
mu4hat = mean(x[inds] * (1 - Fn(x[inds]))^4)
if (type == 0) {
if (censoring[1] == 0 & censoring[2] == Inf) {
xihat <- (mu0hat - 4 * mu1hat)/(mu0hat - 2 * mu1hat)
sigmahat <- mu0hat * (1 - xihat)
thetahat <- c(sigmahat, xihat)
} else {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:1), sigma = theta[1], xi = theta[2], type = 0, censoring = censoring, empiric = empiric,
unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat)
return(matrix(pwm.theor - pwm.empir, ncol = 2))
}
theta0 <- c(sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 10^(-6)), upper = c(Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
}
names(thetahat) <- c("sigma", "xi")
return(thetahat)
} else if (type == 1) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:2), kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat)
return(matrix(pwm.theor - pwm.empir, ncol = 3))
}
theta0 <- c(kappa0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 10^(-6)), upper = c(Inf, Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("kappa", "sigma", "xi")
return(thetahat)
} else if (type == 2) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:2), delta = theta[1], sigma = theta[2], xi = theta[3], type = 2, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat)
return(matrix(pwm.theor - pwm.empir, ncol = 3))
}
theta0 <- c(delta0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 10^(-6)), upper = c(100, Inf, 0.99), vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("delta", "sigma", "xi")
return(thetahat)
} else if (type == 3) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:3), kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4], type = 3,
censoring = censoring, empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat, mu3hat)
return(matrix(pwm.theor - pwm.empir, ncol = 4))
}
theta0 <- c(kappa0, delta0, sigma0, xi0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(1e-04, 1e-04, 1e-04, 10^(-6)), upper = c(Inf, 100, Inf, 0.99),
vcov = "iid")
thetahat <- res$coefficients
names(thetahat) <- c("kappa", "delta", "sigma", "xi")
return(thetahat)
} else if (type == 4) {
fct <- function(theta, x) {
pwm.theor <- extgp.pwm(orders = c(0:4), prob = theta[1], kappa = theta[2], delta = theta[2] + theta[3], sigma = theta[4],
xi = theta[5], type = 4, censoring = censoring, empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples)
pwm.empir <- c(mu0hat, mu1hat, mu2hat, mu3hat, mu4hat)
return(matrix(pwm.theor - pwm.empir, ncol = 5))
}
res0 <- extgp.fit.pwm(x[x > quantile(x, 0.9)] - quantile(x, 0.9), type = 0)
if (is.na(prob0)) {
prob0 <- 0.5
}
if (is.na(kappa0)) {
kappa0 <- mean(x[x < quantile(x, 0.1)])/(quantile(x, 0.1) - mean(x[x < quantile(x, 0.1)]))
}
if (is.na(delta0)) {
delta0 <- kappa0 + 0.01
}
Ddelta0 <- delta0 - kappa0
if (is.na(sigma0)) {
sigma0 <- res0[1]
}
if (is.na(xi0)) {
xi0 <- res0[2]
}
theta0 <- c(prob0, kappa0, Ddelta0, sigma0, xi0)
names(theta0) <- c("prob", "kappa", "Ddelta", "sigma", "xi")
# print(theta0)
res <- gmm::gmm(fct, x, theta0, optfct = "nlminb", lower = c(0, 1e-04, 0, 1e-04, 10^(-6)), upper = c(1, Inf, Inf, Inf, 0.99),
vcov = "iid")
thetahat <- res$coefficients
thetahat[3] <- thetahat[2] + thetahat[3]
names(thetahat) <- c("prob", "kappa", "delta", "sigma", "xi")
return(thetahat)
}
}
extgp.fit.pwm.boot <- function(data, i, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf),
empiric = FALSE, unifsamp = NULL, NbSamples = 10^4) {
return(extgp.fit.pwm(data[i], type = type, prob0 = prob0, kappa0 = kappa0, delta0 = delta0, sigma0 = sigma0, xi0 = xi0, censoring = censoring,
empiric = empiric, unifsamp = unifsamp, NbSamples = NbSamples))
}
######################## ### Inference via ML ### ###
extgp.nll <- function(theta, x, censoring, rounded, type) {
x.cens1 <- x[x < censoring[1]]
x.cens2 <- x[x > censoring[2]]
x.not.cens <- x[x >= censoring[1] & x <= censoring[2]]
if (type == 0) {
if (theta[1] <= 0 | theta[2] <= 10^(-6) | theta[2] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], sigma = theta[1], xi = theta[2], type = 0), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], sigma = theta[1], xi = theta[2], type = 0), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, sigma = theta[1], xi = theta[2], type = 0, log = TRUE),
na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, sigma = theta[1], xi = theta[2], type = 0) - pextgp(x.not.cens,
sigma = theta[1], xi = theta[2], type = 0))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 1) {
if (theta[1] <= 0 | theta[2] <= 0 | theta[3] <= 10^(-6) | theta[3] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1),
0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1),
1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1,
log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, kappa = theta[1], sigma = theta[2], xi = theta[3],
type = 1) - pextgp(x.not.cens, kappa = theta[1], sigma = theta[2], xi = theta[3], type = 1))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 2) {
if (theta[1] <= 0 | theta[1] > 100 | theta[2] <= 0 | theta[3] <= 10^(-6) | theta[3] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], delta = theta[1], sigma = theta[2], xi = theta[3], type = 2),
0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], delta = theta[1], sigma = theta[2], xi = theta[3], type = 2),
1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, delta = theta[1], sigma = theta[2], xi = theta[3], type = 2,
log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, delta = theta[1], sigma = theta[2], xi = theta[3],
type = 2) - pextgp(x.not.cens, delta = theta[1], sigma = theta[2], xi = theta[3], type = 2))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 3) {
if (theta[1] <= 0 | theta[2] <= 0 | theta[2] > 100 | theta[3] <= 0 | theta[4] <= 10^(-6) | theta[4] > 0.99) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3, log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, kappa = theta[1], delta = theta[2], sigma = theta[3],
xi = theta[4], type = 3) - pextgp(x.not.cens, kappa = theta[1], delta = theta[2], sigma = theta[3], xi = theta[4],
type = 3))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
} else if (type == 4) {
if (theta[1] < 0 | theta[1] > 1 | theta[2] <= 0 | theta[3] <= 0 | theta[4] <= 0 | theta[5] <= 10^(-6) | theta[5] > 0.99 |
theta[3] < theta[2]) {
return(Inf)
} else {
censor1 <- ifelse(censoring[1] > 0, pextgp(censoring[1], prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4), 0)
censor2 <- ifelse(censoring[2] < Inf, pextgp(censoring[2], prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4), 1)
contrib.cens1 <- ifelse(length(x.cens1) > 0, length(x.cens1) * log(censor1), 0)
contrib.cens2 <- ifelse(length(x.cens2) > 0, length(x.cens2) * log(1 - censor2), 0)
contrib.not.cens <- ifelse(rounded == 0, sum(dextgp(x.not.cens, prob = theta[1], kappa = theta[2], delta = theta[3], sigma = theta[4],
xi = theta[5], type = 4, log = TRUE), na.rm = TRUE), sum(log(pextgp(x.not.cens + rounded, prob = theta[1], kappa = theta[2],
delta = theta[3], sigma = theta[4], xi = theta[5], type = 4) - pextgp(x.not.cens, prob = theta[1], kappa = theta[2],
delta = theta[3], sigma = theta[4], xi = theta[5], type = 4))))
return(-(contrib.cens1 + contrib.not.cens + contrib.cens2))
}
}
}
extgp.fit.ml <- function(x, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf), rounded = 0.1,
print = FALSE) {
if (type == 0) {
theta0 <- c(sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("sigma", "xi")
return(thetahat)
} else if (type == 1) {
theta0 <- c(kappa0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("kappa", "sigma", "xi")
return(thetahat)
} else if (type == 2) {
theta0 <- c(delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("delta", "sigma", "xi")
return(thetahat)
} else if (type == 3) {
theta0 <- c(kappa0, delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("kappa", "delta", "sigma", "xi")
return(thetahat)
} else if (type == 4) {
theta0 <- c(prob0, kappa0, delta0, sigma0, xi0)
opt <- optim(par = theta0, fn = extgp.nll, x = x, censoring = censoring, rounded = rounded, type = type, method = "Nelder-Mead",
control = list(maxit = 1000), hessian = FALSE)
if (print) {
print(opt)
}
thetahat <- opt$par
names(thetahat) <- c("prob", "kappa", "delta", "sigma", "xi")
return(thetahat)
}
}
extgp.fit.ml.boot <- function(data, i, type = 1, prob0 = NA, kappa0 = NA, delta0 = NA, sigma0 = NA, xi0 = NA, censoring = c(0, Inf),
rounded = 0.1, print = FALSE) {
return(extgp.fit.ml(data[i], type = type, prob0 = prob0, kappa0 = kappa0, delta0 = delta0, sigma0 = sigma0, xi0 = xi0, censoring = censoring,
rounded = rounded, print = print))
}
Try the mev package in your browser
Any scripts or data that you put into this service are public.
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.