{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# 3: Mathematical expectation\n", "\n", "\n", " #### [Back to main page](https://petrosyan.page/fall2020math3215)\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGbCAYAAAD3MIVlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAATy0lEQVR4nO3db4xl9X3f8c93/2F3F2OgLsJ/VFDr4hTctHhwLNkJTlssIElJZUdxikjSVqKJsJT0QSyjPEnxozRRJbcmRVh1W/+JrbQpEjV/bMcCu6nihl0Hm2BDICwR60WGxbJZNmTZWX59MLPyZJll7zL7nTO79/WSrmbuPefMfM+9l9n3nHPnUmOMAABwcm2aegAAgNORyAIAaCCyAAAaiCwAgAYiCwCgwZapB1jNlVdeOe6+++6pxwAAmEWtduOGPJK1b9++qUcAAFiTDRlZAACnOpEFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANJgpsqrqyqp6uKoeraoPrbL8mqr6RlXdX1U7q+pdK5Y9XlUPHFl2MocHANiothxvharanOTmJFck2ZPkvqq6fYzxzRWrfSnJ7WOMUVX/IMnvJXnLiuU/PsbYdxLnBgDY0GY5kvX2JI+OMR4bY7yQ5LNJrlm5whjjuTHGWL66PckIAMAcmyWy3pDkiRXX9yzf9tdU1T+vqoeS3JHkX61YNJJ8oap2VdX1axkWAOBUMUtk1Sq3veRI1RjjtjHGW5L8dJIPr1j0zjHGpUmuSnJDVf3Yqt+k6vrl13PtfPrpp2cYCzhi8cX5Png87/sPbEzHfU1Wlo5cvWnF9Tcm2XuslccYX6mqv1NVf3OMsW+MsXf59qeq6rYsnX78yirb3Zrk1iRZWFjwExNOwJZNlavveia79y9OPcq6u/DMLbnzqnOnHgPgJWaJrPuSvLmqLkzy7STvT/IvVq5QVX83yZ8vv/D90iTbkjxTVduTbBpj7F/+/D1JbjqpewAkSXbvX8wjzx6eegwAlh03ssYYi1X1gSSfT7I5ycfHGA9W1S8tL78lyXuT/HxVHUryfJKfXQ6u85LcVlVHvtfvjjHubtoXAIANY5YjWRlj3JnkzqNuu2XF57+Z5DdX2e6xJD+8xhkBAE453vEdAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYzRVZV/XpVvbB8uWuV5TdX1fPLl2er6mdm3RYA4HR03Miqqq1JfiPJFUnOTvLuqvqpo1b7epIfGmO8OslvJfmvJ7AtAMBpZ5YjWb+Y5PtjjC+PMQ4k+XKSG1auMMa4dYzx+PLVTyZ59azbAgCcjmaJrIuS7Ftx/fEkr3+Z9f9TkkdOdNuq+mRVHaiqA7t3755hLPiBxRfH1CPAZOb9+T/v+8/GtWWGdWqV21Z9RlfVv03yniQXn+i2Y4zrklyXJAsLC/6L4YRs2VS5+q5nsnv/4tSjrLt3nbctH7v87KnHYELz/Py/8MwtufOqc6ceA1Y1S2Q9lOQXVly/IMmTR69UVe9N8u+T/MQY49ET2RZOht37F/PIs4enHmPdXbBj/vaZl5rX5z9sZLOcLvxEkrOq6keranuSy5P8zsoVquodST6T5IYxxhdOZFsAgNPRcSNrjHEwyYeTfCnJ95L8nzHG7VX16ar69PJq/y1LR8U+svw2DgdebtuTvxsAABvLLKcLM8a4KclNR9127YrP33Ii2wIAnO684zsAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1miqyq+vWqemH5ctcqy6+qqv1VNarqf5/ItgAAp6PjRlZVbU3yG0muSHJ2kndX1U8dtdpfJPlAkv/7CrYFADjtzHIk6xeTfH+M8eUxxoEkX05yw8oVxhjfHGP89ySLJ7otAMDpaJbIuijJvhXXH0/y+hm//szbVtUnq+pAVR3YvXv3jF+eIxZfHFOPAMAE/PzfuPfBlhnWqVVum3VvZt52jHFdkuuSZGFhYWPeWxvYlk2Vq+96Jrv3H30w8fT3rvO25WOXnz31GACTmOef/0ly4ZlbcudV5049xqpmiayHkvzCiusXJHlyxq+/lm05Qbv3L+aRZw9PPca6u2DH/O0zwErz+vN/o5vldOEnkpxVVT9aVduTXJ7kd2b8+mvZFgDglHXcI1ljjINV9eEkX8rS6b97xhi3V9Wnl5dfW1VvTfInSTYnSVUtJvnbY4xvr7Zt074AAGwYs5wuzBjjpiQ3HXXbtSs+f+BYX2u1bQEATnfe8R0AoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGswUWVV1ZVU9XFWPVtWHVlleVfUfl5d/o6ouXbHs8ap6oKrur6qdJ3N4AICNasvxVqiqzUluTnJFkj1J7quq28cY31yx2lVJ3rx8+ZEk/3n54xE/PsbYd9KmBgDY4GY5kvX2JI+OMR4bY7yQ5LNJrjlqnWuSfGIs+WqS11bV+Sd5VgCAU8YskfWGJE+suL5n+bZZ1xlJvlBVu6rq+mN9k6q6vqp2VtXOp59+eoaxfmDxxXFC6wMAdDvu6cIktcptR1fNy63zzjHG3qr6W0m+WFUPjTG+8pKVx7g1ya1JsrCwcELVtGVT5eq7nsnu/Ysnstlp413nbcvHLj976jEAgBVmiaw9Sd604vobk+yddZ0xxpGPT1XVbVk6/fiSyFqr3fsX88izh0/2lz0lXLBjPvcbADayWU4X3pfkzVV1YVVtS/L+JLcftc7tSX5++a8M35Hk+2OMJ6tqe1WdmSRVtT3Je5L86UmcHwBgQzrukawxxmJVfSDJ55NsTvLxMcaDVfVLy8tvSXJnkquTPJrkL5P8y+XNz0tyW1Ud+V6/O8a4+6TvBQDABjPL6cKMMe7MUkitvO2WFZ+PJDesst1jSX54jTMCAJxyvOM7AEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQAORBQDQQGQBADQQWQAADUQWAEADkQUA0EBkAQA0EFkAAA1EFgBAA5EFANBAZAGT2nx4Ma/5y2ez6cXDU48CcFJtmXoAYP5sO3Qw79v1uXzwrpvz95/8sxzatCVbX1zMg+f/vfzWVTfkf77tJ/PC1jOmHhNgTRzJAtbVZbv/JE/82qX56KduzCV7H86mMXLG4UPZNEbeuvfhfPRTN+aJX7s0C7vvn3pUgDURWcC6Wdh9f/7gt38m5x74Xl5z8MCq67zm4IGce+B7+YPffp/QAk5pIgtYF9sOHcwdH7k22194fqb1d7zwfO74yLXZduhg82QAPUQWsC7et+tz2bp46IS22bZ4KO/ddUfTRAC9RBawLj54183HPEV4LGcePJAP3vXRpokAetUYY+oZXuLMM88cb3vb205omwe+eyjPH954+7Ieztq6KRe9dsvc3gf2f+Pv/6Yx8uy3vvqKfqt7MclrfugdebFq1eWv3lx56zlb1zTf6WAjP/6dPP5L5vXxTzbGc+Dee+9d9QfUafEWDiOZ/A6e2rzfB/Z/Y+//9sXFHK7KplfwS93hqlx+zuYc2HLsH1cjyeoJNh82+uPfzeM/349/snGfAxvySNbCwsLYuXPn1GO8vO9+N3n3u5O/+qsf3PaqVyX33pucc85UU62veb8P7P/s+3/4cLJ1a/JKft5UJYcOJZs3r2VaTrZ5f/7PO4//0VZtPK/JeqXOOSe58cZkz57kO99Z+njjjfP15Jr3+8D+z77/mzcnF1/8yr7PxRcLrI1o3p//887jPxNHstbq6aeTxx9PLrgged3rpp5mGvN+H9j/2fb/U59KfvmXk+eem/1r79iR3HJLcu21a52SLvP+/J93Hv8jVj2SJbKA9XHwYPL61y+dZpjVOecke/cmZ/hf7AAbmtOFwITOOCO5++5k+/bZ1t++fWl9gQWcokQWsH4uuyy5556lI1Q7dqy+zo4dS8vvuWdpfYBTlMgC1tdlly2dArzlluSSS5b+enDr1qWPl1yydPvevQILOOV5TRYwrcOHl14Mv2OHvyIETlWn75uRAqewzZuTs86aegqAk87pQgCABiILAKCByAIAaCCyAAAaiCwAgAYiCwCggcgCAGggsgAAGogsAIAGIgsAoIHIAgBoILIAABqILACABiILAKCByAIAaCCyAAAaiCwAgAYbMrJ27dr1+SR1qlyq6t9MPcPUl3m/D+z/fO//vF88/vN98fincgw1xjjWMmZUVTvHGAtTzzGleb8P7P987/+88/jPN4//sW3II1kAAKc6kQUA0EBknRy3Tj3ABjDv94H9Z555/Oebx/8YvCYLAKCBI1kAAA1EFgBAA5G1RlV1ZVU9XFWPVtWHpp5nvVXVx6vqqar606lnWW9V9aaquqeqvlVVD1bVr0w903qqqldV1R9X1deX9//fTT0T66+qHq+qB6rq/qraOfU8rJ+qumj5cT9yebaqfnXquTYSr8lag6ranOTPklyRZE+S+5L83Bjjm5MOto6q6seSPJfkE2OMS6aeZz1V1flJzh9jfK2qzkyyK8lPz8vjX1WVZPsY47mq2prkD5P8yhjjqxOPxjqqqseTLIwx9k09C9NZ/vfw20l+ZIzxF1PPs1E4krU2b0/y6BjjsTHGC0k+m+SaiWdaV2OMryT57tRzTGGM8eQY42vLn+9P8q0kb5h2qvUzljy3fHXr8sVvbTCf/kmSPxdYf53IWps3JHlixfU9maN/ZPmBqrogyT9K8v+mnWR9VdXmqro/yVNJvjjGmKv9J8lSWH+hqnZV1fVTD8Nk3p/kM1MPsdGIrLVZ7f9X5Df5OVNVO5L8fpJfHWM8O/U862mMcXiM8Q+TvDHJ26tqrk4ZkyR55xjj0iRXJblh+SUEzJGq2pbknyX5H1PPstGIrLXZk+RNK66/McneiWZhAsuvRfr9JJ8eY/yvqeeZyhjje0nuTXLlxKOwzsYYe5c/PpXktiy9jIL5clWSr40xvjP1IBuNyFqb+5K8uaouXC759ye5feKZWCfLL/z+L0m+Ncb4D1PPs96q6nVV9drlz1+d5J8meWjaqVhPVbV9+Y8+UlXbk7wnydz9pTH5uThVuCqRtQZjjMUkH0jy+Sy96Pn3xhgPTjvV+qqqzyT5oyQXVdWeqvrXU8+0jt6Z5Lok/3jFnzBfPfVQ6+j8JPdU1Tey9AvHF8cYn5t4JtbXeUn+sKq+nuSPk9wxxrh74plYR1X1N7L0F/ZzeyT/5XgLBwCABo5kAQA0EFkAAA1EFgBAA5EFANBAZAEANBBZAAANRBYAQIP/D8Dqqsud+4Q3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "# library\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "plt.rcParams[\"figure.figsize\"] = (10, 7)\n", "\n", "\n", "range_x = np.array([0, 1, 2, 3, 5, 7])\n", "pmf_values = np.array([1/20, 1/10, 2/10, 7/20, 2/10, 1/10])\n", "\n", "fig, ax2 = plt.subplots(num=1, clear=True)\n", "\n", "\n", "ax2.set_ylim(-0.03, 0.4) \n", "ax2.set_xlim(-0.7, 8)\n", "ax2.axhline(y=0, color='k')\n", "ax2.set_xticks(range_x)\n", "ax2.set_yticks(pmf_values)\n", "ax2.spines[\"top\"].set_visible(False) \n", "ax2.spines[\"right\"].set_visible(False)\n", "ax2.spines[\"bottom\"].set_visible(False)\n", "\n", "\n", "# PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n", "ax2.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n", "ax2.bar(range_x, pmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, label=\"Histogran\")\n", " \n", "mean = np.average(range_x, weights=pmf_values)\n", "ax2.scatter(np.array([mean]),np.zeros(1), color =\"red\", s=200, label=\"Mean\", zorder=3)\n", "\n", "plt.show();\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mean and Variance under linear transformation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Theorem**\n", "\n", "
\n", "\n", "Suppose $X$ is a random variable with mean $\\mu_X$ and standard deviation $\\sigma_Y$. Let $Y=aX+b$ where $a$ and $b$ are any two numbers. Then mean and the standard deviation of $Y$ are\n", "$$\\mu_Y=a\\mu_X+b,\\quad \\sigma_Y=|a|\\sigma_X.$$\n", "\n", "
\n", " \n", "**Proof**\n", "$\\mu_Y=E[aX+b]=aE[X]+b=a\\mu_X+b.$ Notice that\n", "\n", "$$\\text{Var}(Y)=E[(Y-E[Y])^2]=E[(aX+b-aE[X]-b)^2]=E[(aX-aE[X])^2]=a^2\\text{Var}(X).$$\n", "\n", "Taking square roots, we have $\\sigma_Y=|a|\\sigma_X$. $\\blacksquare$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example**\n", "\n", "Let $X$ be a random variable with $\\text{range}(X)=\\{-2,-1,0,1,2\\}$ and \n", "\n", "$$f_X(-2)=0.4, \\quad f_X(-1)=0.25, \\quad f_X(0)=0.15, \\quad f_X(1)=0.1, \\quad f_X(2)=0.1.$$\n", "\n", "Take the random variable $Y=2X+1$.Notice that $\\text{range}(Y)=\\{-3,-1,1,3,5\\}$ and \n", "\n", "$$_Y(-3)=0.4, \\quad f_Y(-1)=0.25, \\quad f_Y(1)=0.15, \\quad f_Y(3)=0.1, \\quad f_Y(5)=0.1.$$\n", "\n", "It can be checked that \n", "\n", "$$\\mu_X=-0.75, \\quad \\sigma_X\\approx 1.414$$\n", "\n", "and thus \n", "\n", "$$\\mu_Y=-0.5, \\quad \\sigma_Y\\approx 2.828. $$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABIcAAAJTCAYAAACSDnaYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5Skd13n8c83mQnCkDkECCgJlxwPqJhA8IyAoosIsklwRfBChEU4ilkv8baygqLogu4uoqx75LZRWC4HibAkLkIuXFZh1aBEDBIugUDCIQmQC2AGAkkm+e0f1WF6Oj3T1Zmequr+vl7nzElX1fNU/fqhZ+rLu5+qqjFGAAAAAOjpsHkvAAAAAID5EYcAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoTh4BNo6q+q6ouqKr3VNUbq2r7vNcEALDVmcFg6xOHgM3k00m+f4zx6CSfSvLEOa8HAKADMxhscdvmvQCAaY0xrlp2cU+SW+e1FgCALsxgsPU5cwhIVd29qs6uqq9U1aer6qlrbP/lFX9uqao/WXb731TV15bdfskGr/e4JCcnedsG3NfpVXVhVd1YVa+ZYvv9fu9V9YCqOqeqvlhVn6uql1bVtmX7HvB2AKCf9cxha81Y653p7sBaN2wGW7q/A86UK7Y94Mw272MDm504BCTJy5LclOTeSZ6W5BVV9e3723iMcdfb/izt89Ukb16x2enLtvuWjVpoVe1M8tokTx9j3LQBd3lVkt9L8uppNl7je395kquTfFOSE5M8OsnPL9t9rdsBgH7WNYflwDPWeu9raodgBpt2przNNDPbXI4NbAXiECyAqnpeVb1i2eWjqurmqvqGGTz2jiQ/kuS3xxhfHmP8bZK3Jnn6lHfxo5kEj/+3QevZXlW/X1WXLx2DsfTng0tn2bwxye+OMTbkbKQxxlljjL9Mct0d2H3l935ckjeNMb42xvhckvOSLB861rodAJixTT6Hbdh9zXoGW8UBZ8qDmdk28jjDViUOwWI4IclFyy6fmOSSMcbXptm5qt5WVV/az5+1Tvt9UJJbxhgfX3bdBzN9tHhGkteNMcaK6/9rVV1bVX9XVd835X0lk98IPTbJ9ya5W5J3Jzk7yZOS/ESSRyR5/tKpw09ZuXNV3aeq/riq/m9VvbKqHltVd6mqh1XVf17HOqax8nv/H0lOXXq8YzI57fq8ZduvdTsAMHt3eA47yBksuWNz2P5mrIOd6eY9g+1vplyPQ3VsYMvzXhewGE5I8t+XXT4xkyesVNXrk/zRGOOiqvrVJEePMX5z+c5jjB88iMe+a5J/XXHdvyY5cq0dq+p+mbw06qdX3PScJB/J5NTdU5P8VVWdOMb45Br3d2SSX0rykDHGZ5aue0uSp4wxPpXJp2O8fo1l/V6S92VySvJ3Jvn9JA9O8tEkGxaH9vO9vyfJzyS5PsnhmZx6/ZfruB0AmL1V57CqemySp40xfir5+kz2F2OMr0efg5zBkvXPYQeasQ5mppvrDHaAmXI9DsmxgS6cOQRzVlVHJPnmJB9advVDs/c3WC9P8vNV9YgkpyT57YN8vKcte6O+c5N8OcnOFZvtTLJ7irv7ySR/O8a4bPmVY4x/GGPsHmPcOMZ4bZK/W1r7Wv5Nkk+NMT6x7Lqjknxuin1v8+wkR2QyFGzP5D197pHkqUnuu477Wcs+33tVHZbk/CRnJdmR5J6ZrP1F09wOAMzegeawMca7kzykqnZW1WlJrloehu7g4x3UHLbGjHUwM928Z7BVZ8r1OITHBloQh2D+HpzkyjHGDUlSVZXk+7J05tAY44Ik35rkT5L85BjjlpV3UFXn1u0/7WH54PF1Y4w3LHujvpOTfDzJtqp64LLNHprkw1Os/SczOftlLSNJTbHd0Um+eNuFpWPxpCx9IkZVvb6qTlz6+ler6r+sch9/kMlHrP7vpcd9ZZIvJXlTJm9kuFFWfu93z2TweenSUHJdkv+VvUPJWrcDALN3wDksk+f6F2XyfjjPW7nzemawZMPnsGTfGetg7mutGeyxVfXqZbe/vqpWnjV1MDPYtDPlemzUsYEWvKwM5u+EJPeqqm/O5InzeUnun+Ty5OtPzl9L8jdjjM+udgdLw8UdMsb4SlWdleQFVfWsTE6lfmKS7z7QflX13UmOyYpPlKiqu2XymvT3ZDIgPCWT30b9ytLtr1l63GeucrcXJ/mOpQB0SZLfyeSJ/S+Wbr/tLKpXZRJVTlrlPv7DsoD2niR/uMb3sS2TfwsPT3J4Td58cs8YY88B9rnd9z7GuLaqLkvyc1X1h5mcvvyM7I18B7wdAJiLA85hmQSLy5OcsNpscDAz2NL+U89ha81Y09zXAeawA85gY4x3V9WLavKJZadm9bOo1jWDLVvTqjPlKtvtd2bbiGMD3TlzCObvhExebnRukkuTfD6T13Xf9tup/5Tkn5M8fikUHQo/n+TOmXxCxBuT/NwY4+u/SVn6rdhvrtjnGUnOGmOsPB13eyavOb8mybVJfjHJDy/7ZIv7ZnKa7+2MMS7M5PXp52RyDL4xySljjJuXbl/zLKrVrlvDb2XysanPTfLvl77+rdtuXOf3/uRMgtU1mfxvuSfJr67jdgBgttaaw45I8sExxpWHcA37ncNWzCFrzVgHvK8lq85ha81gSw54FtUdmMFus+pctcoMdqCZbSOODbRWB/dm8MDBWjrl+M/GGG9Z5bbvyeTJ9wlJ/ijJuWOMd8x4iRtm6XX9H8zkzQ5vXmv7VfavTAa4940xnr/R6wMAejnQHLZ0+6OT/NgY4/TZrmzjbcActjN7z6I6lLEMmANnDsH8nZDJpzjso6qOzuSTM54xxrg1ySsy+Y3HpjXGuGmM8W13ZCBZMouzqACAPladw5b59myR96XZgDlsFmdRAXPizCGYo6o6KpPTl3ccxBN1C1vtLCoAYL7MYeuzlc6iAm7PmUMwR2OML44xjjCQHNhWPIsKAJgvc9i6bZmzqIDbc+YQAAAAQGPOHAIAAABoTBwCAAAAaEwcAgAAAGhs27wXsJqTTjppnHfeefNeBgBw6NS8F8C+zF8A0MKqM9hCnjl07bXXznsJAACtmL8AoK+FjEMAAAAAzIY4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANDYVHGoqk6qqkuq6tKqeu4BtvvOqrqlqn50vfsCALAvMxgAMAtrxqGqOjzJy5KcnOTBSX6iqh68n+1elOT89e4LAMC+zGAAwKxMc+bQw5NcOsb41BjjpiRnJnniKtv9YpK3JLn6DuwLAMC+zGAAwExME4eOSfKZZZevWLru66rqmCRPSvLK9e4LAMCqzGAAwExME4dqlevGist/nOQ5Y4xb7sC+kw2rTquqC6vqwmuuuWaKZTGtPbeueshnblHWAQCbxCGfwbbq/LUoM8eirAMA1rJtim2uSHLfZZePTXLVim12JTmzqpLknklOqao9U+6bJBljnJHkjCTZtWuXZ9INtO2wyinnXpfLdu+Z2xqOO3Jbzjn5HnN7fADYhA75DLZV5y+zDwCszzRx6P1JHlhVxyW5MsmpSZ66fIMxxnG3fV1Vr0nytjHGX1bVtrX2ZTYu270nn7h+5S8VAYAFZgY7CGYfAJjemnFojLGnqk7P5BMwDk/y6jHGh6vqZ5duX/ka9zX33ZilAwBsXWYwAGBWpjlzKGOMc5Kcs+K6VQeSMcYz19oXAIC1mcEAgFmY5g2pAQAAANiixCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGxCEAAACAxsQhAAAAgMbEIQAAAIDGpopDVXVSVV1SVZdW1XNXuf1pVfUvS3/+vqoeuuy2y6vqQ1V1UVVduJGLBwDYysxgAMAsbFtrg6o6PMnLkvxAkiuSvL+q3jrG+MiyzS5L8ugxxher6uQkZyR5xLLbHzPGuHYD1w0AsKWZwQCAWZnmzKGHJ7l0jPGpMcZNSc5M8sTlG4wx/n6M8cWli+9LcuzGLhMAoB0zGAAwE9PEoWOSfGbZ5SuWrtufn05y7rLLI8k7quqfquq09S8RAKAlMxgAMBNrvqwsSa1y3Vh1w6rHZDKYfM+yqx81xriqqu6V5J1V9bExxntX2fe0JKclyf3ud78plsVmcvc7VfbcOrLtsNV+nGZrUdYBAGs45DOY+auHRZh9FmENAOzfNHHoiiT3XXb52CRXrdyoqh6S5M+SnDzGuO6268cYVy399+qqOjuTU6RvF4fGGGdk8jr57Nq1a9XBh81r5/bDsu2wyinnXpfLdu+Z2zqOO3Jbzjn5HnN7fABYh0M+g5m/epj3DGb+Alh808Sh9yd5YFUdl+TKJKcmeeryDarqfknOSvL0McbHl12/I8lhY4zdS18/PskLNmrxbD6X7d6TT1x/y7yXAQCbgRmMDWMGA+BA1oxDY4w9VXV6kvOTHJ7k1WOMD1fVzy7d/sokz09yjyQvr6ok2TPG2JXk3knOXrpuW5I/H2Ocd0i+EwCALcQMBgDMyjRnDmWMcU6Sc1Zc98plXz8rybNW2e9TSR56kGsEAGjJDAYAzMI0n1YGAAAAwBYlDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQ2VRyqqpOq6pKqurSqnrvK7d9aVRdU1Y1V9ewVt11eVR+qqouq6sKNWjgAwFZnBgMAZmHbWhtU1eFJXpbkB5JckeT9VfXWMcZHlm32hSS/lOSH93M3jxljXHuwiwUA6MIMBgDMyjRnDj08yaVjjE+NMW5KcmaSJy7fYIxx9Rjj/UluPgRrBADoyAwGAMzENHHomCSfWXb5iqXrpjWSvKOq/qmqTlvP4gAAGjODAQAzsebLypLUKteNdTzGo8YYV1XVvZK8s6o+NsZ47+0eZDK0nJYk97vf/dZx9zC9u9+psufWkW2HrfZjPVuLsg4AFtYhn8HMX3SzKPPXoqwD4DbTxKErktx32eVjk1w17QOMMa5a+u/VVXV2JqdI3y4OjTHOSHJGkuzatWs9gw9Mbef2w7LtsMop516Xy3bvmds6jjtyW845+R5ze3wANoVDPoOZv+jGHAiwumni0PuTPLCqjktyZZJTkzx1mjuvqh1JDhtj7F76+vFJXnBHFwsb5bLde/KJ62+Z9zIA4EDMYHAImAMBbm/NODTG2FNVpyc5P8nhSV49xvhwVf3s0u2vrKpvTHJhkp1Jbq2qX0ny4CT3THJ2Vd32WH8+xjjv0HwrAABbhxkMAJiVac4cyhjjnCTnrLjulcu+/lwmpzqvdH2Shx7MAgEAujKDAQCzMM2nlQEAAACwRYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNiUMAAAAAjYlDAAAAAI2JQwAAAACNTRWHqup5VXXT0p9zV7n95KraXVWjqv5qPfsCALA6MxgAMAtrxqGq2p7kd5P8QJKjknxfVf27FZt9OsnpSf7uDuwLAMAKZjAAYFamOXPomUn+dYzxnjHGV5K8J8kvLN9gjPGRMcZrk+xZ774AAKzqmTGDAQAzME0c+pYk1y67fHmS+0x5/wezLwBAZ2YwAGAmpolDtcp1Y8r7n3rfqnp9VX2lqr5y2WWXTXn3sDnd/U6VPbdO+9fo0FmENSSLsw6ABXPIZzDzF/S2KDOYdSzWGuhp2xTbfCzJM5ZdfkCSz055/1PvO8Z4epKnJ8muXbv8jWBL27n9sGw7rHLKudflst0rXwkwG99z7yPyp48+aq5rSJLjjtyWc06+x9weH2CBHfIZzPwFvc17Hk0Waxac9/FYpGNBP9PEodcleXlVfW+SDyR5dJKnTHn/B7MvbHmX7d6TT1x/y1we+wF3vWXuawDggMxgwCFnFtyX40FXa8ahMcaNVfXCJO/O5BTlvx5jvLWq3rB0+9Oq6oQk/5zk8CSpqj1J7j/GuHK1fQ/R9wIAsGWYwQCAWZnmzKGMMV6Q5AUrrnvasq8/tL/7Wm1fAADWZgYDAGZhmjekBgAAAGCLEocAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4cAAAAAGhOHAAAAABoThwAAAAAaE4fYMo786u4c84WrcuRXd897KQAALNlxw+7kiiuS66+f91IA2I9t814AHIxte27Okz9wTn793Jfm26/6eG7atj1H7Lk5H77Pg/IHJ5+es77jlOzZtn3eywQAaGWfGe2zH0/udERy003J8ccnz3lO8iM/kmw3owEsCnGITetuX/lSzn/JqXnQ5z+ZI2+8IUmy/aY9SZITr/hI/ufrnp1fO/8V+bf/8cx8acfd5rlUAIA2VpvRcsNkRstFFyU/8zPJi1+cvOtdyVFHzW+hAHydl5WxKW3bc3POf8mpOf7Kj+0dOlY48sYbcvyVH8v5Lzk12/bcPOMVAgD0M82Mli9/Obn44uRxj0tuNqMBLAJxiDXdc/d12XXZRbnn7usWZh1P/sA5edDnP5k73XLggeJOt9ycB33+k3nSP59zSNYxL4uwhkVaBwBsJM9v+1rP8Zh2RstNNyWXXJKcddYGrXI2/Gzsy/HYy7Fgs/OyMg7oKf9wdv70tc/++nv5POsZf5Q3PeKH576Oq4+85/5/G7XCkTfekF8/92V583c+ccPXMY/jsQhrWKR1AMBG8vy2r/Uej18/96VTz2j5yleSF70oecpTNmi1h5afjX05Hns5FmwF4lATxx25/v+pj7r+2vzZ656dO9/8tdzl5q8lSV71ul/L5Y94dL64857ruq9jdxy+oeu4/xeuWNd9HH/VJTlx2w05dsedN3Qdsz4ei7CGjV7HHXl8AFiL57d9zeJ47Lhh9+TNp9fj4osnn2K2c+e613dH+dnYl+Oxr3nO54t2LOilxhjzXsPtHHPMMeMJT3jCvBim4S8AAAbCSURBVJexZdw6Rg6rWv+O11yTvP3t+74WfPv25AlPSI4+em7ruFuSFya50zru4sbDD89v//iP50t3vevmPh6LsIZFWgewaZ1xxhn+0i+YrTR/eX7boHWs83jc7ctfzgvf/Obcac+eqR9i+Yw2C342NmgdjsdeW/RYsHXtbwZbyDi0a9euceGFF857GVxzTXL/+ydf/ere6+585+TTn75D/9Bt6DrWa/v25NprD+63UotwPBZhDYu0DmAzM/kuGPNXPL+ttN7jcf31yT3ukawjDm3IjDYLfjb25Xjs5Viw+aw6g3lDavbv6KOTV71q8o/bzp2T/77qVbP/R261ddz//uu7j+OPP/ihYxGOxyKsYZHWAQAbyfPbvtZ7PHbunMxc67ERM9os+NnYl+Oxl2PBFuHMIdZ2zTXJ5ZcnD3jAfP+RW76Od787edazJm9kuJYdOyb/QG/Umx0uwvFYhDUs0jqAzciZQwvG/LWM57d9red4nHnm/Ga0WfCzsS/HYy/Hgs3Dy8rYQm6+OXnkIydvYnjTTfvf7ogjkhNOSC64YHLaMgCLQhxaMOYvNoQZDWDReVkZW8j27cm73jU5FXnHjtW32bFjMnS8852GDgCAWTCjAWxK4hCb11FHJe973+R05Ic9bDJc3OUuk/8+7GGT6y+4YLIdAACzYUYD2HS8rIyt4/rrJ3927twcb2wI0JuXlS0Y8xeHjBkNYJGsOoNtm/Uq4JAxcAAALB4zGsDC87IyAAAAgMbEIQAAAIDGFvI9h6rqvDHGSfNeB3tV1WljjDOsY3HWsQhrWKR1AHBwzF/78vy2L8djL8diX47HXo4Fm9lCxiEWT1VdOMbYZR2Ls45FWMMirQMANpLnt305Hns5FvtyPPZyLNjMvKwMAAAAoDFxCAAAAKAxcYhpLcprZ61jr0VYQ7I46wCAjeT5bV+Ox16Oxb4cj70cCzYt7zkEAAAA0JgzhwAAAAAaE4eYWlW9uKo+VlX/UlVnV9Xd5rSOH6uqD1fVrVU1008DqKqTquqSqrq0qp47y8detoZXV9XVVXXxPB5/2TruW1V/XVUfXfrf45fnuR4A2GhV9cKlueeiqnpHVd1n3mual3nOX4tiUWawRVBV31BV/1hVH1z6ufjP817TvFXV5VX1oaV/Ly6c93pgvcQh1uOdSY4fYzwkyceT/Mac1nFxkicnee8sH7SqDk/ysiQnJ3lwkp+oqgfPcg1LXpPkpDk87kp7kvzaGOPbkjwyyS/M6XgAwKHy4jHGQ8YYJyZ5W5Lnz3tBczSX+WvBvCaLMYMtghuTfP8Y46FJTkxyUlU9cs5rWgSPGWOc6OPs2YzEIaY2xnjHGGPP0sX3JTl2Tuv46Bjjkjk89MOTXDrG+NQY46YkZyZ54qwXMcZ4b5IvzPpxV1nHZ8cYH1j6eneSjyY5Zr6rAoCNM8a4ftnFHUnavlnnHOevhbEoM9giGBNfXrq4felP278fsBWIQ9xRP5Xk3HkvYsaOSfKZZZeviBiSJKmqByR5WJJ/mO9KAGBjVdXvV9Vnkjwtvc8cgn1U1eFVdVGSq5O8c4zRfQ4cSd5RVf9UVafNezGwXtvmvQAWS1W9K8k3rnLT88YY/2dpm+dl8pKiN8xzHXNQq1zX/jckVXXXJG9J8isrfsMKAAtvrZljjPG8JM+rqt9IcnqS35npAmdoQecvFtQY45YkJy69D+nZVXX8GKPz+zE9aoxxVVXdK8k7q+pjS2ebwaYgDrGPMcbjDnR7VT0jyQ8meewY45CFkbXWMSdXJLnvssvHJrlqTmtZCFW1PZMw9IYxxlnzXg8ArNc6Zo4/T/L2bOE4tKDzFwtujPGlqvqbTN6PqW0cGmNctfTfq6vq7EzekkIcYtPwsjKmVlUnJXlOkh8aY9ww7/XMwfuTPLCqjquqI5KcmuStc17T3FRVJXlVko+OMV4y7/UAwEarqgcuu/hDST42r7XAIqmqo2/75OKqunOSx6Xx34+q2lFVR972dZLHp3EoY3MSh1iPlyY5MpPTJC+qqlfOYxFV9aSquiLJdyV5e1WdP4vHXXoz7tOTnJ/Jmy+/aYzx4Vk89nJV9cYkFyT5lqq6oqp+etZrWPKoJE9P8v1LPw8XVdUpc1oLABwK/62qLq6qf8nk/+z98rwXNC/zmr8WyQLNYIvgm5L89dLfjfdn8p5Db5vzmubp3kn+tqo+mOQfk7x9jHHenNcE61KH8JVBAAAAACw4Zw4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADQmDgEAAAA0Jg4BAAAANCYOAQAAADT2/wHC62yUuXnvowAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "plt.figure(figsize=(20,10)) \n", "\n", "# values\n", "range_x = np.array([-2,-1,0,1,2])\n", "range_y = np.array([-3,-1,1,3,5])\n", "xpmf_values = np.array([0.4, 0.25, 0.15, 0.1, 0.1])\n", "ypmf_values = np.array([0.4, 0.25, 0.15, 0.1, 0.1])\n", "\n", "# mean \n", "mean_x = np.average(range_x, weights = xpmf_values)\n", "mean_y = np.average(range_y, weights = ypmf_values)\n", "\n", "# variance\n", "mean_x2 = np.average(np.power(range_x,2), weights = xpmf_values)\n", "mean_y2 = np.average(np.power(range_y,2), weights = ypmf_values)\n", "var_x = mean_x2 - mean_x**2\n", "var_y = mean_y2 - mean_y**2\n", "\n", "# set up the figure\n", "fig, [ax1, ax2] = plt.subplots(1,2, num=1, clear=True)\n", "\n", "ax1.set_ylim(-0.01,0.5) \n", "ax1.set_xlim(-5, 10)\n", "ax1.axhline(y=0, color='k', linewidth = 0.6)\n", "ax1.set_xticks(range_x)\n", "ax1.set_yticks(xpmf_values)\n", "ax1.spines[\"top\"].set_visible(False) \n", "ax1.spines[\"right\"].set_visible(False)\n", "ax1.spines[\"bottom\"].set_visible(False)\n", "\n", "ax2.set_ylim(-0.01, 0.5) \n", "ax2.set_xlim(-5, 10)\n", "ax2.axhline(y=0, color='k', linewidth = 0.6)\n", "ax2.set_xticks(range_y)\n", "ax2.set_yticks(ypmf_values)\n", "ax2.spines[\"top\"].set_visible(False) \n", "ax2.spines[\"right\"].set_visible(False)\n", "ax2.spines[\"bottom\"].set_visible(False)\n", "\n", "# \n", "ax1.scatter(np.array([mean_x]),np.zeros(1), color =\"red\", s=150, zorder=3)\n", "ax1.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20, zorder=2)\n", "ax1.bar(range_x, xpmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, zorder=1)\n", "ax1.set_title(r\"$\\mu_X=${:.2f}, $\\sigma^2_X\\approx${:.3f}\".format(mean_x, var_x))\n", "\n", "\n", "# \n", "ax2.scatter(np.array([mean_y]),np.zeros(1), color =\"red\", s=150, zorder=3)\n", "ax2.scatter(range_y,np.zeros(range_y.shape), color =\"red\", s=20, zorder=2)\n", "ax2.bar(range_y, ypmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, zorder=1)\n", "ax2.set_title(r\"$\\mu_Y=${:.2f}, $\\sigma^2_Y\\approx${:.3f}\".format(mean_y, var_y))\n", "\n", "plt.show();\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Negative binomial distribution\n", "\n", "\n", "**Definition**\n", "\n", "
\n", "\n", "Let $X$ be the number of independent Bernoulli trials with success probability $p$ needed until exactly r successes occur. Then we say $X$ has the **negative binomial distribution** with parameters $(r,p)$.\n", "
\n", "\n", "\n", "* $\\text{range}(X)=\\{r,r+1,\\dots\\}.$\n", "*If $X(s)=x$, means there are $r-1$ successes and $x-r$ fails in the first $x-1$ trials and the $x-th $ trial is a success. Therefore\n", "\n", "$$f(x)={x-1\\choose r-1}p^{r}q^{x-r}.$$\n", "\n", "* When $r=1$, the negative binomial distribution is the same as the geometric distribution. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "29a0363f73f840da9fd6de906a741846", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=30, description='r', max=50, min=1, readout_format=''), FloatSlider(valu…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from ipywidgets import interact, FloatSlider, IntSlider\n", "from scipy.special import comb\n", "\n", "def negbinomial_pmf(r, p):\n", " q = 1-p\n", " N=100\n", " range_x = np.arange(r, N, 1)\n", "\n", " def negbin_pmf(x):\n", " pmf_val = comb(x-1, r-1, exact=True) * p**r * q**(x-r)\n", " return pmf_val\n", " mean = r/p\n", "\n", " pmf_values = np.array([negbin_pmf(x) for x in range_x])\n", "\n", " # plot setup\n", " plt.figure(figsize=(14,7)) \n", " plt.axhline(y=0, color='k', linewidth = 0.6)\n", " plt.ylim(-0.01, max(np.max(pmf_values)+0.05, 0.1))\n", " plt.xlim(np.min(range_x)-2, N+1)\n", " plt.gca().spines['top'].set_visible(False)\n", " plt.gca().spines['right'].set_visible(False)\n", " plt.gca().spines['bottom'].set_visible(False)\n", "\n", " # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n", " plt.scatter(np.array([mean]),np.zeros(1), color =\"red\", s=200,zorder=3)\n", " plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20, zorder=2)\n", " plt.bar(range_x, pmf_values,width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, zorder=1)\n", " plt.title(\"Negative binomial distribution\")\n", " plt.figtext(0.8,0.8, \" r={} \\n p={} \".format(r,p), ha=\"left\", va=\"top\",\n", " backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n", " plt.plot();\n", "\n", "# create interactive variables\n", "r = IntSlider(min=1, max=50, step=1, value=30, readout_format='')\n", "p = FloatSlider(min=0.01, max=1, step=0.01, value=0.5, readout_format='')\n", "\n", "# display the interactive plot\n", "interact(negbinomial_pmf, r=r, p=p);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Poisson distribution\n", "\n", "\n", "**Definition**\n", "\n", "
\n", "\n", "A discrete random variable has Poisson distribution with parameter $\\lambda>0$ if \n", "\n", "$$\\text{range}(X)=\\{0,1,2,3,\\dots\\}$$\n", "\n", "and \n", "\n", "$$f(x)=\\frac{e^{-\\lambda}\\lambda^x}{x!},\\quad x\\in\\text{range}(X).$$\n", "\n", "$\\lambda$ is called the **expected rate**.\n", "
\n", "\n", " Notice that $f(x)$ is indeed a pmf:\n", " \n", "$$\\sum_{x\\in \\text{range}(X)}f(x)=\\sum_{x=0}^\\infty \\frac{e^{-\\lambda}\\lambda^x}{x!}=e^{-\\lambda} \\sum_{x=0}^\\infty \\frac{\\lambda^x}{x!}=e^{-\\lambda}e^{\\lambda}=1.$$\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ce4cad0bfde54e178a59f4800324c50a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(FloatSlider(value=6.0, description='$\\\\lambda$', readout_format='', step=0.5), Output())…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy.special import factorial\n", "from ipywidgets import interact, FloatSlider\n", "\n", "\n", "def poisson_pmf(lam):\n", " N=50\n", " range_x = np.arange(0, N, 1)\n", "\n", " def poiss_pmf(x):\n", " pmf_val = np.exp(-lam) * np.power(lam, x) / factorial(x)\n", " return pmf_val\n", " mean = lam\n", "\n", " pmf_values = np.array([poiss_pmf(x) for x in range_x])\n", "\n", " # plot setup\n", " plt.figure(figsize=(14,7)) \n", " plt.axhline(y=0, color='k', linewidth = 0.6)\n", " plt.ylim(-0.01, max(np.max(pmf_values)+0.05, 0.1))\n", " plt.xlim(0, N+1)\n", " plt.xticks(np.arange(0, N+1, 5))\n", " plt.gca().spines['top'].set_visible(False)\n", " plt.gca().spines['right'].set_visible(False)\n", " plt.gca().spines['bottom'].set_visible(False)\n", "\n", " # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n", " plt.scatter(np.array([mean]),np.zeros(1), color =\"red\", s=200, zorder=3)\n", " plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20, zorder=2)\n", " plt.bar(range_x, pmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, zorder=1)\n", " plt.title(\"Poisson distribution\")\n", " plt.figtext(0.8,0.8, r\"$\\lambda$={}\".format(lam), ha=\"left\", va=\"top\",\n", " backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n", " plt.plot();\n", "\n", "# create interactive variables\n", "lam = FloatSlider(min=0.0, max=100, step=0.5, value=6, readout_format='', description=r\"$\\lambda$\")\n", "\n", "# display the interactive plot\n", "interact(poisson_pmf, lam=lam);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Motivation \n", "\n", "The Poisson distribution arises in the following way (made up scenario): grocery route bus can arrive every $\\frac{1}{n}$ hour. It is a very unpredictable bus but we know that there are on average $\\lambda$ buses arriving every hour. Assuming whether a bus arrives or not at given time is independent from other times, this becomes a binomial distribution with parameter $p=\\frac{\\lambda}{n}$. We know that the expect number of successes in a binomial distribution with parameters $(n,p)$ is $np$ and so $p=\\frac{\\lambda}{n}$. Therefore, the probability of $x$ number of buses arriving in an hour will be\n", "\n", "$$f(x)={n\\choose x}p^x(1-p)^{n-x}={n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x\\left(1-\\frac{\\lambda}{n}\\right)^{n-x},\\quad x=0,1,2,\\dots, n. $$\n", "\n", "* Now imagine the bus starts arriving at more irregular times. To accommodate that we can increase the $n$. \n", "* When $n\\to \\infty$ that will correspond to the buses potentially arriving at any moment in the continuous interval. \n", "* In that case\n", "\n", "$$f(x)=\\lim_{n\\to \\infty}{n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x\\left(1-\\frac{\\lambda}{n}\\right)^{n-x}, \\quad x=0,1,2,\\dots.$$\n", "\n", "To compute this limit notice that \n", "\n", "$$\\lim_{n\\to \\infty}{n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x=\\frac{\\lambda^x}{x!}\\lim_{n\\to \\infty}\\frac{n(n-1)\\cdots (n-x+1)}{•n^x}=\\frac{\\lambda^x}{x!}$$\n", "\n", "and (you should know this from calculus)\n", "\n", "$$\\lim_{n\\to \\infty}\\left(1-\\frac{\\lambda}{n}\\right)^{n-x}=e^{-\\lambda}.$$\n", "\n", "Consequently, if the buses can potentially arrive any moment within the hour, with $\\lambda$ buses arriving in average at every hour, then the probability that $x$ buses will arrive in a given hour is \n", "$$f(x)=\\frac{e^{-\\lambda}\\lambda^x}{x!},\\quad x=0,1,2,\\dots.$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f3e9a8ae4e5d42de9a755b2ffbca4fce", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=100, description='n', max=200, min=6, readout_format=''), FloatSlider(va…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# nbi:hide_in\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy.special import comb, factorial\n", "from ipywidgets import interact, IntSlider, FloatSlider\n", "\n", "def poiss_binom_pmf(n, lam):\n", " p = lam/n\n", " q = 1-p\n", " N = 50\n", " brange_x = np.arange(0, n, 1)\n", " prange_x = np.arange(0, N, 1)\n", "\n", " def poiss_pmf(x):\n", " pmf_val = np.exp(-lam) * np.power(lam, x) / factorial(x)\n", " return pmf_val\n", " mean = lam\n", "\n", " ppmf_values = np.array([poiss_pmf(x) for x in prange_x])\n", " \n", " def binom_pmf(x):\n", " pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n", " return pmf_val\n", " mean = n*p\n", "\n", " bpmf_values = np.array([binom_pmf(x) for x in brange_x])\n", "\n", " # plot setup\n", " plt.figure(figsize=(14,7)) \n", " plt.axhline(y=0, color='k', linewidth = 0.6)\n", " plt.xlim(-2,20)\n", "\n", " # Plotting hypergeometric\n", " plt.bar(prange_x, ppmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', \n", " linewidth=1.3, label=\"Poisson\")\n", " \n", " # Plotting binomial\n", " plt.bar(brange_x, bpmf_values, width=1, color=(0.3, 0.5, 0.3, 0.1), edgecolor='green',\n", " linewidth=1.3, label=\"Binomial\")\n", " plt.figtext(0.83,0.75, r\" $\\lambda$={}\".format(lam)+\"\\n n = {}\".format(n), ha=\"left\", va=\"top\",\n", " backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n", " plt.gca().spines['top'].set_visible(False)\n", " plt.gca().spines['right'].set_visible(False)\n", " plt.legend()\n", " plt.plot();\n", "\n", "# create interactive variables\n", "lam = FloatSlider(min=0.0, max=100, step=0.5, value=5, readout_format='', description=r\"$\\lambda$\")\n", "n = IntSlider(min=np.floor(lam.value)+1, max=200, step=1, value=100, readout_format='')\n", "\n", "# enforce K<=N and n<=N\n", "def update_range(*args):\n", " n.min = np.floor(lam.value)+1\n", " \n", "lam.observe(update_range, 'value')\n", "\n", "# display the interactive plot\n", "interact(poiss_binom_pmf, lam=lam, n=n);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us see that the Poisson distribution with parameter $\\lambda$ and the binomial distribution with parameters $(n,p=\\frac{\\lambda}{n})$ are close to each other when $n$ is large by looking at their histograms.\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }