R/egp3.R In texmex: Statistical Modelling of Extreme Values

```#' @include texmexFamily.R
#' @export egp3
NULL

egp3 <- texmexFamily(name="EGP3",
param = c(lambda=0, phi=0, xi=0), # Standard exponential
density = function(n, param, model){
degp3(n, exp(c(param[, 1])), exp(c(param[, 2])), c(param[, 3]), u=model\$threshold)
},
rng = function(n, param, model){
regp3(n, exp(c(param[, 1])), exp(c(param[, 2])), c(param[, 3]), u=model\$threshold)
},
prob = function(n, param, model){
pegp3(n, exp(c(param[, 1])), exp(c(param[, 2])), c(param[, 3]), u=model\$threshold)
},
quant = function(n, param, model){
qegp3(n, exp(c(param[, 1])), exp(c(param[, 2])), c(param[, 3]), u=model\$threshold)
},
log.lik = function(data, th, ...){
y <- data\$y

X.lambda <- data\$D\$lambda
X.phi <- data\$D\$phi
X.xi <- data\$D\$xi

n.lambda <- ncol(X.lambda)
n.phi <- n.lambda + ncol(X.phi)
n.end <- n.phi + ncol(X.xi)

function(param) {
stopifnot(length(param) == n.end)
lambda <- X.lambda %*% param[1:n.lambda]
phi <- X.phi %*% param[(1 + n.lambda):n.phi]
xi <- X.xi %*% param[(1 + n.phi):n.end]
sum(degp3(y, exp(lambda), exp(phi), xi, u=th, log.d=TRUE))
}
}, # Close log.lik
start = function(data){
c(0, log(sd(data\$y)), rep(1e-05, sum(sapply(data\$D, ncol)) - 2))
}, # Close start
rl = function(m, param, model){
qegp3(1/(m * model\$rate),
exp(param[, 1]), exp(param[, 2]), param[, 3], u=model\$threshold,
lower.tail=FALSE)
},
endpoint = function(param, model){
res <- model\$threshold - exp(param[, 2]) / param[, 3]
res[param[, 3] >= 0] <- Inf
res
},
delta = function(param, m, model){
kappa <- exp(param[1])
sigma <- exp(param[2])
xi <- param[3]

out <- matrix(0, nrow=3, ncol=length(m))

z <- 1 - 1/(m*model\$rate)

dk <- -(z^(1/kappa) * sigma * (1 - z^(1/kappa))^(-xi - 1) * log(z)) / kappa^2
ds <- ((1 - z^(1/kappa))^(-xi) -1 ) / xi

dz1 <- -sigma * (1 - z^(1/kappa))^(-xi) * log(1 - z^(1/kappa)) / xi
dz2 <- -sigma * ((1 - z^(1/kappa))^(-xi) -1) / xi^2
dz <- dz1 + dz2

out[1, ] <- dk * kappa # To account for working with exp(lambda)
out[2, ] <- ds * sigma # To account for working with exp(phi)
out[3, ] <- dz
out
},
resid = function(o){
p <- texmexMakeParams(coef(o), o\$data\$D)
delta <- (o\$data\$y - o\$threshold) / exp(p[, 2]) * p[, 3]
#(texmex:::.log1prel(delta * p[, 3]) * delta)^{exp(p[, 1])} # Standard exponential
r <- (1 - (1 + delta)^(-1/p[, 3]))^exp(p[, 1])
-log(1-r)
} # Close resid
)
```

Try the texmex package in your browser

Any scripts or data that you put into this service are public.

texmex documentation built on May 2, 2019, 5:41 a.m.