#' @import foreach
NULL
## Set parameters
hazDWParUpdate <- function(parameters) {
# set threshold dependent on tau
# and compute also where the GPD is pure
parameters$loc <- hazDWSetThreshold(parameters)
parameters$loc2 <- parameters$loc + parameters$epsilon
return(parameters)
}
hazDWSetThreshold <- function(parameters) {
# help function for hazDWParUpdate
tau <- parameters$tau
beta <- parameters$beta
kappa <- parameters$kappa
return((-log(1- tau) * beta^kappa)^(1/kappa))
}
## transition function
linearTransition <- function(x) {
y <- rep(1, length(x))
isBetweenZeroAndOne <- (x > 0 & x < 1)
y[x >= 1] <- 0
y[isBetweenZeroAndOne] <- 1 - x[isBetweenZeroAndOne]
return(y)
}
polynomialTransition <- function(x) {
y <- rep(1, length(x))
isBetweenZeroAndOne <- (x > 0 & x < 1)
y[x >= 1] <- 0
y[isBetweenZeroAndOne] <- 2 * x[isBetweenZeroAndOne]^3 - 3 * x[isBetweenZeroAndOne]^2 + 1
return(y)
}
#smoothTransition <- function(x) {
# # the scaling has to be done before
# a <- smoothTransitionHelpFunction(x)
# b <- smoothTransitionHelpFunction(1 - x)
# return(1 - a/(a + b))
#}
#
#smoothTransitionHelpFunction <- function(x) {
# y <- rep(0, length(x))
# index <- x > 0
# y[index] <- exp(-1/x[index])
# return(y)
#}
#' @export
#' @rdname DistortedWeibull
hazDW <- function(x, parameters) {
a <- hazW(x, parameters)
b <- hazGpd(x, parameters)
c <- polynomialTransition((x-parameters$loc)/parameters$epsilon)
return(c * a + (1 - c) * b)
}
hazW <- function(x, parameters) {
kappa <- parameters$kappa
beta <- parameters$beta
return(kappa/beta * (x / beta)^(kappa - 1))
#return(hazardRateWeibull(x, parameters))
}
hazGpd <- function(x, parameters) {
xi <- parameters$xi
gamma <- parameters$gamma
loc <- parameters$loc
isGreaterLoc <- x >= loc
y <- rep(0, length(x))
y[isGreaterLoc] <- (gamma * loc + xi * (x[isGreaterLoc] - loc))^(-1)
return(y)
}
## Cumulative hazard rate
HazDW <- function(x, parameters) {
xi <- parameters$xi
gamma <- parameters$gamma
loc <- parameters$loc
loc2 <- parameters$loc2
y <- rep(NA, length(x))
isSmaller <- x <= loc
isGreater <- x >+ loc2
isBetween <- !(isSmaller | isGreater)
y[isSmaller] <- HazW(x[isSmaller], parameters)
y[isBetween] <- HazW(loc, parameters) + HazTransitionWGpd(x[isBetween], parameters)
y[isGreater] <- HazW(loc, parameters) + HazTransitionWGpd(loc2, parameters) +
HazGpd(x[isGreater], parameters) - HazGpd(loc2, parameters)
return(y)
}
HazW <- function(x, parameters) {
kappa <- parameters$kappa
beta <- parameters$beta
return((x / beta)^kappa)
}
HazGpd <- function(x, parameters) {
xi <- parameters$xi
gamma <- parameters$gamma
loc <- parameters$loc
scale <- gamma * loc
y <- 1/xi * log(1 + xi * (x - loc)/scale)
return(y)
}
HazTransitionWGpd <- function(x, parameters) {
a <- HazTransitionW(x, parameters)
b <- HazTransitionGpd(x, parameters)
return(a + b)
}
HazTransitionW <- function(x, parameters) {
y <- HazTransitionWHelp(x, parameters) - HazTransitionWHelp(parameters$loc, parameters)
return(y)
}
HazTransitionGpd <- function(x, parameters) {
y <- HazTransitionGpdHelp(x - parameters$loc, parameters) - HazTransitionGpdHelp(0, parameters)
return(y)
}
HazTransitionWHelp <- function(x, parameters) {
kappa <- parameters$kappa
beta <- parameters$beta
loc <- parameters$loc
epsilon <- parameters$epsilon
b <- loc / epsilon
c <- 2 * kappa / (epsilon^3 * (kappa + 3)) * x^(kappa + 3)
d <- -3 * kappa * (2 * b + 1) / (epsilon^2 * (kappa+2)) * x^(kappa + 2)
e <- 6 * loc * kappa * (b + 1) / (epsilon^2 * (kappa + 1)) * x^(kappa + 1)
f <- -(b^2 * (2 * b + 3) - 1) * x^kappa
y <- (c + d + e + f) / beta^kappa
return(y)
}
HazTransitionGpdHelp <- function(x, parameters) {
gamma <- parameters$gamma
xi <- parameters$xi
loc <- parameters$loc
epsilon <- parameters$epsilon
sigma <- gamma * loc
stopifnot(xi > -sigma/epsilon)
a <- (epsilon * xi)^(-1)
X <- xi * x + sigma
b <- -2/3 * a * X^3
c <- (3 * a * sigma + 3/2) * X^2
d <- -6 * sigma * (a * sigma + 1) * X
e <- sigma^2 * ( 2 * a * sigma + 3) * log(X)
y <- a^2 / xi * (b + c + d + e)
return(y)
}
#' Get marginal parameters of the distorted Weibull distribution
#' @param tau Quantile where GPD behaviors sets in
#' @param epsilon Length of transition interval
#' @param betas Weibull scale parameter
#' @param kappas Weibull shape parameter
#' @param gamma GPD dispersion coefficient
#' @param xi GPD shape parameter
#' @param xpar For the spatial and temporal dimension
#' @export
setMarginalParameters <- function(tau, epsilon, betas, kappas, gamma, xi, xpar) {
mParams <- list()
for (site in 1 : xpar$S) {
mParams[[site]] <- list(tau = tau, epsilon = epsilon
,beta = betas[site], kappa = kappas[site]
,gamma = gamma, xi = xi)
}
mParams <- lapply(mParams, hazDWParUpdate)
return(mParams)
}
#' Calculates the Weibull shape parameter
#'
#' that corresponds to a GPD shape parameter xi at quantile tau
#' @inheritParams setMarginalParameters
#' @export
calculateCorrespondingWeibullShape <- function(xi, tau) {
kappa <- (xi * (- log(1 - tau)) + 1)^(-1)
return(kappa)
}
calculatePenUltimateShapeParameter <- function(kappa, tau) {
xi <- (1 - kappa)/kappa * (-log(1 - tau))^(-1)
return(xi)
}
#' Generates correlated, distorted weibull data
#' @param marginalParameters list of marginal parameters
#' @param copula copula object
#' @param xpar named list with `S` the number of sites
#' and `T` the number of days
#' @importFrom copula rCopula
#' @export
generateDWData <- function(marginalParameters, copula, xpar) {
s <- NULL
uniform <- rCopula(xpar$T, copula)
data <- foreach(s = 1 : xpar$S, .combine = "cbind") %dopar% {
return(qDW(uniform[, s], parameters = marginalParameters[[s]]))
}
return(data)
}
#' @importFrom evd qgpd
#' @note 2016-03-10: We switched from POT to evd (which has no lambda)
getTrueReturnLevels <- function(params, xpar, returnPeriods, blockLength) {
i <- NULL
tmp <- foreach(i = 1 : xpar$S, .combine = "cbind") %do% {
loc = params[[i]]$loc
scale = params[[i]]$gamma * loc
shape = params[[i]]$xi
tau = params[[i]]$tau
p = 1 - 1/(returnPeriods * (1-tau) * blockLength)
return(qgpd(p = p, loc = loc, scale = scale, shape = shape))
}
colnames(tmp) = 1 : xpar$S
rownames(tmp) = returnPeriods
return(tmp)
}
#' Distorted Weibull distribution
#' @param parameters Weibull parameters
#' @note TODO: Switch to named parameters
#' @param x,q quantile
#' @param tau Belongs to parameters
#' @param n number of observations
#' @name DistortedWeibull
NULL
#' @rdname DistortedWeibull
#' @export
rDW <- function(n, parameters) {
f <- function(x, parameters, tau) {
return(1 - exp(-HazDW(x, parameters)) - tau)
}
unif <- runif(n)
result <- numeric(n)
for (i in 1 : n) {
result[i] <- uniroot(f, c(0, 1000), parameters, tau = unif[i])$root
}
return(result)
}
#' @rdname DistortedWeibull
#' @export
qDW <- function(tau, parameters) {
f <- function(x, parameters, tau) {
return(1 - exp(-HazDW(x, parameters)) - tau)
}
result <- numeric()
for (i in 1 : length(tau)) {
result[i] <- uniroot(f, c(0, 1000), parameters, tau = tau[i])$root
}
return(result)
}
#' @rdname DistortedWeibull
#' @export
pDW <- function(q, parameters) {
return(1 - exp(-HazDW(q, parameters)))
}
#' @rdname DistortedWeibull
#' @export
dDW <- function(x, parameters) {
return(hazDW(x, parameters) * exp(-HazDW(x, parameters)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.