GEV.ipynb

https://github.com/pymc-devs/pymc-examples/blob/main/examples/case_studies/GEV.ipynb

Generalized Extreme Value Distribution

:::{post} Sept 27, 2022 :tags: extreme, inference, posterior predictive :category: beginner :author: Colin Caprani :::

Introduction

The Generalized Extreme Value (GEV) distribution is a meta-distribution containing the Weibull, Gumbel, and Frechet families of extreme value distributions. It is used for modelling the distribution of extremes (maxima or minima) of stationary processes, such as the annual maximum wind speed, annual maximum truck weight on a bridge, and so on, without needing a priori decision on the tail behaviour.

Following the parametrization used in {cite:t}coles2001gev, the GEV distribution for maxima is given by:

G(x)=exp⁡{[1−ξ(x−μσ)]−1ξ}G(x) = \exp \left\{ \left[ 1 - \xi \left( \frac{x - \mu}{\sigma} \right) \right]^{-\frac{1}{\xi}} \right\}

when:

  • $\xi < 0$ we get the Weibull distribution with a bounded upper tail;
  • $\xi = 0$, in the limit, we get the Gumbel distribution, unbonded in both tails;
  • $\xi > 0$, we get the Frechet distribution which is bounded in the lower tail.

Note that this parametrization of the shape parameter $\xi$ is opposite in sign to that used in SciPy (where it is denoted c). Further, the distribution for minima is readily examined by studying the distribution of maxima of the negative of the data.

We will use the example of the Port Pirie annual maximum sea-level data used in {cite:t}coles2001gev, and compare with the frequentist results presented there.

Data

The Port Pirie data is provided by {cite:t}coles2001gev, and repeated here:

Modelling & Prediction

In the modelling we wish to do two thing:

  • parameter inference on the GEV parameters, based on some fairly non-informative priors, and;
  • prediction of the 10-year return level.

Predictions of extreme values considering parameter uncertainty are easily accomplished in the Bayesian setting. It is interesting to compare this ease to the difficulties encountered by {cite:t}caprani2010gev in doing this in a frequentist setting. In any case, the predicted value at a probability of exceedance $p$ is given by:

xp=μ−σξ{1−[−log⁡(1−p)]}x_p = \mu - \frac{\sigma}{\xi} \left\{1 - \left[-\log\left(1-p\right)\right] \right\}

This is a deterministic function of the parameter values, and so is accomplished using pm.Deterministic within the model context.

Consider then, the 10-year return period, for which $p = 1/10$:

And now set up the model using priors estimated from a quick review of the histogram above:

  • $\mu$: there is no real basis for considering anything other than a Normal distribution with a standard deviation limiting negative outcomes;
  • $\sigma$: this must be positive, and has a small value, so use HalfNormal with a unit standard deviation;
  • $\xi$: we are agnostic to the tail behaviour so centre this at zero, but limit to physically reasonable bounds of $\pm 0.6$, and keep it somewhat tight near zero.

Prior Predictive Checks

Let's get a feel for how well our selected priors cover the range of the data:

And we can look at the sampled values of the parameters, using the plot_posterior function, but passing in the idata object and specifying the group to be "prior":

Inference

Press the magic Inference Button$^{\mathrm{TM}}$:

Divergences

The trace exhibits divergences (usually). The HMC/NUTS sampler can have problems when the bounds of support for parameters are approached. And since the bounds of the GEV change with the sign of $\xi$, it is difficult to offer a transformation that resolves this problem. One possible transformation - the Box-Cox - has been proposed by {cite:t}bali2003gev, but {cite:t}caprani2009gev find it numerically unstable, even for just maximum likelihood estimation. In any case, recommendations to alleviate divergence problems are:

  • Increase the target acceptance ratio;
  • Use more informative priors, especially limit the shape parameter to physically reasonable values, typically $\xi \in [-0.5,0.5]$;
  • Decide upon the domain of attraction of the tail (i.e. Weibull, Gumbel, or Frechet), and use that distribution directly.

Inferences

The 95% credible interval range of the parameter estimates is:

And examine the prediction distribution, considering parameter variability (and without needing to assume normality):

And let's compare the prior and posterior predictions of $z_p$ to see how the data has influenced things:

Comparison

To compare with the results given in {cite:t}coles2001gev, we approximate the maximum likelihood estimates (MLE) using the mode of the posterior distributions (the maximum a posteriori or MAP estimate). These are close when the prior is reasonably flat around the posterior estimate.

The MLE results given in {cite:t}coles2001gev are:

(μ^,σ^,ξ^)=(3.87,0.198,−0.050)\left(\hat{\mu}, \hat{\sigma}, \hat{\xi} \right) = \left( 3.87, 0.198, -0.050 \right)

And the variance-covariance matrix of the estimates is:

0.000197 & 0.000410 & -0.000778 \\ -0.00107 & -0.000778 & 0.00965 \end{array} \right] $$

Note that extracting the MLE estimates from our inference involves accessing some of the Arviz back end functions to bash the xarray into something examinable:

The results are a good match, but the benefit of doing this in a Bayesian setting is we get the full posterior joint distribution of the parameters and the return level, essentially for free. Compare this to the loose normality assumption and computational effort to get even the variance-covarince matrix, as done in {cite:t}coles2001gev.

Finally, we examine the pairs plots and see where any difficulties in inference lie using the divergences

Authors

References

:::{bibliography} :filter: docname in docnames :::

Watermark