{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "Lab 2 - Revenue Management", "version": "0.3.2", "provenance": [], "collapsed_sections": [], "include_colab_link": true }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "metadata": { "id": "RoDTAKe3eWDE", "colab_type": "text" }, "cell_type": "markdown", "source": [ "# Lab 2 - Revenue Management\n", "\n", "In this lab, we will study implementations of the algorithms we have discussed in the Revenue Management lectures. Since we will be using a lot of numerical methods, we would need to program a large set of statistical distributions, mathematical functions, and other mathematical quantities. Luckily for us, this has already been done.\n", "\n", "A **package** is a set of functions, constants, and other data that comes pre-packaged and can be used once installed locally. Many packages are available (exactly [168,180](https://pypi.python.org/pypi)), but we will focus only on a set of them of course.\n", "\n", "The fist one is [Numpy](http://www.numpy.org/), self-described as \"the fundamental package for scientific computing with Python\". It comes with a very large number of scientific functions. These range from simply implementing mathetical constants (such as $\\pi$ or $e$), to mathematical functions (such as the logistic functions), random number generators, and much, much, more. As we move along with the activities of the module, we will use many of Numpy's functions, but the packages that we will use will most certainly be using Numpy under the hood.\n", "\n", "## Loading packages\n", "\n", "To load a package so it can be used in your terminal, write the following line:\n", "\n", "```\n", "import PACKAGE as SHORT_NAME\n", "```\n", "\n", "For example, to load Numpy and assign it the (well-known) alias \"np\" we run" ] }, { "metadata": { "id": "FxPqqILaeWva", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "import numpy as np" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "3WV0FYJReqOH", "colab_type": "text" }, "cell_type": "markdown", "source": [ "With this, all of numpy is available by calling ```np.NAME```. For example, to check the value of $\\pi$." ] }, { "metadata": { "id": "DhBKYRSAencB", "colab_type": "code", "outputId": "66df7113-d46e-405c-eab2-6f45cc685f73", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "\n", "np.pi" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "3.141592653589793" ] }, "metadata": { "tags": [] }, "execution_count": 2 } ] }, { "metadata": { "id": "geCgr2iGezDN", "colab_type": "text" }, "cell_type": "markdown", "source": [ "... or to calculate the value of $e^4$." ] }, { "metadata": { "id": "FEFFkAOveyVG", "colab_type": "code", "outputId": "a25a120e-af7f-4247-f145-f67c3a621cc1", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "np.e ** 4" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "54.59815003314423" ] }, "metadata": { "tags": [] }, "execution_count": 3 } ] }, { "metadata": { "id": "mEQ9jjbve4dV", "colab_type": "text" }, "cell_type": "markdown", "source": [ "Note that the power function is two asterisks, so for example $4^2$ would require this line of code:" ] }, { "metadata": { "id": "2Dhr4EPYe15k", "colab_type": "code", "outputId": "a69e3e62-f529-4425-f769-3a1aa9211b69", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "4 ** 2" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "16" ] }, "metadata": { "tags": [] }, "execution_count": 4 } ] }, { "metadata": { "id": "oO40Wo37fBEi", "colab_type": "text" }, "cell_type": "markdown", "source": [ "Numpy is a **very** extensive package. You can explore the details of it [here](https://docs.scipy.org/doc/numpy-dev/user/quickstart.html). We will now apply it in order to calculate the functions we need.\n", "\n", "We first start by studying elasticity.\n", "\n", "## Constant elasticity functions\n", "\n", "This means the demand function is $𝑑(𝑝)=𝐢𝑝^{βˆ’πœ–}$. What is the shape of this function?\n", "\n", "### Calculating the demand function.\n", "\n", "We will first create a function using Python. As you studied last week (right????), a function is defined following this convention:\n", "\n", "```\n", "def NAME(PARAMS):\n", " FUNCTION CODE\n", " MORE CODE\n", "```\n", "\n", "So, let's define the function \"d_const_elast\" which will receive the price, the elasticity, and an optional constant with the number of cases that we have." ] }, { "metadata": { "id": "2AFXpYXVfAUb", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "def d_const_elast(p, elast, C = 100):\n", " out = C * p ** (-1 * elast)\n", " return(out)" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "aligxYJEfU52", "colab_type": "text" }, "cell_type": "markdown", "source": [ "Now our function is ready to use! See how we give a value to ```C```? This is an **optional variable**. By saying ```C = 100```, we are allowing the function to receive two or three parameters.\n", "\n", "1. Calling ```d_const_elast(10, 1)``` will calculate the value of the demand using ```p = 10, elast = 1, C = 100```.\n", "2. Calling ```d_const_elast(10, 1, 400)``` will calculate the value of the demand using ```C = 400```.\n", "3. Calling ```d_const_elast(10)``` will fail.\n", "4. Calling ```d_const_elast(elast = 1, p = 10, C = 400)``` (or any other order) also works as long as all the inputs are named!" ] }, { "metadata": { "id": "rVAbSsKWfO6H", "colab_type": "code", "outputId": "1075a704-c0ea-4374-92cb-bd75469fb1be", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "d_const_elast(10, 1, 400)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "40.0" ] }, "metadata": { "tags": [] }, "execution_count": 7 } ] }, { "metadata": { "id": "xb48cR4ffkgy", "colab_type": "code", "outputId": "f7491799-e8da-420c-b230-676b95d31ef0", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "d_const_elast(elast = 1, p = 10, C = 400)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "40.0" ] }, "metadata": { "tags": [] }, "execution_count": 8 } ] }, { "metadata": { "id": "YGz38n4aflwD", "colab_type": "code", "outputId": "e985c463-a470-4181-c8a3-e0b77aa6f9d5", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "d_const_elast(C = 400, elast = 1, p = 10)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "40.0" ] }, "metadata": { "tags": [] }, "execution_count": 9 } ] }, { "metadata": { "id": "oyJJB7Aefm3q", "colab_type": "code", "outputId": "e2fc100d-af01-4438-912e-f3b299a26a9d", "colab": { "base_uri": "https://localhost:8080/", "height": 172 } }, "cell_type": "code", "source": [ "d_const_elast(C = 400, elast = 1)" ], "execution_count": 0, "outputs": [ { "output_type": "error", "ename": "TypeError", "evalue": "ignored", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0md_const_elast\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m400\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0melast\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mTypeError\u001b[0m: d_const_elast() missing 1 required positional argument: 'p'" ] } ] }, { "metadata": { "id": "tIw9LnGDfpRS", "colab_type": "text" }, "cell_type": "markdown", "source": [ "The last commands fails because the argument p is **not** an optional argument.\n", "\n", "### Plotting the function\n", "\n", "We can now plot the function. We will do so using the extremely powerful package ```matplotlib```, in particular the python implementation ```pyplot```. We will use the common alias ```plt```." ] }, { "metadata": { "id": "VlK7Afshfofi", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "WeSkUb3Of24_", "colab_type": "text" }, "cell_type": "markdown", "source": [ "The line ```%matplotlib inline``` is called a \"magic\", and is a command that tells Jupyter to do something specific. This commands tells Jupyter to plot inline (i.e. in this notebook directly), instead of saving the image elsewhere. A list of all magic commands is available [here](http://ipython.readthedocs.io/en/stable/interactive/magics.html).\n", "\n", "We can now plot a function with a constant elasticity over the range [1, 100] with the following functions:" ] }, { "metadata": { "id": "Mm4KTO16f0w_", "colab_type": "code", "outputId": "d967bf73-93f3-48fe-ba13-5359eda83f65", "colab": { "base_uri": "https://localhost:8080/", "height": 265 } }, "cell_type": "code", "source": [ "p = np.arange(1, 100.1, 0.1) # Calculates evenly spaced points in the interval 1 to 100.\n", "\n", "plt.plot(p, d_const_elast(p, 1)) # Applies the vector p to the d_const_elast function using elast = 1\n", "plt.show()" ], "execution_count": 0, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGq9JREFUeJzt3XmQnPV95/H3c/Q1p2ZGrQNZQuL6\ncYMPltsWjoi5UqwtvF4XIXgh5azLSfmobDZVSWxMkvJWUi4765AsVOK1TZJypVxxgmIHKFjHJoBB\nJjgcln8gQEL3jKTRaO7pa/94np7pGXpGmp4etZ7n+byqpvrp53n6eX5fRnz6N7/nciqVCiIiEg9u\nqxsgIiLNo1AXEYkRhbqISIwo1EVEYkShLiISI34rdz4wMLyoU296etoYHBxbruactpJYdxJrhmTW\nncSaYWl15/OdznzLItVT932v1U1oiSTWncSaIZl1J7FmWL66IxXqIiKyMIW6iEiMKNRFRGJEoS4i\nEiMKdRGRGDmpUxqNMRcD/wR81Vr758aY9cDDgAccAO6y1k4aY+4EPguUgYestX+9TO0WEZE6TthT\nN8a0A18HnqyZfT/wgLX2emAncE+43heALcBm4HPGmN6mt1hEROZ1MsMvk8AtwP6aeZuBR8LpbQRB\nfiWw3Vo7ZK0dB54Grm1eU2fs7R/hH596k7JuGywiMssJh1+stUWgaIypnd1urZ0Mp/uBtcAaYKBm\nner8efX0tC36BPx8vpPvPb2LbU/v4kPXbGLDms5FfT6q8vlk1FkriTVDMutOYs2wPHU34zYB812u\nOu9lrFWLvUQ2n+9kYGCY8fEpAPoHRsh5J9xN5FXrTpIk1gzJrDuJNcPS6l7oy6DRs19GjDG5cHod\nwdDMfoLeOnPmN53vBs0ulTX8IiJSq9FQfwLYGk5vBR4FngOuMMasMMZ0EIynP7X0Jr6T6wa987JC\nXURklhMOvxhj3gt8BdgIFIwxdwB3At80xvwGsBv4lrW2YIz5XeAxoAJ8yVo7tByN9sJQL5XLy7F5\nEZHIOpkDpS8QnO0y14111v0u8N2lN2thXjiOXlRPXURklkheUepp+EVEpK6Ihnp4oLSkUBcRqRXR\nUNeYuohIPdEMda8a6uqpi4jUimaouwp1EZF6Ihnq1fPUNaYuIjJbJEN95opSjamLiNSKZKhXx9R1\nSqOIyGzRDHVXFx+JiNQTyVDXmLqISH2RDHVPY+oiInVFM9Q1pi4iUlckQ93XeeoiInVFMtRdhbqI\nSF2RDHXd0EtEpL6Ihrp66iIi9UQz1D3dpVFEpJ5ohrp66iIidSnURURiJKKhrgOlIiL1RDPUNaYu\nIlJXJENd56mLiNQXyVCvXlGq2wSIiMwWyVCfuaGXQl1EpFY0Q706pl7SmLqISK1IhrrG1EVE6otm\nqDsOjqNQFxGZK5KhDsG4ukJdRGS2yIZ6yncoakxdRGSWyIa657oUdUWpiMgskQ31lO9SLKqnLiJS\ny2/kQ8aYDuDbQA+QAb4EHAT+EqgAL1lrP9WsRtbjuQ5F3SZARGSWRnvqnwCstfYG4A7gz4CvAZ+x\n1l4LdBtjbm5OE+tTT11E5J0aDfXDQF843QMcBTZZa7eH87YBW5bYtgX5nktBY+oiIrM0NPxirf2O\nMeYTxpidBKH+K8ADNav0A2tPtJ2enjZ831vUvvP5TgCyGZ9SuTL9Pu6SUmetJNYMyaw7iTXD8tTd\n6Jj6rwJvW2tvMsZcBnwPGKpZxTmZ7QwOji1qv/l8JwMDw8GbSoVCsTTzPsZm1Z0QSawZkll3EmuG\npdW90JdBo8Mv1wKPAVhr/wPIAStrlq8D9je47ZPiey6Viu6pLiJSq9FQ3wlcCWCMORMYBnYYY64L\nl38EeHTpzZuf7wVN17nqIiIzGhp+AR4EvmGM+VG4jf9OcErjg8YYF3jOWvtEk9pYlx/eqbFYKpNJ\nLW5cXkQkrho9UDoC/Jc6i65fWnNO3nRPXac1iohMi+wVpRp+ERF5pwiH+szwi4iIBKIb6n7Q9IJC\nXURkWmRDPRUOv5Q0/CIiMi2yoV59Tql66iIiMyIb6jM9dYW6iEhVZEPd8zSmLiIyV2RDPTV9nrrG\n1EVEqiIb6jqlUUTknSIc6tWLjxTqIiJVkQ91jamLiMyIbqj71eEXjamLiFRFNtTT4ROTCoVSi1si\nInL6iG6op4KmT+oujSIi06Ib6tWeelE9dRGRqsiGeiq8oddUQT11EZGqyIZ6Onza0ZSGX0REpkU3\n1Ku33tWBUhGRaZEPdR0oFRGZEd1QT+mURhGRuSIb6tMHStVTFxGZFtlQ9z0Xz3WY0imNIiLTIhvq\nEPTWCzqlUURkWqRDPe27OlAqIlIj2qGe8nRFqYhIjUiHesp3dUWpiEiNSId62vd0oFREpEa0Qz0V\nHCitVHRPdRERiHqo+y4V9KAMEZGqaIf69E29NAQjIgIRD3XdfldEZDa/0Q8aY+4EfgcoAl8AXgIe\nBjzgAHCXtXayGY2cT/VBGeqpi4gEGuqpG2P6gC8C1wG3AbcD9wMPWGuvB3YC9zSrkfOpPtJOV5WK\niAQaHX7ZAjxhrR221h6w1n4S2Aw8Ei7fFq6zrGZ66gp1ERFofPhlI9BmjHkE6AHuA9prhlv6gbUn\n2khPTxt+GMwnK5/vnJ7u7s4C0NaemTU/juJeXz1JrBmSWXcSa4blqbvRUHeAPuDDwJnAD8N5tctP\naHBwbFE7zec7GRgYnn5fnCoC0H94hIHuzKK2FSVz606CJNYMyaw7iTXD0upe6Mug0eGXQ8Az1tqi\ntfYNYBgYNsbkwuXrgP0NbvukTZ/SqAdliIgAjYf648AHjTFueNC0A3gC2Bou3wo82oT2LSibDkJ9\nYkqhLiICDYa6tXYf8F3gJ8C/AL9FcDbM3caYp4Be4FvNauR8culg9Gg8HIYREUm6hs9Tt9Y+CDw4\nZ/aNS2vO4mQz6qmLiNSK9BWl2bCnPqGeuogIEPlQD3vqk+qpi4hAxEM9p566iMgskQ716pj6uHrq\nIiJA1EN9+pRG9dRFRCDioe65LmnfZVxnv4iIABEPdYBsxtcpjSIioeiHetrT8IuISCjyoZ5L+zql\nUUQkFPlQz6Y9JgslymU9fFpEJBahDrpVgIgIxCDUcxldgCQiUhX5UFdPXURkRvRDPaPb74qIVEU+\n1KvDL+MTCnURkciHekc2CPVRhbqISPRDvT2XAmB0otDiloiItF70Qz0bhvq4Ql1EJPqhntPwi4hI\nVeRDvSPsqY+opy4iEv1Qnx5TV6iLiEQ/1LNpD9dxNPwiIkIMQt1xHNpzvs5+EREhBqEOwRkwGn4R\nEYlLqOd8RieKVCq6/a6IJFs8Qj2bolSu6KZeIpJ4sQl10FWlIiKxCPXOtiDUh8cU6iKSbLEI9e72\nNADHR6da3BIRkdaKRah3haE+pFAXkYSLRah3K9RFRADwl/JhY0wOeAX4Q+BJ4GHAAw4Ad1lrJ5fc\nwpPQpeEXERFg6T313weOhtP3Aw9Ya68HdgL3LHHbJ009dRGRQMOhbow5H7gQ+H44azPwSDi9Ddiy\npJYtQmdbGsdRT11EZCnDL18BfhO4O3zfXjPc0g+sPdEGenra8H1vUTvN5zvrzu9uzzA6UZh3edTF\nta6FJLFmSGbdSawZlqfuhkLdGPNrwLPW2reMMfVWcU5mO4ODY4vabz7fycDAcN1lHbkUR46Pz7s8\nyhaqO66SWDMks+4k1gxLq3uhL4NGe+q3AmcZY24D3gVMAiPGmJy1dhxYB+xvcNsN6e5Is3dghKlC\niXRqcb1/EZG4aCjUrbUfq04bY+4DdgHXAFuBvwlfH116807eivBg6bGRSVb1tJ3KXYuInDaaeZ76\nF4G7jTFPAb3At5q47RPq7coCcPT4KTmLUkTktLSk89QBrLX31by9canba1RvVwaAI8cnWtUEEZGW\ni8UVpQB90z11hbqIJFdsQr0nDPUjGn4RkQSLTaj3dgbDL+qpi0iSxSbUcxmf9qzP0WH11EUkuWIT\n6gA9nVmODE3oWaUiklixCvX8iiyThRLD43oCkogkU6xCfXV40dGho4u7/YCISFzEK9R7cwAcVKiL\nSELFKtTX9AY99f7B8Ra3RESkNWIV6qvDUFdPXUSSKlah3t2eJpP2NKYuIokVq1B3HIfVPTn6B8cp\n67RGEUmgWIU6BOPqU8Uyx3QRkogkUOxCvXpa4/4joy1uiYjIqRe7UF+/qgOAPf0jLW6JiMipF7tQ\n37A6DPVDCnURSZ7YhfrKFTlyGY+31VMXkQSKXai7jsP6fAcHjowyVSi1ujkiIqdU7EIdYP3qTioV\n2HdYB0tFJFliGeobwoOluw8Nt7glIiKnVixDfdPaLgDe3He8xS0RETm1YhnqZ+TbyWV8Xt97rNVN\nERE5pWIZ6q7jcM66bg4NjjM0OtXq5oiInDKxDHWAc9/VDcBO9dZFJEFiH+qv7RlqcUtERE6d2Ib6\nprVd+J7Ljt2DrW6KiMgpE9tQT6c8zIYV7B0YYVB3bBSRhIhtqANcclYfAK+8eaTFLREROTViHuq9\nALysUBeRhIh1qK/pbWNld5ZXdx2lUCy3ujkiIssu1qHuOA7vNXnGJ0u8+tbRVjdHRGTZ+Y1+0Bjz\nJ8D14Ta+DGwHHgY84ABwl7W25Uco/9MFq3ns+T08/4tDXH7uylY3R0RkWTXUUzfG3ABcbK29GrgJ\n+BpwP/CAtfZ6YCdwT9NauQQb13SSX5HlxdcP61a8IhJ7jQ6//Bj4aDh9DGgHNgOPhPO2AVuW1LIm\ncRyHKy9czeRUiZ/a/lY3R0RkWTU0/GKtLQHVm5XfC/wA+FDNcEs/sPZE2+npacP3vUXtO5/vXNT6\nALdvPpfvP7ubp185xO03nLfoz58OGqk76pJYMySz7iTWDMtTd8Nj6gDGmNsJQv2XgddrFjkn8/nB\nwbFF7S+f72RgYPH3SPeAizf18fKbR/j3Vw9MP5w6KhqtO8qSWDMks+4k1gxLq3uhL4OGz34xxnwI\n+D3gZmvtEDBijMmFi9cB+xvd9nLY/O4zAPh//763xS0REVk+jR4o7Qb+FLjNWls9V/AJYGs4vRV4\ndOnNa57Lzl5JfkWWp18+oNsGiEhsNdpT/xiwEvh7Y8y/GmP+Ffhj4G5jzFNAL/Ct5jSxOVzX4dar\nN1IsVfiX53a3ujkiIsui0QOlDwEP1Vl049Kas7yuuXgN255+ix/9bD+3XHUmKzoyrW6SiEhTxfqK\n0rl8z+XWazZSKJb5hx+92ermiIg0XaJCHeD9l57Bu/Id/NvLB3jrgB5MLSLxkrhQd12HO288F4Bv\nP2YplnSjLxGJj8SFOoDZ0MO1F69h98Fhvv+sDpqKSHwkMtQBPr7lPHq7Mmx7epeGYUQkNhIb6m1Z\nn3tvuYBypcJffO9ljo9OtbpJIiJLlthQB7hgYy8fvn4TR45P8hffe1nj6yISeYkOdYDbrtnI+85f\nxWt7h/irf/455XKl1U0SEWnYkm7oFQeO43DvrRdwbGSS53f0k/Y9PnHL+bjOSd2TTETktJL4njpA\nJuXx2TsuY+OaTv7t5QM89MireqapiESSQj3UlvX5/Mcu55x13Ty/o5+v/v3PGJsotLpZIiKLolCv\n0ZFL8dv/9XLefe5KfvH2Mb70ze3sPpi8+zyLSHQp1OdIpzw+/eFLuPXqMxk4NsEfP/xTHt++RwdQ\nRSQSFOp1uK7D1g+czWc/ehnZtM93nnydL//tC+w/PHriD4uItJBCfQGXnt3HH/36lVxx/ire2Hec\nL37jef7uidcYGddYu4icnhTqJ9DVnuZT//lifusjl9DbleGJn+7lf/6fZ/nnZ3YxNlFsdfNERGZJ\n/HnqJ+vd5+W5+Kw+fvjiPrY9/Rb/8OM3+cFPdrP53eu48X3r6enUAzdEpPUU6ouQ8l1++Yr1XHfJ\nWn70s308vn0Pjz73No8/v4fLzunj/ZedwSVn9eG6unBJRFpDod6AtqzPzVedyZb3vYtnXjnID1/c\nx4uvH+bF1w/T05nhivNX8T6zirPWdenKVBE5pRTqS5DyPT5w+To+cPk6dh08zo9/tp/ndvTz+PY9\nPL59Dz2dGd5zXp5Lz+7jvPUryKS8VjdZRGJOod4kG9d0sfGmLj6+5Tx27D7K9l/08+Jrh3nyhb08\n+cJefM/lvPXdXLSpl/M39LB+VQe+p+PUItJcCvUmS/kul569kkvPXknxpjKv7TnGq28d5dW3jvLz\nXYP8fNcgAGnfZePaLs5Z180567rZuLaT7vY0joZrRGQJFOrLyPdcLtzYy4Ube/noDTA0OsWOXUd5\nfe8QO/cN8fqeY7y259j0+l1tKdav7mTDqg7Wr+5g/apOVvfkWliBiESNQv0U6m5Pc9VFa7jqojUA\njE8WeXP/cd7YN8TuQ8Ps6R+Z7tVXuY7D6r428t1Z1vS2saavjbW9bazqaaO7I60DsSIyi0K9hXIZ\nn4s29XLRpt7peWMTBfb0j/D2oRH2DIxw8OgY/YPjHDg8yktvHJn1ed9z6O3K0teVZWV3lr7umene\nriwrOtKkfB2cFUkShfpppi2bwmzowWzomZ6Xz3fy1ttHOXh0jINHxsKgH+PI8QmODE2wY/fg/NvL\n+HR3pFnRkQle24PX7vY03R0ZOnMpOtpSdORSOnArEgMK9YjoyKWmD6rONVkocTQM+MNDExw5PsHR\n4xMcG5ni+OgUx0YmOXBk7IT7yKY9OnKpmZ+21Kz3bRmfXPjTlvWn32fSnoaBRE4TCvUYyKQ81va1\ns7avfd51CsUyQ6OTDI1OMTQyxdBIMD08XmB0vMDwWPg6XmDf4dFFPfnJAbIZn7aMRy6TCl99clmf\nXNonk/LIpL3p12zKI53yyNbMy6RcMmmfTMolrfP5RRqmUE+IlO+ysjvHyu6TO5tmslBiZKzAyPjM\nz9hkkfHwpzo9NjEzb3yyyJHj4+ybLLHUu89n0x5pPwj4dMoj5bmkUi5p3w2ng+Wp8Cfte+Grix++\nVudV56fC977v4nsOvjt72vMcPNfRaaUSaQp1qSuT8sh0e/R1Zxf92XKlwsRkKQj6qSKThRJTUyUm\nCiUmp0pM1rxOFEpMTZWZKBSZLJSn55cqFUbHCkwWShwfnWKqWKJQKC/5y+JEHMDzwqCf9RpMe17w\npVI77VXXcZ3gS6LmC8LzHFwnWNdzg3mu6+CHr17tuq5Dz6ERRkYmZ607s9yd89l3brc6rS+m5FKo\nS9O5jhOMuWcb/+eVz3cyMDD7UYKVSoVSucJUoUyhVKZQKFEolYP3xTKFYompYnW6HHwRTE/PzCuV\nKhRKZUqlMsVShWLNa6lUplCqBMvKFYrFMsVymbHJ4vR0qRS043TmOg6uGzzwxXXCH9fBdcCZNQ9c\n18V15q4bbKO6rlf3s87sz7n1P7tQW1zXobMjy9jYVLD9cJkTTjtOuC3mvJ+73HGmP19dNt/7+p+f\n/7PunP3U/Swz71139raqbT9Vmh7qxpivAlcBFeAz1trtzd6HJJPjONM951YrVyp1vxRqp8vlIPyD\nn/B9aWZeuVyhWJ69Xi6X5vjwBKVSedZ6pepnw/3O3vacbYXrlcs1P5UK5Qoz0+UKlUqFYglKheL0\n+2AZ0+uUKxUqp/f3VyTUfik5jkPKd/ncx9/DOWs6mr6vpoa6MeYDwLnW2quNMRcA3wCubuY+RE4H\nruPg+h6pJneL6v2F0mr1wr5S50ti/i+Pmc+Vqp8th+tVKnR2Zhk8Ngbh+0rNa6Xee2rWC7dXXVau\nQIU572ctP4nt13w2mL/wuuW5+69U6u8bqJSD9ruus2zPYGh2T/2XgH8EsNbuMMb0GGO6rLXHm7wf\nETlFHMfBcxyW6w+k0/GL7FRYrrqbHeprgBdq3g+E8+qGek9PG/4ir3jM5zsbblyUJbHuJNYMyaw7\niTXD8tS93AdKFzw6MDh44gtiaukbPTmSWDMks+4k1gxLq3uhL4Nm/0G1n6BnXnUGcKDJ+xARkXk0\nO9QfB+4AMMa8B9hvrU3eV7CISIs0NdSttc8ALxhjngH+N/DpZm5fREQW1vQxdWvt7zZ7myIicnJa\nfxWHiIg0jUJdRCRGnIquARYRiQ311EVEYkShLiISIwp1EZEYUaiLiMSIQl1EJEYU6iIiMaJQFxGJ\nkcg8ozRJj8kzxvwJcD3B7+fLwHbgYcAjuOvlXdbayda1cHkYY3LAK8AfAk+SjJrvBH4HKAJfAF4i\nxnUbYzqAbwM9QAb4EnAQ+EuC/7dfstZ+qnUtbC5jzMXAPwFftdb+uTFmPXV+v+G/g88CZeAha+1f\nN7rPSPTUax+TB9xLcLOwWDLG3ABcHNZ6E/A14H7gAWvt9cBO4J4WNnE5/T5wNJyOfc3GmD7gi8B1\nwG3A7cS/7k8A1lp7A8EdXf+M4N/4Z6y11wLdxpibW9i+pjHGtANfJ+igVL3j9xuu9wVgC7AZ+Jwx\nprfR/UYi1JnzmDygxxjT1domLZsfAx8Np48B7QS/6EfCedsIfvmxYow5H7gQ+H44azMxr5mgpies\ntcPW2gPW2k8S/7oPA33hdA/Bl/immr+841TzJHALwXMmqjbzzt/vlcB2a+2QtXYceBq4ttGdRiXU\n1xA8Gq+q+pi82LHWlqy1o+Hbe4EfAO01f4L3A2tb0rjl9RXg8zXvk1DzRqDNGPOIMeYpY8wvEfO6\nrbXfATYYY3YSdGB+GxisWSU2NVtri2FI16r3+52bb0v6bxCVUJ9rwcfkxYEx5naCUP/NOYtiV7sx\n5teAZ621b82zSuxqDjkEvdaPEAxL/F9m1xq7uo0xvwq8ba09B/gg8DdzVoldzQuYr9Yl/TeISqgn\n6jF5xpgPAb8H3GytHQJGwoOIAOuY/edcHNwK3G6M+Qnw68AfEP+aAQ4Bz4Q9ujeAYWA45nVfCzwG\nYK39DyAHrKxZHseaa9X7dz0335b03yAqoZ6Yx+QZY7qBPwVus9ZWDxo+AWwNp7cCj7aibcvFWvsx\na+0V1tqrgL8iOPsl1jWHHgc+aIxxw4OmHcS/7p0EY8gYY84k+CLbYYy5Llz+EeJXc616v9/ngCuM\nMSvCs4OuBZ5qdAeRufWuMeZ/Ae8nOOXn0+G3fOwYYz4J3Ae8VjP7boKwywK7gf9mrS2c+tYtP2PM\nfcAugt7ct4l5zcaY3yAYZgP4I4LTV2Nbdxha3wBWE5yy+wcEpzQ+SNDJfM5a+/n5txAdxpj3Ehwr\n2ggUgH3AncA3mfP7NcbcAfwPgtM6v26t/dtG9xuZUBcRkROLyvCLiIicBIW6iEiMKNRFRGJEoS4i\nEiMKdRGRGFGoi4jEiEJdRCRG/j+bBqplQUeepwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "tags": [] } } ] }, { "metadata": { "id": "30iAu5rEgAGl", "colab_type": "text" }, "cell_type": "markdown", "source": [ "We can even compare multiple functions. For example, the following code compares multiple elasticities." ] }, { "metadata": { "id": "C_mL5D03f7VW", "colab_type": "code", "outputId": "259d1d01-8640-49d3-f95b-8e3ca2b49eb7", "colab": { "base_uri": "https://localhost:8080/", "height": 265 } }, "cell_type": "code", "source": [ "plt.plot(p, d_const_elast(p, 0.5))\n", "plt.plot(p, d_const_elast(p, 1))\n", "plt.plot(p, d_const_elast(p, 1.5))\n", "\n", "plt.legend(['$\\epsilon = 0.5$', '$\\epsilon = 1$', '$\\epsilon = 1.5$'])\n", "\n", "plt.show()" ], "execution_count": 0, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8JHWd//FXHX13J+kkncxk7gNq\nGGZgBIYbRORQwEUBXUWQAVFXUXFZd9efrqKuP3/7c9f1XsVrBdyHrrgqeAE/FFQcBgaEmQFmijmY\neybJ5OwkfVf9/qjqpJNJJlcnne7+PH20XVdXfb8T8q5vvvXtKsW2bYQQQlQGtdQFEEIIUTwS6kII\nUUEk1IUQooJIqAshRAWRUBdCiAqil/Lg7e3xSQ29iUaDdHUNzFRx5qxqrHc11hmqs97VWGeYXr1j\nsYgy1rqyaqnrulbqIpRENda7GusM1VnvaqwzzFy9yyrUhRBCnJiEuhBCVBAJdSGEqCAS6kIIUUEk\n1IUQooJMaEijYRhrgAeBL5mm+XXDMBYB9wMacAS42TTNlGEY7wQ+AljAt03T/N4MlVsIIcQoxm2p\nG4YRAr4G/K5g8WeBb5imeRGwC7jN3e5TwGXAJcDfGoZRX/QSCyGEGNNEul9SwFXA4YJllwAPudO/\nxAnyc4DNpmn2mKaZAP4MXFC8og55YfNGfvmVz5HJZmZi90IIUbbG7X4xTTMLZA3DKFwcMk0z5U63\nAfOBeUB7wTb55WOKRoOTHoAfi0U4/IffYOw4yIE92znnghk5b8w5sVik1EWYddVYZ6jOeldjnWFm\n6l2M2wSM9XXVMb/GmjfZr8jGYhHa2+PYljN/rKOb9vb4pPZRjvL1ribVWGeoznpXY51hevU+0clg\nqqHeZxhGwO1mWYDTNXMYp7WetwDYNMX9n5CiOr1GmUx2JnYvhBBTtmnTRr7ylX/DsiyuuebN3Hzz\nhlG3u/TSS/H5/KiqhqZpfO979xfl+FMN9ceA64Efuu8PA08D3zUMow7I4vSnf6QYhRwpH+rpjPSp\nCyHmjlwux7//+//lS1/6Bk1Nzdx++7u48MKLWbZs+ajbf/Wr91BXV1fUMowb6oZhnAl8EVgKZAzD\nuAF4J/ADwzDeB+wD7jVNM2MYxseARwAb+Ixpmj1FLa1LUZ1++JyEuhCiiDZt2si3vvV1ADweD/fc\n85+o6sS/zrN9+0ssXLiIBQsWAnDZZVfw5JN/GDPUZ8JELpQ+hzPaZaTLR9n2p8BPp1+sE1M0J9Rl\n9IsQlecnv9/F5h1tRd3n+lVNvO3SleNu9+Uv/ytf//p3aGxsPG7dBz5wOwMDx18HvOOOO1m//hwA\n2tvbaGpqHlwXizXx8ssvjnm8u+66A1C49trruPba6yZQk/GV9H7qU6Vqzpkzm5U+dSFE8Zx77gVs\n2PB2Lr/8jdx5598NW/cf//Hdoh7rRz/6EaoapKurk4985A6WLFnKunVnTHu/ZRrqTktdQl2IyvO2\nS1dOqFVdbNu2bQFsfvGLh9H146NxIi31WKyJtrbWwXXt7W3EYk2jHq+5uZn29jjRaD0XX3wJL7/8\nUvWGer77JSehLoQokscff4xFixaj6zq2bTMw0E8oFB5cP5GW+qpVqzlw4ACHDx8iFmviscce5e67\nP3fcdolEgr4+ZXB68+an2bDh9qLUoyxDXZVQF0IU2WWXXcm//Ms/89BDP8fr9fF3f/cxVq06ZVL7\n0HWdu+76e+6660NYVo6rr/4rli9fMbj+ox/9MB/72CdJpVJ88IO3k81a5HI5Lr/8Ss499/yi1KMs\nQ113Q92yciUuiRCiUqxevYb77vvvae/nvPMu5LzzLhx13b/921cHpx966KEZ+dJVWd56V9Wcc5GV\nk5a6EEIUKstQ19yLGLmstNSFEKJQWYa67pHuFyGEGE1Zhrrm9qnbEupCCDFMWYa67vEAYOUk1IUQ\nolBZhnr+QqktoS6EEMOUZ6i7N/SybQl1IYQoVJ6h7vapY1mlLYgQQswxZR3qcqFUCDHXfP7zn+Ga\nay7n5pvfVpLjl2eoq+4XYaWlLoSYY6666k188YtfK9nxy/I2AfkLpYMPKxVCiCKY7kMyANatO4Mj\nRw7PRPEmpExD3el+UWybnGWhTfIfXQgxd/1s1694vm1bUff5mqa1XLfymnG3m+5DMuaC8gx1d/SL\nikU6YxHwSagLIaZvNh+SMVPKMtTz935RbJt01iLgK3GBhBBFc93KaybUqi62YjwkYy4oy1DPXyhV\nsUmms9SGvCUukRCi3BXjIRlzQVn2W+QvlCq2RTIlwxqFENN32WVX8uCDP+OWW97Oe9+7gQMHDkxp\nP3ff/XH+5m9uZf/+fbzlLVfxq1/9osglPbHybKnrQy31REruqS6EmL5iPSTjM5/5fBFKM3Vl2VLX\n8hdKbZtkWlrqQgiRV56hrjl3aVSwSKSlpS6EEHnlGeruhVIFSEr3ixBCDCrLUB/sU7ctEtL9IoQQ\ng8oy1BV16BulcqFUCCGGlGmoK0B+nLq01IUQIq8sQ518Sx1b+tSFEKJAWYa64t7AS7Vt6VMXQogC\nZfnlIwpDXVrqQog55POf/wwbNz5JNBrl/vt/MuZ2l156KT6fH1XV0DSN733v/qIcf0qhbhhGGLgP\niAI+4DPAUeCbgA1sNU3z/UUp4SgUzQl1xb33ixBCzBVXXfUmrr/+r/nc5z417rZf/eo91NXVFfX4\nU+1+2QCYpmm+DrgB+ArwZeBO0zQvAGoNw3hjcYo4CmWopS4XSoUQxbJp00Y2bLiRDRtu5D3vuQVr\nCk9XW7fuDGpqamagdBMz1e6XY8Bp7nQU6ASWmaa52V32S+Ay4LfTK97oBvvU5d4vQlSc9gd+TPzZ\nzeNvOAmRs9YTe+vbx91uth+ScddddwAK1157Hddee92U9jHSlELdNM0fG4axwTCMXTih/ibgGwWb\ntAHzx9tPNBpE17VJHTsWi5ANqOzGGaeeyuSIxSKT2kc5qoY6jlSNdYbqrHdhnfsCXga04o7hCAS8\nE/p3fd3rLuHWW9/Bm970Jj7xiU8MW/fAAxO/2VcqFULXtRMe80c/+hHNzc10dHRw6623cvrpq1m/\nfv2EjzGWqfap3wTsN03zDYZhnA78HOgp2ESZyH66uo4/651ILBahvT2OlUwAQxdKW1t7UdUJHbIs\n5etdTaqxzlCd9R5Z5/A11xG+pjit1kLj/btu27aFRCLNz3/+W3RdP277ybTUOzv7yWZzJzxmc3Oz\nu97L+edfzMaNm1m6dNWE6nKik8VUu18uAB4BME1zi2EYAcBTsH4BMHNPXnX71BXb+b+BVJZwwHPi\nzwghxAnM5kMyEokEfX3K4PTmzU+zYcPtRdn3VP/G2QWcA2AYxhIgDmw3DONCd/11wMPTL97olIIH\nT6NY9CcyM3UoIUSVmI2HZHz0ox/m2LF2Ojs7uPHGG7nllnfwnvfcwnnnXcC5555flHootm1P+kPu\nkMbvA804rf1P4gxpvAfnRPG0aZp3jbef9vb4pA6e/zPNtix2vvc2DjR7+K/o9XzixvNY0VI76XqU\nC/mTvHpUY72rsc4wvXrHYpEx+5uneqG0D3jbKKsumsr+Jk1x6uN0v0hLXQgh8srzNgGKgq0oqDYo\nikV/QoY1CiEElGmoA9iqgmLZoFr0JaWlLoQQUOahrlqAYkv3ixBCuMo21FFVNNtpqUv3ixBCOMo2\n1G1NRbXcPnXpfhFCCKCMQx3VCXVUiz7pfhFCCKCcQ11TUS0bTbOlpS6EEK6yDXVb01Bt8PsV6VMX\nQghX2Ya6ojotdZ9PIS7dL0IIAZRxqKNpaBYEfCqJVJZMdvI3sxdCiEpT1qGuWjZ+n3PLgPhAusQF\nEkKI0ivbUFfcPnWvz5nv6ZdQF0KIqd5PveQUTUO1wOtxbvQooS6EEGXcUld153zk0Z2+9F4JdSGE\nKN9QV/KhruQACXUhhIAyDnVVcx5fp0moCyHEoPINdd15pJ2qOqEufepCCFHWoe601BU7g4K01IUQ\nAso51N3uFyubIRTwSEtdCCEo41DX3JZ6NpumPuKjK55iKg/RFkKISlK2oZ4f0pjNZKiv8ZPK5BhI\nyY29hBDVrWxDXdGcULeyaeprnK+VdvQkS1kkIYQouTIOdWf0Sy6Xpb7GD0BnPFXKIgkhRMmVf6hn\nMtRHnJZ6V6+01IUQ1a3sQ93KpQdb6h290lIXQlS38g1190KplckN9ql3xqWlLoSobmUb6gy21DPU\nhX0oCnTKhVIhRJUr21DPd78oloWi2NSFfXKhVAhR9co+1FUL0lZm8AtIliVfQBJCVK8KCHWbjJWh\nodZPzrLp7pPWuhCiek35yUeGYbwT+AcgC3wK2ArcD2jAEeBm0zRnLmHdLx+pFqRzGZqiAQBauxKD\no2GEEKLaTKmlbhhGA3A3cCFwDXAt8FngG6ZpXgTsAm4rViFHM9hSt52WenM0CEBr18BMHlYIIea0\nqXa/XAY8Zppm3DTNI6Zpvhe4BHjIXf9Ld5sZkw91zYJ0Lj0Y6m2diZk8rBBCzGlT7X5ZCgQNw3gI\niAKfBkIF3S1twPzxdhKNBtHdh11MVCwWAUBpiHAUUHM2oRoPLc0xALr604PbVJJKrNN4qrHOUJ31\nrsY6w8zUe6qhrgANwFuAJcDj7rLC9ePqmmRXSSwWob09DkDfgHNHRs2Cto4eGupTBH06+4/2Dm5T\nKQrrXS2qsc5QnfWuxjrD9Op9opPBVLtfWoGNpmlmTdPcDcSBuGEYAXf9AuDwFPc9IYrHuZ+6nrNJ\nWxkURaG5PkB7d0KGNQohqtZUQ/1R4FLDMFT3omkYeAy43l1/PfBwEco3pvxtAjTLJpV1en2ao0Gy\nOZtOubGXEKJKTSnUTdM8BPwU2AT8FvgQzmiYWwzD+BNQD9xbrEKOZjDUc5DKOaGeH9Z4VEbACCGq\n1JTHqZumeQ9wz4jFl0+vOBOX737RLJuk21JvaQwBcLi9nzXLGmarKEIIMWeU7zdK3Za6nrNJui31\nhbEwAAfb+0tWLiGEKKXyDfXBljokc04fenN9AF1TOdDeV8qiCSFEyZRvqOtuqOeGul80VaWlMcjh\nY/0yAkYIUZXKNtTVwdEvDHa/gNMFk8lacrsAIURVKttQVzz50S82yezQEMZ8v/oh6VcXQlSh8g11\nt/tFH9lSb3JGwOxvk351IUT1KdtQR9NAUfBYymCfOsDiZufrs/uOVt/XjoUQomxDXVEUFF1Ht5TB\nLx8B1AS9NNb6efVIL7YtF0uFENWlbEMdnLHqHothfeoAy1tq6EtkaO+W2/AKIapLeYe6x+PeJiCN\nZVuDy5fPrwFgz5HeUhVNCCFKorxDXfegWTY2NulcenD5shYn1F89LP3qQojqUt6h7tHRck6/eeEI\nmCXNETRVYffhnlIVTQghSqK8Q133oLrfHC0cAeP1aCydF2HvkTjJdLZUxRNCiFlX3qHu8aBmnb70\nwhEwAMbiKJZts/OgtNaFENWjvENd11FyTqgPZIePdFm1pA6AHfu7Zr1cQghRKmUd6qrHg2LbqJbN\nQGb4vV5OWlCHpirs2NddotIJIcTsK+tQV7xewLmnen9meEvd59VYNr+GfUfjJFLSry6EqA5lHepq\nPtSzNgPZ4+/KuGpJHZZtY+6X1roQojqUdagrXh8Aeg76M8eHev6Rdlt3H5vVcgkhRKmUeag7LXVP\n1h411FcsqCEc8LBld4fcB0YIURXKOtTH637RVJW1yxvoiqfY3yq34hVCVL6yDvV8S91rcdyF0rx1\nJzUC8PzO9lkrlxBClEpZh7rq9qmHbO9xQxrz1iyrR1MVXtgp/epCiMpX1qGu+JyWegjPqH3qAAGf\nzqnL6tnf1seRDnnEnRCispV1qKse55F2QctDf3ZgzIuh56xuBuDpl1tnrWxCCFEKZR3q+SGNQVvD\nsq3j7v+S95qTGvHqKk+/3CqjYIQQFa2sQz0/+sVva8DYF0v9Xp11JzXS2pVgrzy7VAhRwco61POj\nX/yWDkBfZuxhi+edOg+AP205PPMFE0KIEinrUM+PfglYTku9Nz12K3zt8gbqa3w89XKr3AtGCFGx\nyjrU86NffJZTjROFuqoqXHx6C6l0Ti6YCiEqVlmHer5P3Ztz5ntTJ/7W6EWntaAqCr//yyG5YCqE\nqEj6dD5sGEYAeBH4Z+B3wP2ABhwBbjZNc/ThKEWSH/3iyToBfaKWOkA04uMMI8azO9p4eW8Xpy6r\nn8niCSHErJtuS/2fgE53+rPAN0zTvAjYBdw2zX2PSw0EANDTTlM9Pk6oA1x17mIAfrNp38wVTAgh\nSmTKoW4YxipgNfBrd9ElwEPu9C+By6ZVsglQfU5LXUmnURV13JY6wNJ5NZy6NMr2fV28eqR3poso\nhBCzajrdL18EPgjc4s6HCrpb2oD54+0gGg2i69qkDhqLRYbN7wkE0LIZanxh+nP9x60fzTuuPIV/\numcjj/3lEB/fsGBSxy+VidSr0lRjnaE6612NdYaZqfeUQt0wjHcBT5mm+aphGKNtokxkP11do9+v\nZSyxWIT29uGtccXvJx3vJ6y30J44dtz60cyv87FiQQ1PbTvCpi0HWdFSO6lyzLbR6l3pqrHOUJ31\nrsY6w/TqfaKTwVS7X64GrjUMYxNwO/BJoM+9cAqwAJiVb/mo/gBWMkGNN0IqlyaZHf/arKIo3PDa\nFQD89PHdMhJGCFExptRSN03zr/PThmF8GtgLnA9cD/zQfX94+sUbnxrwk2lvo9ZXA0BPuhe/Hhv3\nc8biKKetaGDr7g627engtBWNM11UIYSYccUcp343cIthGH8C6oF7i7jvMan+AHY2S73u/DnSlZz4\nQ6ZveO0KFAV+9LtdZLLWTBVRCCFmzbTGqQOYpvnpgtnLp7u/ycoPa6wnCEDnJEJ9YVOY15+xkMee\nO8hvN+3jry5cNiNlFEKI2VLW3ygFp6UOUIcfgK5k16Q+/5aLl1Mb9vKrp/bR2jm5C7dCCDHXlH+o\nB5wwr7GdWwZ0pibeUgfnyUg3XnYy2ZzFvQ/vwJKLpkKIMlb+oe621MM55ylIk+lTzzvLiLFuZSM7\n9nfz6DMHilo+IYSYTeUf6m6fuprOEPGEpxTqiqKw4apV1IS8/M8fdrNPHqQhhChT5R/qbkvdSiaI\n+uvoTHVj2ZMfyVIT9PLuq08hZ9nc89BLcs91IURZKv9Qd/vUrUSCxkA9WStLT2pq93RZu7yBK9Yv\n4mjnAN/91cvSvy6EKDvlH+r+fKgnaQo6XzpqGzg25f299XUrOGVJlOd3HuOXf95bjCIKIcSsqYBQ\nH+p+aQo43wptS7RPeX+aqvI3155KY62fB598lY0vHilKOYUQYjaUfahrQedLR7n+/qK01AEiQS93\n3nAaQZ/O93+9g627p7c/IYSYLeUf6mHn9gC5vj6agm5LfWDqLfW8BbEwd771NDRN4T9+/iI7D05+\nVI0QQsy2Cgj1MAC5vjghT5CwJzTtlnreSQvreP+b15DN2fz7T7Zg7p/ct1WFEGK2lX2oK7qOGgiQ\n63MeOt0UbORYspOclSvK/tetbHSCPWvxpQe28PLezvE/JIQQJVL2oQ5OF0yuz/nCUFMghmVbHEsW\nL3zPNGLccd1aLMvmyw9s5ZntrUXbtxBCFFNlhHokjNXXh23bzAs1AXCk72hRj7FuZSN3vvV0PLrC\ntx58iV8/tVceriGEmHMqI9TDEexsFjuVZGG4BYCDfcV/8NKpS+v5X+88k/oaH//zhz18/zfbSWeK\n080jhBDFUCGh7l4sjfexMOKE+oH4zDxNb2FTmE/cfBZL5kX487aj/O/7n6N1ks9aFUKImVIhoZ4f\n1hgn4g1T662ZkZZ6XjTi4+M3ncEl61o40NbHZ3+wmc072mbseEIIMVEVEur5YY3OCJiFkRa6Uz30\npftn7JgeXeNdb1jF7dc4NwH75i9e5NsPvURfIjNjxxRCiPFUSKgPtdQBFs1gv/pI56+Zz90b1rNs\nfg2bXm7lk997mhd2yTdQhRClURmhHsn3qTuhvjCyAID98YOzcvz5DSE+fvMZXP/a5fQNZPjqT7fy\njZ9vo6MnOSvHF0KIvGk/eHou0GpqAcj2OF/lX1a7GIA9PXtnrwyqytXnLeX0FY3c94jJc2Y723Z3\ncPX5S3nD2Yvx6BVx/hRCzHEVkTR6NApAtssJ9TpfLQ3+evZ075vSAzOmY2FTmP910xm8++pT8Pt0\nfv7HPXziO5vY+OIRLEvGtQshZlZlhHpNLSgK2e6he7OsqFtKf3aAo/2zPypFURQuWDufz7/nXK5Y\nv4juvhTf/dV27v7+M/zllXb50pIQYsZURKgruo5WU0u2a+jWACtqlwKwexa7YEYK+nXe/vqT+Px7\nz+XCtfM53NHP13+2jc/e+yzP7miTlrsQougqItTB6YLJdnUNtoJX1i0DYFf3nlIWC4DG2gC3XX0K\nn7v9HM5a1cT+o3H+4xcv8vHvbOKJ5w+Rycq3UoUQxVERF0rBCfXU3lex+vvRwmGag03UeCPs6NyJ\nZVuoSunPX/MbQnzgzWs40tHPI88cYOOLR7jvEZNf/GkPF53ewmvXtdBYGyh1MYUQZaz0SVckel3+\nYqnTBaMoCqsbDPoy/RyIHypl0Y4zvyHEhjeu4gvvP583nruYnGXz66f28Y/feoqvPLCFrbuPSdeM\nEGJKKqal7nFHwGS6uvAtcoY0ntqwik1HnuWljh0sqVlUyuKNqi7s462XrOTaC5bxzPY2Hn/+EFt2\nd7Bldwf1NT7OXT2P805tJhaLlLqoQogyUTGhPjSscehi6aroSaiKyksdJlctu7xURRuX16Nx4Wnz\nufC0+ew7Gufx5w+xeUcrv9m0j99s2sfyllrWr4px9inNRCO+UhdXCDGHVUyoe2LOfdQzbUNDGIOe\nACtql7Kr+1W6kt1E/XWlKt6ELZkXYcMbV3HjZSfxwq5jbHqplW17OthzuIef/H4XKxbUcsbJMc4w\nYjTVSf+7EGK4KYe6YRhfAC5y9/F/gM3A/YAGHAFuNk0zVYxCToSnqRkYHuoAZzSdzs7uPTzftpVL\nF188W8WZNq9H4+xTmjn7lGa8AS+/fXIPm3e0sfNAN7sO9fCTx3exqCnMmSfHOH1lI4uaw6iKUupi\nCyFKbEqhbhjG64A1pmmeZxhGA/A88DvgG6ZpPmAYxueB24BvFq+oJ6ZFIqiBAOm24Y+ae03TWh7Y\n+SDPlVmoF6oN+3j9mQt5/ZkL6e1P88KuYzxntvPy3k4OtPXxiydfpSbo4dRlDaxdXs+py+qJBL2l\nLrYQogSm2lL/I/CMO90NhIBLgL9xl/0S+CizGOqKouCJNZE+egTbslBUZ2BPxBvGiK5ke+crtA0c\noynYOFtFmhE1IS8Xn97Cxae3MJDMsm1PBy/u6WDbq5089dJRnnrpKAqwdH6E1UvrMRbXsXJBLX5v\nxfS0CSFOYEq/6aZp5oD8zcrfDfwGuLKgu6UNmD/94k2Op6mZ1P59ZHt6BkfDAJw770y2d77Ck4c3\ncd3Ka2a7WDMm6Nc5Z3Uz56xuxrJtDrb1uSHfya5DPbx6JM6vn9qHpiosmRfBWFSHsTjKSQtrCfgk\n5IWoRNP6zTYM41qcUL8C2FmwakKdu9FoEF3XJnXMEw3vG1i6kL5nnyGU6qU2tnhw+WX15/E/u3/J\n00ef49b11+PVy69rYiLDGpubajhzjXMv+YFkhu17O9m26xgv7ulg14Fu9hzu5bdP70dVYMn8Gk5e\nHGXVkijGknoWxMKo6tzqk6/WoZzVWO9qrDPMTL2nc6H0SuATwBtM0+wxDKPPMIyAaZoJYAEw7hMq\nuib5bM9YLEJ7e3zM9dnaBgDatu8iPW/JsHXnzlvPo/se59GXN3LO/DMnddxSG6/eY1ncEGRxw2Ku\nPmcxyXSWXYd6MPd388qBbvYejfPq4V4e2bQPgIBPZ/n8CMtaalneUsOS5gh1YS9KiS6+TrXO5a4a\n612NdYbp1ftEJ4OpXiitBf4VuMw0zfzA8MeA64Efuu8PT2Xf0+Fb6HzBKHXgwHHrLmw5h/+37wl+\nf+BPnD3vjJKFVan4vTprljWwZplz4svmLA6297HncC+7D/Wy50gvL+3t4qW9Q3e6jAQ9LG4Ks6g5\nwuLmMIubIsyrD865Fr0QYshUW+p/DTQCPzEMI7/sFuC7hmG8D9gH3Dv94k2Od34LaBqpg8eHekOg\nnjObT+fZ1hfYduxlToudOtvFm1N0TWXpvBqWzqvh0jOcZX2JDK8e6eXVI70caO1jX2v8uKD36ioL\nm8IsjIVpaQyxoDFES2OopK16IcSQqV4o/Tbw7VFWlfRrm4qu4503n9Shg8NGwOS9Yenrea51C7/Z\n+xhrG1dLCI0QDnhYu7yBtcsbBpcNJDMcaOtjf2sf+1vj7G/rY9/ROHsO9w77bMCn0dLgBPzgqyFE\ntMYn4+eFmEUVNwTCt2gR6UMHybS34W2eN2zd/FAzr2lay1/atvJ8+zbOaDqtRKUsH0G/B2NxFGPx\n0GiiTNaitWuAw8f6B1+HjvWz92ic3SPC3qOrNNUFaIoGaK4POu/RIM3RAHURCXwhiq3iQt2/eAnx\nTU+R3PvqcaEO8KblV7Kl/SV+tvNXrGlYhVcrv5EwpebRVRbGnC6YQtmcRWtXgiNuyB/p6Ke1M0Fr\n1wCHjvUftx+vrhJzQ74pGqCx1k9jrZ+G2gCRGrkFghBTUXmhvmIlAIldu6g557zj1jcFY1y66CL+\n3/4neHTfE1yz/IrZLmLF0jWVBW4/+1kFy23bJj6QobVrYDDk27qc99auBIfajw98cC7UNtT43bAP\n0DAY+n4aavwy1l6IUVTcb4V/yVIUXSe5a+eY27xh6aVsbn2eR/b9ntNiq1kcWTiLJaw+iqJQE/JS\nE/Jy0sLhN1WzbZvegQxtXQN09CQ51pOkozdJ70CGw8f6OdjudOuMJuDTqAv7qI/4qIv4iEb8RCM+\n5xX2Ea3xEQl45NqJqCoVF+qKruNftpzErp1YyQSq//g/4/26n5tOeStff+G73PvSj/nH9Xfi1Twl\nKK1QFIXakJfakJeTCs6t+TG8lm3T2592wr4nybGehPPem6Q7nqIrnuJIx9jfd9A1hbqwbzDs68I+\natzj1bonmtqQl0jQK0M1RUV+yttTAAATz0lEQVSouFAHCJx0MomdrzBgmoRPXzfqNqfUn8xrF17A\nHw7+mZ/ufIgbV10/y6UUE6EqTijXhX2sXFA76japTI7uvhTd8RSd8VHe+1LsOtSDfYKHSSkKRAIe\nakI+asNeaoLeYe+1IWc6EvQQCnjQtYp5aJioMBUZ6qG1p9H5m1/Rv3XLmKEO8OYVV7Grew9/Pvw0\nC8PzuXjh+bNYSlEsPo/mjqgJjrlNzrLo7c/Q05+ipy9NT7/z6s2/96XoGcjQ0ZvgYHvfuMcM+HQi\nQQ+RgIdwwEMk6CWcnw96iATceXdZwKdLN5CYFRUZ6v7lK1CDQfq3bcW27TF/mbyah/et3cAXnv0q\nD+x8iMZAA6sbjFG3FeVNU9XBLpjxpDI5egsCv6c/TU9fivhAhngiQ99A2n3P0NGTJDeB58lqqkIo\n4IR8yO8h5NcJ+vXB6VDAQ9Cv09I8QDaVHVrm06VbSExKRYa6ommETl1DfPMzpA8fwrdg7AuhDYEo\n71n7Lr72wnf49rb7+MDpt3FydMUsllbMNT6PRqwuQGwCT5aybZtEKkc8kabPDf34QJo+N/TjAxn6\nEpnB9V29KQ639zOZx4oHfBohv2fYSSDo9xAK6AR9zitw3Esj6NPxe+WkUG0qMtQBQqetI775Gfqe\ne/aEoQ6wsm4Z71lzM9/edh/f3PqfvG/tLayqP2mWSirKmaIoBN1Wd3N0/O0BLMtmIJVlIJmhP5ml\nP5lhIJmlP5EBTaO9o5++gmX9ySwDqQytnQlSmfG7hkbyebWC4Necd68zH3SX+UecHJwTgua+dLwe\nVbqPykTFhnr4Na9B8Xrp3fhn6q/5q+NuGTDSmsZTePeam/j+iz/kG1u+x02r3lp2d3MU5UFVFcJu\nX/xI496JNGc5Ye+eEAaSGQZSWZKpHIlUloFUlkQqS8Kdz78GUll6+9O0dmYn1F00koJzcvC5Ie/3\nOIHvKwj+/Elg9G10fF6NQME22ji/k2JqKjbUVX+A8JlnEX9qI4ldOwmePH5f+emxU/nguvdwz7Z7\nuW/7f3N0oI1rll2Bpk7unu9CzBRdUwfH/E+Fbduks1ZB4A8P/6ETQ45kOksynSOVyZFMZUlmciTT\nznR3PEUqk5tWXTy6is+jEQx40FUFn0fFqzuB73XXeT2a+z40XzjtG2V9fj/V2u1UsaEOUHv+hcSf\n2kjPH5+YUKgDnBRdzt+d+QG+tfUHPLrvcfb07GXD6ncQ9deN/2Eh5jhFUQaDsC48/kXjE7Fsm1Q6\nNxT86SypdI5EOucuzw6uHzwhFC5z57OWTW9/hnQmRzprFammQyeNwcDX3cB36+/xqHh1FY+uue/O\nOo/uLPfq7rRHxaOpeDwF2+lDn59rJ5CKDvWAsQpvSwvxzc/Q+Jbr8TRM7Pmk80PNfGz9h/nh9p/y\nQvs2Pvf0v/PmlVdxQcvZqIr8ySgEON8hyPfBT0dhl5Nl22QyFqmMc6JIZ3KkMpb7nl/mrE9nnZNH\nOuvOp3Okss62Q9s7832JDOnM9P+6GIumKk7462MHv2dwuUrAp3PD5cbEHhE3SRUd6oqqUv+Gqzn6\n/e/Q9cjDNN1404Q/G9AD3L7mJjYefoaf7fo1PzZ/xuajz3P9SdewpGbRDJZaiOqlKspg3/1MyHc/\npTI5MhmLdDZHJmuRyVqksxaZrHPScOZz7jLnxDBsm6zlfr5gvmB9fCBDJpsinc2N+aW3JQvqOPvk\niTU0J6OiQx0gcvY5HHvwZ/T88QnqLr8Cb6xpwp9VFIULFpzDqY2r+MkrD7Kl/UW+8OzXOKt5Hdcs\nu5JYsGH8nQgh5ozC7qfZYNs2Oct2TxRO+KezFrZtc/qqeXR0TH4003gqvi9B0XUar38rdjZL+09+\nPKV91Plqee/ad/Hhde9lUWQBz7a+wGc2fYHvv/hfHIgfKnKJhRCVQlEUdE0l6NepDfuI1QVY0Bhi\n4Qw+6L3iW+oAkfXn0PP47+l//i/En91M5Kz1U9qPUb+Sf4h+iL+0beXRfY/zXNsWnmvbwsnRlVzY\ncjanxdbgUavin1QIMUdVRQIpikLzLbey77N303r/D/AvX46nfmpdJ6qiclbzOs5sOp2XO1/hsX1P\n8ErXLl7p2kVID3L2vDM4o/l0ltYskouqQohZVxWhDuCdN5/Y295O2w/v4/DXvsyif/z4qLflnShF\nUTi1weDUBoPW/jaeOvIsm448y+MHn+Txg09S56tlXWwN62JrWV67RMa6CyFmhWKf6H6kM6y9PT6p\ng4/3bbvx2LZN2w/vpecPTxA8dQ0td3wY1Vu8x9nlrBwvd5q80PYiW469RCKbAMCv+Tk5uoJT6k9i\nVf3JxAINk/rK9XTrXY6qsc5QnfWuxjrD9Oodi0XGDJCqaamD07puuvFmsl1d9G/dwqGvfomWOz6M\nFijO8zA1VWNt42rWNq7mHVaWV7p2s/XYy+zofIWtx15i67GXAKj11rCsdgnL3dfCyALpixdCFEVV\ntdTzrEyGo9/+Fn3PP4dn3jxa3v8hfAsWTHu/J3Is0cmOzlfY0bmTPT176UkP1UNXNFrC81gQbmFh\nuIWFkRYWhOcR0J2TTTW2ZKqxzlCd9a7GOsPMtdSrMtQB7FyOYz97gK5HHkbxeml883XUvf5yFG3m\n+75t26Yz2c2rPXvZ07uPV3v2cbi/layVHbZdgz9Kc7CJpQ0tRJRamoIxmoMx6ny1FX/HPPlFrx7V\nWGeQUAdm5ocff+5ZWu//AVZfH75Fi2h4yw2E1p4266GZs3K0DrRzsO8wB/sOcyh+hEP9R4inj/9y\nglf1EAs2Uu+vo94fLXg58xFPuOxDX37Rq0c11hkk1IGZ++Hn4nHaH/hvejc+CYB/+XLqXn8F4TPO\nRPWU9oHUA5kEGd8A5uF9tA600zbQTutAO+2JDtK59Kif8ag6UV8dNb4Itd4aarwRanwRarzuvDsd\n8gTn7LBL+UWvHtVYZ5ALpTNKi0SYd9vtRK+4ko4Hf0Hf889xdM+30CIRIuecR+TMs/CvWDnuPdln\nQtATINbQRK01fFy9bdv0ZwfoTHbRmeymK9ntTjuvrmQP7YkO7BM8Y0dVVEJ6kJA3REgPEvYECXlC\nhDxBwvllXmc+pAcJeAIEND8erbQnOiHE2CTUC/gWLqLljg+Rbj1Kzx+eoOfPf6L7sUfpfuxRtNpa\nQqvXEFi1iuCqUyZ8x8eZoigKYU+IsCfE4sjoT3bKWTnimT56U3F603F60r0F03F6U3H6s/3EU3Fa\n+9tOeAIopCsaft1PwH359YAzreXnC9f58aoefJoPn+bFq3nxDb58Mn5fiCKT7pcTsLNZBra/TPy5\nZ+nf8jy5eMGIlfp6fIuX4F+yFN+SJfgWLEKPRmekNT8b9bZsi4FMgv5MP32ZAfoz/fRnBuhz3/sz\nAySyCRLZJIlckmQ26Uxnk2SszJSPqynaiKB3gj8SCEJOw6d68WgePKqOR/U4Ly0/XfBesI2uevCq\nuvOuedBVHa/qQVXm/iPZqrErohrrDNL9UhKKrhNaexqhtadhWxbpw4cY2LGdgR3bSb66h/4Xnqf/\nheeHtvd48DQ1421uxtM8D09jDD1ahydaj14XRQ3P3QuYqqIS9oYIe0M0T/KzOSt3XNAnsu58Lkk6\nlyadS5MafKVI5zLustTg8oFMgq5kN2krA93Fr6OCMhj+uqKhue+6qqOpmrtMQ1d0dHVofX7d0Hb6\nKMuO/6ymqGiKhqqozrSanx75PrQukFJJZJNoioamqGVxIhJzi4T6BCmqim/hInwLFxG97AoAst3d\nJPfvJbVvH+kjh0kfPUq6tZX0oYOj70PX0aP1aLW1aOGw+4qghcJokaFpNRhA9QdQ/X5Uv382qzkl\nmqoRVp2uoGKwbIvaqI9DbZ2ksmkyVsZ9ZclYGbJWlnTOfXfnM7nh2wxO5wqXDb1bVo6MlSVppchl\nBsjaWXJWjqw9Mw9RmA7VDXcn5DX3ZDE0raojThAjThrOS0HJT6MMniyc+RHrFQWVgvUFn3E+527j\nLj/+cydYr6goBftSFZVjhOjtSaIMlkMZmiY/P/Q++vITbe/WZcS6SlX07hfDML4EnAvYwJ2maW4e\na9u53v0yFbZtk+vpId16lGxnB9muLjJdXWS7u8h2dZHt6iTX28uYd84fher1ovj8qIGhoFf9fhSv\nF9XjRfF6UDxeZzuPx3l5vajuu7POg6K76zTdGY+vaSi6M63omrNcd5drekkuDOeV6mdt2zY5O0fW\nyg2+Z60sOTs7yrLC9+HT+XWWbZGzc+57wbRlYdk5crY1uE3OtvB4VAaSKawRy4c+lzvBOmef+WNN\n9BpJtcqH+8iwVxj9JDDaiSV/whi5jTpiXyM/q6s6t5x5HZFc/ZTKPmvdL4ZhvBY4yTTN8wzDOAX4\nPnBeMY8x1ymKgl5Xh1439jNNbcvCSiTI9cXJ9fW5r6FpK5nASiSwkkmsZBItmyYV78dKJsn2dGOn\nUrNVmcHQJx/+buCjqU7oK+67qqJoI+bddxR3nTryM4r7rhV8RgFVpS/kJ5HMOssUBVQF3OGXgycb\nRSmYdvbnLM5PO59zWmX5fTi/WPlpUJxjjrFeRcFXcGxnfyoo7jHdFi6ooHiH/dtR0Bocahkqzmd1\nZfi27oPNovUhursHBued44zcjzK4uvCzzpsyuLmFjYWNbdtYtoUNWIoT9ZbtPKjBAmzFnbbd7RWw\nbWtoHoucbYO7P2e5uz8KtrNtbMV2TzrO1sO2t/Oft9ztnZNPIKjT1592y+j+L182d592ft2wY7nr\nj1tuD69P4TKG7yu/j8HtRrwPLXfKn19n5ZdZmQl89vh3VVE52neMSGBqoX4ixe5+eT3wCwDTNLcb\nhhE1DKPGNM3eIh+nrCmqihYKoYVCTKQDe2Sr1bYsrGQSO5PGzmSw0hlnOp3Bzmaw0unBeSuTX5d2\nprNZyOWws1nsXA4757yTnx9cnhtzWyuThqSFbVtgOS97xPt0zUCXelk4UOoCjMM9fVW+0bpnlJEn\n4tE2mdjnFE3j5MVZrBXTKeToih3q84DnCubb3WWjhno0GkTXJzekLRaLTLlw5ez4eteWpBwTZY8I\nejuXGz4/YtrO5UbMF5wk3BYStu3Mw+CJ4/j1TosSy2mRYQO2Nbjctmxwlzvr7aHPW0P7Gdpnwbaj\nri8oEwxfV7is4N0eMX+i5c5kwfKCeXtw3l02uP3Q/NBH3X+X/Hp33h4xf9z2BccdqpO7XcG8PWJ+\n+PYFdRgxyRjbDC/XaNsWLh5nG9sefdsxthltH2MeY9imoxxnjPWKpuGpq6NmBvJspi+UnvBqRFfX\nwKR2Vg596jOhsuo9oq2Xnx3xX+JM1VkZY3quqKyf9cRUY50BaqY3pHHMdcX+S+owTss8rwU4UuRj\nCCGEGEOxQ/1R4AYAwzDOAA6bpll9p2AhhCiRooa6aZobgecMw9gIfBW4o5j7F0IIcWJF71M3TfNj\nxd6nEEKIiamK0UlCCFEtJNSFEKKCSKgLIUQFkVAXQogKUtL7qQshhCguaakLIUQFkVAXQogKIqEu\nhBAVREJdCCEqiIS6EEJUEAl1IYSoIBLqQghRQWb6IRlFM5kHWpc7wzC+AFyE8/P5P8Bm4H5Aw7k/\n/c2mac7Sg0pnj2EYAeBF4J+B31EddX4n8A9AFvgUsJUKrrdhGGHgPiAK+IDPAEeBb+L8bm81TfP9\npSthcRmGsQZ4EPiSaZpfNwxjEaP8fN3/Dj4CWMC3TdP83lSPWRYt9cIHWgPvxrmtb0UyDON1wBq3\nrm8Avgx8FviGaZoXAbuA20pYxJn0T0CnO13xdTYMowG4G7gQuAa4lsqv9wbANE3zdTjPXvgKzn/j\nd5qmeQFQaxjGG0tYvqIxDCMEfA2ngZJ33M/X3e5TwGXAJcDfGoYx5SdSl0WoM+KB1kDUMIya0hZp\nxvwReKs73Q2EcH7QD7nLfonzw68ohmGsAlYDv3YXXUKF1xmnTo+Zphk3TfOIaZrvpfLrfQxocKej\nOCfxZQV/eVdSnVPAVThPhMu7hON/vucAm03T7DFNMwH8Gbhgqgctl1Cfh/MQ67z8A60rjmmaOdM0\n+93ZdwO/AUIFf4K3AfNLUriZ9UXgroL5aqjzUiBoGMZDhmH8yTCM11Ph9TZN88fAYsMwduE0YD4K\ndBVsUjF1Nk0z64Z0odF+viPzbVr/BuUS6iPNxWcGF5VhGNfihPoHR6yquLobhvEu4CnTNF8dY5OK\nq7NLwWm1XofTLfGfzP1nY0+LYRg3AftN01wJXAr8cMQmFVfnExirrtP6NyiXUK+qB1obhnEl8Ang\njaZp9gB97kVEgAUM/3OuElwNXGsYxibgduCTVH6dAVqBjW6LbjcQB+IVXu8LgEcATNPcAgSAxoL1\nlVjnQqP9dz0y36b1b1AuoV41D7Q2DKMW+FfgGtM08xcNHwOud6evBx4uRdlmimmaf22a5nrTNM8F\nvosz+qWi6+x6FLjUMAzVvWgapvLrvQunDxnDMJbgnMi2G4Zxobv+OiqvzoVG+/k+Daw3DKPOHR10\nAfCnqR6gbG69axjGvwAX4wz5ucM9y1ccwzDeC3waeKVg8S04YecH9gG3mqaZmf3SzTzDMD4N7MVp\nzd1HhdfZMIz34XSzAXwOZ/hqxdbbDa3vA804Q3Y/iTOk8R6cRubTpmneNfYeyodhGGfiXCtaCmSA\nQ8A7gR8w4udrGMYNwN/jDOv8mmma/zXV45ZNqAshhBhfuXS/CCGEmAAJdSGEqCAS6kIIUUEk1IUQ\nooJIqAshRAWRUBdCiAoioS6EEBXk/wNoIngxHsxqFAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "tags": [] } } ] }, { "metadata": { "id": "v-PkqmeUgFFU", "colab_type": "text" }, "cell_type": "markdown", "source": [ "## Littlewood's rule and protection levels.\n", "\n", "We will now calculate the Littlewood rule for a continuous function. The idea is to calculate the limit $y_1$ such that:\n", "\n", "$$\n", "p_2 = p_1 P(D_1 \\ge y_1)\n", "$$\n", "\n", "Let's assume the following:\n", "\n", "- $p_1 = 300$\n", "- $p_2 = 100$\n", "- The first class demand distributes as a [gamma](http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm) distribution of parameter 5.\n", "- The plane has 1,000 seats (i.e. the scale of the distribution is 50).\n", "\n", "For this, we will use the \"sister\" package of Numpy, Scipy. Scipy comes with the useful \"stats\" subpackage, with includes a large number of statistical distributions. To use one, we will first import the package, create our distribution, and then plot the CDF and PDF." ] }, { "metadata": { "id": "QUjV0OJygBns", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "import scipy.stats as stats" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "HZ1mb1vLgHZb", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "fp = stats.gamma(a = 5, scale = 50) # Creates the distribution with the parameters we got." ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "69EGfrKngJLD", "colab_type": "code", "outputId": "70f88b9e-8eea-4528-ea9e-c20c9b95acb6", "colab": { "base_uri": "https://localhost:8080/", "height": 305 } }, "cell_type": "code", "source": [ "# Sequence over which to plot.\n", "t = np.arange(0, 1000, 0.01)\n", "\n", "# Create an overall figure\n", "plt.figure(1)\n", "\n", "# CDF plot\n", "plt.subplot(121)\n", "plt.plot(t, fp.cdf(t))\n", "plt.title('CDF')\n", "\n", "# PDF plot\n", "plt.subplot(122)\n", "plt.plot(t, fp.pdf(t))\n", "plt.title('PDF')\n", "\n", "# Minor adjustments to size and shape\n", "plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,\n", " wspace=0.35)\n", "\n", "\n", "plt.show()" ], "execution_count": 0, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEgCAYAAACXa1X+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUXPV16PtvVVfPc7d6ULcmJKQt\nCTA2o2RiA4YYO8FJ/Ix9k0dwSEjuDcEvdnydhHedZRsTkyx7+XJDhuVwn82L7cAjxgaDIxtiYixA\nDGIwYCG20NBSq+d5Hqqr6v1xTrVKTVd3dXcNp6r2Z9lL3eec36ldTZ/e9Zt9kUgEY4wxJln8mQ7A\nGGNMbrHEYowxJqkssRhjjEkqSyzGGGOSyhKLMcaYpLLEYowxJqkCmQ4g34mID/gz4A+AQpz/Jo8D\n/zfwUeAfgXb3HMCjwB2qOuyWfwrYAYwuuPWVqtqT6viN8QIRiQDHgDmcD8wjwO2q+uSC8yGgHPgF\n8BVVfc49fzNnnrVY/0tVv5GO95BLLLFk3t8CVwHXqWqHiJQDfwf8CPgW8JyqXgsgItXu9U+JyB5V\nnXbv8Req+t30h26Mp1ylqqcBROQK4DEREVXtiz3vfpi7AfihiNygqvvd8/PPmlkbawrLIBGpA/4U\n+D1V7QBQ1QngU8BXAV/s9ao6oqq3AmPAJ9McrjFZQ1WfBY4Cexc5F1HV7wH/A+eDmkkySyyZtQc4\nrapvxR5U1WlVfQwIxyn3GHB1qoMzJssVAjNLnH8UuFxEStMUT96wprDMqgNW0w8yClTHfP9VEfmr\nmO9nVPXCNUVmTBYTkQ8DzcCzS1w2ivPhutL9fq+IvLXgms+q6r4UhJjTLLFkVj/QuopyW4DemO+t\nj8UYp+8x2nnfBnxYVceXuH4LEASG3e+tjyVJLLFk1vNAk4hcpKqvRA+KSCHwJaBjYQERKQB+C/hy\nuoI0JkvMd94n6AbgKVWdFZFUxZSXrI8lg9whw18Fvi0i5wKISBlwL/AeYDL2enfE2L3AEPBv6Y3W\nmNwgIj4RuQH4DE4Hvkkyq7FkmKp+SUQGgUfd2kgY+CFwK/DbnGn3LQBK3XPXqepcpmI2JktFm8qq\ngTeBX1fVlzIcU07y2X4sxhhjksmawowxxiSVJRZjjDFJZYnFGGNMUlliMcYYk1SWWIwxxiRV2ocb\n9/WNxR2GVltbxtDQZLzTaeeleLwUC3grnqViaWio9C16IgfYs7Q6XooFvBVPsp4lT9VYAoGCTIdw\nFi/F46VYwFvxeCkWr/Daz8RL8XgpFvBWPMmKxSZIGpMgEbkbZ0XqCPBpVT0Yc+5a4C6cjaT2qeqd\nCZS5DviJqvrc72/EmQ0eBu5V1W+m5Y0Zk2SeqrEY41UiciWwXVX3ArcA9yy45B7gY8AVwAdFZPdS\nZUSkBGeX0C73+3LgC8C1OBu//Zm7X48xWccSizGJuQZ4BEBVDwO1IlIFICJbgUFVbVfVMLDPvT5u\nGZw1qv4RmHW/vxw46G7mNoWz3PsVaXlnxiSZJRZjEtMM9MV83+ceW+xcL7A+XhkR2QFc6O5iGO/+\n0XsYk3Wsj8WY1VlqhEy8c9Hjd+NsSb3a+wPOCJ6lOlsbGirjnssEL8XjpVjAW/EkI5aEEouInI+z\nqu7dqvoPC84t2mlpTI7p5EwNBaAFt39kkXOt7rHZRcrMADuBf3X3AFkvIj8HvrjIPZ5fKqClhqg2\nNFTS1ze2VPG08lI8XooFvBXPUrGsJOEs2xTmdir+PfBknEve0WmZ8Ksbkz2ewNkYChG5COhU1TEA\nVW0DqkRki4gEgOvd6xcrc1JVt6nqHlXdA3Sp6pXAC8ClIlIjIhU4z9PT6X2LxiRHIn0sM8Cv4XwC\nO8sSnZbG5BRVPQC8LCIHcD5M3SYiN4vIR91LbgUewEkGD6rqkcXKLHH/KeB24HHgp8AdqjqSundk\nTOos2xTmbig1F2frzsU6HLclJzSTKpFIhKmZECMTM4xOzDI+FWRqJsT07BxTsyGmZ+aYng0xPRti\nLhRmLhQmGAozNxdmLhxhbi4MPh9TM3NEwhHCkQiRCERw/g2HI0QiEcIR57Wi58763t0HaLHtgKLH\nIsScjJz1z1mqK4r40u9fRnV5UXJ/UAuo6u0LDr0Wc24/sDeBMgvPb4n5+iHgobVFmV7hSITv/ewo\nP/9FJ5ubKvmjj+ymrqok02GZDEt25711OKbQSmIJhyN09I1zsnuUrv4J5/8DE/QOTjI0NkNwLryq\nGPx+H4UBP4ECP4UFfgoKfPh8Pgr84PP58ft8+Hzg8/nw+334o1/7fPj80a/Bh2/+tyX6S+PzvfPX\nJ/bQwvPRb6sritnQUk1JkY1FSbefvnSax19sp6SoAG0f5u9/8AZ/9cmLKfDbgNN8ttYnMV6nZVzW\n4bg6y8UyMxtC24c5fHKQE52jnOwdZ2Y2dNY1PqCmspgNDeVUlRVRXVFEVXkRFaVFlBYXUFoUoLQ4\nQElRASXFAUoKCwgE/BQW+ChwE4nf70sonnSKxrJYNF76YJBrpmfn+OEzxykvCfCV/7qHB598m+cO\n9fDsG928/8KWTIdnMmhNiUVV20SkSkS2AKdxOi1vTEZgZnnjU0Fe0l4OHu7lSPswobDTUOTzwfr6\ncjY3VbKpqYKmujKaaktZV11KYcA+SZrkeP5QD1MzIX7rV86hqqyIj199Lgff6uXxF0/xvnetX7QG\navLDsolFRC4Gvg5sAYIicgPwKHBCVR/mTKcluJ2WKYrV4PRN6Klh/uOldl4/NjCfTDY3VbL7nFp2\nb6nj3JZqiou8s7CdyU3PHerG54P3ubWTmopiLtnZyPOHetBTw+zcXJvhCE2mJNJ5/zLO2kXxzi/a\naWmSKxKJ8MqRPn74zAnae8cB2NRYwZ7zmrl0ZyP11dZhatJnfCrI0Y4RtrVUU1tZPH/8fe9q4flD\nPRx8q9cSSx6z3s4scKxzhK8+8CpvnRzC54NLdzbyq5dsZFtrlTU3mIz45YkBIhG4YFv9Wcd3bKym\norSQV470ceMHd+C338+8ZInFw2ZmQ/xg/3F++lI7EeDiHQ38H1duZX19eaZDM3nurZPDAFyw9ewF\nmAv8ft69fR3PvN7F8c5Rzm2tzkR4JsMssXhUe+84//TwG/QMTdFUW8pnfucimqqKly9oTBoc7xyh\nKOBnY2PFO869a2s9z7zexeG2QUssecqGCHnQc4e6+cq3X6JnaIrrLtvIHX9wGedvW5fpsIwBYGpm\njo6+Cbasr1p0vsrOzbX4gMMnh9IfnPEEq7F4SCQS4ccvnOKhp45RWlzAp37jAi7a0ZDpsIw5S1vX\nKBFgW0vVoucrSgvZ2FTB0Y5RZoMhigpthGK+sRqLR0QiEb731DEeeuoYtZXFfP6mSyypGE863jUK\nwNY4iQVg1+Za5kJhjnbYcmf5yBKLRzzy9Al+8sIp1teX8fmbLqZlnXXQG2/q6JsAYGNT/FUNdmyo\nAeCYJZa8ZE1hHvDEwXYeO9BGY00pf/E776G6wjrpjXd19E9QVOhn3RJzp6K1meOdo+kKy3iI1Vgy\n7BdH+/n/nnybmooi/vtvv9uSivG0UDhM18AELfXlS85Rqa4opr6qmONdo/MrWZv8YYklg7oHJ/nf\njx2iMODn0zdcSENNaaZDMmZJvUNTzIUitCbQVHtOSzVjk0H6R6bTEJnxEkssGTIbDPEPP3iDqZkQ\nN39oJ5ubbRVe433R/pXWhnfOX1lo63prDstXllgy5KGnjtHZP8E1F21g7/nNyxcwxgM6+53Eksjg\nkmg/y4kuSyz5xhJLBhxqG+SnL59mfX0ZH7/aNtw02aPH3U+pub5s2Ws3NlbgA071eGPfHpM+lljS\nbGY2xH37DlPg9/GH1++2yWMmq/QOT1Hg91GfwPJCpcUBGmtLae8dtw78PGPDjdPsR8+1MTg6w6/v\n3cw56+NPMDPeIyJ3A3uACPBpVT0Yc+5a4C4gBOxT1TvjlRGRvcDXgCAwA9ykqn0iEgSejXnJa1T1\n7G1AM6xvaIr6qpKEtx7e2FjBS9pH/7B14OcTq7GkUffgJD954RR1VcVcv3dLpsMxKyAiVwLbVXUv\ncAtwz4JL7gE+BlwBfFBEdi9R5rPAJ1X1auA54I/c4yOqelXM/z2VVKZn5xidDNJQm/joxegkyhNd\nNlEyn1hiSaMHfvo2oXCE37lmu+3wmH2uAR4BUNXDQK2IVAGIyFZgUFXbVTUM7HOvX7SMqn5cVY+L\niA9oxdnW2/P63FpH4wqGxUdXPz5hM/DziiWWNHnr5BBvHB9g1+ZaWwMsOzUDfTHf97nHFjvXC6xf\nqoyIfAhQoAn4rnu+RETuF5FnReSzSX8Ha9Q7NAWwovlWm6KJxYYc5xXrY0mDSCTCD/YfB+BjV26z\nXR9zw1L/EeOdmz+uqj8REQH+Frgdp3/mczhJJgLsF5H9qvpSvBeprS0jEIhf821oSO7cqMlf9gBw\n7ubahO+9bl0FlWWFHO8cSXo8a+GlWMBb8SQjFkssafD6sQGOdozwnu3rllwR1nhaJ2dqKAAtQFec\nc63usdnFyojIR1X1YVWNiMj3gS8BqOo3oheKyJPABUDcxDLkDv1dTENDJX19yR3m29bh7BpZ7Pet\n6N6t68p569Qw7R1DlBRl/k9OKn42a+GleJaKZSUJx5rCUiwSifDI0yfwAR99/9ZMh2NW7wngBgAR\nuQjoVNUxAFVtA6pEZIuIBIDr3evjlfmSiLzbve/lgIrjfhHxufe4AjiUvre3vL7haFNY/MUnFxOd\npd89GD8RmtyS+Y8POe7wySFO9oxxyc5GNiSwDIbxJlU9ICIvi8gBIAzcJiI344zkehi4FXjAvfxB\nVT0CHFlYxj1/C/BPIjIHTOEMN+4VkXbgRffaR1X1xbS9wQQMjE5TXhJYca0jOku/s3+CLc1WY88H\nllhS7McvnALgw5dvynAkZq1U9fYFh16LObcf2JtAGdx+k/cucvwvkxBmygyNzbCueuULpba4s/Q7\n+63Gki+sKSyFTvWMcejEIDs31dhkSJPVJqfnmJ4NUZfAjPuFYmssJj9YYkmhx190aisfstqKyXJD\nY84clrrKlSeWyrIiqiuKLLHkEUssKTI2OcvBt3pZX1/GBVvrMx2OMWsyODYDQG3VyjruozY2VdI3\nPMVs0FOLCZgUscSSIs++0c1cKMKV7261eSsm6w2Orr7GAk5iiWAjw/KFJZYUiEQi/PwXHQQK/LzX\n9loxOWDIrbGsNrFsctcMs+aw/GCJJQXeOjlEz9AUl+5spKK0MNPhGLNmg6NrbwoD6BywxJIPLLGk\nwM9f6wTg6ve0ZjgSY5Ij2nlfu+YaizWF5QNLLEk2OT3HK0f6WV9fxrZWG2JscsPg2AzlJQGKV7kx\nXU1lMeUlAWsKyxOWWJLs5SO9zIXC7NndZJ32JidEIhEGR2eoW2UzGIDP56O5roy+4SlC4XASozNe\nZIklyZ4/5KwAe/l51mlvcsPUzBwzwdCqm8GimurKCIUj9I/YbpK5zhJLEg2NzfDWySG2tVataDMk\nY7xsZGIWgJqKojXdp7nOWdqlx4Yc5zxLLEn04uEeIsCe3VZbMbljZNxJLFXla6uxRBNL94AlllyX\n0CKUInI3sAdnA6JPq+rBmHO3Ab8LhICXVPUzqQg0G7x4uAe/z8eluxozHYoxSROtsVSXr63G0hRN\nLO5OlCZ3LVtjEZErge2quhdnue97Ys5VAX8OvE9VfwXYLSJ7UhWslw2OTnOia4ydm2uoKlvbA2iM\nlyQrsTTWOs3D1hSW+xJpCrsGeARAVQ8DtW5CAWeHvFmgwt2cqAwYTEWgXvfq2/0Atp+9yTkjE87k\nyOo19rEUFxZQX1Vsy7rkgUQSSzPQF/N9n3sMVZ0G7gCOAyeBF9wNjvLOK0ecH9F7tltiMblldDw5\nNRZwmsOGxmaYmbXFKHPZajb6mp+c4dZc/gewAxgF/lNELlTV1+IVrq0tIxCIP8lqJfsqp0Mi8YxO\nzKLtw+zYVMOOresyGks6eSkeL8WSa6JNYVVJSixvtg3RMzQ5Pxvf5J5EEksnbg3F1QJ0uV/vAo6r\naj+AiDwNXEzMznoLDQ3FrwY3NFTS1zeWQEjpkWg8z77RRTgc4YJz6lIWf7b+bNJhqVgs4azdyMQs\nxUUFK96SeDHNtW4H/qAlllyWSFPYE8ANACJyEdCpqtGnuA3YJSLRSRuXAG8nO0ivizaDWf+KyUUj\nE7NJaQaDMyPDrAM/ty37EURVD4jIyyJyAAgDt4nIzcCIqj4sIl8DfiYic8ABVX06tSF7S3AuzJtt\nQzTXlbG+vjzT4ZgUWmbY/bXAXTjD7vep6p3xyojIXuBrQBCYAW5S1T4RuRH4DM5zdq+qfjN9725x\n4XCEsclZmmqrk3K/5vozNRaTuxKq26rq7QsOvRZz7p+Bf05mUNnk6OlhZoIh2yUyx8UOuxeRXcC3\ngL0xl9wDXAd0AD8Xke8DDXHKfBb4pKoeF5EvAn8kIn8HfAG4DGek5UEReVhVMzrKcmxylkgkOR33\nAOuqSijw++getLksucxm3q/RGyec5/6CrXUZjsSkWNxh9yKyFRhU1XZVDQP73OsXLaOqH3eTig9o\nBU4DlwMHVXVEVaeAZ4Er0vsW3+nMHJa1zbqP8vt9NNaW0jM4SSQSSco9jfdYYlmjXx4foDDgZ8fG\nmkyHYlIr7rD7Rc71AuuXKiMiHwIUaAK+u8Q9Mmp+RNga57DEaq4rY3JmjrGpYNLuabxl7cM88tjQ\n2Ayn+yY4f2sdRavcp8JkraX2RIh3bv64qv5ERAT4W+B2nIEwid4fSM/Q/fCJIQA2NFet+X7R8ls3\n1PDq2/1Mh2BbhkbteW20oJfiSUYslljW4JfHBwC44BzrX8kDSw27X3iu1T02u1gZEfmoqj6sqhG3\nL+ZLwIFF7vH8UgGlY+j+6e4RAPzh8JruFxtPZYnzZ0eP99NYmf7lj7w0VB68FU+yhu5bU9gaRPtX\nzrf+lXwQd9i9qrYBVSKyxV3a6Hr3+nhlviQi73bvezlOk9gLwKUiUiMiFTj9KxkfYZnMyZFRTdE1\nw2wxypxlNZZVCocjvHlikPqqkvnlwE3uWm7YPXAr8IB7+YPu0kZHFpZxz98C/JM7RH8KZ7jxlIjc\nDjyOMzT5DlUdSdsbjGN80ukHSebCqjaXJfdZYlmlU71jTM7McbE02BbEeWKZYff7OXv4cbwyqOpL\nwHsXOf4Q8NDaI02e8WknsVSUFibtntXlRRQXFdCzRFOeyW7WFLZKb50cBmDn5toMR2JM6kxMBQkU\n+CkqTN6fCp/PR1NtKb1DU4RtyHFOssSySm+dckbL7NxkicXkrvGpIBWlgaTXypvrypidCzM8NpPU\n+xpvsMSyCqFwmLdPD9NUW0ptZXImjhnjReNTc0ltBotqrLV+llxmiWUVTvWMMzUTQqy2YnJYKBxm\namaO8pLkJ5boyDDbpjg3WWJZhTPNYDbb3uSuiak5ILkd91HNNjIsp1liWQU95XTcW43F5LJxd8mV\n8hQkFhtynNsssaxQKBzmSPswTXVl1r9iclo0saSixlJRWkh5ScAmSeYoSywrdKpnnOnZkDWDmZw3\nkcLEAk6tpW94ilA4nJL7m8yxxLJCR087k6G3b0jOxkfGeNWZprDUzKNuqi0jFI4wMDKdkvubzLHE\nskJHO5zEcm6rJRaT2yamU9d5D9BUZ2uG5SpLLCt0tGOEqrJCGmpKMx2KMSmVyj4WcGosYNsU5yJL\nLCswODrN0NgM21qrbX0wk/NSnViiQ457bZvinGOJZQXmm8Gsf8XkgYkUDjcGaJyfJGk1llxjiWUF\noh331r9i8sF8531JajrvS4sDVJcX2VyWHGSJZQWOdoxQ4Pexpdk724gakyrj00FKiwMU+FP3Z6Kp\ntpSB0WmCczbkOJdYYknQzGyIUz3jbGmupHCJfcaNyRXRlY1TqamujEgE+oatnyWXWGJJUFv3KOFI\nhG3WDGbyQCQSYWIqmLKO+6j5pV2snyWnWGJJkM1fMflkJhhiLhRJWcd9VHSV4x4bGZZTLLEk6ETX\nGABbW6oyHIkxqZfqocZRVmPJTbbnfYJOdI1SXV5kC0/mMRG5G9gDRIBPq+rBmHPXAncBIWCfqt4Z\nr4yIbATuAwqBIPC7qtotIkHg2ZiXvEZVQ2l4a+8wv2R+CvZiidVYE62xWGLJJZZYEjAyPsPQ2AwX\nbqu3iZF5SkSuBLar6l4R2QV8C9gbc8k9wHVAB/BzEfk+0BCnzF8D96rqv4nIbcBngb8ARlT1qrS9\nqSWMT6enxlJUWEB9VbEt65JjrCksASe6nWawc9ZbM1geuwZ4BEBVDwO1IlIFICJbgUFVbVfVMLDP\nvT5emT8Bvu/etw+oT+cbSUSqJ0fGaqwtY2hshpnZjFTOTApYYklAW9coAFssseSzZpwkENXnHlvs\nXC+wPl4ZVZ1Q1ZCIFAC3Afe750tE5H4ReVZEPpuKN5GodPWxQMxuktbPkjOsKSwBbW6NZct6mxhp\n5i3VJhrv3PxxN6l8B/hPVX3SPfw54Ls4/TH7RWS/qr4U70Vqa8sILDGnqqFh9b+vEXdSZGtz1Zru\nk0g8WzfW8rNXO5gOrS3mZMSSKV6KJxmxWGJZRiQSoa1rlPqqEqrKijIdjsmcTs7UUABagK4451rd\nY7NLlLkPeFtV74ieVNVvRL8WkSeBC4C4iWVoiU/4DQ2V9PWNxX83y+jpHwdgbja4pvskEk95kZPE\njrQNsKMl9X9g1/qzSTYvxbNULCtJONYUtoyhsRlGJ4NWWzFPADcAiMhFQKeqjgGoahtQJSJbRCQA\nXO9ev2gZEbkRmFXVL0ZvLo77RcTn3uMK4FD63t7Z5nePTPGoMLCmsFxkNZZlnHD7V6zjPr+p6gER\neVlEDgBh4DYRuRlnJNfDwK3AA+7lD6rqEeDIwjLu+dtw+lOecr9/U1X/RETagRfdax9V1RfT8uYW\nMe4ON05H5/266hL8Pp9NkswhlliWMd+/YgtP5j1VvX3Boddizu3n7OHH8cqgqu+Nc/+/XGuMyTI+\nFaTA76OkKPXr4gUK/KyrKbEaSw5JKLEsMzFsI84ntSLgFVX941QEminRGoslFpNPouuEpWveVlNt\nGW8cH2ByOkhZGprfTGot28cSOzEMuAVnIlisrwNfV9XLgJCIbEp+mJnhdNyP0VRbar/sJq9MTKd+\nAcpYTXXuDHybKJkTEum8X2pimB94H/Coe/42VT2VoljTrmdwksmZOTZbbcXkkXA4wuT0XFr6V6Ka\nap0O/G5b2iUnJJJYlpoY1gCMAXeLyDMi8jdJji+jTnQ6KxpvbrLEYvLHxHSQCOmZHBk1X2OxxJIT\nVtN571vwdSvwd0Ab8O8i8uuq+u/xCqdyUleyPfFyBwDn72j0RFxeiCGWl+LxUizZ7sys+/SN7Wl2\nayy91hSWExL5zVlqYlg/cFJVj8H8pK7zgLiJJZWTupItWmOpKi7IeFxe+9l4KZ5kTeoyjok0DjWO\nqqsqIVDgs6awHJFIU9hSE8PmgOMist299mJAUxFoJhzvHKG6oojqcptxb/LHeBonR0b5/T4aa8vo\nGZoiEomk7XVNaiybWFT1ABCd5HUP7sQwEfmoe8lngPvc8yPAYymLNo3Gp4L0DU2xqdE+8Zr8Mp7G\nlY1jNdWWMjUzx9hkMK2va5IvoUbUZSaGHQV+JZlBeUF7r7NW0qamigxHYkx6pXNl41ixu0lWWStB\nVrO1wuJo73Ha7Dc2WmIx+WUiTZt8LdRU64wMs36W7GeJJY5T8zUWawoz+SVzTWE2MixXWGKJ41TP\nOCVFBfN7chuTLyYy3BRmNZbsZ4llEcG5MF0DE2xZX4Xfb3vcm/wyX2MpSe8atTUVRRQXFtgqxznA\nEssiOvsnCIUjnNNanelQjEm78ak5SosLCBSk98+Dz+ejqbaU3uFJwjbkOKtZYlnEqV6n435riyUW\nk38mpoOUZ2jR1ca6MmaDYYbHZjLy+iY5LLEsor3H6bjfajUWk4fGp9K7snGsZlvlOCdYYlnEqd5x\nfD7YZKsamzwzEwwRnAunfURYVHRkmC1Gmd0ssSwQiURo7x2nqbaMkiLbYNPkl0yNCIuaTyy2m2RW\ns8SywNDYDFMzc2xoKM90KMakXSbWCYt1Zvl8awrLZpZYFujsnwCgZZ0lFpN/zkyOzExtvaK0kLLi\ngNVYspy19Sxwus9JLK0NtpSLOZuI3A3sASLAp1X1YMy5a4G7gBCwT1XvjFdGRDYC9wGFQBD4XVXt\nFpEbcRZ1DQP3quo30/fuHJlaJyzK5/PRVFfGqZ4xwuGIzSPLUlZjWSBaY2m1GouJISJXAttVdS9w\nC85K37HuAT4GXAF8UER2L1Hmr3ESx5XAw8BnRaQc+AJwLXAV8GciUpfit/UOE9POXiyZSizgNIeF\nwhH6R6czFoNZG0ssC3T0T1Dg99FYa0u5mLNcAzwCoKqHgVoRqQIQka3AoKq2q2oY2OdeH6/MnwDf\nd+/bB9QDlwMHVXVEVaeAZ3GSVFplusYCMbtJ2siwrGVNYTHCkQidAxM015elfdax8bxm4OWY7/vc\nY6Puv30x53qBbcC6xcqo6hEAESkAbgO+HOce65cKKBXbfIfdncc3tFQnfffNRO937uY6eOYE47Ph\nlO0A6rWdRb0UTzJiscQSY3BkmpnZkDWDmUQs1fgf79z8cTepfAf4T1V9UkT+zxXcH0jNNt99g05T\ncHA6mNStp1cST2nAeevH2odSsv21l7bVBm/Fk6xtvu1jeYwO618x8XXi1CqiWoCuOOda3WNLlbkP\neFtV71jmHmnlhaYwmySZ/SyxxDgz1NhGhJl3eAK4AUBELgI6VXUMQFXbgCoR2SIiAeB69/pFy7ij\nv2ZV9Ysx938BuFREakSkAqd/5en0vLUzJqaCFPh9lBTFb2JLtbKSAFVlhTbkOItZU1iMM0ONrcZi\nzqaqB0TkZRE5gDMc+DYRuRkYUdWHgVuBB9zLH3T7UY4sLOOevw0oEZGn3O/fVNU/EZHbgcdxhibf\noaojaXlzMcangpSXBPD5MjvMt6mujKMdI8yFwtbfmYUsscTo7J8gUOC3zb3MolT19gWHXos5tx/Y\nm0AZVPW9ce7/EPDQGsNck/HIlFvAAAAefklEQVSpoCf2m2+qLePt0yP0DU+xvt4+6GUb+yjgCkci\ndA1M0FJfZpOyTF4KhyNMTs9ltH8lypZ2yW6WWFz9w1PMzoVpsWYwk6cmZ+aIkNmO+6hoB75tU5yd\nLLG4OvpsRJjJbxPz64RlPrE019kqx9nMEovrzFBjGxFm8pMXhhpHNdWV4vOdGalpsoslFtf8UGNr\nCjN5ykuJpTBQQGNNKZ39E0QikUyHY1bIEovrdN8ERYV+1lWXZDoUYzLCS4kFnK0rJqbnGJmYzXQo\nZoUssQChcJjuwQla6svxZ3j8vjGZMt/HkqFNvhaKziez5rDsY4kF6B2aYi4UsY57k9fGp6M1Fm9M\nb4tuttdhiSXrWGLhzIgw618x+Wx8ytmLxQujwuDMQBqrsWQfSyzEbu5lI8JM/vJaH0tzXRl+n89q\nLFnIEgu2qrExcKaPxSuJpTDgp7G2lM4+GxmWbSyx4NRYSooKqKsqznQoxmTM+FSQ4qICTy362Lqu\nnMmZOYbHbWRYNvHOb1CGzIXCdA9O0rKuPOMruhqTSRPTQSo8MiIsKtqBb/0s2SXvE0vP4CShsI0I\nM2Z8KuiZZrAoGxmWnfI+sVj/ijEQnAsxGwx7ZqhxVOt8jWU8w5GYlUjot0hE7gb24GxA9GlVPbjI\nNX8D7FXVq5IaYYrZUGNjvDfUOKrJHRnW2W+LUWaTZWssInIlsF1V9wK3APcscs1u4P3JDy/1bKix\nMd4bahxVGPDTVFdKh60ZllUSaQq7BngEQFUPA7UiUrXgmq8Dn09ybGnR0T9BWXGAmorM75pnTKZ4\nNbGA088yZSPDskoiiaUZ6Iv5vs89BoC77/fPgbZkBpYOwbkQvUNTtDTYiDCT37y0F8tCGxqc1oT2\nXutnyRar6amb/wssInXA7wPXAq2JFK6tLSMQKIh7vqGhchUhrc6JzhHCkQjbNtTEfd10xrMcL8UC\n3oonHbEs1dcoItcCdwEhYJ+q3rlUGRH5U5yafq2qjrvHgsCzMS95jaqGUv7GiKmxeGy4McDGxmhi\nGeNd2+ozHI1JRCKJpZOYGgrQAnS5X38AaACeBoqBbSJyt6r+WbybDS2xI1xDQyV9fWMJhJQcvzzS\nC0B9RdGir5vueJbipVjAW/EsFUuyEk5sX6OI7AK+BeyNueQe4DqgA/i5iHwf59l4RxkR+STQhPNs\nxRrJ1OCXcQ/XWDY1Wo0l2yTSFPYEcAOAiFwEdKrqGICqPqSqu1V1D/BR4JWlkorX2FBjswJx+xpF\nZCswqKrtqhoG9rnXxyvzsKp+HqcW4wle7mOpry6htDhgiSWLLFtjUdUDIvKyiBwAwsBtbr/KiKo+\nnOoAU+nMUGMbEWaW1Qy8HPN9tK9xlHf2Q/YC24B1i5VR1SNxXqNERO4HNgPfV9X/uVRAyWxWnnNT\n3KbWGhpS9EFrLbXHra3VHD4xQGV1KSVFa59r46VmXPBWPMmIJaH/Qqp6+4JDry1yTRtw1ZojSqPO\n/gkqSgupKvPepzTjeUuN9oh3brkRIp8DvotTk9kvIvtV9aV4FyezWXlweAqA2akZ+vrCCZdL1Fqb\nTptrSjkUgdcO97C1ZeGg1PTGkmxeiidZzcremmabRjPBEH3DU+zYWGMjwkwiluprXHiu1T02u0SZ\nd1DVb0S/FpEngQuAuIklmcangvh9PkqLvfknYWPTmQ78tSYWk3p5u6RL18AEEWzGvUnYUn2NbUCV\niGwRkQBwvXt93DILieN+EfG597gCOJTi9zRvfCpIeWnAsx+yNloHflbx5seTNIj2r1jHvUlEAn2N\ntwIPuJc/6PajHFlYBkBEPg/8Kk5t5sci8pyq/oWItAMvutc+qqovpuv9jU8FqfRwk3DrunJ8Pjhl\niSUr5G1i6bQRYWaFluprVNX9nD38OF4ZVPUrwFcWOf6XSQhzxcKRCBPTQZrryzLx8gkpKixgfX05\np3vHCUci+D1aszKOvG0Kmx9qbCPCTJ6bmpkjEvHm5MhYGxsrmJ4N0T8ynelQzDLyN7H0TVBVXuTJ\ncfvGpNOZyZHebsCY72fp8cYIKhNfXiaW6dk5BkanrRnMGLw9OTJWdAb+yR7rZ/G6vEws0b0dLLEY\nc2YBSq8nls3NzjyKk91WY/G6vEwsHX3OJx4bamyMt9cJi1VZVsS66hJOdI3a3iwel5+JxUaEGTMv\nuntkpccTC8CW5krGp4IMWAe+p+VlYrGhxsacMT7lbKDl9aYwgHPWO7Pu26w5zNPyMrF09E9QW1lM\nmceHVxqTDtEaSzYkli1uYjnRNZrhSMxS8i6xTE7PMTQ2Q4vVVowBYkaFlXl/e+4tbge+JRZvy7vE\nYs1gxpxtfNJpCisv8fY8FoDS4gDNdWWc7BkjbB34npV3ieV0vzMirNVGhBkDOE1hpcUBAgXZ8efg\nnPWVTM2E6BmMv22Ayazs+E1KoujikxtsKRdjAKfzvsLjs+5jRftZ2rqsA9+r8jCxuHNY6q3GYkwk\nEmF8ao6KUu/3r0Sd0+x24HdbP4tX5V1i6eyfoKGmhOKi+Fu6GpMvZoIh5kLhrBgRFrWpqYICv4/j\nnZZYvCqvEsvoxCyjk0Fa11kzmDGQPeuExSoqLGBTUwUnu8eYDYYyHY5ZRF4lljNL5VszmDGQnYkF\n4NzWGkLhiA079qj8Sixu/4oNNTbGcWYOS3Yllu0bqgE42jGS4UjMYvIrsdjmXsacZXwyS2ssbmJ5\n+7QlFi/Ku8Ti9/lorvPuFqzGpFO0xpINC1DGqqkoZl11Ccc6RmyipAdlz+D1NYpEInT0TdBUV0ph\nIK/yqUkSEbkb2ANEgE+r6sGYc9cCdwEhYJ+q3rlUGRH5U+DrQK2qjrvHbgQ+A4SBe1X1m6l+T9my\nZP5itm+o5rlDPXQPTNoSTR6TN39hh8ZmmJqZs/4VsyoiciWwXVX3ArcA9yy45B7gY8AVwAdFZHe8\nMiLySaAJ6Iy5fznwBeBa4Crgz0SkLqVviuytsQCcu6EGsH4WL8qbxGL9K2aNrgEeAVDVw0CtiFQB\niMhWYFBV21U1DOxzr49X5mFV/TxOLSbqcuCgqo6o6hTwLE6SSqmsrrG0RvtZhjMciVkofxJLny0+\nadakGeiL+b7PPbbYuV5gfbwyqrrYWiTx7pFS2TrcGJwdYMuKAxxpt8TiNXnTx9Jhi0+a5PKt4txS\nZVZ8bW1tGYFA/BUkGhoql32R6WCY0uICWtZXryC01UkknpW64Nx1vHCom0hBAY0rGJSTiljWwkvx\nJCOW/EksfRMECnw01pZmOhSTnTo5U0MBaAG64pxrdY/NLlFmufu3As8vFdDQUPzVfRsaKunrW36R\nxuGxacpLChO6di0SjWeltq6v5IVD3Tzzajvve1dLRmNZLS/Fs1QsK0k4edEUFo5E6ByYYH19OQX+\nvHjLJvmeAG4AEJGLgM5ok5aqtgFVIrJFRALA9e71ccss4gXgUhGpEZEKnP6Vp1P4fgCnKSwb+1ei\ndm2uBeCtk0MZjsTEyosaS//wFLPBsPWvmFVT1QMi8rKIHMAZDnybiNwMjKjqw8CtwAPu5Q+q6hHg\nyMIyACLyeeBXcWooPxaR51T1L0TkduBxnE79O1Q1pcOdZoMhZoPZtQDlQq3ryqksK+TwySEikQg+\n30paG02q5EViae91+lc2NtqIMLN6qnr7gkOvxZzbD+xNoAyq+hXgK4scfwh4aO2RJmbMnXVflWXL\nucTy+Xzs2lzLi4d76R6cZL1th+EJedEudKrHEosxC426WxJXZsFe90vZac1hnpMXiWW+xtLknZEX\nxmTa6ISTWKrKszuxRPtZ3rTE4hl5k1iqyouozvIHyJhkOlNjyd6mMIDGmlLqq0o43DZEKBzOdDiG\nBPtYllkj6Wrgb3DWSFLgD93Zx54wOR1kYHSa885J+eoYxmSVaB9Ltn/g8vl8vGtbPT97tYNjHaPs\n2FiT6ZDy3rI1lgTWSLoXuEFVrwAqgQ8lPco1sI57YxYXbQrL9j4WgAu21QPw2rH+DEdiILGmsLhr\nJLkuVtXT7td9QH1yQ1ybU5ZYjFlUtCmsKgcSy67NtRQG/Lx+bCDToRgSSyxLrZGEqo4CiMh64IM4\nC/B5RrTGsskSizFnGZvIjT4WgOLCAnZuqqWjb4KBkelMh5P3VjOP5R0zkESkEXgM+BNVXfIjQzLW\nN1qJ7sFJCgN+LpAmCgpWPlYh19bwSSYvxeOlWLLF6GSQkqICigrjP4/Z5F3b6nnj+ACvHx/g6ve0\nZjqcvJZIYllqjSTcZrEfA59X1SeWu1ky1jdKVCgcpq1rjNZ15QwOTqy4fLas4ZMJXoonWesb5ZvR\nydmsH2oc613b6vnX/4DXjvZbYsmwRD7CL7fe0deBu1X1JymIb026ByaZC4Wtf8WYBcKRCOOTwZzo\nX4lqqCmltaGcN9uGmJqZy3Q4eW3ZGstSayThrGv0SWC7iPyhW+R+Vb03VQGvhI0IM2Zxk9NzhMKR\nnOhfiXWpNPLIMyd47Wg/e85rXr6ASYmE+liWWiMJKE5eOMl1ssepWG1qssRiTKxcmXW/0MU7ncTy\nkvZZYsmgnJ5539Y1hg/Y3Gzt7MbEGsuRdcIWal1XTsu6ct44PsD0rDWHZUrOJpZwJEJbzxjr15VT\nUpQXizgbk7DRHJl1v5hLpIHgXNjmtGRQziaW7oFJZmZDbLHaijHvMJpDc1gWumRnIwAHD/dmOJL8\nlbOJpa17FMASizGLGB6fAaCmwrNdpKvWuq6c1oZyXjvWz/hUMNPh5KXcTSxdTsf9lvVVy1xpTP4Z\nHnMTS2XuJRafz8cV569nLhThhTd7Mh1OXsrdxNI9ht/ns6VcjFnEfI0lB/tYAPae14Tf5+OZN7qW\nv9gkXU4mllA4zKmeMVobynNmuQpjkml4fJbykkDOPh/VFcVcsLWOk91jnHbns5n0ycnhUl39k8zO\nha1/xSTVMvsSXQvchbMv0T5VvTNeGRHZCHwHKMBZHukmVZ0RkSDwbMxLXqOqoVS8l+HxmZzsX4l1\nxQXree3YAM+80cVvX7M90+HklZyssRzvcjvurX/FJEkC+xLdA3wMuAL4oIjsXqLMl4F/VNX3AUeB\nP3CPj6jqVTH/T0lSmQ2GmJieo6YiN5vBoi48dx1VZYU883oXM7Mp+VGaOHIysbx9ehiA7a3VGY7E\n5JC4+xKJyFZgUFXb3d1T97nXxytzFfCoe9/HgGvT+D4Ydoca53qNpTDg56r3tDI5M8dzh7ozHU5e\nycnEcrRjlNLiAlrWlWc6FJM7ltqXaOG5XmD9EmXKVXVmwbUAJSJyv4g8KyKfTXL883J5RNhCV767\nlQK/jydfPk0kEsl0OHkj5/pYRidn6Rmc5Pxz6vD737F1jDHJstQvV7xzix2PPfY54Ls4/TH7RWS/\nqr4U70VWu7fRWx1OU/GG5qq0biuQiS0MGhoqueLCFva/2kHX8AwX7mjIWCxL8VI8yYgl5xLLsdMj\nAJy7wZrBTFIttS/RwnOt7rHZOGXGRaRUVadirkVVvxG9UESeBC4A4iaW1e5tdKrTeUYCRNK2p04m\n9+953/nN7H+1gweeeIuW2hJP7SUEubm3Uc41hb3d4SYW618xyRV3XyJVbQOqRGSLiASA693r45X5\nKU5HP+6/PxHH/SLic+9xBXAoFW9kviksx/tYora1VrNrcy2HTgxyzE2qJrVyLrEc7RjB7/OxtcVG\nhJnkUdUDQHRfontw9yUSkY+6l9wKPAA8DTyoqkcWK+Ne+0Xg90TkaaAO+BdVVaAdeBFnyPE+VX0x\nFe+lf9TZE76+uiQVt/ekj7x3CwCPPduW0TjyRU41hc0GQ7R1jbGxscJWNDZJt9S+RKq6H9ibQBlU\ntQv41UWO/2USwlzWwMgUgQJfzu3FshTZVMOODdW8fmyAo+3DVJfk5sRQr8ipGsvbHSPMhcLs3FyT\n6VCM8ayBkWnqq0rw+/JncIvP5+M3f+UcAO770SEbIZZiOZVYDrcNAbBrc12GIzHGm2aCIUYng3nV\nDBa1a0sd79pWz+tH+3ntqO3Vkkq5lVhODlLg97Fjo3XcG7OYgRGnf2VdHiYWgE9cfS5+v48Hf3aU\nuVA40+HkrJxJLJPTQdq6x9jaUmX9K8bEMRDtuK/Kz8TSsq6cD+/dQs/gJI+/eCrT4eSsnEksb50a\nJhKBXZtrMx2KMZ7VP19jKc1wJJnzux/aSXV5ET98po2ugYlMh5OTciaxvHHcaTM97xzrXzEmnv7h\nKSC/hhovVFFWxE3XCXOhMPfte4tw2Dryky0nEks4EuEXR/upKC1kW4v1rxgTT/egM1u/qa4sw5Fk\n1kU7GrhkZyNHO0b40YG2TIeTc3IisbR1jTEyPsuF59bb+mDGLKF7cJLS4gBVZYWZDiXjPnmdUF9V\nwg+fOcHhtsFMh5NTciKx/OKos4Dsu89tyHAkxnhXKBymd2iK5royfHk0hyWeitJC/vi3zsPv9/HP\njx6abyY0a5f1iSUSifCy9hEo8HPeOdZxb0w8/cPThMIRmvO8GSzWtpZqfufa7YxOBrn7e68xMR3M\ndEg5IesTy6mecboGJrnw3HobZmzMErrc/pXmeksssT5w0Qauu2wjXQOT3PPQ60zPzmU6pKyX9Ykl\nujPce89rXuZKY/JbZ78ztLbFEss7fPzqc7lsVyNvnx7hfz74GpPTllzWIqsTSygc5oU3eygvCXDB\ntvpMh2OMp53sdvbZ2NzknU2lvMLv8/FHH9nNnt1NHO0Y4WsPvMrQ2MzyBc2isjqxvHKkn5GJWfac\n10ygIKvfijEpd7JnjPKSQF7PYVlKgd/PH16/m/dfuJ6TPWN8+V8O2v4tq5TVf43/46V2AK65eEOG\nIzHG2yan5+gdmmJTU6WNCFuC3+/j9z60k9/+wLmMTszyt999hR8daCMUtnXFViJrE8vRjhGOnh7h\nXdvqbZSLMcs43uV88t7SbM1gy/H5fHzwsk189r+8m8qyQn6w/zh3fecVTnSNZjq0rJGViSUSifDQ\nU8cA+LU9mzMcjTHe99bJYQB22lp6CTtvSx1fvuVyLt/dxImuUe78l5f434+9Sc/QZKZD87ysHJ/7\nypE+jrQPc+G2enZstE29jFnO4ZNDFPh9bN9gSx6tREVpIf/tN87j/Re28OCTb/PcoW6eP9TNxdLA\nr166kXNbq61pcRFZl1hGJ2f5zuNKoMDPJz5wbqbDMcbzBkenOdE1ys5NNTbXa5V2ba7lCzdfykva\ny49fOMVL2sdL2kdjTSl7z2/mPdvXsbGxwpKMK6HfMhG5G9gDRIBPq+rBmHPXAncBIWCfqt6ZikDB\n2dP+H37wBqOTQT5x9bmsry9P1UsZ8w6reQ4WKyMiG4HvAAVAF3CTqs6IyI3AZ4AwcK+qfjMZcT//\nZg8Al+1qSsbt8pbf7+OyXU1curORt04N88zrXbx8pJcfPnOCHz5zguqKIs4/p47tG2o4Z30VrevK\n83btwmUTi4hcCWxX1b0isgv4FrA35pJ7gOuADuDnIvJ9VX0z2YEOjEzzjUd/ybGOUS7b1cgHL9uY\n7JcwJq7VPAdAQ5wyXwb+UVW/JyJ3AX8gIt8GvgBcBswCB0XkYVVd0+qI07Nz/MfBdoqLCrhkZ+Na\nbmVcPp+PXZtr2bW5lptmd/CLt/t54/gAvzwxyLNvdPPsG86k7eLCAlrWldFcV0ZTnfNvXVUJNRVF\nVJcXUxjIyi7uhCRSY7kGeARAVQ+LSK2IVKnqqIhsBQZVtR1ARPa5168qsYxPzjI0NsNcKMxcKExw\nLkzP0BSHTgzw/KEeZufC7Dmvid//8C78VuU06bWa56BhsTLAVcAfu/d9DPgcoMBBVR1x7/EscIV7\nfsWGx2Zo7x3n4f3HGZmY5fr3bqGi1FY0TraSogB7zmtmz3nNhCMRTveOc7xrlOOdo5zoGqW9d5wT\nXWOLlq0oLaSqvIjqimICfh9lJQFKiwOUFhdQUhSgsMBPYcD5f6DAR2GgYP5Ygd+Hz+fUovw+H/7o\n9z7ne5/fh9/9Pvq1D+eaqGiznc8H0cOVSdoALpHE0gy8HPN9n3ts1P23L+ZcL7BtNYE89WoH335c\n456vqyrmo+/bynvPb7Z2TJMJq3kO1sUpU66qMzHXro9zj/WrCfQnL5zi3352dP77nZtq+I0rtqzm\nVmYF/D4fm5oq2dRUyVXvbgWc1UEGRqbpHpyiZ3CSofEZRsZnGB6fZdj9untgAq/sNVZaHOCrt+6l\nvGRtH0JW05O31F/1Zf/i19aWEQgUvOP4RbubaesZBx8Eopm6wE9DbRnbN9Ww+5x6CjLQXtnQ4J1x\n/16KBbwVTwZiWc1zsNjxlVx7lnjP0qXnr6dneJqS4gJ2bKrl6os3eqLZJV9/X5qbqjlvifORSITp\n2RCT00HGp4JMTAWZng0xNxdmdi7EbNBpvQnOhQjOhZkNhgiFI4TDEcKR6L8s+N79N+br+dwVgUgE\nIkRw/zcfR311KRtba9f8tzaRxNKJ82kqqgWnw3Gxc63usbiG4owBryou4M9vuoS+vsWrjYMD4wmE\nmlwNDZVx40k3L8UC3opnqViS+AdkNc/BbJwy4yJSqqpTMdcudo/nlwoo3rNUX17I53734vmfyfBQ\n5vd1z5bfl0xoaKhkfNTZC6aswEdZRVFGY0nGs5TIx5gngBsAROQioFNVxwBUtQ2oEpEtIhIArnev\nNybXrOY5iFfmp8DH3Pt+DPgJ8AJwqYjUiEgFTv/K02l6b8Yk1bI1FlU9ICIvi8gBnGGQt4nIzcCI\nqj4M3Ao84F7+oKoeSVm0xmTIKp+DIwvLuOe/CHxbRP4bcBL4F1UNisjtwOM4rRN3RDvyjck2vkgk\nvb1GfX1jcV/Qi1VUr8TjpVjAW/EsU33P2ZEe9iytjpdiAW/Fk6xnKfM9esYYY3KKJRZjjDFJZYnF\nGGNMUlliMcYYk1SWWIwxxiSVJRZjjDFJlfbhxsYYY3Kb1ViMMcYklSUWY4wxSWWJxRhjTFJZYjHG\nGJNUlliMMcYklSUWY4wxSbWaHSRTQkTuBvbgLBn+aVU9mKbX/SrwPpyfxd8AvwFcDAy4l3xNVf9d\nRG4EPoOz/Pm9qvrNJMdxFfA94JB76A3gq8B3gAKcDaJuUtWZVMfixnMLcFPMoUuAl4ByILpz1H9X\n1ZdF5M+Bj3Nmufd9SYzjfOCHwN2q+g8ispEEfyYiUgj8v8BmIAT8vqoeT1ZsXmXPkj1LceJI27Pk\niXksInIl8Oeqer2I7AK+pap70/C6V7uv+2siUg+8Cvwn8JCq/ijmunLgFeAynF0BDwLvV9XBJMZy\nFfApVb0h5th9wD5V/Z6I3AW0A99OdSyLxHYl8AngPDfGX8acOwd4CNgLVONsTnWeqoaS8LrlwI+A\nt4HX3Ych4Z8J8BHgMlW9TUQ+CNyiqv9lrXF5mT1L9izFed20PkteaQq7BngEQFUPA7UiUpWG192P\n8+kAYBjnE8Q7NxGHy4GDqjribif7LM4Of6l2FfCo+/VjwLUZiuULwJ1xzl0N/FhVZ1W1D2fjqt1J\net0Z4Nc4e7vrq0j8Z3IN8LB77U9Jz3+zTLNnaXFXYc9S2p4lrzSFNQMvx3zf5x4bTeWLup8EolXR\nW4B9ONW8T4nIZ4Fe4FNuLH0xRXuB9SkIabeIPArUAXcA5ao6s+A10xULACJyKdCuqt0iAvBlEVkH\nHMapLseL5421vraqzgFz7utGreRnMn9cVcMiEhGRIlWdXWtsHmbPksOepRjpfpa8UmNZKK27/onI\nb+I8DJ/CaXO8XVU/APwC+NIiRVIR39s4D8BvAr8HfJOzE3+810z1z+oPcdpWAf4Op7nj/Zy91W46\n40nktTL1s/Iie5bsWUpEUp8lrySWTpyMGNWC05mUciJyHfB54MNu9e9JVf2Fe/pR4IJF4mvl7Crl\nmqlqh6o+qKoRVT0GdOM0Y5QueM2Ux7LAVcABN8aH3djAqTqn5WezwPgKfibzx93OR1+O11bAniV7\nlhKXsmfJK4nlCeAGABG5COhU1ZRvAi0i1cDXgOujHXYi8n0R2epechXwS+AF4FIRqRGRCpz2xaeT\nHMuNIvI59+tmoAm4D/iYe8nHgJ+kI5aYmFqAcVWdFRGfiPxURGrc01fh/Gz+E/h1ESlyr28F3kxF\nPK6fkvjP5AnOtPt/BPhZCuPyCnuW7FlKVMqeJU+MCgMQkb/FGX0QBm5T1dfS8Jr/Fad6fiTm8H04\n1fhJYBxnWF2viNwA/DnOMMC/V9V/TXIslcD9QA1QhFOVfxVnlEYJTkfe76tqMNWxxMR0MfDXqvph\n9/tPAH+J05begTMyZFJE/i/gRjeev1LVJ5P4+l8HtgBB9zVvxGlOWPZnIiIFwP8DbMfpvLxZVduT\nEZuX2bNkz1Kc10/bs+SZxGKMMSY3eKUpzBhjTI6wxGKMMSapLLEYY4xJKkssxhhjksoSizHGmKSy\nxGKMMSapLLEYY4xJKkssxhhjkur/BxZ9oIaPCwaTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "tags": [] } } ] }, { "metadata": { "id": "KiGyuN9hgNuV", "colab_type": "text" }, "cell_type": "markdown", "source": [ "Now we can solve the equation. For this, we need to calculate the point where the CDF is equal to the proportion of prices. We can use the function ```fsolve``` from the subpackage ```scipy.optimize```. We don't need the whole package, just one function, so we can import it as follows:" ] }, { "metadata": { "id": "U1V63ZMtgLdT", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "from scipy.optimize import fsolve" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "zq-yZHg-gWDZ", "colab_type": "text" }, "cell_type": "markdown", "source": [ "Now the function is available without calling any other package before it (so no ```x.fsolve```).\n", "\n", "The function ```fsolve``` finds the solution to equations of the form $f(x)=0$. It takes as mandatory arguments the function, and a first guess (which must be sufficiently close to the real solution!).\n", "\n", "Let's program the function to optimize." ] }, { "metadata": { "id": "Xy0XSuppgWX3", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "def flittle(y, p1 = 300, p2 = 100):\n", " out = 1- fp.cdf(y) - p2/p1\n", " return(out)" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "fr-xq-hCgXx4", "colab_type": "code", "outputId": "45b3c7f0-6640-4c87-a1e4-be0e0434fda9", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "flittle(100)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "0.6140136493229555" ] }, "metadata": { "tags": [] }, "execution_count": 20 } ] }, { "metadata": { "id": "But0Tom1gZIP", "colab_type": "code", "outputId": "7402580e-ddb8-45b0-d514-7b7eba133a02", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "fsolve(flittle, 100)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "array([282.93393485])" ] }, "metadata": { "tags": [] }, "execution_count": 21 } ] }, { "metadata": { "id": "s2YzMrjkgeVC", "colab_type": "text" }, "cell_type": "markdown", "source": [ "This is telling us that the optimal reservation quantity is 282 seats, i.e., we should only sell $1000 - 282 = 718$ seats to class 2.\n", "\n", "## EMSR-b\n", "\n", "Finally, we will estimate EMSR-b for continuous distributions. The pooling of the demand makes this a bit more of a complicated measure, but this can easily be solved with a good estimation of the demand.\n", "\n", "For example, if we assume that the demands are normal, such that each demand has a mean $\\mu_j$ and a variance $\\sigma^2_j$, then the remaining demand will simply be the sum of the remaining elements, i.e.\n", "\n", "$$\n", "\\mu_j = \\sum_{i=1}^j \\mu_i\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\sigma^2_j = \\sum_{i=1}^j \\sigma^2_i\n", "$$\n", "\n", "As an example, let us consider three classes such that:\n", "\n", "- Class one has $p_1 = 500, \\mu_1 = 50, \\sigma^2_1 = 100$.\n", "- Class two has $p_2 = 250, \\mu_1 = 300, \\sigma^2_2 = 1000$.\n", "- Class three has $p_3 = 100, \\mu_1 = 1000, \\sigma^2_2 = 5000$.\n", "\n", "This means that setting the price from class three requires pooling demand and applying Littlewood's rule. We start by defining this problem using vectors." ] }, { "metadata": { "id": "cj7knotpgalw", "colab_type": "code", "colab": {} }, "cell_type": "code", "source": [ "p = [500, 250, 100]\n", "mu = [50, 300, 1000]\n", "sigmasq = [100, 1000, 5000]" ], "execution_count": 0, "outputs": [] }, { "metadata": { "id": "mSkOnnD_gnQc", "colab_type": "text" }, "cell_type": "markdown", "source": [ "We can now calculate first the weighted average price observed by class three. The argument \"weights\" of ```np.average``` allows to quickly calculate weighted averages." ] }, { "metadata": { "id": "C31KAQyFglWN", "colab_type": "code", "outputId": "4003db4a-0353-4057-d253-147b5ca082ad", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "pavg = np.average(p[1:2], weights=mu[1:2])\n", "pavg" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "250.0" ] }, "metadata": { "tags": [] }, "execution_count": 23 } ] }, { "metadata": { "id": "1T-pPnK4grO8", "colab_type": "text" }, "cell_type": "markdown", "source": [ "And the parameters for the normal distribution will be:" ] }, { "metadata": { "id": "z76HD-HGgpDl", "colab_type": "code", "outputId": "d91b6189-daac-4a6b-9a1c-a677b0db1e02", "colab": { "base_uri": "https://localhost:8080/", "height": 54 } }, "cell_type": "code", "source": [ "muavg = np.sum(mu[0:2])\n", "ssqavg = np.sum(sigmasq[0:2])\n", "print('mu = ', muavg)\n", "print('sigmasq = ', ssqavg)" ], "execution_count": 0, "outputs": [ { "output_type": "stream", "text": [ "mu = 350\n", "sigmasq = 1100\n" ], "name": "stdout" } ] }, { "metadata": { "id": "y_LuffXrgugc", "colab_type": "text" }, "cell_type": "markdown", "source": [ "We can now apply Littlewood's rule iteratively using the normal function as a base. Starting from the first class to the second." ] }, { "metadata": { "id": "L0lo1dq5gsn8", "colab_type": "code", "outputId": "edb667b5-83ac-46e7-94ad-e78373acff9c", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "normal = stats.norm(loc = mu[0], scale = np.sqrt(sigmasq[0]))\n", "def flittle(y, p1 = p[0], p2 = p[1]):\n", " out = 1 - normal.cdf(y) - p2/p1\n", " return(out)\n", "\n", "fsolve(flittle, 70)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "array([50.])" ] }, "metadata": { "tags": [] }, "execution_count": 25 } ] }, { "metadata": { "id": "2nAtGK_QgzHi", "colab_type": "text" }, "cell_type": "markdown", "source": [ "This means, 50 seats should be saved for class 1. \n", "\n", "We can now calculate the prices for the aggregated class to class three." ] }, { "metadata": { "id": "hIH697nZgw3j", "colab_type": "code", "outputId": "4dcb05e0-0874-49d5-9eee-710d7288cd75", "colab": { "base_uri": "https://localhost:8080/", "height": 35 } }, "cell_type": "code", "source": [ "normal = stats.norm(loc = muavg, scale = np.sqrt(ssqavg))\n", "def flittle(y, p1 = pavg, p2 = p[2]):\n", " out = 1 - normal.cdf(y) - p2/p1\n", " return(out)\n", "\n", "fsolve(flittle, 300)" ], "execution_count": 0, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "array([358.40257283])" ] }, "metadata": { "tags": [] }, "execution_count": 26 } ] }, { "metadata": { "id": "dxoz9E4vg3Qa", "colab_type": "text" }, "cell_type": "markdown", "source": [ "So 358 seats should be saved for the remaining classes. That means that $358 - 50 = 308$ are saved for class 2. The remaining $1000 - 358 = 642$ can be freely sold to class three.\n", "\n", "**Note: An earlier version of this notebook incorrectly showed 692 seats for class 3. This has been corrected (thanks Hana!)**\n", "\n", "## Self-study\n", "\n", "Implement EMSR-a for the example above, and compare the results. Which one gives better results? Does that satisfy what you expected.\n", "\n", "Now you are ready to answer question 2 of the first coursework!" ] } ] }