{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to predict catastrophic events ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the end of that short but preciser tutorial, you'll understand what is a GEV (Generalized extreme value distribution) and how to fit it on your data in python. \n", "\n", "What we need to think about before speaking of extreme values statistics, is that in general statistics focus on average values and observable phenomenoms, for which we **do have datas**. Extreme values statistics is different from \"average\" statistics in this regard, indeed as we will see, we here try to model something for which **we have no data!** You will see soon what we mean by that. Before getting there I want to show you why extreme values modelling is an actual problem, of high interest for insurance companies and public services. Let's start with two examples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###
Montpellier, France, 2014\n", "**3 hours of rain = 50% of total annual average.** *Montpellier is not known to be a city with a high annual rain level, located in the dry location of Provence, therefore if we look at the average value we have a very common one. But it is not excluded that sometimes, rarely, very high values appear, as we can see:*\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###
Saltina river in Brigue, Switzerland, 1993\n", "**Unprecedented measured flows and precipitations in the twentieth century**.\n", "*Same problem here, this rivers has reasonable flows in average, but sometimes the flow can be exceptional.*\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###
Number of natural catastrophes since 1980 is growing \n", "**Those \"extreme\" events tends to become more and numerous those days, and make extreme value analysis \"popular\".**\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###
Some could wonder what statistical model shall we use to predict a Nadal loss at Roland Garros \n", "\n", "\n", "
\n", "
(Almost) just kidding" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question: What if we would like to quantify the probability of such events ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we said the particularity of the extreme value analysis is not to model an average phenomenom but an extreme phenomenom." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to understand what we mean by that, let's work with the Nidd river dataset here : https://www.stat.auckland.ac.nz/~wild/data/Rdatasets/csv/evir/nidd.thresh.csv, it contains $n=154$ trimestrial measures of the river's flows (in $m³/s$), over a period of 38.5 years (source: Natural Environment Research, 1975). \n", "\n", "Let's plot a histogram of the data.\n" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAH59JREFUeJzt3Xm8HHWZ7/HPlyQQIEBYIwQ0AeKCRhnIIIviCaCIiAQvzEURgYsT2UFxrlEZFhFFr6A4OiIjXEAZAkQEFBe4kQPCIJIIGJLAECEQMLJJyIKAIc/9o35NiqZPd52lunNOfd+vV7+69np+VdX9VP1qU0RgZmbVtVanAzAzs85yIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4J4JekjRXUlen4+gkSQdJWiRpuaR/6MP4b0zjDuuh/5mSftxk/IWS9unF/C6U9K+9jbNMko6V9GRaDptKCknbdzouqyYngpxGfzCSjpR0e609It4eEd0tpjMu/bCHlxRqp30TOCEiRkXEPfU9U9nnSFor1+0rki4FiIjH0rivtCPYiDgmIs4uY9qSunu7YyBpBHA+8IG0HJ4tI7aqkPQpSQtSUv2VpK1y/X6Zutc+L0uaM0Dz/aKkKyXd0eB/4yZJHxiI+bSDE8EgtAYkmDcBc1sMsxVwaBtiaaqno44BmK7yia6XxgAjab0MrQVJ7wO+ChwIbAI8AlxZ6x8R+6VkOyoiRgH/BVwzEPOOiK9GxMeAC4C9cjGtD+wM3DoQ82kHJ4Jeyh81SNpF0ixJS9Nh/vlpsNvS95K0F7KbpLUknSbpUUlPSbpc0ka56X4y9XtW0r/WzedMSTMk/VjSUuDINO87JS2RtFjSdyWtnZteSDpO0kOSlkk6W9J2aZylkq7OD19XxoaxSlpH0nJgGHCfpD81WVTfAM5qlLTqj5gkjZd0a4rzZmCzuuEPzy2bL7VYP5dK+r6kX0haAUxO3b6S+s+X9OHc8MMlPSNpp9S+q6T/Ssv1vvzeftr7P0fSHcALwLZ18+5pe8gP82bgwdS6RNJvGgyzUVrmT6dyn1ZLOql959T8ibQcd0jtn5J0XdFYelh+Z0q6Jm1ry5Qd2b1Z0hfStrBIuT3dFOvFaRt8QtmR37DUbztJv0nr7RlJV0ganRt3oaTPSfqjpOclXSVpZJE4cw4AromIuRHxMnA2sKek7RqUbRzwXuBHPZS9tl0elcr5nKRjJP1jinGJpO/WjXMJcA5wfa7z3sAdEfFSX9dD20WEP+kDLAT2qet2JHB7o2GAO4HDU/MoYNfUPA4IYHhuvP8FLCD78xgFXAv8KPXbAVgOvAdYm6zq5e+5+ZyZ2qeQJe91yfY4dgWGp/nNB07JzS+AG4ANgbcDLwEz0/w3AuYBR/SwHHqMNTft7ZssxwAmALOBT6VuXwEubbR80nI8H1gH2BNYBvy4btnsmfqfD6ysX0+5eV8KPA/skZbVyNTtK6n/6cAVueH3Bx5IzWOBZ4EPpXHfn9o3T/27gcfS8hwOjKibd8PtoUGMjbaPV5cpcDnZH8sGadj/Bo7O9Ts1NV8E/Ak4NtfvM72JpUFsZwIvAvumMl5Otpf9JWAE8M/AI7nhrwN+AKwPbAH8Hvh06rd9WobrAJuT7SB9u+639Huyo8dNyLbhY1K/NwJLmnw+noY7D/j33DTHpmV5YIOynQ50Nyl7bb1cmLabD6RlcV0q21jgKeB9dePtBvws135hbhn0aT20+9PxANakT9owl9dtcC/QcyK4DTgL2KyHDSr/Q58JHJdrfwvZn/vwtIFemeu3HvAyr00Et7WI/RTgp7n2APbItc8GPp9rPy//o6ybVo+x5qbdKhFsT/aH+lj6I2iYCNIPfiWwfm78/2R1IjgdmJ7rt35+2TSY96XA5Q261RLB9mSJZr3UfgVwemr+PLmEl7r9mpQwyRLBl5uUu+H20GC4RttHbZkNI0vaO+T6fZr0BwYcDdyQmucDn6otH+BRYKfexNIgtjOBm3PtB5D9Joal9g1SrKPJqrheAtbNDf8x4JYepj0FuKfut/SJXPs3gAt7Ge/ewDPAO8l2kH4ArAI+1mDYBcCRBdbL2Fy3Z4H/mWv/CWmHC/guWTXTTOCDuWEeBbbpz3po98dVQ683JSJG1z7AcU2GPRp4M/CApLvzVQ4NbEW2gdQ8SvZHOCb1W1TrEREvkG2AeYvyLelw/eeS/qKsuuir1FWpAE/mmv/WoH1UH2ItLCJ+QZYIpjYZbCvguYhYUTe/fP/8slnB65dNvUU99YiIBWR/oAdIWg/4CFnigezcxyGpCmCJpCVkR2lbFpk2vdseerIZ2VFh/fIfm5pvBd4r6Q1kSeMqYI9U7bERcO8AxFK/nTwTq0/s/y19jyJbXiOAxbnl9QOyvWckbSFpeqoyWgr8mNdvo3/JNb9Az9tkQxExEziD7A/6UbLksgx4PD+cpPcAbwBmFJhsod9NRJwQEYdExN4R8as0n4nA0oiobScDsU2UzomgHyLiochOFm0BfB2YoexEUaNHuv6Z7IdTU9sTfhJYDGxd6yFpXWDT+tnVtX8feACYEBEbAl8E1PfSFI61t04jq1ZYr4f+i4GN03LLzy/ff5taS/rzrl829Vo9UvdKsj3XA4F5KTlA9if/o/yOQESsHxHnFpl2k+2hN54hO/qqX/5PpHksIPvDPInsKHEZ2Z/pVLIj11UDGEsri8iOCDbLLa8NI+Ltqf/XyJbXO9M2+gkKbqNafYlxT5/DasNGxPciYkJEbEGWEIYD99dN8gjg2ohY3r8it/Qh4MZcbO1YD/3mRNAP6WTd5unHtyR1fgV4muzwNH8y8UrgM8pOjI4i24O/KiJWku2lHCBpd2UncM+i9Q9mA2ApsFzSW4FjB6xgzWPtlcgutZ1D9kNs1P9RYBbZieW1057bAblBZgAflvSetGy+TP+32+lk9b/HsvpoALI91gMk7StpmKSRkrokbd1wKnWabA+FpT3vq4FzJG0g6U3AZ1NsNbcCJ7D6qpTuuvamsaSTtEf2Jq4eYl0M3AScJ2lDZRcZbKfsSh7IttHlZCfFxwL/0otp1y4x7ulzRSrLSEnvUOaNZOdNLoiI52rTSjtWh5BVEZZtf+AXuXn3e5toByeC/vkgMFfZlTQXAIdGxIupaucc4I50yLwrcAnZ1Qq3kZ18exE4ESAi5qbm6WR7wMvITkq91GTenwM+nob9D7IqgoHSY6x9dBrZycCefBx4N/BXssP8y2s90rI5nuwPezHwHHWH/b2V/sDuBHYnt9zS4fyBZEdXT5Pt8f4LxX8nDbeHPoR4IrACeBi4nazsl+T630r2J3tbD+09xpKS6abA7/oQVyOfJKvKmke2bmawuirtLGAnspP3N5JddDDQRpItn+VkJ57vBOpvHpySYrilhPm/StlVgG8ju0S1ZqC2iVIpndCwNUjaC19CVu3zSKfjsaEjHXEdn6orbABJ+ifg4Ij4p07H0ltOBGsISQeQXX0gsit63k12BYhXkNkgoOz+imURcWenY+mtTt+haqsdSFYdI7I680OdBMwGj4i4qdMx9JWPCMzMKs4ni83MKm5QVA1tttlmMW7cuE6H0RYrVqxg/fXXuMuM26LKZYdql99lL6fss2fPfiYiNm813KBIBOPGjWPWrFmdDqMturu76erq6nQYHVHlskO1y++yd5UybUmPth7KVUNmZpXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFTco7izuj3HTbmw9UAkWnrt/R+ZrZtZbPiIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOziis1EUj6jKS5ku6XdKWkkZLGS7pL0kOSrpK0dpkxmJlZc6UlAkljgZOASRHxDmAYcCjwdeBbETEBeA44uqwYzMystbKrhoYD60oaDqwHLAb2Amak/pcBU0qOwczMmlBElDdx6WTgHOBvwE3AycDvImL71H8b4JfpiKF+3KnAVIAxY8bsPH369D7FMOeJ5/sWfD9NHLtRn8Zbvnw5o0aNGuBoBocqlx2qXX6XvZyyT548eXZETGo13PBS5g5I2hg4EBgPLAGuAfZrMGjDTBQRFwEXAUyaNCm6urr6FMeR027s03j9tfCwrj6N193dTV/LOthVuexQ7fK77F0djaHMqqF9gEci4umI+DtwLbA7MDpVFQFsDfy5xBjMzKyFMhPBY8CuktaTJGBvYB5wC3BwGuYI4PoSYzAzsxZKSwQRcRfZSeE/AHPSvC4CPg98VtICYFPg4rJiMDOz1ko7RwAQEWcAZ9R1fhjYpcz5mplZcb6z2Mys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKziWiYCSetLWis1v1nSRySNKD80MzNrhyJHBLcBI9Mbx2YCRwGXlhmUmZm1T5FEoIh4Afgo8G8RcRCwQ7lhmZlZuxRKBJJ2Aw4Dam95KfVhdWZm1j5FEsEpwBeAn0bEXEnbkr1TwMzMhoCWe/YRcStwq6T1U/vDwEllB2ZmZu1R5Kqh3STNA+an9ndJ+vfSIzMzs7YoUjX0bWBf4FmAiLgP2LPMoMzMrH0K3VAWEYvqOr1SQixmZtYBRa7+WSRpdyAkrU12fmB+uWGZmVm7FDkiOAY4HhgLPA7smNrNzGwIaHpEIGkYcHhEHNameMzMrM2aHhFExCvAgW2KxczMOqDIOYI7JH0XuApYUesYEX8oLSozM2ubIolg9/T95Vy3APYa+HDMzKzditxZPLkdgZiZWWe0TASSTm/UPSK+3Ki7mZkNLkWqhlbkmkcCH8b3EZiZDRlFqobOy7dL+iZwQ2kRmZlZW/XlncXrAdsOdCBmZtYZRc4RzCG7SghgGLA5r72CyMzMBrEi5wg+nGteCTwZEStLisfMzNqsSNXQcOAvEfEoMAE4TtLocsMyM7N2KZIIfgK8Iml74GJgPPCfpUZlZmZtUyQRrEpVQR8Fvh0RnwG2LDcsMzNrlyKJ4O+SPgZ8Evh56jaivJDMzKydiiSCo4DdgHMi4hFJ44EfF5m4pNGSZkh6QNL89P7jTSTdLOmh9L1xfwpgZmb90zIRRMS8iDgpIq5M7Y9ExLkFp38B8KuIeCvwLrI7kqcBMyNiAjAztZuZWYe0TASSJqS9+nmSHq59Coy3IdlL7i8GiIiXI2IJ2fsNLkuDXQZM6Xv4ZmbWX4qI5gNItwNnAN8CDiCrKlJEnNFivB2Bi4B5ZEcDs4GTgSciYnRuuOci4nXVQ5KmAlMBxowZs/P06dN7UazV5jzxfJ/G66+JYzfq03jLly9n1KhRAxzN4FDlskO1y++yl1P2yZMnz46ISa2GK5IIZkfEzpLmRMTE1O23EfHeFuNNAn4H7BERd0m6AFgKnFgkEeRNmjQpZs2a1aosDY2bdmOfxuuvhefu36fxuru76erqGthgBokqlx2qXX6XvauUaaf/75aJoMjJ4hclrQU8JOkESQcBWxQY73Hg8Yi4K7XPAHYCnpS0ZQpyS+CpAtMyM7OSFEkEp5A9aO4kYGfgE8ARrUaKiL8AiyS9JXXam6ya6Ibc+EcA1/cyZjMzG0BFHkN9N4CkiIijejn9E4ErJK0NPEx2fmEt4GpJRwOPAYf0cppmZjaAijx9dDeyK39GAW+U9C7g0xFxXKtxI+JeoFH91N69DdTMzMpRpGro28C+wLMAEXEf2WWhZmY2BBR6MU1ELKrr9EoJsZiZWQcUeR/BIkm7A5Hq+k/C7yw2MxsyihwRHAMcD4wluyR0x9RuZmZDQJGrhp4BDmtDLGZm1gFFrhoaT3YZ6Lj88BHxkfLCMjOzdilyjuA6sstHfwasKjccMzNrtyKJ4MWI+E7pkZiZWUcUSQQXSDoDuAl4qdYxIv5QWlRmZtY2RRLBROBwYC9WVw1Fajczs0GuSCI4CNg2Il4uOxgzM2u/IvcR3AeMbjmUmZkNSkWOCMYAD0i6m9eeI/Dlo2ZmQ0CRRND0lZRmZja4Fbmz+NZ2BGJmZp1R6OmjZmY2dDkRmJlVXI+JQNLM9P319oVjZmbt1uwcwZaS3gd8RNJ0QPmevrPYzGxoaJYITgemAVsD59f1853FZmZDRI+JICJmADMk/WtEnN3GmMzMrI2KXD56tqSPsPqF9d0R8fNywzIzs3ZpedWQpK8BJwPz0ufk1M3MzIaAIncW7w/sGBGrACRdBtwDfKHMwMzMrD2K3keQf+jcRmUEYmZmnVHkiOBrwD2SbiG7hHRPfDRgZjZkFDlZfKWkbuAfyRLB5yPiL2UHZmZm7VHkiICIWAzcUHIsZmbWAX7WkJlZxTkRmJlVXNNEIGktSfe3KxgzM2u/pokg3Ttwn6Q3tikeMzNrsyIni7cE5kr6PbCi1tHvLDYzGxqKJIKzSo/CzMw6puXJ4vTO4oXAiNR8N1D4XQSShkm6R9LPU/t4SXdJekjSVZLW7mPsZmY2AIo8dO6fgRnAD1KnscB1vZjHycD8XPvXgW9FxATgOeDoXkzLzMwGWJHLR48H9gCWAkTEQ8AWRSYuaWuyh9b9MLWL7IU2M9IglwFTeheymZkNpCKJ4KWIeLnWImk42RvKivg28L+BVal9U2BJRKxM7Y+THWGYmVmHFDlZfKukLwLrSno/cBzws1YjSfow8FREzJbUVevcYNCGSUXSVGAqwJgxY+ju7i4Q6uudOnFl64FK0Nd4ly9f3udxB7sqlx2qXX6XvbujMSii+c69pLXI6vE/QPZH/mvgh9FixPTymsOBlcBIYEPgp8C+wBsiYqWk3YAzI2LfZtOaNGlSzJo1q1iJ6oybdmOfxuuvhefu36fxuru76erqGthgBokqlx2qXX6XvauUaUuaHRGTWg1X5Omjq9LLaO4i23t/sFUSSON9gfS46nRE8LmIOEzSNcDBwHTgCOD6VtMyM7PyFLlqaH/gT8B3gO8CCyTt1495fh74rKQFZOcMLu7HtMzMrJ+KnCM4D5gcEQsAJG0H3Aj8suhMIqIb6E7NDwO79DZQMzMrR5Grhp6qJYHkYeCpkuIxM7M26/GIQNJHU+NcSb8AriY7R3AI2d3FZmY2BDSrGjog1/wk8L7U/DSwcWkRmZlZW/WYCCLiqHYGYmZmndHyZLGk8cCJwLj88H4MtZnZ0FDkqqHryC7x/BmrHxVhZmZDRJFE8GJEfKf0SMzMrCOKJIILJJ0B3AS8VOsYEYXfSWBmZmuuIolgItkzg/ZiddVQpHYzMxvkiiSCg4Bt84+iNjOzoaPIncX3AaPLDsTMzDqjyBHBGOABSXfz2nMEvnzUzGwIKJIIzig9iiGor+9BOHXiSo7s5zsU+vouBDOrpiLvI7i1HYGYmVlnFLmzeBmrXye5NjACWBERG5YZmJmZtUeRI4IN8u2SpuD3CZiZDRlFrhp6jYi4Dt9DYGY2ZBSpGvpornUtYBKrq4rMzGyQK3LVUP69BCuBhcCBpURjZmZtV+Qcgd9LYGY2hDV7VeXpTcaLiDi7hHjMzKzNmh0RrGjQbX3gaGBTwInAzGwIaPaqyvNqzZI2AE4GjgKmA+f1NJ6ZmQ0uTc8RSNoE+CxwGHAZsFNEPNeOwMzMrD2anSP4P8BHgYuAiRGxvG1RmZlZ2zS7oexUYCvgNODPkpamzzJJS9sTnpmZla3ZOYJe33VsZmaDj//szcwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKziSksEkraRdIuk+ZLmSjo5dd9E0s2SHkrfG5cVg5mZtVbmEcFK4NSIeBuwK3C8pB2AacDMiJgAzEztZmbWIaUlgohYHBF/SM3LgPnAWLLXXF6WBrsMmFJWDGZm1poiyn8PvaRxwG3AO4DHImJ0rt9zEfG66iFJU4GpAGPGjNl5+vTpfZr3nCee79N4nTJmXXjyb/2bxsSxGw1MMG22fPlyRo0a1ekwOqbK5XfZyyn75MmTZ0fEpFbDlZ4IJI0CbgXOiYhrJS0pkgjyJk2aFLNmzerT/MdNu7FP43XKqRNXct6clq+SbmrhufsPUDTt1d3dTVdXV6fD6Jgql99l7ypl2pIKJYJSrxqSNAL4CXBFRFybOj8pacvUf0vgqTJjMDOz5sq8akjAxcD8iDg/1+sG4IjUfARwfVkxmJlZa/2rg2huD+BwYI6ke1O3LwLnAldLOhp4DDikxBjMzKyF0hJBRNwOqIfee5c1XzMz6x3fWWxmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFVfay+utc8ZNu7HTIfTJqRNXcmQfY1947v4DHI1ZdfiIwMys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzjeU2ZDQqZvofCObDQU+IjAzqzgnAjOzinMiMDOrOJ8jMBukqnhepIplbgcfEZiZVVxHEoGkD0p6UNICSdM6EYOZmWXaXjUkaRjwPeD9wOPA3ZJuiIh57Y7FrL8GsqqiP+9jaKcyqmfW9LKXWSXVrOztqpLqxBHBLsCCiHg4Il4GpgMHdiAOMzMDFBHtnaF0MPDBiPhUaj8ceHdEnFA33FRgamp9C/BgWwPtnM2AZzodRIdUuexQ7fK77OV4U0Rs3mqgTlw1pAbdXpeNIuIi4KLyw1mzSJoVEZM6HUcnVLnsUO3yu+ydLXsnqoYeB7bJtW8N/LkDcZiZGZ1JBHcDEySNl7Q2cChwQwfiMDMzOlA1FBErJZ0A/BoYBlwSEXPbHccarHLVYTlVLjtUu/wuewe1/WSxmZmtWXxnsZlZxTkRmJlVnBNBm0m6RNJTku7PddtE0s2SHkrfG6fukvSd9CiOP0raqXOR918PZT9T0hOS7k2fD+X6fSGV/UFJ+3Ym6oEhaRtJt0iaL2mupJNT9yG/7puUfcive0kjJf1e0n2p7Gel7uMl3ZXW+1XpwhkkrZPaF6T+49oSaET408YPsCewE3B/rts3gGmpeRrw9dT8IeCXZPde7Arc1en4Syj7mcDnGgy7A3AfsA4wHvgTMKzTZehH2bcEdkrNGwD/nco45Nd9k7IP+XWf1t+o1DwCuCutz6uBQ1P3C4FjU/NxwIWp+VDgqnbE6SOCNouI24C/1nU+ELgsNV8GTMl1vzwyvwNGS9qyPZEOvB7K3pMDgekR8VJEPAIsIHs8yaAUEYsj4g+peRkwHxhLBdZ9k7L3ZMis+7T+lqfWEekTwF7AjNS9fr3XtocZwN6SGt2EO6CcCNYMYyJiMWQ/GmCL1H0ssCg33OM0/wENViek6o9LalUjDOGyp8P9fyDbO6zUuq8rO1Rg3UsaJule4CngZrIjnCURsTINki/fq2VP/Z8HNi07RieCNVuhx3EMct8HtgN2BBYD56XuQ7LskkYBPwFOiYilzQZt0G1Ql79B2Sux7iPilYjYkewpCrsAb2s0WPruSNmdCNYMT9YO+9P3U6n7kH8cR0Q8mX4oq4D/YHUVwJAru6QRZH+EV0TEtalzJdZ9o7JXad0DRMQSoJvsHMFoSbUbevPle7Xsqf9GFK9O7TMngjXDDcARqfkI4Ppc90+mK0h2BZ6vVSMMFXX13gcBtSuKbgAOTVdRjAcmAL9vd3wDJdXzXgzMj4jzc72G/LrvqexVWPeSNpc0OjWvC+xDdo7kFuDgNFj9eq9tDwcDv4l05rhUnT6rXrUPcCXZYfDfybL/0WR1gDOBh9L3JrH6ioPvkdUpzgEmdTr+Esr+o1S2P5L9CLbMDf+lVPYHgf06HX8/y/4eskP8PwL3ps+HqrDum5R9yK974J3APamM9wOnp+7bkiW3BcA1wDqp+8jUviD137YdcfoRE2ZmFeeqITOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIrAhQ9IruSdZ3itpnKQuST8fwHlsJWlG6yF7Pd190tNGr8tfXy/pV5IG7eMVbHBo+6sqzUr0t8hu5X/VQD/GNyL+zOobgQZyuv9PUhewFbASXr0BaZOIeGKg52eW5yMCq4z07P/r0kPOfifpnan7HEmj0128z0r6ZOr+I0n71E1jnNL7FCQdKenatNf+kKRv9DDfhZK+KulOSbMk7STp15L+JOmY2nARcRpwBbB96tRF9kgCJJ0raV6K/ZsDu2Ss6nxEYEPJuukpjwCPRMRBdf3PAu6JiCmS9gIuJ3vg2R3AHsCjwMPAe1O/XYFjW8xzR7Knab4EPCjp3yJiUYPhFkXEbpK+BVya5jcSmAtcKOnTwLuAN5PdZQuwH3CdpE3IHsHw1oiI2iMLzAaKE4ENJa+rGqrzHuB/AETEbyRtKmkj4LdkL815lOyJmFNTvfxfY/Wz5HsyMyKeB5A0D3gTr32Eck3tz30O2YtKlgHLJL0oaXRE/KDBOHsAnwNWAS8CP5R0IzBg5zzMwFVDVi09PeL3NrKjgPeSVcU8TXYe4LcFpvlSrvkVet65qg23qm6cVY3GkbQt2VHEy5E9l34Xsqd3TgF+VSAus8KcCKxKbgMOA0gnZp+JiKWpKmczYEJEPAzcTrYnXiQRlGU/0h9+eo7/RhHxC+AUsuooswHjqiGrkjOB/yvpj8ALrH7cL2RvzBqWmn8LfI0sIXTKB4ETU/MGwPWSRpId1XymY1HZkOSnj5qtYSStA9wREZM6HYtVgxOBmVnF+RyBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlV3P8HFmCh6woIA3kAAAAASUVORK5CYII=\n", "text/plain": [ "