dgev <- function(x, loc = 0, scale = 1, shape = 0){
t1 <- if(shape == 0) {
exp(-(x - loc)/scale)
} else {
(1 + shape * (x - loc)/scale )^(-1/shape)
}
d <- 1/scale * t1^(shape + 1) * exp(-t1)
# setting desity to NA outside of range
outside <- numeric()
if(shape < 0) outside <- x > loc - scale/shape
if(shape > 0) outside <- x < loc - scale/shape
d[outside] <- NA
return(d)
}
draw_support <- function(y, loc = 0, scale = 1, shape = 0, ...) {
if (shape == 0) segments(x0 = par("xaxp")[1], x1 = par("xaxp")[2], y0 = y, ...)
end <- loc - scale/shape
if (shape < 0) segments(x0 = par("xaxp")[1], x1 = end, y0 = y, ...)
if (shape > 0) segments(x0 = end, x1 = par("xaxp")[2], y0 = y, ...)
points(x = end, y = y, pch = 21, ...)
}
x <- seq(-4, 4, 0.01)
shape <- c("3" = -0.5, "1" = 0, "2" = 0.5)
col <- c("blue", "black", "green")
prob <- seq(0, 1, 0.001)
svg("GevDensity.svg")
par(mar = c(5, 4, 2, 0) + 0.1)
plot(NA, xlim=range(x), ylim=c(-0.05, 0.5), type = "n",
xlab = expression(italic(x)), ylab = "Density",
main = "Generalized extreme value densities",
sub = "horizontal lines denote the support of the distribution")
for(i in seq_along(shape)) {
lines(x, dgev(x, shape = shape[i]), col=col[i], lwd = 2)
draw_support(shape = shape[i], col=col[i], y = -i/100 -0.02 )
}
legend("topleft", lwd = 2, col = col, inset = 0.02,
title = "all with \u03bc = 0, \u03c3 = 1",
legend = paste0("GEV Type ", names(shape), ": \u03be = ", shape))
dev.off()