{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# Generalized Extreme Value Distribution\n", "\n", ":::{post} Sept 27, 2022\n", ":tags: extreme, inference, posterior predictive\n", ":category: beginner \n", ":author: Colin Caprani \n", ":::" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Introduction\n", "\n", "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.\n", "\n", "Following the parametrization used in {cite:t}coles2001gev, the GEV distribution for maxima is given by:\n", "\n", "$$G(x) = \\exp \\left\\{ \\left[ 1 - \\xi \\left( \\frac{x - \\mu}{\\sigma} \\right) \\right]^{-\\frac{1}{\\xi}} \\right\\}$$\n", "\n", "when:\n", "- $\\xi < 0$ we get the Weibull distribution with a bounded upper tail;\n", "- $\\xi = 0$, in the limit, we get the Gumbel distribution, unbonded in both tails;\n", "- $\\xi > 0$, we get the Frechet distribution which is bounded in the lower tail.\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pymc as pm\n", "import pymc_experimental.distributions as pmx\n", "import pytensor.tensor as pt\n", "\n", "from arviz.plots import plot_utils as azpu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data\n", "The Port Pirie data is provided by {cite:t}coles2001gev, and repeated here:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbdUlEQVR4nO3df5DUdf3A8dfK6eIP7hSKHxcg6CikKDaSheZXCH8MItqUP0oDopwoKTBKBS31bPTUiiElNc3UGQUcTczyRzKWgPkjfkg6OCOpKJdKjGZ3gLkKfL5/mDceIHKw+95deDxm9o/97Hv38+I9N7PP+ewdm8uyLAsAgER2KfcAAMDORXwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBSNeUeYGMbNmyI1157LTp16hS5XK7c4wAAWyHLsli9enXU19fHLrts+dpGxcXHa6+9Fr169Sr3GADANmhqaoqePXtucU3FxUenTp0i4v3ha2tryzwNALA1WlpaolevXq3v41tScfHxwUcttbW14gMAqszW/MqEXzgFAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLtjo958+bFyJEjo76+PnK5XNx7770fuXbcuHGRy+Vi2rRp2zEiALAjaXd8rF27NgYOHBjTp0/f4rp77703nnrqqaivr9/m4QCAHU+7v1hu+PDhMXz48C2uefXVV+N73/te/OlPf4oRI0Zs83AAwI6n6N9qu2HDhhg1alScd955cfDBB3/s+kKhEIVCofV+S0tLsUcCACpI0ePjqquuipqampgwYcJWrW9sbIyGhoZij0GZ9Zl8f7lHaLeXr3SVDiCFov61y6JFi+KXv/xl3HrrrZHL5bbqOVOmTInm5ubWW1NTUzFHAgAqTFHjY/78+bFq1aro3bt31NTURE1NTbzyyivxwx/+MPr06bPZ5+Tz+aitrW1zAwB2XEX92GXUqFFx7LHHtjl2wgknxKhRo2Ls2LHFPBUAUKXaHR9r1qyJF154ofX+8uXLY8mSJdG5c+fo3bt3dOnSpc36XXfdNbp37x79+vXb/mkBgKrX7vhYuHBhDB06tPX+pEmTIiJizJgxceuttxZtMABgx9Tu+BgyZEhkWbbV619++eX2ngIA2IH5bhcAICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBIqt3xMW/evBg5cmTU19dHLpeLe++9t/Wx9957Ly644II45JBDYs8994z6+voYPXp0vPbaa8WcGQCoYu2Oj7Vr18bAgQNj+vTpmzz29ttvx+LFi+MnP/lJLF68OO65555YtmxZnHzyyUUZFgCofjXtfcLw4cNj+PDhm32srq4u5syZ0+bYtddeG0cccUSsWLEievfuvW1TAgA7jJL/zkdzc3PkcrnYe++9S30qAKAKtPvKR3u88847MXny5DjzzDOjtrZ2s2sKhUIUCoXW+y0tLaUcCQAos5LFx3vvvRdf/epXY8OGDXHdddd95LrGxsZoaGgo1Riw1fpMvr/cI7Tby1eOKPcIAO1Wko9d3nvvvTj99NNj+fLlMWfOnI+86hERMWXKlGhubm69NTU1lWIkAKBCFP3Kxwfh8Y9//CP+8pe/RJcuXba4Pp/PRz6fL/YYAECFand8rFmzJl544YXW+8uXL48lS5ZE586do76+Pk499dRYvHhx/PGPf4z169fHypUrIyKic+fOsdtuuxVvcgCgKrU7PhYuXBhDhw5tvT9p0qSIiBgzZkxceumlcd9990VExGGHHdbmeX/5y19iyJAh2z4pALBDaHd8DBkyJLIs+8jHt/QYAIDvdgEAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkFS742PevHkxcuTIqK+vj1wuF/fee2+bx7Msi0svvTTq6+tj9913jyFDhsTSpUuLNS8AUOXaHR9r166NgQMHxvTp0zf7+NVXXx1Tp06N6dOnx4IFC6J79+5x3HHHxerVq7d7WACg+tW09wnDhw+P4cOHb/axLMti2rRpcdFFF8WXv/zliIi47bbbolu3bjFjxowYN27c9k0LAFS9ov7Ox/Lly2PlypVx/PHHtx7L5/NxzDHHxOOPP77Z5xQKhWhpaWlzAwB2XEWNj5UrV0ZERLdu3doc79atW+tjG2tsbIy6urrWW69evYo5EgBQYUry1y65XK7N/SzLNjn2gSlTpkRzc3PrrampqRQjAQAVot2/87El3bt3j4j3r4D06NGj9fiqVas2uRrygXw+H/l8vphjAAAVrKhXPvr27Rvdu3ePOXPmtB579913Y+7cuXHkkUcW81QAQJVq95WPNWvWxAsvvNB6f/ny5bFkyZLo3Llz9O7dO84999y44oor4oADDogDDjggrrjiithjjz3izDPPLOrgAEB1and8LFy4MIYOHdp6f9KkSRERMWbMmLj11lvj/PPPj//+979xzjnnxFtvvRWf+9zn4uGHH45OnToVb2oAoGrlsizLyj3Eh7W0tERdXV00NzdHbW1tucdhG/WZfH+5R9gpvHzliHKPABAR7Xv/9t0uAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBIqqbcA/Dx+ky+v9wjUKGq8Wfj5StHlHsEoMxc+QAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEkVPT7WrVsXP/7xj6Nv376x++67x3777ReXXXZZbNiwodinAgCqUE2xX/Cqq66KG264IW677bY4+OCDY+HChTF27Nioq6uLiRMnFvt0AECVKXp8PPHEE3HKKafEiBEjIiKiT58+MXPmzFi4cGGxTwUAVKGif+zyhS98IR555JFYtmxZRET8/e9/j8ceeyxOPPHEza4vFArR0tLS5gYA7LiKfuXjggsuiObm5ujfv3906NAh1q9fH5dffnl87Wtf2+z6xsbGaGhoKPYYAECFKvqVjzvvvDNuv/32mDFjRixevDhuu+22+PnPfx633XbbZtdPmTIlmpubW29NTU3FHgkAqCBFv/Jx3nnnxeTJk+OrX/1qREQccsgh8corr0RjY2OMGTNmk/X5fD7y+XyxxwAAKlTRr3y8/fbbscsubV+2Q4cO/tQWAIiIElz5GDlyZFx++eXRu3fvOPjgg+Ppp5+OqVOnxje/+c1inwoAqEJFj49rr702fvKTn8Q555wTq1ativr6+hg3blxcfPHFxT4VAFCFih4fnTp1imnTpsW0adOK/dIAwA7Ad7sAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACCpmnIPAOxc+ky+v9wjbJOXrxxR7hFgh+HKBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASKok8fHqq6/G17/+9ejSpUvssccecdhhh8WiRYtKcSoAoMrUFPsF33rrrTjqqKNi6NCh8eCDD0bXrl3jxRdfjL333rvYpwIAqlDR4+Oqq66KXr16xS233NJ6rE+fPsU+DQBQpYr+sct9990XgwYNitNOOy26du0an/nMZ+Kmm276yPWFQiFaWlra3ACAHVfRr3y89NJLcf3118ekSZPiwgsvjL/97W8xYcKEyOfzMXr06E3WNzY2RkNDQ7HH+Eh9Jt+f7FwAwKZyWZZlxXzB3XbbLQYNGhSPP/5467EJEybEggUL4oknnthkfaFQiEKh0Hq/paUlevXqFc3NzVFbW1vM0SJCfADb5uUrR5R7BKhoLS0tUVdXt1Xv30X/2KVHjx5x0EEHtTn26U9/OlasWLHZ9fl8Pmpra9vcAIAdV9Hj46ijjornn3++zbFly5bFvvvuW+xTAQBVqOjx8YMf/CCefPLJuOKKK+KFF16IGTNmxI033hjjx48v9qkAgCpU9Pj47Gc/G7Nnz46ZM2fGgAED4qc//WlMmzYtzjrrrGKfCgCoQkX/a5eIiJNOOilOOumkUrw0AFDlfLcLAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkSh4fjY2Nkcvl4txzzy31qQCAKlDS+FiwYEHceOONceihh5byNABAFSlZfKxZsybOOuusuOmmm2KfffYp1WkAgCpTsvgYP358jBgxIo499tgtrisUCtHS0tLmBgDsuGpK8aKzZs2KxYsXx4IFCz52bWNjYzQ0NJRiDACgAhX9ykdTU1NMnDgxbr/99ujYsePHrp8yZUo0Nze33pqamoo9EgBQQYp+5WPRokWxatWqOPzww1uPrV+/PubNmxfTp0+PQqEQHTp0aH0sn89HPp8v9hgAQIUqenwMGzYsnn322TbHxo4dG/37948LLrigTXgAADufosdHp06dYsCAAW2O7bnnntGlS5dNjgMAOx//wykAkFRJ/tplY48++miK0wAAVcCVDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJFVT7gEAqkGfyfeXe4R2e/nKEeUeATbLlQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJBU0eOjsbExPvvZz0anTp2ia9eu8aUvfSmef/75Yp8GAKhSRY+PuXPnxvjx4+PJJ5+MOXPmxLp16+L444+PtWvXFvtUAEAVqin2Cz700ENt7t9yyy3RtWvXWLRoUfzf//1fsU8HAFSZosfHxpqbmyMionPnzpt9vFAoRKFQaL3f0tJS6pEAgDIqaXxkWRaTJk2KL3zhCzFgwIDNrmlsbIyGhoZSjgGwU+oz+f5yj7BTePnKEeUeoeqU9K9dvve978UzzzwTM2fO/Mg1U6ZMiebm5tZbU1NTKUcCAMqsZFc+vv/978d9990X8+bNi549e37kunw+H/l8vlRjAAAVpujxkWVZfP/734/Zs2fHo48+Gn379i32KQCAKlb0+Bg/fnzMmDEjfv/730enTp1i5cqVERFRV1cXu+++e7FPBwBUmaL/zsf1118fzc3NMWTIkOjRo0fr7c477yz2qQCAKlSSj10AAD6K73YBAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBSNeUeAACqWZ/J95d7hHZ7+coRZT2/Kx8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACCpksXHddddF3379o2OHTvG4YcfHvPnzy/VqQCAKlKS+Ljzzjvj3HPPjYsuuiiefvrpOProo2P48OGxYsWKUpwOAKgiJYmPqVOnxre+9a04++yz49Of/nRMmzYtevXqFddff30pTgcAVJGaYr/gu+++G4sWLYrJkye3OX788cfH448/vsn6QqEQhUKh9X5zc3NERLS0tBR7tIiI2FB4uySvCwDVohTvsR+8ZpZlH7u26PHxxhtvxPr166Nbt25tjnfr1i1Wrly5yfrGxsZoaGjY5HivXr2KPRoAEBF100r32qtXr466urotril6fHwgl8u1uZ9l2SbHIiKmTJkSkyZNar2/YcOG+Pe//x1dunTZ7PqUWlpaolevXtHU1BS1tbVlnaVa2cPtY/+2nz3cPvZv++0se5hlWaxevTrq6+s/dm3R4+MTn/hEdOjQYZOrHKtWrdrkakhERD6fj3w+3+bY3nvvXeyxtkttbe0O/QOTgj3cPvZv+9nD7WP/tt/OsIcfd8XjA0X/hdPddtstDj/88JgzZ06b43PmzIkjjzyy2KcDAKpMST52mTRpUowaNSoGDRoUgwcPjhtvvDFWrFgR3/nOd0pxOgCgipQkPs4444x4880347LLLovXX389BgwYEA888EDsu+++pThdyeTz+bjkkks2+ViIrWcPt4/92372cPvYv+1nDzeVy7bmb2IAAIrEd7sAAEmJDwAgKfEBACQlPgCApHba+Lj++uvj0EMPbf1PXwYPHhwPPvjgFp9TKBTioosuin333Tfy+Xzsv//+8dvf/jbRxJVnW/bwjjvuiIEDB8Yee+wRPXr0iLFjx8abb76ZaOLK1tjYGLlcLs4999wtrps7d24cfvjh0bFjx9hvv/3ihhtuSDNghdua/bvnnnviuOOOi09+8pOtP7N/+tOf0g1Z4bb2Z/ADf/3rX6OmpiYOO+ywks5VLbZ2/7yX7MTx0bNnz7jyyitj4cKFsXDhwvjiF78Yp5xySixduvQjn3P66afHI488EjfffHM8//zzMXPmzOjfv3/CqStLe/fwsccei9GjR8e3vvWtWLp0adx1112xYMGCOPvssxNPXnkWLFgQN954Yxx66KFbXLd8+fI48cQT4+ijj46nn346LrzwwpgwYUL87ne/SzRpZdra/Zs3b14cd9xx8cADD8SiRYti6NChMXLkyHj66acTTVq5tnYPP9Dc3ByjR4+OYcOGlXiy6tCe/fNeEhEZrfbZZ5/sN7/5zWYfe/DBB7O6urrszTffTDxVddnSHv7sZz/L9ttvvzbHrrnmmqxnz54pRqtYq1evzg444IBszpw52THHHJNNnDjxI9eef/75Wf/+/dscGzduXPb5z3++xFNWrvbs3+YcdNBBWUNDQ2mGqxLbsodnnHFG9uMf/zi75JJLsoEDB5Z8xkrWnv3zXvK+nfbKx4etX78+Zs2aFWvXro3Bgwdvds19990XgwYNiquvvjo+9alPxYEHHhg/+tGP4r///W/iaSvT1uzhkUceGf/85z/jgQceiCzL4l//+lfcfffdMWLEiMTTVpbx48fHiBEj4thjj/3YtU888UQcf/zxbY6dcMIJsXDhwnjvvfdKNWJFa8/+bWzDhg2xevXq6Ny5cwkmqx7t3cNbbrklXnzxxbjkkktKPFl1aM/+eS95X8m+1bYaPPvsszF48OB45513Yq+99orZs2fHQQcdtNm1L730Ujz22GPRsWPHmD17drzxxhtxzjnnxL///e+d7rO6D2vPHh555JFxxx13xBlnnBHvvPNOrFu3Lk4++eS49tprE09dOWbNmhWLFy+OBQsWbNX6lStXbvIFjd26dYt169bFG2+8ET169CjFmBWrvfu3sV/84hexdu3aOP3004s8WfVo7x7+4x//iMmTJ8f8+fOjpmanfguJiPbvn/eS9+3UVz769esXS5YsiSeffDK++93vxpgxY+K5557b7NoNGzZELpeLO+64I4444og48cQTY+rUqXHrrbfudMX6Ye3Zw+eeey4mTJgQF198cSxatCgeeuihWL58+U77nT9NTU0xceLEuP3226Njx45b/bxcLtfmfva//6R44+M7um3dvw/MnDkzLr300rjzzjuja9euJZiw8rV3D9evXx9nnnlmNDQ0xIEHHphgwsq2LT+D3kv+p9yf+1SSYcOGZd/+9rc3+9jo0aOz/fffv82x5557LouIbNmyZSnGqwpb2sOvf/3r2amnntrm2Pz587OIyF577bUU41WU2bNnZxGRdejQofUWEVkul8s6dOiQrVu3bpPnHH300dmECRPaHLvnnnuympqa7N133001ekXYlv37wKxZs7Ldd989++Mf/5hw4srT3j186623Nlmfy+Vajz3yyCNl+peUx7b8DHoveZ9rZh+SZVkUCoXNPnbUUUfFXXfdFWvWrIm99torIiKWLVsWu+yyS/Ts2TPlmBVtS3v49ttvb3KZtkOHDq3P29kMGzYsnn322TbHxo4dG/37948LLrigdW8+bPDgwfGHP/yhzbGHH344Bg0aFLvuumtJ560027J/Ee9f8fjmN78ZM2fO3Ol/36i9e1hbW7vJ+uuuuy7+/Oc/x9133x19+/Yt+cyVZFt+Br2X/E+Z46dspkyZks2bNy9bvnx59swzz2QXXnhhtssuu2QPP/xwlmVZNnny5GzUqFGt61evXp317NkzO/XUU7OlS5dmc+fOzQ444IDs7LPPLtc/oezau4e33HJLVlNTk1133XXZiy++mD322GPZoEGDsiOOOKJc/4SKs/Fvym+8hy+99FK2xx57ZD/4wQ+y5557Lrv55puzXXfdNbv77rvLMG3l+bj9mzFjRlZTU5P96le/yl5//fXW23/+858yTFuZPm4PN+avXdr6uP3zXvK+nfZ3Pv71r3/FqFGjol+/fjFs2LB46qmn4qGHHorjjjsuIiJef/31WLFiRev6vfbaK+bMmRP/+c9/YtCgQXHWWWfFyJEj45prrinXP6Hs2ruH3/jGN2Lq1Kkxffr0GDBgQJx22mnRr1+/uOeee8r1T6h4G+9h375944EHHohHH300DjvssPjpT38a11xzTXzlK18p45SVa+P9+/Wvfx3r1q2L8ePHR48ePVpvEydOLOOUlW3jPaR9vJdsXi7LdsLr3QBA2ey0Vz4AgPIQHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEn9P3BFymEv1oWFAAAAAElFTkSuQmCC", "text/plain": [ "