{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('fivethirtyeight')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Probability: Mathematical/Theoretical and Computational Approaches\n",
"\n",
"*Dan L. Nicolae* \n",
"\n",
"We will illustrate throughout this textbook how some problems can be solved either mathematically, or computationally (using simulations), or through a combination of mathematical/analytical and computational approaches. Each approach has its strengths and, in this chapter, we use the calculation versus estimation of probabilities to highlight some of them.\n",
"\n",
"We illustrate the concepts we want to introduce with a classic probability exercise called **the birthday problem**. Suppose you and a friend go to a party where there are 30 people (all unknown to both of you) and your friend wants to bet you that there are two people at that party who share their birthday. Would you be willing to take that bet? \n",
"\n",
"Your willingness to take the bet should be related to the chance of winning the bet. What do you think it is more likely to happen: finding a pair with shared birthdays or having 30 distinct birthdays? \n",
"\n",
"We will answer this question using the language of probability; we will calculate the probability of the event that at least two people share birthdays in a group of 30 people. In the following sections, we will introduce the rules we need for deriving this probability and then show how to estimate them using simulations.\n",
"\n",
"Let's start with a simpler problem: what is the probability that two people share their birthday. Can you justify the following result?\n",
"\n",
"$$P(\\mbox{two random people have the same birthday}) ~=~ \\frac{1}{365}$$\n",
"\n",
"Think about the assumptions you implicitly or explicitly made in your justification. \n",
"\n",
"We will show in the next section that, given a number of people, **n** (with $2\\leq n\\leq 365$), the probability, $P_n$, that at least two share a birthday is given by:\n",
"\n",
"$$P_n ~=~ 1-\\frac{365\\times364\\times ...\\times (365-n+1)}{365^n}$$\n",
"\n",
"which can also be written as\n",
"\n",
"$$P_n ~=~ 1-\\frac{365}{365}\\times\\frac{364}{365}\\times\\frac{363}{365}\\times ...\\times \\frac{(365-n+1)}{365}$$\n",
"\n",
"The asssumptions used to obtain the above formula for $P_n$ are:\n",
"\n",
"a. 365 days in a year;\n",
"\n",
"b. All days are equally likely to be a birthday;\n",
"\n",
"c. Subjects' birthdays are independent of each other (for example, no twins in the room).\n",
"\n",
"The function below calculates these probabilities. Note that, for computational reasons, we implement the second formula for $P_n$: consider how large $365^n$ is for n= 30 or 50."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# a function that calculates the probability for 1\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
"

Number of people

Probability

0

2

0.0027

1

3

0.0082

2

4

0.0164

3

5

0.0271

4

6

0.0405

5

7

0.0562

6

8

0.0743

\n",
""
],
"text/plain": [
" Number of people Probability\n",
"0 2 0.0027\n",
"1 3 0.0082\n",
"2 4 0.0164\n",
"3 5 0.0271\n",
"4 6 0.0405\n",
"5 7 0.0562\n",
"6 8 0.0743"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Construct a data frame with the probabilities for a range of n's\n",
"number_people=np.arange(2,61,1)\n",
"b_probs= np.array([]) # an empty array\n",
"\n",
"for i in number_people: \n",
" b_probs= np.append(b_probs,birthday_prob(i))\n",
"\n",
"Birthday_df=pd.DataFrame(\n",
" {\"Number of people\":number_people,\n",
" \"Probability\":b_probs})\n",
"Birthday_df.head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We construct below a line graph of these probabilities that shows the trend. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAEfCAYAAADGLVhVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvcElEQVR4nO3de1SU1f4/8PeAInhrEmG8AN4gBZJQULwejvSNTt5SFBXtpChqQh4zMS9ZhJZoqOkx5KdOePSoHW9YpJZ2lIMXBE1FLU1HAUWTQSACCoSY+f3hmtGRGWBghmGeeb/Wcq2e/TzPzGcPKz7svT+zH1FRUZESREREAmBl6gCIiIgMhUmNiIgEg0mNiIgEg0mNiIgEg0mNiIgEg0mNiIgEg0mNiIgEg0mNiIgEw+KSmkwmM3UIjYZ9FSb2VZjYV8OwuKRGRETCxaRGRESCYdKkdubMGUyaNAnu7u4Qi8XYtWtXrff89NNPGD58ODp06AB3d3esXr0aSiW3ryQiIhMntd9//x0eHh5YtWoV7Ozsar2+uLgYY8eOhaOjI06cOIFVq1Zh48aN+PzzzxshWiIiauqamfLNAwMDERgYCAAIDw+v9fp9+/ahrKwM8fHxsLOzg4eHB27evIlNmzbh7bffhkgkMnbIRETUhJnVmtq5c+cwcOBAjVHdyy+/jAcPHuDOnTsmjIyIiJoCk47U9JWXl4dOnTpptDk4OKjPde3aVet9z5aPsnRWmNhXYWJfTe9emQgf3rRBUaUI4uZKrHihAp3tlA1ot4X48l11u77c3Nx0njOrpAag2hSjqkikpqnHpz8AmUxW4wciJOyrMLGvwtSYfc0ursSsk0XIf1SF9i2ssdVfjC5tmutsjzj0EFdLKgAAOeXAyrvP4dhIB4O1G5JZTT86OjoiLy9Poy0/Px/AkxEbEZGlyS6uROChh+h7IBeBhx7iTkllje2zThbh3MMKZBZX4dzDCsxMKaqxPf9Rlcb7qY4N1W5IZpXU+vfvj7Nnz6K8vFzdlpycjI4dO6JLly4mjIyIqHFoS1TGTlLtW1hrtKuODdVuSCZNaqWlpbhy5QquXLkChUKBe/fu4cqVK8jJyQEAREdHY/To0errx48fDzs7O4SHh+PatWtISkrC+vXrER4ezspHIhIUfUZZxk5SW/3F6O9gg+5trdHfwQZb/cUNane2VWi0G5JJ19QuXbqEUaNGqY9jYmIQExODkJAQxMfHIzc3F1lZWerzzz33HA4ePIjIyEgMGzYMYrEYERERePvtt00RPhFRg6nWsR6U2KLjjYfqdSxV8gKATFRhZkoRjo100Jqo2rewRiaetD+dpLS1b/UXq5Ohau2spvYubZprXfuqb/vj9UNnPT6lujNpUhs6dCiKiop0no+Pj6/W5unpiW+//daIURERGZ6uIownycsKOeUVNSYvQHuiMnaSMidmV/1IRNSU1Z68ah95AfqNsoScpPTFpEZEVA+mSF6AZSYqfTCpERHVoLGS14OScnRsY8vk1UBMakRENWiskZcxiycsCZMaERF0j8g4bWhemNSIiKB7RMbkZV6Y1IjIoug7ImPyMi9MakRkUfQdkTF5mRcmNSISLG2jMn1HZGRemNSISLC0jco4IhM2JjUiMnv6rJMdDLTniEzAmNSIyOzps07GEZmwMakRkdkwVOUiCReTGhGZDVYuUm3M6snXRGTZahqRaXsoJVkejtSIqMnR9eBMjsioNhypEVGTo5pmzCm3wrmHjx+cCXBERrXjSI2ITEbfwg+OyKg2HKkRkcmoRmSZxVUaIzLVtKLKs8dEujCpEZHJ1Fb44Wyr4DQj6YXTj0RkdLqmGWsr/OCDM0lfHKkRkdHpmmZk4QcZGkdqRGR0LPygxsKRGhEZHQs/qLFwpEZEBqNr7Yx7MFJjYVIjIoPRtTcjpxmpsXD6kYgMRtfaGVFj4UiNiPSmb4k+UWPhSI2I9MYSfWqqOFIjIr2xRJ+aKo7UiEhvLNGnpoojNSLSiSX6ZG5MPlKTSqXw8vKCRCKBv78/UlNTa7z++PHjeOWVV+Dk5ITu3bsjJCQEt27daqRoiSyLrrUz1TTjxXEd1CX7RE2BSZNaYmIiFi9ejAULFuDkyZPo378/goODkZOTo/X67OxsTJ48GQMHDsTJkyfx1Vdfoby8HMHBwY0cOZFlYIk+mRuTJrW4uDhMnjwZU6dORc+ePREbGwuJRIKEhASt11++fBmVlZWIiopC9+7d4eXlhfnz5yMrKwsFBQWNHD2R8HHtjMyNyZJaRUUFMjIyEBAQoNEeEBCA9PR0rfd4e3ujefPm2LFjB6qqqlBSUoIvv/wSffv2hb29fWOETSRY2cWVCDz0EH0P5CLw0EPcKalkiT6ZHVFRUZHSFG/84MEDuLu74/Dhwxg8eLC6ffXq1di3bx9++OEHrfelpqZi2rRpKCgogEKhgJeXF/bv3w8HB91lxDKZzODxEwnN9MstcLXkyUjMq00VvnjpkQkjItLOzc1N5zmTVz+KRCKNY6VSWa1NRS6XY+7cuZg0aRLGjRuH0tJSrFy5EtOmTcM333wDKyvtA8+nP4DHDx3U/YEICfsqTMbq6+9XcoGndgMpFdnAzc3F4O+jD/5chcmYfTVZUrO3t4e1tTXy8vI02vPz83WOurZu3YqWLVti+fLl6rYtW7bA09MT6enpGDhwoFFjJhIybnFFQmCyNTUbGxt4e3sjOTlZoz05ORl+fn5a7ykrK4O1teb/aKpjhUJhnECJBEbb2hnALa5IGExa/RgREYHdu3djx44duHHjBhYtWoTc3FyEhoYCAKKjozF69Gj19YGBgbh8+TJWrVqF27dvIyMjAxEREXBycoK3t7eJekFkXvjdMxIyk66pBQUFobCwELGxsZDL5XB3d8fevXvh4vJ4Hj83NxdZWVnq6/39/SGVSrFhwwZs3LgRtra28PX1xf79+9GqVStTdYPIrPC7ZyRkJi8UCQsLQ1hYmNZz8fHx1drGjRuHcePGGTssIsHi2hkJmcm3ySKixsW1MxIyk4/UiMg4dG1GzMfDkJBxpEYkULoKQoiEjEmNSKBYEEKWiEmNSKC4GTFZIiY1IoFiQQhZIhaKEJk5FoQQPcGRGpGZY0EI0RNMakRmjgUhRE8wqRGZORaEED3BpEZk5lgQQvQEC0WIzISqIORBiS063njIghAiLThSIzITqoKQnHIrFoQQ6cCkRmQmWBBCVDsmNSIzwYIQotoxqRGZCVVBiLOtggUhRDrondSkUil+/fVXY8RCRDVQFYQk+pbj2EgHdGnT3NQhETU5eie1hQsXolevXpgyZQqSkpJQUVFhjLiILFZ2cSUCDz1E3wO5CDz0EHdKKk0dEpHZ0DuppaWlITw8HFeuXMHUqVPxwgsvYP78+UhLSzNGfEQWh9teEdWf3kmtZ8+eiIqKwtWrV5GUlIRRo0YhMTERw4cPh7e3N2JiYnD79m1jxEpkEVjlSFR/DSoUGTp0KDZu3AiZTIZt27bB1dUVsbGx6NevH1555RVs27YNf/zxh6FiJbIIrHIkqj+DVD/+8MMPOHHiBM6fPw+lUgkvLy+IRCK8++67eOmll3Dy5ElDvA2RReC2V0T1V+9tsm7evIk9e/Zg3759uHfvHiQSCaZOnYqQkBC4u7sDAK5du4Y5c+Zg4cKFSE9PN1jQRELA56ARGZ7eSS0+Ph579uzBlStX0KJFCwwfPhzr1q1DQEAArKw0B34eHh6YPXs25s6da7CAiYRCVRACAJmowsyUIiYzogbSO6ktXboUfn5++OyzzzB27Fi0bdu2xuv79OmDhQsX1jtAIqFiQQiR4emd1C5cuIDu3bvX+Xp3d3f1dCQRPdG+hTUyUaVxTEQNo3ehyLx585CSkqLz/MmTJzFq1KgGBUVkCVgQQmR4eo/UTp8+jTfffFPn+fz8fJw5c6ZBQRFZAhaEEBmewR8Sev/+fbRq1crQL0tktnRVORKR4dUpqR0+fBhHjhxRH//rX//C//73v2rXFRUVISUlBT4+PgYLkMjcscqRqPHUKaldv34dBw4cAACIRCKcP38eFy5c0LhGJBKhZcuWGDBgAFatWmX4SInMFKsciRpPnZJaZGQkIiMjAQDPP/884uLiEBwcbNTAiISCVY5EjUfv6sdff/3VoAlNKpXCy8sLEokE/v7+SE1NrfF6pVKJTZs2oV+/fnB0dETPnj3x0UcfGSweIkNjlSNR4zF4oYg+EhMTsXjxYqxduxYDBgyAVCpFcHAw0tLS4OzsrPWe999/H0ePHsXy5cvh6emJ3377DXK5vJEjJ6o7VjkSNZ5ak5qXlxesrKxw/vx5NG/eXL1ZcU1EIhEyMjJqffO4uDhMnjwZU6dOBQDExsbi+PHjSEhIQFRUVLXrZTIZtmzZgjNnzqBnz561vj5RY2KVI5Hp1ZrUBg8eDJFIpN7XUXXcUBUVFcjIyKi2L2RAQIDOzY+PHDmCrl274r///S8mTJgAhUKBwYMHY8WKFXBw4F/CZFqsciQyvVqTWnx8fI3H9VVQUICqqqpqycjBwQF5eXla78nOzkZOTg4SExOxadMmiEQifPDBB5g0aRK+//77ahsqq8hkshqPhYx9bTwPSmzx9DL1g5Jyo8Vk6r42JvZVmBrSVzc3N53nTLqmBqDaqE+pVOocCSoUCjx69AibN2+Gq6srAGDz5s3w9fXFxYsX4evrq/W+pz8AmUxW4wciJOxr4+p44yFyyiueHLexhZub9rXhhmgKfW0s7KswGbOvtSa1nJycer2wrkIPFXt7e1hbW1cbleXn5+ucSpRIJGjWrJk6oQFAjx490KxZM9y7d09nUiNqDFv9xZiZormmRkSNq06FIvVZQyssLKzxvI2NDby9vZGcnIwxY8ao25OTkzF69Git9wwYMAB//vknsrKy0K1bNwCPpyT//PPPWpMokbGxypHI9GpNap9//rlBCkO0iYiIwOzZs+Hj4wM/Pz8kJCQgNzcXoaGhAIDo6GhcuHABSUlJAIC//vWveOmllxAREYGYmBgAwJIlS+Dr64s+ffoYJUYiIjIftSa1KVOmGO3Ng4KCUFhYiNjYWMjlcri7u2Pv3r1wcXEBAOTm5iIrK0t9vZWVFfbs2YNFixZhxIgRsLW1xbBhw/DJJ5/oLBIhMgaW7xM1TSYvFAkLC0NYWJjWc9oqLTt06IDt27cbOyyiGrF8n6hpqjWpffnllwCASZMmQSQSqY9rExIS0rDIiJowblJM1DTVmtTCw8MhEokwbtw42NjYIDw8vNYXFYlETGokaNykmKhpqjWpXb58GcDjasWnj4ksGcv3iZqmWpOaqmhD1zGRJWL5PlHT1KBCkR9//FH95WxnZ2d4enoarfyfyBRY5UhkXuqV1A4cOICoqCj88ssvGu2dOnVCVFQUHyBKgsEqRyLzondS27VrF95++224ubkhOjoarq6uUCqVuH37Nnbs2IHZs2ejoqLCqN9vI2osrHIkMi96J7V169bBx8cHhw4dgq2trca5mTNnYvjw4Vi3bh2TGgkCqxyJzIve23Dcv38fwcHB1RIaANja2mLixInVpiWJzNVWfzH6O9ige1tr9HewYZUjUROn90itV69eePDggc7zv/zyC59KTYLBKkci86L3SG358uXYvn07Dh48WO3cgQMHsGPHDqxYscIgwREREemj1pGatkpGe3t7zJgxA4sXL0a3bt0gEomQmZmJhw8fokePHti4cSOGDh1qlICJjIGl+0TCUGtS+/nnn7V+98zJyQkA1OtnLVq0gJOTEx49eoQbN24YOEwi42LpPpEw1JrUrl692hhxEJkUS/eJhIEPISNC9VJ9lu4TmacGbZNVUlKC4uJiKBSKauecnZ0b8tJEjYobFBMJQ72S2o4dO/DPf/4TmZmZOq8pLCysd1BEjY2l+0TCoPf047///W/MmzcPzs7OWLZsGZRKJebMmYP58+fD0dERvXv3xsaNG40RKxERUY30Tmrx8fEYOnQoDh48iGnTpgEAAgMD8cEHHyAtLQ1FRUUoLi42dJxEBpFdXInAQw/R90AuAg89xJ2SSlOHREQGpHdSy8zMxMiRIx/fbPX49srKx78YxGIx3nzzTUilUgOGSGQ4qtL9zOIqnHtYgZkpRaYOiYgMSO+k1qpVKyiVSgBA69atYW1tjdzcXPX5du3ace9HarJYuk8kbHonNTc3N1y7dg0A0KxZM/Tu3Rv/+c9/UFlZifLycuzZswddunQxeKBEhsDSfSJh0zupjRgxAt9//z3Ky8sBAJGRkUhNTUXXrl3h6uqK9PR0zJ8/3+CBEhkCd90nEja9S/rnzp2LuXPnqo9HjBiBI0eO4KuvvkKzZs3wt7/9DUOGDDFokESGwtJ9ImFr0JevVQYMGIABAwYY4qWIDIIbFBNZpnontZKSEpw6dQo5OTkAABcXFwwePBht27Y1WHBE9cUNioksU72S2meffYY1a9agrKxMXQkJAHZ2dliwYAEWLFhgsACJ6oNVjkSWSe+ktmHDBixfvhxDhgxBWFgYXF1doVQqcfv2bUilUnzyySdo1qwZ5s2bZ4x4ieqkfQtrZKJK45iIhE/vpLZ161YMGzYMiYmJGu0vvvgiRo8ejbFjx2Lr1q1MamRS3KCYyDLpndQKCwsxfPhwredEIhFGjhyJDz/8sMGBETUEqxyJLJPe31N76aWX8PPPP+s8f/36dXh7ezckJiIionrRO6nFxsbim2++wYYNG1BaWqpuLy0txfr163H48GHExsbW+fWkUim8vLwgkUjg7++P1NTUOt13+/ZtODk5oXPnzvp2gQREtUFx0A+23KCYiGqffvTz86vWJhKJEB0djRUrVsDR0REikQhyuRwKhQKOjo6YMWMG0tLSan3zxMRELF68GGvXrsWAAQMglUoRHByMtLS0Gh8yWlFRgenTp2PQoEE4c+ZMre9DwvWkdN8KOeUVLN0nsnC1JrX27dtDJBJptDk4OMDV1VWjrVu3bnq/eVxcHCZPnoypU6cCeDwKPH78OBISEhAVFaXzvqioKHh6emLw4MFMahaOpftE9LRak9rhw4eN8sYVFRXIyMjQ2HILAAICApCenq7zvqNHj+Lo0aNISUlBUlKSUWIj88HSfSJ6mt5raoZSUFCAqqoqODhoThU5ODggLy9P6z25ubmYN28eNm/ejDZt2jRGmNTEqTYodrZVcINiIqrfjiJVVVXYvXs3jh07hrt37wJ4vE3Wq6++ipCQEFhb1/2v5WenNpVKZbU2lVmzZmH69Ono16+fXvHKZLIaj4XMEvoa11P1X+WoyC2CLLemq4XBEn6uKuyrMDWkr25ubjrPiYqKipQ6z2pRXFyMoKAgXLx4Ea1bt0bXrl2hVCpx584dlJaWwsfHB4mJibWOpCoqKtCxY0d88cUXGDNmjLo9MjIS165dw5EjR6rdIxaLNRKmUqmEQqGAtbU11q5di2nTptUav0wmq/EDERL2VZjYV2FiXw1D75Haxx9/jEuXLmHlypWYPn06bGxsAACVlZVISEjA0qVL8fHHH2P16tU1vo6NjQ28vb2RnJyskdSSk5MxevRorfc8W+5/5MgRrF27FsePH0enTp307QqZEe66T0R1oXdSO3ToEEJDQ/HWW29ptDdv3hyzZ8/GzZs38c0339Sa1AAgIiICs2fPho+PD/z8/JCQkIDc3FyEhoYCAKKjo3HhwgV1QYiHh4fG/ZcuXYKVlVW1dhIe7rpPRHWhd1IrKCiAu7u7zvMeHh7YuXNnnV4rKCgIhYWFiI2NhVwuh7u7O/bu3QsXFxcAjwtDsrKy9A2RBIil+0RUF3pXPzo7OyM5OVnn+eTk5Bq/OP2ssLAwXL16FXl5eUhJScHgwYPV5+Lj43H16lWd906ZMgX379+v83uR+Xq2VJ+l+0Skjd5J7Y033sDhw4cxZ84cXL9+HZWVlaisrMS1a9cQERGBI0eO4M033zRGrGTBVKX73dtas3SfiHTSe/px3rx5uHPnDv71r39hz5496vJ7pVIJpVKJ0NBQ/OMf/zB4oGTZuOs+EdWF3klNJBLhs88+w6xZs3D06FGN76kFBgayaIOIiExGr6RWVlaGCRMmYOLEiXjjjTdqLBghqg+W7hNRQ+i1pmZnZ4fLly+jqoqVZ2QcqtL9zOIqnHv4eNd9IqK60rtQZMiQIXV+5hmRvli6T0QNoXdSW716NS5evIgPPvgA2dnZUCgUxoiLLBRL94moIfQuFOnXrx8UCgXi4uIQFxcHKysrNG+uueYhEonwyy+/GCxIshxb/cWYmaK5pkZEVFd6J7WgoCBjxEEEgKX7RNQwdU5qjx49wpEjR+Dm5oZ27drh1VdfRYcOHYwZGxERkV7qlNTkcjmGDx+OrKws9fPOWrZsiT179mhsa0VUVyzdJyJjqFOhyMcff4zs7GyEh4djz549WLlyJVq0aIH33nvP2PGRQLF0n4iMoU4jtRMnTiAkJAQff/yxus3R0RFhYWG4f/8+OnfubLQASZhYuk9ExlCnkZpcLoefn59G24ABA6BUKnHv3j2jBEbCxtJ9IjKGOiW1qqoq2NraarSpjsvLyw0fFQked90nImOoc/VjdnY2Lly4oD4uLi4GAMhkMrRu3bra9T4+PgYIj4SKpftEZAx1TmoxMTGIiYmp1v5ssYiqOrKwsLDh0REREemhTkktLi7O2HGQgLF8n4gaS52S2uTJk40dBwmYqnwfADJRhZkpRZx6JCKj0HtDYyJ9sXyfiBoLkxoZHcv3iaixMKmR0bF8n4gai9679BPpi+X7RNRYmNTIYFjlSESmxulHMhhuUkxEpsakRgbDKkciMjUmNTIYVjkSkakxqZHBsMqRiEyNhSJkMKxyJCJT40iNiIgEgyM10htL94moqeJIjfTG0n0iaqpMntSkUim8vLwgkUjg7++P1NRUndeeOnUKISEh6NmzJzp27IhBgwbh3//+dyNGSwBL94mo6TJpUktMTMTixYuxYMECnDx5Ev3790dwcDBycnK0Xn/u3Dl4enpi+/btOHv2LGbMmIF33nkH+/bta+TILRtL94moqTLpmlpcXBwmT56MqVOnAgBiY2Nx/PhxJCQkICoqqtr1CxYs0DieMWMGTp06haSkJAQHBzdKzPS4dH9miuaaGhFRU2CypFZRUYGMjAzMnTtXoz0gIADp6el1fp2SkhJ06tTJ0OFRDVi6T0RNlcmSWkFBAaqqquDgoPnL0cHBAXl5eXV6je+++w4pKSk4evRojdfJZLIaj4WsIX29VybChzdtUFQpgri5EiteqEBnO6UBozMs/lyFiX0Vpob01c3NTec5k5f0i0QijWOlUlmtTZu0tDTMnDkTq1evho+PT43XPv0ByGSyGj8QIWloXyMOPcTVkgoAQE45sPLuc012hMafqzCxr8JkzL6arFDE3t4e1tbW1UZl+fn51UZvzzp79iyCg4OxZMkSzJgxw5hhWjRWORKRuTFZUrOxsYG3tzeSk5M12pOTk+Hn56fzvjNnziA4OBjvvfcewsPDjR2mRWOVIxGZG5OW9EdERGD37t3YsWMHbty4gUWLFiE3NxehoaEAgOjoaIwePVp9/alTpxAcHIzQ0FBMmDABcrkccrkc+fn5puqCoHGDYiIyNyZdUwsKCkJhYSFiY2Mhl8vh7u6OvXv3wsXFBQCQm5uLrKws9fW7d+/GH3/8gY0bN2Ljxo3qdmdnZ1y9erXR4xc6VjkSkbkxeaFIWFgYwsLCtJ6Lj4+vdvxsGzUc93IkIqEw+TZZZHrcy5GIhIJJjVjlSESCwaRGrHIkIsFgUiNWORKRYJi8UIQaj66CEFY5EpFQcKRmQVgQQkRCx6RmQVgQQkRCx6RmQVgQQkRCx6RmQVgQQkRCx0IRAVIVhDwosUXHGw9ZEEJEFoMjNQFSFYTklFuxIISILAqTmgCxIISILBWTmgCxIISILBWTmgCpCkKcbRUsCCEii8JCETNW2w4hMpkMbm7Opg6TiKjRcKRmxrhDCBGRJiY1M8aCECIiTUxqZowFIUREmrimZgZ0rZ1t9RdjZopmOxGRJWNSMwOqtTMAyEQVZqYU4dhIB+4QQkT0DE4/mgGunRER1Q2Tmhng2hkRUd1w+rEJ4doZEVHDMKk1IVw7IyJqGE4/NiFcOyMiahgmtSaEa2dERA3D6UcT0bZ+xrUzIqKGYVIzEV3rZ1w7IyKqPyY1I9NV0cj1MyIiw+OampHp2kmf62dERIbHpGZkukZkqgd5dm9rzQd5EhEZCKcfDUTXNGP7FtbIxJPEphqR8btnRESGZ/KRmlQqhZeXFyQSCfz9/ZGamlrj9T/99BOGDx+ODh06wN3dHatXr4ZSqWykaHXTNc3IERkRUeMx6UgtMTERixcvxtq1azFgwABIpVIEBwcjLS0Nzs7O1a4vLi7G2LFjMWjQIJw4cQIymQwRERFo2bIl5s6d2ygx61v4wREZEVHjMelILS4uDpMnT8bUqVPRs2dPxMbGQiKRICEhQev1+/btQ1lZGeLj4+Hh4YHXX38d8+bNw6ZNmww+WssurkTgoYfoeyAXgYce4k5JJQAWfhARNWUmS2oVFRXIyMhAQECARntAQADS09O13nPu3DkMHDgQdnZ26raXX34ZDx48wJ07dwwan67kxcIPIqKmy2TTjwUFBaiqqoKDg+bUnIODA/Ly8rTek5eXh06dOlW7XnWua9euWu+TyWQ1HmvzoMQWT+f8ByXlkMlkaKVsAeDJKKyVskL9enE9n9xfkVsEWW6tb2N0demrULCvwsS+ClND+urm5qbznMmrH0UikcaxUqms1lbb9dran/b0ByCTyWr8QFQ63niInPKKJ8dtbOHm5oydHSqrbWXVpU3zWl/PFOraVyFgX4WJfRUmY/bVZEnN3t4e1tbW1UZl+fn51UZvKo6OjlqvB6DznvrStQ8jCz+IiJouk62p2djYwNvbG8nJyRrtycnJ8PPz03pP//79cfbsWZSXl2tc37FjR3Tp0sWg8amS18VxHdTPNCMioqbNpNWPERER2L17N3bs2IEbN25g0aJFyM3NRWhoKAAgOjoao0ePVl8/fvx42NnZITw8HNeuXUNSUhLWr1+P8PDwGqcfiYjIMph0TS0oKAiFhYWIjY2FXC6Hu7s79u7dCxcXFwBAbm4usrKy1Nc/99xzOHjwICIjIzFs2DCIxWJERETg7bffNlUXiIioCTF5oUhYWBjCwsK0nouPj6/W5unpiW+//dbYYRERkRky+TZZREREhsKkRkREgiEqKioy/W7AREREBsCRGhERCQaTGhERCQaTGhERCQaTGhERCQaTGhERCYZFJTWpVAovLy9IJBL4+/sjNTXV1CE12JkzZzBp0iS4u7tDLBZj165dGueVSiViYmLQq1cvdOjQASNGjMD169dNFG39rVu3DsOGDYOzszN69OiBiRMn4tq1axrXCKWvW7duxaBBg+Ds7AxnZ2e88sorOHr0qPq8UPqpzdq1ayEWi7Fw4UJ1m1D6GxMTA7FYrPHvhRdeUJ8XSj9VcnNz8dZbb6FHjx6QSCTw8/PD6dOn1eeN1V+LSWqJiYlYvHgxFixYgJMnT6J///4IDg5GTk6OqUNrkN9//x0eHh5YtWqVxsNTVTZs2IC4uDisXr0aJ06cgIODA8aOHYuSkhITRFt/p0+fxowZM3D06FEkJSWhWbNmGDNmDH799Vf1NULpa6dOnRAdHY2UlBQkJyfjL3/5C6ZMmYIff/wRgHD6+azz589j+/bt8PT01GgXUn/d3Nxw48YN9b+n/7AWUj+Liorw6quvQqlUYu/evUhPT8enn36q8TQVY/XXYr6n9vLLL8PT0xP//Oc/1W19+/bF66+/jqioKBNGZjidO3fGp59+iilTpgB4/JdQr169MHPmTERGRgIAysrK4ObmhhUrVqg3jjZHpaWlcHFxwa5du/Daa68Juq8A0LVrV0RFRWHatGmC7Odvv/0Gf39/bNiwAZ9++ik8PDwQGxsrqJ9rTEwMkpKScPbs2WrnhNRPAFi+fDnOnDmjMcPwNGP21yJGahUVFcjIyEBAQIBGe0BAANLT000UlfHduXMHcrlco992dnYYNGiQ2fe7tLQUCoUCYrEYgHD7WlVVhQMHDuD3339H//79BdvPd955B6+//jr8/f012oXW3+zsbLi7u8PLywvTp09HdnY2AOH18/Dhw/Dx8UFoaChcXV0xZMgQbNmyRf1QZ2P21+QbGjeGgoICVFVVVXuQqIODQ7WHjgqJXC4HUP0Bqg4ODnjw4IEpQjKYxYsXo3fv3ujfvz8A4fX1p59+QmBgIMrLy9GqVSvs3LkTnp6e6v/hhdJPANi+fTsyMzOxefPmaueE9HP19fXFpk2b4Obmhvz8fMTGxiIwMBBpaWmC6ifwOHl/8cUXCA8PxzvvvIOrV69i0aJFAIBZs2YZtb8WkdRUnn3mmlKptIjnsAmt30uXLkVaWhq+++47WFtba5wTSl/d3Nxw6tQp/Pbbb0hKSsKcOXNw6NAh9Xmh9FMmk2H58uX49ttvYWNjo/M6IfT3lVde0Tj29fWFt7c3du/ejX79+gEQRj8BQKFQoE+fPuqlnZdeegmZmZmQSqWYNWuW+jpj9Nciph/t7e1hbW1dbVSWn59f7S8FIZFIJAAgqH4vWbIEBw4cQFJSErp27apuF1pfbWxs0L17d/Uvht69e2PTpk2C6+e5c+dQUFCAgQMHwt7eHvb29jhz5gykUins7e3Rrl07AMLp79Nat26NXr16ITMzU3A/V4lEgp49e2q0vfDCC7h37576PGCc/lpEUrOxsYG3tzeSk5M12pOTk+Hn52eiqIyvS5cukEgkGv0uLy/H2bNnzbLfixYtwv79+5GUlKRRCg0Ir6/PUigUqKioEFw/R4wYgdTUVJw6dUr9r0+fPhg3bhxOnToFV1dXQfX3aeXl5ZDJZJBIJIL7uQ4YMAC3bt3SaLt16xacnZ0BGPf/V4uZfoyIiMDs2bPh4+MDPz8/JCQkIDc31+yqip5VWlqKzMxMAI9/8d27dw9XrlzB888/D2dnZ8yZMwdr166Fm5sbXF1dsWbNGrRq1Qrjx483ceT6iYyMxJ49e7Bz506IxWL1nHyrVq3QunVriEQiwfT1o48+QmBgIDp37ozS0lLs378fp0+fxt69ewXVTwDq72s9rWXLlnj++efh4eEBAILp77Jly/C3v/0NTk5O6jW1P/74AyEhIYL7uYaHhyMwMBBr1qxBUFAQrly5gi1btuCDDz4AAKP212KSWlBQEAoLCxEbGwu5XA53d3fs3bsXLi4upg6tQS5duoRRo0apj2NiYhATE4OQkBDEx8dj3rx5KCsrw8KFC1FUVAQfHx8kJiaiTZs2Joxaf1KpFADw+uuva7QvWrQIS5YsAQDB9FUul2PWrFnIy8tD27Zt4enpif379+Pll18GIJx+1pVQ+vvLL78gLCwMBQUFaN++PXx9ffH999+rfwcJpZ/A469L7dq1C8uXL0dsbCycnJywdOlShIWFqa8xVn8t5ntqREQkfBaxpkZERJaBSY2IiASDSY2IiASDSY2IiASDSY2IiASDSY2IiASDSY2oAcRiMebPn2/qMOosKysL48ePR5cuXbQ+VNYcicVixMTEmDoMaiKY1KhJ27VrF8RiMRwdHdX7xj1t3Lhx6N27twkiM09z587FxYsXsXjxYmzevBmDBw82dUhEBsWkRmahoqIC69atM3UYZq2qqgpnz57FhAkTMGfOHEycOFFjU2giIWBSI7PQu3dv7Ny5U+toTeiUSiXKy8sb/DqFhYWoqqrCc889Z4CoiJomJjUyC++++y4A1Dpau3Pnjs61ot69e2POnDnqY9XU5unTp7F06VK4urrCxcUFERERKC8vx++//4533nkH3bt3h4uLCyIjI/Hnn39qfd/ExET4+flBIpFg0KBBWh9jX1xcjGXLlqF3795wdHTEiy++iI8++giPHj3SuE61TvfVV19h0KBBcHR0xIEDB2rs99mzZzFq1Ch07twZTk5OGDNmDH744Qf1+ZiYGLi5uQEAVq9erXUj4Wep4jBk3xQKBdavXw8fHx84OjrC3d0dCxcuxG+//aZx3YgRI9CvXz9cvXoVr732Gjp27AhPT0+sX7++xpj1jYeEx2I2NCbz5uTkhMmTJ2Pnzp1499134eTkZLDXXrJkCdq3b49FixYhIyMDu3btQsuWLZGdnQ07Ozu8//77OHnyJKRSKbp3747w8HCN+9PT03Hw4EHMnj0brVu3xvbt2zFlyhR8/fXX6jWrsrIyjBw5Enfu3MG0adPQrVs3XL16FZ9//jlu3ryJ3bt3a7zm2bNn8fXXX2PmzJmQSCTVHrXztDNnzmDs2LHo1KkTIiMjoVAosG3bNowYMQKHDx+Gr68vRo0ahfbt22PhwoUYOXKkxibYNTF03xYsWIBt27bhtddew1tvvYXr16/jiy++wIULF3D06FE0b95cfW1xcTHGjRuHkSNHYuzYsThy5Ag++ugjVFVVYcGCBTpj1vezJmFhUiOzsWDBAuzevRvr1q0z6Pqavb09EhMT1U/cvXv3LqRSKYKDg7FlyxYAwIwZM+Dn54edO3dWS2rXrl3D0aNH1c+BmjJlCvr27Yvo6GgcO3YMALBp0ybIZDL873//03h4oru7OyIjI5GamopBgwap22/cuIGUlBR4eXnVGv/777+PVq1a4b///S/at28PAAgJCUH//v2xbNkyfPfdd3jxxRfh4OCAhQsXwtPTExMnTqzTZ2PIvl27dg3btm3DhAkT1J8r8Pgp30uWLMGXX36JN998U90ul8vx4YcfqkfpYWFhGD16NNasWYOwsDCd06j6ftYkLJx+JLPh7OysHq0Zcm3tjTfe0HiEvK+vL5RKJf7+979rXOfj44OsrKxq9/fp00fjwYbt2rVDcHAwzp07h6KiIgDAwYMH4efnh/bt26OgoED9769//SsA4OTJkxqv6efnV6eEJpfLkZGRgZCQEHVCA4BOnTph/PjxSE9PV8dQH4bsm2ra8h//+IfGe0yfPh1t27atNq1pZWWl8agSKysrzJw5E2VlZTh16pTOmPX9rElYOFIjs2KM0dqzU5lt27bV2V5WVoZHjx6hRYsW6vYePXpUe01VW05ODsRiMW7fvo0ff/xR67XA48fYP62uVYl3794FAK3Tkz179oRSqVTHUB+G7Nvdu3chEonUa3sqLVq0QJcuXdR9UXF0dFT/LLS9ty76ftYkLExqZFaeHq2ppqWe9vSI61kKhUJru7W1tdZ2KyvtExlKpeYjCLW957PXKBQK/OUvf9EaM/B4ZPU0Ozs7rdfp49kY6sMYfdNGqVRWe6+6vLc2hoiHzBeTGpmdp0drz3r++ecBoFo13aNHj5Cbm2uUeG7dulWtLTMzE8DjJAwA3bp1Q2lpqXoKzFBUT02+efNmtXMymQwikUgdQ30Ysm8uLi5QKpWQyWR48cUX1e0VFRW4e/cuhg4dqnG9XC5HcXGxxmjt2ffWxlifNZkHrqmR2Xl6tHb//n2Nc23atEH79u2rrbkkJCSgqqrKKPFcunQJ586dUx8XFhZi37596Nevn3raLygoCBcvXsSRI0eq3V9WVobS0tJ6vbdEIoG3tzf+85//oKCgQN3+4MED7Nu3D35+fvWeegQM27fAwEAAQFxcnMY127ZtQ3FxMV599VWNdoVCAalUWu3Y1tYWQ4YM0RmzsT5rMg8cqZFZUo3Wfv7552p/tU+bNg1r1qxBeHg4+vXrh0uXLiElJQX29vZGicXDwwMTJ07ErFmz1GXvJSUl+PDDD9XXzJ07F8eOHcPf//53TJgwAT4+Pnj06BFu3bqFgwcPqhNFfXzyyScYM2YM/u///g9Tp06FUqnEF198gcrKSqxYsaLJ9M3T0xOhoaHqJDZs2DBcv34d27ZtQ9++fRESEqLx3hKJBP/v//0/3Lt3D+7u7jh8+DBOnTqFpUuX1piojflZU9PHpEZmydnZGVOmTMG2bduqnYuMjERhYSESExPx1VdfYciQIfj666/r/N0sffn5+WHo0KFYtWoVsrOz0aNHD+zcuVNjOs3Ozg5JSUnYsGEDEhMTceDAAbRq1Qpdu3bFnDlzqhVP6GPw4MH4+uuvsXLlSnz66acQiUTw9fXFtm3bGvzL29B9W7t2Lbp06YIdO3bg2LFjsLe3x4wZM7Bs2TKN76gBjwtzEhIS8N5772H37t1o164dPvzww1o3kDbmZ01Nn6ioqKjhq8lEJDhisRihoaH47LPPGv29R4wYgby8PJw/f77R35vMG9fUiIhIMJjUiIhIMJjUiIhIMLimRkREgsGRGhERCQaTGhERCQaTGhERCQaTGhERCQaTGhERCQaTGhERCcb/B5l2ru+C5SJsAAAAAElFTkSuQmCC",
"text/plain": [
"