2

I have fitted my data to a GEV distribution, and I wonder how to find the probability of P(x<=40). Thanks for any help.

library(extRemes)
ams <- c(44.5,43.2,38.1,39.1,32.3,25.4,33.0,32.5,48.5,34.3,45.7,35.3,76.7,34.0,86.6,48.5,59.4,53.3,30.5,42.7,83.3,59.2,37.3,67.3,38.4,47.0,38.1,72.4,40.9,47.0,36.3,85.3,35.6,55.9,44.2,45.2,51.6,59.4,47.8,55.4,42.4,40.1,36.6,47.0,48.8,51.3,39.4,45.7)
fit_mle <- fevd(x=ams, method = "MLE", type="GEV",period.basis = "year")
Yang Yang
  • 858
  • 3
  • 26
  • 49

2 Answers2

2

According to the help page for fevd, section Details:

The GEV df is given by

Pr(X <= x) = G(x) = exp[-(1 + shape*(x - location)/scale)^(-1/shape)]

So you can do the following.

location <- fit_mle$results$par[1]
scale <- fit_mle$results$par[2]
shape <- fit_mle$results$par[3]
x <- 40
exp(-(1 + shape*(x - location)/scale)^(-1/shape))
#    shape 
#0.3381735

Or you can simply use the built-in cummulative distribution function.

pevd(x, location, scale, shape)
#[1] 0.3381735
Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • Thank you so much for your help. However, when I check the arguments of `pevd`, I find that it requires `q`, which is a numeric vector of quantiles. So, can we just use `40` here (as 40 is a value, not a quantile)? – Yang Yang Jan 25 '18 at 19:37
  • @YangYang Yes, in the documentation it calls the argument `q`, but I already had an `x` from the code above, that's just it. `40` *is* a quantile. – Rui Barradas Jan 25 '18 at 19:43
  • I am sorry but I am a bit confused. In my case, 40 is a data point, not quantile. For example, 40 stands for 40 mm rainfall. – Yang Yang Jan 25 '18 at 19:45
  • @YangYang That is a matter of the name you wish to call it. You want `P(X <= 40)` and that is what the code in my answer computes, given the (correct) fit by `fevd`. – Rui Barradas Jan 25 '18 at 19:52
  • OK, thanks for your explanation. I have posted another question at `https://stackoverflow.com/questions/48434236/how-to-plot-a-fevd-object-using-ggplot2-in-r`, and I wonder if you could take a look? Thank you for your time. – Yang Yang Jan 25 '18 at 19:54
2
library(EnvStats)
ams <- c(44.5,43.2,38.1,39.1,32.3,25.4,33.0,32.5,48.5,34.3,45.7,35.3,76.7,34.0,86.6,48.5,59.4,53.3,30.5,42.7,83.3,59.2,37.3,67.3,38.4,47.0,38.1,72.4,40.9,47.0,36.3,85.3,35.6,55.9,44.2,45.2,51.6,59.4,47.8,55.4,42.4,40.1,36.6,47.0,48.8,51.3,39.4,45.7)
fit_gev <- egevd(ams, method = "mle")# Parameters estimation
pgevd(40, location = fit_gev$parameters[[1]], scale = fit_gev$parameters[[2]],
  shape = fit_gev$parameters[[3]])

0.3381751
Marcel
  • 147
  • 4