Nothing
## ----shape--------------------------------------------------------------------
# load package
library(eva)
# A naive implementation of the GEV cumulative density function
pgev_naive <- function(q, loc = 0, scale = 1, shape = 1) {
exp(-(1 + (shape * (q - loc))/scale)^(-1/shape))
}
curve(pgev_naive(1, 0, 1, x), 1e-20, .01, log = "x", n = 1025,
xlab = "Shape", ylab = "GEV CDF", col = "red", lty = 1, lwd = 1)
curve(pgev(1, 0, 1, x), 1e-20, .01, log = "x", n = 1025, col = "blue", lty = 6, lwd = 3, add = TRUE)
# Similarly for the GPD cdf
pgpd_naive <- function(q, loc = 0, scale = 1, shape = 1) {
(1 - (1 + (shape * (q - loc))/scale)^(-1/shape))
}
curve(pgpd_naive(1, 0, 1, x), 1e-20, .01, log = "x", n = 1025,
xlab = "Shape", ylab = "GPD CDF", col = "red", lty = 1, lwd = 1)
curve(pgpd(1, 0, 1, x), 1e-20, .01, log = "x", n = 1025, col = "blue", lty = 6, lwd = 3, add = TRUE)
## ----seqtesting---------------------------------------------------------------
data(lowestoft)
gevrSeqTests(lowestoft, method = "ed")
## ----returnlevel, fig.height = 12, fig.width = 8------------------------------
# Make 250 year return level plot using gevr for r = 1 to 10 with the LoweStoft data
data(lowestoft)
result <- matrix(0, 20, 4)
period <- 250
for(i in 1:10) {
z <- gevrFit(as.matrix(lowestoft[, 1:i]))
y1 <- gevrRl(z, period, conf = 0.95, method = "delta")
y2 <- gevrRl(z, period, conf = 0.95, method = "profile", plot = FALSE)
result[i, 1] <- i
result[i, 2] <- y1$Estimate
result[i, 3:4] <- y1$CI
result[(i + 10), 1] <- i
result[(i + 10), 2] <- y2$Estimate
result[(i + 10), 3:4] <- y2$CI
}
result <- cbind.data.frame(result, c(rep("Delta", 10), rep("Profile", 10)))
colnames(result) <- c("r", "Est", "Lower", "Upper", "Method")
result <- as.data.frame(result)
prof <- subset(result, Method == "Profile")
del <- subset(result, Method == "Delta")
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 1))
plot(prof$r, prof$Est, main = "Profile Likelihood",
xlab = "r", ylab = "250 Year Return Level",
xlim = c(1, 10), ylim = c(4, 7))
polygon(c(rev(prof$r), prof$r), c(rev(prof$Lower), prof$Upper), col = 'grey80', border = NA)
points(prof$r, prof$Est, pch = 19, col = 'black')
lines(prof$r, prof$Est, lty = 'solid', col = 'black')
lines(prof$r, prof$Lower, lty = 'dashed', col = 'red')
lines(prof$r, prof$Upper, lty = 'dashed', col = 'red')
plot(del$r, del$Est, main = "Delta Method",
xlab = "r", ylab = "250 Year Return Level",
xlim = c(1, 10), ylim = c(4, 7))
polygon(c(rev(del$r), del$r), c(rev(del$Lower), del$Upper), col = 'grey80', border = NA)
points(del$r, del$Est, pch = 19, col = 'black')
lines(del$r, del$Est, lty = 'solid', col = 'black')
lines(del$r, del$Lower, lty = 'dashed', col = 'red')
lines(del$r, del$Upper, lty = 'dashed', col = 'red')
par(oldpar)
## ----nonstatfit1--------------------------------------------------------------
set.seed(7)
n <- 100
r <- 10
x <- rgevr(n, r, loc = 100 + 1:n / 50, scale = 1 + 1:n / 100, shape = 0)
## Plot the largest order statistic
plot(x[, 1])
## Creating covariates (linear trend first)
covs <- as.data.frame(seq(1, n, 1))
names(covs) <- c("Trend1")
## Create some unrelated covariates
covs$Trend2 <- rnorm(n)
covs$Trend3 <- 30 * runif(n)
## Use full data
fit_full <- gevrFit(data = x, method = "mle", locvars = covs, locform = ~ Trend1 + Trend2*Trend3,
scalevars = covs, scaleform = ~ Trend1)
## Only use r = 1
fit_top_only <- gevrFit(data = x[, 1], method = "mle", locvars = covs, locform = ~ Trend1 + Trend2*Trend3,
scalevars = covs, scaleform = ~ Trend1)
## ----nonstatfit2--------------------------------------------------------------
## Show summary of estimates
fit_full
fit_top_only
## ----nonstatfit3--------------------------------------------------------------
## Compare AIC of three models
fit_reduced1 <- gevrFit(data = x, method = "mle", locvars = covs, locform = ~ Trend1,
scalevars = covs, scaleform = ~ Trend1)
fit_reduced2 <- gevrFit(data = x, method = "mle", locvars = covs, locform = ~ Trend1,
scalevars = covs, scaleform = ~ Trend1, gumbel = TRUE)
AIC(fit_full)
AIC(fit_reduced1)
AIC(fit_reduced2)
LRT <- as.numeric(2 * (logLik(fit_reduced1) - logLik(fit_reduced2)))
pval <- pchisq(LRT, df = 1, ncp = 0, lower.tail = FALSE, log.p = FALSE)
round(pval, digits = 3)
## ----rfa_example1-------------------------------------------------------------
set.seed(7)
require(SpatialExtremes)
n.site <- 4
n.obs <- 15
## Simulate a max-stable random field
locations <- matrix(runif(2 * n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")
## Smith model
U <- rmaxstab(n.obs, locations, "gauss", cov11 = 16, cov12 = 0, cov22 = 16)
cor(U, method = "spearman")
## ----rfa_example3-------------------------------------------------------------
## Transform to GEV margins
locations <- c(8, 10, 12, 9)
out <- frech2gev(U, loc = 0, scale = 1, shape = 0.2)
out <- out + t(matrix(rep(locations, nrow(out)), ncol = nrow(out)))
out <- as.vector(out)
## Create design matrix for the location parameters
loc <- cbind.data.frame(as.factor(sort(rep(seq(1, n.site, 1), n.obs))))
colnames(loc) <- c("Site")
z <- gevrFit(out, locvars = loc, locform = ~ Site)
z
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.