{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Thinking probabilistically - Discrete variables\n", "> A Summary of lecture \"Statistical Thinking in Python (Part 1)\", via datacamp\n", "\n", "- toc: true \n", "- badges: true\n", "- comments: true\n", "- author: Chanseok Kang\n", "- categories: [Python, Datacamp, Data_Science, Statistics]\n", "- image: images/bin-dist.png" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "%matplotlib inline\n", "\n", "sns.set()\n", "\n", "df = pd.read_csv('./dataset/iris.csv')\n", "renamed_columns = ['sepal length (cm)', 'sepal width (cm)', \n", " 'petal length (cm)', 'petal width (cm)', 'species']\n", "df.columns = renamed_columns\n", "versicolor_petal_length = df[df['species'] == 'Versicolor']['petal length (cm)']\n", "setosa_petal_length = df[df['species'] == 'Setosa']['petal length (cm)']\n", "virginica_petal_length = df[df['species'] == 'Virginica']['petal length (cm)']\n", "versicolor_petal_width = df[df['species'] == 'Versicolor']['petal width (cm)']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Random number generators and hacker statistics\n", "- Hacker statistic\n", " - Uses simulated repeated measurements to compute probabilities\n", " - E.g. Coin Flips\n", " \n", "- np.random module\n", " - Suite of functions based on random number generation\n", " - ```np.random.random()```: draw a number between 0 and 1\n", " \n", "- Bernoulli trial\n", " - An experiment that has two options, \"success\" (True) and \"failure\" (False).\n", " \n", "- Random number seed\n", " - Integer fed into random number generating algorithm\n", " - Manually seed random number generator if you need reproducibility\n", " - Specified using ```np.random.seed()```\n", "\n", "- Hacker stats probabilities\n", " - Determine how to simulate data\n", " - Simulate many many times\n", " - Probability is approximately fraction of trials with the outcome of interest\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generating random numbers using the np.random module\n", "We will be hammering the np.random module for the rest of this course and its sequel. Actually, you will probably call functions from this module more than any other while wearing your hacker statistician hat. Let's start by taking its simplest function, ```np.random.random()``` for a test spin. The function returns a random number between zero and one. Call ```np.random.random()``` a few times in the IPython shell. You should see numbers jumping around between zero and one.\n", "\n", "In this exercise, we'll generate lots of random numbers between zero and one, and then plot a histogram of the results. If the numbers are truly random, all bars in the histogram should be of (close to) equal height.\n", "\n", "You may have noticed that, in the video, Justin generated 4 random numbers by passing the keyword argument ```size=4``` to ```np.random.random()```. Such an approach is more efficient than a for loop: in this exercise, however, you will write a for loop to experience hacker statistics as the practice of repeating an experiment over and over again." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD7CAYAAACFfIhNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAW/klEQVR4nO3cf0xV9/3H8deFe4v6hW8o9l4xzJhsbSTDTpbd/XBtLrNpuSBg41UzhY4tW0Olm3Vuo2NCIK4j2o6JWTZ0Js2WEfsHtRWswYvbHHaOJlUyZ7Q0bVK11h9wUVYBhXHvPd8/Gj9fka3A5cJVeT4Sg/fcc7ifN+B93sv1HptlWZYAAJAUF+sFAADuHEQBAGAQBQCAQRQAAAZRAAAYRAEAYBAFAIBhj/UCJqu3d0DhcGRvtZg7N1FXrvRHeUV3rpk2r8TMMwUzj19cnE333/8///X6uz4K4bAVcRRuHj+TzLR5JWaeKZg5Ovj1EQDAIAoAAIMoAAAMogAAMIgCAMAgCgAAY1xR6O/vV35+vj766CNJUnt7uwoKCpSdna26ujqzX2dnp3w+n7xeryoqKhQMBiVJFy9eVFFRkXJyclRaWqqBgQFJ0rVr11RSUqLc3FwVFRUpEAhEez4AwASMGYV//vOfWrdunc6ePStJGhwc1ObNm1VfX6+WlhadOnVKR44ckSSVlZWpqqpKra2tsixLjY2NkqQtW7aosLBQfr9fixcvVn19vSRpx44dcrvdOnjwoNasWaOampopGhOYOf49HJLTmRSTP0n/OzvW42OSxnzzWmNjo6qrq/X8889Lkk6ePKmFCxdqwYIFkqSCggL5/X49+OCDGhwcVGZmpiTJ5/Pp17/+tdasWaNjx47pt7/9rdn+1FNPqaysTG1tbdqzZ48kKT8/Xz//+c81PDwsh8MxJcMC0ynpf2drVkJs3h9a8OPmmNzuG796Un0xuN2bIZxug0NB9V27Me23O5XG/Im9/dF7d3e3nE6nuexyudTV1TVqu9PpVFdXl3p7e5WYmCi73T5i++2fy263KzExUVevXtW8efMmP9k4zLQfpFjNK92b/3jGMivBHpM75zd+9eS032as3eeIj9nXOhYRnEoTfhgTDodls9nMZcuyZLPZ/uv2mx9vdfvlW4+Ji5vYa99z5yZOaP/bxeIH6bVt+TG7c47lI8hZMZo5Vl/rmSiWDzxiJZbzTsVtTzgKqampI14QDgQCcrlco7b39PTI5XIpJSVFfX19CoVCio+PN/tLnzzL6OnpUWpqqoLBoAYGBpScnDyh9Vy50h/x+T9i9c2M5aOaWPn3cEj3OeKn/XaHhkNKiMHtzlQz8Wc7EIjNcwWnMymi246Ls33qg+kJR2HJkiU6c+aMzp07p8985jM6cOCAVq1apbS0NCUkJKijo0Nf+tKX1NzcLI/HI4fDIbfbrZaWFhUUFKipqUkej0eSlJWVpaamJq1fv14tLS1yu928nnCPiuWdRSyfHQF3mwlHISEhQdu2bdOGDRs0NDSkrKws5eTkSJJqa2tVWVmp/v5+ZWRkqLi4WJJUXV2t8vJy7dy5U/Pnz9f27dslSRs3blR5ebny8vKUlJSk2traKI4GAFMrlr8u+/dwaEo+77ijcPjwYfP3pUuXav/+/aP2SU9P1969e0dtT0tLU0NDw6jtycnJ2rVr13iXAAB3lFg9A5am7pko72gGABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAQRQAAAZRAAAYRAEAYBAFAIBBFAAABlEAABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAQRQAAAZRAAAYRAEAYBAFAIBBFAAABlEAABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAQRQAAMakotDc3Ky8vDzl5eXpxRdflCR1dnbK5/PJ6/WqoqJCwWBQknTx4kUVFRUpJydHpaWlGhgYkCRdu3ZNJSUlys3NVVFRkQKBwCRHAgBEKuIo3LhxQzU1NWpoaFBzc7OOHz+u9vZ2lZWVqaqqSq2trbIsS42NjZKkLVu2qLCwUH6/X4sXL1Z9fb0kaceOHXK73Tp48KDWrFmjmpqa6EwGAJiwiKMQCoUUDod148YNBYNBBYNB2e12DQ4OKjMzU5Lk8/nk9/s1PDysY8eOyev1jtguSW1tbSooKJAk5efn680339Tw8PBk5wIARMAe6YGJiYnauHGjcnNzNXv2bH35y1+Ww+GQ0+k0+zidTnV1dam3t1eJiYmy2+0jtktSd3e3OcZutysxMVFXr17VvHnzJjMXACACEUfh3Xff1Wuvvaa//vWvSkpK0k9+8hP9/e9/l81mM/tYliWbzWY+3ur2y7ceExc3/icwc+cmRjYAANzlnM6kqH/OiKNw9OhRLV26VHPnzpX0ya+EXn755REvFPf09MjlciklJUV9fX0KhUKKj49XIBCQy+WSJLlcLvX09Cg1NVXBYFADAwNKTk4e9zquXOlXOGxFNMNUfEEBYLoEAn0TPiYuzvapD6Yjfk0hPT1d7e3tun79uizL0uHDh/WVr3xFCQkJ6ujokPTJ/07yeDxyOBxyu91qaWmRJDU1Ncnj8UiSsrKy1NTUJElqaWmR2+2Ww+GIdFkAgEmI+JnCo48+qnfeeUc+n08Oh0MPP/ywSkpK9MQTT6iyslL9/f3KyMhQcXGxJKm6ulrl5eXauXOn5s+fr+3bt0uSNm7cqPLycuXl5SkpKUm1tbXRmQwAMGERR0GSSkpKVFJSMmJbenq69u7dO2rftLQ0NTQ0jNqenJysXbt2TWYZAIAo4R3NAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMCYVBQOHz4sn8+n3Nxc/eIXv5Aktbe3q6CgQNnZ2aqrqzP7dnZ2yufzyev1qqKiQsFgUJJ08eJFFRUVKScnR6WlpRoYGJjMkgAAkxBxFM6fP6/q6mrV19dr//79euedd3TkyBFt3rxZ9fX1amlp0alTp3TkyBFJUllZmaqqqtTa2irLstTY2ChJ2rJliwoLC+X3+7V48WLV19dHZzIAwIRFHIU//elPWr58uVJTU+VwOFRXV6fZs2dr4cKFWrBggex2uwoKCuT3+3XhwgUNDg4qMzNTkuTz+eT3+zU8PKxjx47J6/WO2A4AiA17pAeeO3dODodD69ev16VLl/SNb3xDDz30kJxOp9nH5XKpq6tL3d3dI7Y7nU51dXWpt7dXiYmJstvtI7ZPxNy5iZGOAAB3NaczKeqfM+IohEIhHT9+XA0NDZozZ45KS0s1a9Ys2Ww2s49lWbLZbAqHw/9x+82Pt7r98liuXOlXOGxFNMNUfEEBYLoEAn0TPiYuzvapD6YjjsIDDzygpUuXKiUlRZL0+OOPy+/3Kz4+3uwTCATkcrmUmpqqQCBgtvf09MjlciklJUV9fX0KhUKKj483+wMAYiPi1xSWLVumo0eP6tq1awqFQvrb3/6mnJwcnTlzRufOnVMoFNKBAwfk8XiUlpamhIQEdXR0SJKam5vl8XjkcDjkdrvV0tIiSWpqapLH44nOZACACYv4mcKSJUv09NNPq7CwUMPDw3rkkUe0bt06ffazn9WGDRs0NDSkrKws5eTkSJJqa2tVWVmp/v5+ZWRkqLi4WJJUXV2t8vJy7dy5U/Pnz9f27dujMxkAYMIijoIkrV69WqtXrx6xbenSpdq/f/+ofdPT07V3795R29PS0tTQ0DCZZQAAooR3NAMADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAY9JRePHFF1VeXi5J6uzslM/nk9frVUVFhYLBoCTp4sWLKioqUk5OjkpLSzUwMCBJunbtmkpKSpSbm6uioiIFAoHJLgcAMAmTisJbb72lffv2mctlZWWqqqpSa2urLMtSY2OjJGnLli0qLCyU3+/X4sWLVV9fL0nasWOH3G63Dh48qDVr1qimpmYyywEATFLEUfjXv/6luro6rV+/XpJ04cIFDQ4OKjMzU5Lk8/nk9/s1PDysY8eOyev1jtguSW1tbSooKJAk5efn680339Tw8PCkBgIARM4e6YFVVVXatGmTLl26JEnq7u6W0+k01zudTnV1dam3t1eJiYmy2+0jtt9+jN1uV2Jioq5evap58+aNex1z5yZGOgIA3NWczqSof86IovDqq69q/vz5Wrp0qV5//XVJUjgcls1mM/tYliWbzWY+3ur2y7ceExc3sScvV670Kxy2JjjBJ6biCwoA0yUQ6JvwMXFxtk99MB1RFFpaWhQIBPTkk0/q448/1vXr12Wz2Ua8UNzT0yOXy6WUlBT19fUpFAopPj5egUBALpdLkuRyudTT06PU1FQFg0ENDAwoOTk5kiUBAKIgotcUfv/73+vAgQNqbm7Wc889p8cee0xbt25VQkKCOjo6JEnNzc3yeDxyOBxyu91qaWmRJDU1Ncnj8UiSsrKy1NTUJOmT0LjdbjkcjmjMBQCIQFTfp1BbW6utW7cqJydH169fV3FxsSSpurpajY2NWr58uY4fP64f/vCHkqSNGzfqxIkTysvL0yuvvKKqqqpoLgcAMEERv9B8k8/nk8/nkySlp6dr7969o/ZJS0tTQ0PDqO3JycnatWvXZJcAAIgS3tEMADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAAyiAAAwiAIAwCAKAACDKAAADKIAADCIAgDAIAoAAIMoAAAMogAAMIgCAMAgCgAAgygAAIxJReE3v/mN8vLylJeXp5deekmS1N7eroKCAmVnZ6uurs7s29nZKZ/PJ6/Xq4qKCgWDQUnSxYsXVVRUpJycHJWWlmpgYGAySwIATELEUWhvb9fRo0e1b98+NTU16fTp0zpw4IA2b96s+vp6tbS06NSpUzpy5IgkqaysTFVVVWptbZVlWWpsbJQkbdmyRYWFhfL7/Vq8eLHq6+ujMxkAYMIijoLT6VR5ebnuu+8+ORwOfe5zn9PZs2e1cOFCLViwQHa7XQUFBfL7/bpw4YIGBweVmZkpSfL5fPL7/RoeHtaxY8fk9XpHbAcAxIY90gMfeugh8/ezZ8/q4MGDeuqpp+R0Os12l8ulrq4udXd3j9judDrV1dWl3t5eJSYmym63j9g+EXPnJkY6AgDc1ZzOpKh/zoijcNP777+vZ555Rs8//7zi4+N19uxZc51lWbLZbAqHw7LZbKO23/x4q9svj+XKlX6Fw1ZEa5+KLygATJdAoG/Cx8TF2T71wfSkXmju6OjQd77zHf34xz/WypUrlZqaqkAgYK4PBAJyuVyjtvf09MjlciklJUV9fX0KhUIj9gcAxEbEUbh06ZK+//3vq7a2Vnl5eZKkJUuW6MyZMzp37pxCoZAOHDggj8ejtLQ0JSQkqKOjQ5LU3Nwsj8cjh8Mht9utlpYWSVJTU5M8Hk8UxgIARCLiXx+9/PLLGhoa0rZt28y2tWvXatu2bdqwYYOGhoaUlZWlnJwcSVJtba0qKyvV39+vjIwMFRcXS5Kqq6tVXl6unTt3av78+dq+ffskRwIARCriKFRWVqqysvI/Xrd///5R29LT07V3795R29PS0tTQ0BDpMgAAUcQ7mgEABlEAABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAQRQAAAZRAAAYRAEAYBAFAIBBFAAABlEAABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAQRQAAAZRAAAYRAEAYBAFAIBBFAAABlEAABhEAQBgEAUAgEEUAAAGUQAAGEQBAGAQBQCAcUdE4Y033tDy5cuVnZ2tPXv2xHo5ADBj2WO9gK6uLtXV1en111/Xfffdp7Vr1+qrX/2qHnzwwVgvDQBmnJhHob29XV/72teUnJwsSfJ6vfL7/frBD34wruPj4myTun3X/bMndTy3e+ffNjPPjNueabcrRXb/N9YxNsuyrEgXFA2/+93vdP36dW3atEmS9Oqrr+rkyZN64YUXYrksAJiRYv6aQjgcls32/+WyLGvEZQDA9Il5FFJTUxUIBMzlQCAgl8sVwxUBwMwV8yh8/etf11tvvaWrV6/qxo0bOnTokDweT6yXBQAzUsxfaJ43b542bdqk4uJiDQ8Pa/Xq1frCF74Q62UBwIwU8xeaAQB3jpj/+ggAcOcgCgAAgygAAAyiAAAw7vkojHWyvc7OTvl8Pnm9XlVUVCgYDMZgldE11sx//vOf9eSTT2rFihV69tln9fHHH8dgldE13pMqtrW16bHHHpvGlU2dsWb+4IMP9K1vfUsrVqzQ9773vRnxfT59+rRWrVqlFStW6JlnntG1a9disMro6u/vV35+vj766KNR103J/Zd1D7t8+bK1bNkyq7e31xoYGLAKCgqs999/f8Q+eXl51j/+8Q/LsizrZz/7mbVnz55YLDVqxpq5r6/PeuSRR6zLly9blmVZO3bssF544YVYLTcqxvN9tizLCgQCVk5OjrVs2bIYrDK6xpo5HA5b2dnZ1pEjRyzLsqxf/vKX1ksvvRSr5UbFeL7P69ats9ra2izLsqytW7da27dvj8VSo+bEiRNWfn6+lZGRYZ0/f37U9VNx/3VPP1O49WR7c+bMMSfbu+nChQsaHBxUZmamJMnn8424/m401szDw8Oqrq7WvHnzJEmLFi3SpUuXYrXcqBhr5psqKyvHfaLFO91YM58+fVpz5swxbwRdv369ioqKYrXcqBjP9zkcDmtgYECSdOPGDc2aNSsWS42axsZGVVdX/8ezPEzV/dc9HYXu7m45nU5z2eVyqaur679e73Q6R1x/Nxpr5vvvv19PPPGEJGlwcFC7d+/W448/Pu3rjKaxZpakP/7xj/r85z+vJUuWTPfypsRYM3/44Yd64IEHtHnzZq1cuVLV1dWaM2dOLJYaNeP5PpeXl6uyslKPPvqo2tvbtXbt2uleZlTV1NTI7Xb/x+um6v7rno7CWCfbuxdPxjfemfr6+lRSUqL09HStXLlyOpcYdWPN/N577+nQoUN69tlnY7G8KTHWzMFgUG+//bbWrVunffv2acGCBdq2bVsslho1Y808ODioiooK/eEPf9DRo0dVWFion/70p7FY6rSYqvuvezoKY51s7/bre3p67vqT8Y3nBIPd3d0qLCzUokWLVFNTM91LjLqxZvb7/QoEAlq1apVKSkrM/HezsWZ2Op1auHChHn74YUlSfn6+Tp48Oe3rjKaxZn7vvfeUkJBgTpPzzW9+U2+//fa0r3O6TNX91z0dhbFOtpeWlqaEhAR1dHRIkpqbm+/6k/GNNXMoFNL69euVm5urioqKu/6ZkTT2zM8995xaW1vV3Nys3bt3y+Vy6ZVXXonhiidvrJm/+MUv6urVq3r33XclSYcPH1ZGRkaslhsVY828cOFCXb58WR988IEk6S9/+YuJ4r1oyu6/Jv1S9R1u//79Vl5enpWdnW3t3r3bsizLevrpp62TJ09almVZnZ2d1qpVqyyv12v96Ec/soaGhmK53Kj4tJkPHTpkLVq0yFqxYoX5s3nz5hivePLG+j7fdP78+Xvifx9Z1tgznzhxwlq1apW1fPly67vf/a7V09MTy+VGxVgzt7W1WQUFBVZ+fr717W9/2/rwww9judyoWbZsmfnfR1N9/8UJ8QAAxj396yMAwMQQBQCAQRQAAAZRAAAYRAEAYBAFAIBBFAAABlEAABj/BxcIDm96R/+zAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Seed the random number generator\n", "np.random.seed(42)\n", "\n", "# Initialize random numbers: random_numbers\n", "random_numbers = np.empty(100000)\n", "\n", "# Generate random numbers by looping over range(100000)\n", "for i in range(100000):\n", " random_numbers[i] = np.random.random()\n", "\n", "# Plot a histogram\n", "_ = plt.hist(random_numbers)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The np.random module and Bernoulli trials\n", "You can think of a Bernoulli trial as a flip of a possibly biased coin. Specifically, each coin flip has a probability p of landing heads (success) and probability 1−p of landing tails (failure). In this exercise, you will write a function to perform n Bernoulli trials, ```perform_bernoulli_trials(n, p)```, which returns the number of successes out of n Bernoulli trials, each of which has probability p of success. To perform each Bernoulli trial, use the ```np.random.random()``` function, which returns a random number between zero and one." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def perform_bernoulli_trials(n, p):\n", " \"\"\"Perform n Bernoulli trials with success probability p and return number of successes.\"\"\"\n", " # Initialize number of successes: n_successes\n", " n_success = 0\n", " \n", " # Perform trials\n", " for i in range(n):\n", " # Choose random number between zero and one: random_number\n", " random_number = np.random.random()\n", " \n", " # If less than p, it`s a success so add one to n_success\n", " if random_number < p:\n", " n_success += 1\n", " \n", " return n_success" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How many defaults might we expect?\n", "Let's say a bank made 100 mortgage loans. It is possible that anywhere between 0 and 100 of the loans will be defaulted upon. You would like to know the probability of getting a given number of defaults, given that the probability of a default is p = 0.05. To investigate this, you will do a simulation. You will perform 100 Bernoulli trials using the ```perform_bernoulli_trials()``` function you wrote in the previous exercise and record how many defaults we get. Here, a success is a default. (Remember that the word \"success\" just means that the Bernoulli trial evaluates to True, i.e., did the loan recipient default?) You will do this for another 100 Bernoulli trials. And again and again until we have tried it 1000 times. Then, you will plot a histogram describing the probability of the number of defaults.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEJCAYAAACKWmBmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxU9eI+8GfYJMVCbECvmqWWlkvmkghK1w0EGVFEI1Q01FxS+5obKoWSK6WY+5LX+3oJJu6K1xC9phmgueQuLt3rcl0AQQiQZZbP7w9+npzQD4M6zJjP+y/OnHPmPGdm4OGcM/MZlRBCgIiI6DFsLB2AiIisG4uCiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSdpYOYA737hXAYKj4x0Nq1nRCVla+GRI9O8z49Kw9H2D9Ga09H8CMFWFjo0KNGtUeO/8vWRQGg3iioniwrrVjxqdn7fkA689o7fkAZnxWeOqJiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSLAoiIpL6S36OgkxX/eWX4Fjl2b8M1Orq5S5TVKxD3u+Fz3zbRPRssShecI5V7KAZv8Mi206YH4A8i2yZiCqCp56IiEiKRUFERFIsCiIikmJREBGRFIuCiIikWBRERCTFoiAiIikWBRERSbEoiIhIyqxFkZCQAD8/P3h7eyMuLq7M/H379iEgIAA9e/bEqFGjkJubCwC4desW+vfvj+7du2PkyJEoKCgwZ0wiIpIwW1Gkp6cjJiYG69evx/bt2xEfH48rV64o8/Pz8zF9+nSsWrUKO3fuROPGjbF48WIAwIwZMxASEoLExEQ0a9YMy5YtM1dMIiIqh9mKIiUlBe7u7nB2dkbVqlXh4+ODxMREZb5Wq0VkZCTc3NwAAI0bN8bt27eh1Wpx9OhR+Pj4AAACAwON1iMiospltqLIyMiAWq1Wpl1dXZGenq5M16hRA926dQMAFBUVYdWqVejatSvu3bsHJycn2NmVjleoVquN1iMiosplttFjDQYDVCqVMi2EMJp+IC8vD59++imaNGmC3r17Iz09vcxyj1pPpmZNpycLDdOGx7a05yGjqSy1L8/DY2jtGa09H8CMz4rZiqJWrVo4duyYMp2ZmQlXV1ejZTIyMjBkyBC4u7tj6tSpAAAXFxfk5eVBr9fD1tb2keuVJysrHwaDqHBmtbo6MjOte+DrZ53R0i9SSzzeL+Lz/KxZez6AGSvCxkYl/QfbbKeePDw8kJqaiuzsbBQWFiIpKQleXl7KfL1ejxEjRsDX1xfTpk1Tjhrs7e3Rpk0b7N69GwCwfft2o/WIiKhyme2Iws3NDePGjUNoaCi0Wi2CgoLQokULDBs2DGPHjsWdO3dw/vx56PV67NmzBwDQrFkzzJo1C5GRkQgPD8fy5ctRu3ZtLFiwwFwxiYioHGb9hjuNRgONRmN02+rVqwEAzZs3R1pa2iPXq1OnDtatW2fOaEREZCJ+MpuIiKT4ndn0winR6i12Eb+oWIe83wstsm2iJ8WioBeOg70tNON3WGTbCfMDYPn3uBBVDE89ERGRFIuCiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSLAoiIpJiURARkRSLgoiIpFgUREQkxaIgIiIpFgUREUmxKIiISIpFQUREUiwKIiKSYlEQEZEUi4KIiKRYFEREJMWiICIiKRYFERFJsSiIiEiKRUFERFIsCiIikmJREBGRFIuCiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSLAoiIpJiURARkRSLgoiIpFgUREQkxaIgIiIpsxZFQkIC/Pz84O3tjbi4uMcuN2nSJGzdulWZ3rZtGzp06ICAgAAEBAQgJibGnDGJiEjCzlx3nJ6ejpiYGGzduhUODg4IDg5Gu3bt0KhRI6NlIiMjkZqaCnd3d+X2s2fPIjw8HP7+/uaKR0REJjLbEUVKSgrc3d3h7OyMqlWrwsfHB4mJiUbLJCQkoEuXLvD19TW6/cyZM9i2bRs0Gg0mTJiA3Nxcc8UkIqJymK0oMjIyoFarlWlXV1ekp6cbLTN06FD07du3zLpqtRqjRo3Czp07Ubt2bURFRZkrJhERlcNsp54MBgNUKpUyLYQwmpZZunSp8vPQoUPRrVu3Cm27Zk2nCi3/MLW6+hOvW1meh4ym+ivti6lM3Wdrf2ysPR/AjM+K2YqiVq1aOHbsmDKdmZkJV1fXctfLy8vDli1bMHjwYAClBWNra1uhbWdl5cNgEBVaByh9wjIz8yq8XmV61hkt/SK1xOP9POyztb8WrT0fwIwVYWOjkv6DbbZTTx4eHkhNTUV2djYKCwuRlJQELy+vcterWrUqvvvuO5w6dQoAEBsbW+EjCiIienbMdkTh5uaGcePGITQ0FFqtFkFBQWjRogWGDRuGsWPHonnz5o9cz9bWFgsXLsT06dNRVFSE119/HdHR0eaKSURE5TBbUQCARqOBRqMxum316tVllps7d67RdJs2bbBt2zZzRiMiIhPxk9lERCTFoiAiIikWBRERSbEoiIhIikVBRERSJhXFvXv3zJ2DiIislElF0aNHD4wfP97ok9ZERPRiMKko9u/fDw8PD0RHR0Oj0SAuLg75+fnmzkZERFbApKJwdHREnz59sHHjRkREROAf//gHOnbsiBkzZvC0FBHRX5zJF7N/+uknjBkzBuPGjUPXrl2xYcMG1K5dG6NGjTJnPiIisjCThvDo1KkTnJ2dERISgq+//hqOjo4AgMaNGyM+Pt6sAYmIyLJMKoro6Gi0bdvW6LYrV66gUaNG+Pe//22WYEREZB2kp55ycnKQk5ODr776Crm5ucr03bt3MXr06MrKSEREFiQ9ohg/fjySk5MBAO3atftjJTs7+Pj4mDcZERFZBWlRrFmzBgAwZcoUzJkzp1ICERGRdZEWxW+//YaGDRtiwIABOHfuXJn5TZs2NVswIiKyDtKimDdvHlatWoUxY8aUmadSqXghm4joBSAtilWrVgEo/WQ2ERG9mKRFMXPmTOnKERERzzQMERFZH2lRODs7V1YOIiKyUtKi4GcliIhIWhQfffQRvv/+e7z33ntQqVRl5p84ccJswYiIyDpIi+Lbb78FAOzatatSwhARkfWRFoWrqysAoE6dOjh06BBSUlJgZ2cHLy+vMmM/ERHRX5NJw4yvWLECc+bMgaOjI2xsbBAREYG4uDhzZyMiIitg0uixu3btwsaNG+Hk5AQACAsLQ0hICPr372/WcEREZHkmHVFUqVIF1apVU6ZfeeUVVKlSxWyhiIjIekiPKJKSkgAAb7zxBkaNGoW+ffvC1tYW27dvR7NmzSolIBERWZa0KNatW2c0vXbtWuXnrKws8yQiIiKrUqGiICKiF49JF7OvXr2K2NhY3L9/H0IIGAwGXLt2DRs2bDB3PiIisjCTLmaPHz8eWq0Wv/76K+rUqYMrV67grbfeMnc2IiKyAiYVRUFBAWbMmIEOHTrAy8sLa9euxcmTJ82djYiIrIBJRfFgFNn69evj8uXLePnllx859hMREf31mHSNon79+pg1axZ69+6NadOm4f79+9DpdObORkREVsCkI4rp06ejTZs2eOedd9CvXz8cPnwYUVFR5s5GRERWwKQjipdeegnvv/8+fvzxR9StWxdz5szByy+/bO5sRERkBUwqigMHDmDy5Mlo1KgRDAYDbty4gZiYGI4gS0T0AjCpKL799lvExsbizTffBACcO3cOX3zxBbZu3WrWcEREZHkmXaNQqVRKSQBA06ZNIYQod72EhAT4+fnB29tbOiz5pEmTjErn1q1b6N+/P7p3746RI0eioKDAlJhERGQG0qLIyclBTk4OmjVrhjVr1qCgoACFhYWIi4uDu7u79I7T09MRExOD9evXY/v27YiPj8eVK1fKLDNixAjs2bPH6PYZM2YgJCQEiYmJaNasGZYtW/aEu0dERE9LeurJ3d0dKpVKOXr4+uuvlXkqlQqTJ09+7LopKSlwd3dXPoPh4+ODxMREjB49WlkmISEBXbp0UZYBAK1Wi6NHj2Lp0qUAgMDAQAwYMAATJ058gt0jIqKnJS2KtLS0J77jjIwMqNVqZdrV1RWnT582Wmbo0KEAgOPHjyu33bt3D05OTrCzK42mVquRnp7+xDmIiOjpmHQx22AwYM2aNfjpp5+g0+ng6emJESNGKH/MH7fOw5/eFkKY9GnuRy1X0U+B16zpVKHlH6ZWV3/idSvL85DRVH+lfTGVqfts7Y+NtecDmPFZMako5s+fj7S0NAwaNAgGgwHx8fGYN28epk2b9th1atWqhWPHjinTmZmZcHV1LXdbLi4uyMvLg16vh62trcnrPSwrKx8GQ/kX2/9Mra6OzMy8Cq9XmZ51Rku/SC3xeD8P+2ztr0VrzwcwY0XY2Kik/2Cb9K6nQ4cOYcWKFejatSu8vb2xfPlyHDp0SLqOh4cHUlNTkZ2djcLCQiQlJcHLy6vcbdnb26NNmzbYvXs3AGD79u0mrUdEROZh0hGFEAL29vbKtIODg9H0o7i5uWHcuHEIDQ2FVqtFUFAQWrRogWHDhmHs2LFo3rz5Y9eNjIxEeHg4li9fjtq1a2PBggUm7g49T0q0eov/d09E5TOpKJo0aYLZs2djwIABUKlUiI2NNen7KDQaDTQajdFtq1evLrPc3Llzjabr1KnDb9d7ATjY20IzfkelbzdhfkClb5PoeWbSqafIyEj8/vvvCA4ORr9+/ZCdnY0vvvjC3NmIiMgKmHREsXLlyjL/9RMR0YvBpCOKAwcOmDkGERFZK5OOKOrWrYuwsDC0atUK1apVU27/+OOPzRaMiIisg0lF8WCIjXPnzsHW1hbVq/OdKkRELwqTTj0NHToUly5dwsGDB7F//37cuHEDY8aMMXc2IiKyAiYVxdSpU9GvXz+cOnUKJ0+ehI+Pj/RT2URE9NdhUlEUFhbiww8/hL29PRwcHDBw4EDcvXvX3NmIiMgKmFQUDRo0wIkTJ5TpS5cuoW7dumYLRURE1sOki9m3bt3CwIED0bhxY9jZ2eH8+fNQq9XKp64TEhLMGpKIiCzHpKKYMGGCuXMQEZGVMqko3n//fXPnICIiK2XSNQoiInpxsSiIiEiKRUFERFIsCiIikmJREBGRFIuCiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSLAoiIpJiURARkRSLgoiIpFgUREQkxaIgIiIpFgUREUmxKIiISIpFQUREUiwKIiKSYlEQEZEUi4KIiKRYFEREJMWiICIiKTtLB6BS1V9+CY5VTHs61OrqZk5DRPQHFoWVcKxiB834HZW+3YT5AZW+TSJ6vvDUExERSbEoiIhIyqxFkZCQAD8/P3h7eyMuLq7M/AsXLiAwMBA+Pj6YNm0adDodAGDbtm3o0KEDAgICEBAQgJiYGHPGJCIiCbNdo0hPT0dMTAy2bt0KBwcHBAcHo127dmjUqJGyzMSJEzFz5ky0bNkSU6dOxcaNGxESEoKzZ88iPDwc/v7+5opHREQmMtsRRUpKCtzd3eHs7IyqVavCx8cHiYmJyvybN2+iqKgILVu2BAAEBgYq88+cOYNt27ZBo9FgwoQJyM3NNVdMIiIqh9mOKDIyMqBWq5VpV1dXnD59+rHz1Wo10tPTlZ/DwsLQqlUrLFiwAFFRUZg/f77J265Z0+mJc/Otp2Rupr7GrP21aO35AGZ8VsxWFAaDASqVSpkWQhhNy+YvXbpUuX3o0KHo1q1bhbadlZUPg0FUOLNaXR2ZmXkVXu9ZeB5eLPRsmPIas+Rr0RTWng9gxoqwsVFJ/8E226mnWrVqITMzU5nOzMyEq6vrY+ffvXsXrq6uyMvLwz//+U/ldiEEbG1tzRWTiIjKYbai8PDwQGpqKrKzs1FYWIikpCR4eXkp8+vUqYMqVarg+PHjAIAdO3bAy8sLVatWxXfffYdTp04BAGJjYyt8REFERM+O2U49ubm5Ydy4cQgNDYVWq0VQUBBatGiBYcOGYezYsWjevDm++eYbREREID8/H02bNkVoaChsbW2xcOFCTJ8+HUVFRXj99dcRHR1trphERFQOsw7hodFooNFojG5bvXq18nOTJk2wefPmMuu1adMG27ZtM2c0IiIyET+ZTUREUhwUkKgSlWj1Fnl7bFGxDnm/Fz6z+6MXC4uCqBI52NtabJRgy78Jk55XPPVERERSLAoiIpJiURARkRSLgoiIpFgUREQkxaIgIiIpFgUREUmxKIiISIpFQUREUiwKIiKSYlEQEZEUi4KIiKRYFEREJMWiICIiKQ4z/pCKfFcA0fPEHK9tU++P34Xx/GNRPMRS3xUAlH5fAJG5WPq1ze/CeL7x1BMREUmxKIiISIpFQUREUiwKIiKSYlEQEZEUi4KIiKRYFEREJMWiICIiKRYFERFJsSiIiEiKRUFERFIsCiIikmJREBGRFIuCiIikWBRERCTFoiAiIikWBRERSbEoiIhIikVBRERSLAoiIpIya1EkJCTAz88P3t7eiIuLKzP/woULCAwMhI+PD6ZNmwadTgcAuHXrFvr374/u3btj5MiRKCgoMGdMIiKSsDPXHaenpyMmJgZbt26Fg4MDgoOD0a5dOzRq1EhZZuLEiZg5cyZatmyJqVOnYuPGjQgJCcGMGTMQEhKCHj16YOnSpVi2bBkmTpxorqhEZEYlWj3U6uqVvt1iS223RI8qDrYmL/8sMxYV65D3e+Ezu78HzFYUKSkpcHd3h7OzMwDAx8cHiYmJGD16NADg5s2bKCoqQsuWLQEAgYGBWLRoEfr27YujR49i6dKlyu0DBgyoUFHY2KieOLdrjZeeeN2nZaltc59fjG1barsO9rYYMjOp0re7JsL7hdrug20XPMHfv/L+ZqqEEOJJQ8msXLkS9+/fx7hx4wAAmzZtwunTp/HVV18BAH799VdER0fj+++/BwBcu3YNn3zyCdatW4egoCD89NNPAACdToeWLVvi7Nmz5ohJRETlMNs1CoPBAJXqj5YSQhhNP27+n5cDUGaaiIgqj9mKolatWsjMzFSmMzMz4erq+tj5d+/ehaurK1xcXJCXlwe9Xv/I9YiIqHKZrSg8PDyQmpqK7OxsFBYWIikpCV5eXsr8OnXqoEqVKjh+/DgAYMeOHfDy8oK9vT3atGmD3bt3AwC2b99utB4REVUus12jAErfHrty5UpotVoEBQVh2LBhGDZsGMaOHYvmzZsjLS0NERERyM/PR9OmTTFnzhw4ODjg5s2bCA8PR1ZWFmrXro0FCxbglVdeMVdMIiKSMGtREBHR84+fzCYiIikWBRERSbEoiIhIikVBRERSLIr/r7wBDC1tyZIl6NGjB3r06IHo6GhLx5GaN28ewsPDLR3jkfbv34/AwED4+vpi5syZlo5Txo4dO5Tned68eZaOo8jPz4e/vz/+97//ASgdokej0cDb2xsxMTEWTlfqzxnj4+Ph7+8PjUaDKVOmoKSkxMIJy2Z8IDY2FgMHDrRQKhMIEnfu3BGdOnUS9+7dEwUFBUKj0YjLly9bOpYiOTlZfPjhh6K4uFiUlJSI0NBQkZSUZOlYj5SSkiLatWsnJk+ebOkoZVy/fl106NBB3L59W5SUlIiPPvpIHDhwwNKxFPfv3xdt27YVWVlZQqvViqCgIJGcnGzpWOLkyZPC399fNG3aVNy4cUMUFhaKDz74QFy/fl1otVoRFhZm8cfxzxn/85//iG7duom8vDxhMBjEpEmTxNq1a60q4wOXL18WHTt2FAMGDLBgOjkeUcB4AMOqVasqAxhaC7VajfDwcDg4OMDe3h4NGzbErVu3LB2rjJycHMTExGDEiBGWjvJIe/fuhZ+fH2rVqgV7e3vExMTg3XfftXQshV6vh8FgQGFhIXQ6HXQ6HapUqWLpWNi4cSMiIyOVERJOnz6N+vXro169erCzs4NGo7H478ufMzo4OCAyMhJOTk5QqVR46623LP478+eMAFBSUoIvv/wSY8eOtWCy8plt9NjnSUZGBtRqtTLt6uqK06dPWzCRsTfffFP5+erVq/jhhx+UwRStyZdffolx48bh9u3blo7ySNeuXYO9vT1GjBiB27dv4+9//zv+7//+z9KxFE5OTvjss8/g6+uLl156CW3btkWrVq0sHQuzZs0ymn7U70t6enplxzLy54x16tRBnTp1AADZ2dmIi4vDnDlzLBFN8eeMADB//nz06dMHdevWtUAi0/GIAuUPYGgtLl++jLCwMEyaNAmvv/66peMY2bRpE2rXro327dtbOspj6fV6pKamYvbs2YiPj8fp06exbds2S8dSpKWlYcuWLfjxxx9x6NAh2NjYYM2aNZaOVcbz8vsClH4vzqBBg9CnTx+0a9fO0nGMJCcn4/bt2+jTp4+lo5SLRYHyBzC0BsePH8fgwYMxfvx49O7d29Jxyti9ezeSk5MREBCARYsWYf/+/Zg9e7alYxl59dVX0b59e7i4uMDR0RFdu3a1qiPHn3/+Ge3bt0fNmjXh4OCAwMBA/PLLL5aOVcbz8PsCAL/99huCg4PRu3dvfPrpp5aOU8auXbtw+fJlBAQEICIiAmfPnrWqI1wjlr5IYg0eXMzOysoS9+/fFz179hSnTp2ydCzFrVu3RLt27URKSoqlo5hky5YtVnkx++TJk8LHx0fk5uYKnU4nhg8fLjZu3GjpWIpDhw6Jnj17ioKCAmEwGMQXX3whFi1aZOlYik6dOokbN26IoqIi4eXlJa5evSp0Op0YMmSI2L17t6XjCSH+yJiXlyc++OADsW3bNktHKuNBxocdPnzYqi9m8xoFADc3N4wbNw6hoaHKAIYtWrSwdCzFmjVrUFxcjLlz5yq3BQcH46OPPrJgqufPu+++i6FDhyIkJARarRaenp5WddjfoUMHnD9/HoGBgbC3t0fz5s3xySefWDpWGVWqVMHcuXMxZswYFBcX44MPPkD37t0tHcvI5s2bcffuXaxduxZr164FAHTu3BmfffaZhZM9nzgoIBERSfEaBRERSbEoiIhIikVBRERSLAoiIpJiURARkRSLgsyic+fOOHPmTKVsKz8/H8HBwejRoweSkpKky/r7++PIkSPSZfR6PUaOHAkfHx/ExsY+Uabw8HDlU9VLlizBvn37nuh+KmLTpk0VHvnYlMcuLS0NHTp0MLrt1KlT6NOnD3x9fTFo0CBkZGQo81auXInu3bujW7duWLx4MR71xsrFixcjKiqqQlnJcvg5CnruXbhwAVlZWdi7d+8zub/09HT8/PPPOHnyJGxtbZ/6/o4cOYJGjRo9g2Ryx48fNxoXzBSyx06n0yE2NharV6/G/fv3ldtLSkowduxYLFiwAK1bt8b69esxbdo0rF69GgcPHsQPP/yArVu3wtbWFkOGDEHDhg3h5+f31PtHlsOieEEdOXIEMTExqFevHi5fvgydTocZM2agdevWCA8Px5tvvokhQ4YAgNF0586d4e/vj8OHDyM3NxdDhw7FiRMncO7cOdjZ2WH58uVwc3MDAKxfvx5paWkoKSnBxx9/jKCgIACl3wmxfPlyaLVaODo6YvLkyXjvvfewePFinDx5EhkZGWjcuDG++eYbo8z79u3DkiVLYDAYUK1aNUyZMgVOTk6YOnUq0tPTERAQgPj4eDg6OirrXLlyBVOnTkVhYSEaNGhg9AfvxIkT+Oabb1BYWAgbGxuMHj0abdu2xdChQ6HT6RAYGIjFixfjl19+QXx8PLRaLXJzczFs2DCEhIRg69at2LNnD1auXAkAZaYBIC4uDmfPnkV0dDRsbW1Ro0YNzJ07FwaDAQAwfPhw+Pj4lHl+4uPjsW7dOtjY2ODVV1/FF198gTfeeOOxz81rr72G/fv3Izk5GY6Ojujfv/9TP3bnz5/HxYsXsWTJEoSFhSm3nzlzBk5OTmjdujUAICgoCLNnz8a9e/ewd+9e+Pv7o2rVqgCAwMBA7Ny5U1oUly9fRlRUFHJycqBSqRAWFoZevXrBYDBg9uzZOHXqFAoKCiCEwMyZM5XXqJOTEy5evIg7d+6gcePGmDdvHqpVq4ZFixZh7969sLe3R40aNTBnzhyrHGLkuWLZD4aTpRw+fFi8/fbb4vz580IIIdasWSP69+8vhBBi8uTJ4rvvvlOWfXi6U6dOYvbs2UIIIf71r3+JJk2aiAsXLgghhBg1apRYvny5slxkZKQQonSIlPbt24tLly6J//73v8Lf319kZ2cLIYS4dOmS8PT0FAUFBWLRokXCx8dHaLXaMnmvXLkiPDw8xPXr14UQpd974enpKfLy8sThw4dFjx49HrmfAQEByjAdx44dE40bNxaHDx8WOTk5wtvbWxlK4c6dO8LLy0vcvHlT3LhxQ7Rs2VIIIUR+fr7o16+fkvfXX39V5m3ZskV88sknyrYenn74MRswYID44YcfhBBChIaGil27dgkhhLhw4YKYPn16mcwpKSmia9euIisrS7lfX19fYTAYpM/Nn+c97WP3wMOPhxBC7Nq1S4SFhRkt07FjR3HhwgURFham7J8Qpd+l0qtXrzL3uWjRIjFjxgyh1WpFly5dxJ49e4QQpc9Dx44dxYkTJ8SJEyfEmDFjhF6vF0IIsXLlSjF8+HBlXx/+jpZevXqJzZs3i1u3bolWrVqJ4uJiIUTp63rv3r3S/aPy8YjiBfa3v/0Nb7/9NgDgnXfeMXkkVW9vbwBAvXr18Oqrr6JJkyYAgNdeew25ubnKcsHBwQBKh0jx9PREamoqbG1tkZGRgcGDByvLqVQqXL9+HQDQsmVL2NmVfVkePnwY7u7uqFevHgAog/udPXv2sSOX3rt3DxcvXkSvXr0AAK1bt1ZOzZw8eRKZmZlGg8WpVCpcvHjR6PRNtWrVsGLFChw8eBBXr15FWlqa0VFJRfn6+iIqKgr79++Hh4cHPv/88zLLHDp0CH5+fnBxcQFQ+l/5rFmzynwrmqme5LGT+fPosUDpCLK2trZlRpIVQsDG5vGXQq9evYri4mLlNeXm5gZvb28cOnQIY8eOxSuvvIINGzbgxo0bOHLkCKpVq6as27FjRzg4OAAA3nrrLeTm5sLNzQ1NmjRB79694eXlBS8vL6se0fh5wWnce+IAAAP1SURBVKJ4gT18mkGlUikXHR/+GQC0Wq3Reg9+OQHA3t7+sff/8B8Ig8EAOzs76PV6tG/fHgsXLlTm3b59G66urti7d69yyuLPHvfHSafTSTM8WO6BByWk1+vRsGFDbNq0SZmXnp4OFxcXo+9WuHPnDj788EP069cPrVu3Rvfu3fHjjz8CKP9xepTg4GB06tQJycnJOHToEJYsWYLExESjLyh6cFrqUfv6JNt8msfuUWrXrm108Vqr1SInJwdubm5l5mVkZKBWrVqPvS+9Xv/YbAcOHMCsWbPw8ccfo0uXLmjQoAF27typLPeo16+NjQ1iY2Nx5swZZUj5jh07YtKkSRXeT/oD3/VEZdSoUQNnz54FUPrH80mHun5whHLr1i2kpqaiffv2aN++PZKTk/Hbb78BAA4ePIiePXuiqKhIel/t27fHzz//jBs3bgAAUlNTcfv2bek31NWoUQNNmzZVyuDcuXO4dOkSgNIjl2vXruHo0aMASi/q+vj4lPkCnrNnz8LFxQWjRo1Chw4dlJLQ6/VwcXHB5cuXUVxcDK1Wiz179jwyh62tLXQ6HYDSorhw4QICAwPx1Vdf4ffffzcashso/U959+7dyM7OBgBs2bIFzs7OqF+/vvS5eXg7T/vYybz77rvIycnBiRMnlHwtW7bEyy+/jC5dumDnzp24f/8+SkpKsHXrVnTt2vWx99WgQQPY2dkp77hKT0/Hnj174OHhgeTkZHTq1AkhISFo1qwZ9u3bB71eL82WlpYGf39/NGzYEMOHD8fgwYMr7d13f2U8oqAyBg4ciAkTJsDHxwd169aFu7v7E91PcXExevfuDa1Wi4iICLzxxhsAgKioKHz++ecQQigXwB8+pfAojRo1QmRkJEaPHg29Xg9HR0esWLEC1atXl663YMECTJkyBRs2bMBrr72GBg0aAABcXFywaNEiREdHo7i4GEIIREdHo27dukaneDw9PbF582Z0794dKpUK77//PlxcXHDt2jV4enqibdu28PX1hVqtRrt27XDx4sUyGTp37owFCxZAq9ViwoQJmD17NhYuXAiVSoXRo0eX+XYzT09PDB48GIMGDYLBYICLiwtWrlwJGxsb6XPj5eWljDA8fPjwp37sHsfe3h5LlixBVFQUCgsL4ezsjHnz5in7eunSJfTt2xdarRZdunRRTv097r6WLVuGmTNnYvHixdDr9fj000/h7u4OtVqN8ePHQ6PRQKfTwdPTE0lJSY884nqgSZMm8PX1RZ8+fVC1alU4OjoiIiLiifaT/sDRY4mISIqnnoiISIpFQUREUiwKIiKSYlEQEZEUi4KIiKRYFEREJMWiICIiKRYFERFJ/T+jp1MgSAx6cAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Seed random number generator\n", "np.random.seed(42)\n", "\n", "# Initialize the number of defaults: n_defaults\n", "n_defaults = np.empty(1000)\n", "\n", "# Compute the number of defaults\n", "for i in range(1000):\n", " n_defaults[i] = perform_bernoulli_trials(100, 0.05)\n", " \n", "# Plot the histogram with default number of bins; label your axes\n", "_ = plt.hist(n_defaults, density=True)\n", "_ = plt.xlabel('number of defaults out of 100 loans')\n", "_ = plt.ylabel('probability')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Will the bank fail?\n", "Plot the number of defaults you got from the previous exercise, in your namespace as n_defaults, as a CDF. \n", "\n", "If interest rates are such that the bank will lose money if 10 or more of its loans are defaulted upon, what is the probability that the bank will lose money?" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def ecdf(data):\n", " \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n", " # Number of data points: n\n", " n = len(data)\n", "\n", " # x-data for the ECDF: x\n", " x = np.sort(data)\n", "\n", " # y-data for the ECDF: y\n", " y = np.arange(1, n + 1) / n\n", "\n", " return x, y" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Probability of losing money = 0.022\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAa7UlEQVR4nO3df1BVBf7/8ddF8QerI0b3QoOV21K6y4Ka21d0W8xGIQmwkD6hJbtTYe62UYxf+6GmzbT+qE8uu7XbzAc/js2OWGo/VPZTyG7m99uKrqOr4WKmrF8ySeEKRiKgwD3fP7rc9Yrx8x7OQZ6Pf/TNuffc11wGXpxz7jnHYRiGIQBAvxdkdQAAgD1QCAAASRQCAMCLQgAASKIQAABeFAIAQBKFAADwGmh1gJ44d+6CPJ6un0YRFjZM1dV1JiQKHDL2nN3zSfbPaPd8kv0z2ilfUJBDI0d+7zuX9+lC8HiMbhVC63Ptjow9Z/d8kv0z2j2fZP+Mds/Xil1GAABJFAIAwItCAABI6oVCqKurU3Jysk6dOtVm2Weffaa0tDQlJiZqyZIlam5uNjsOAOA7mFoIn376qebMmaPy8vKrLl+0aJGWLVumHTt2yDAMbd682cw4AIB2mFoImzdv1vLly+Vyudosq6ioUGNjo8aPHy9JSktLU2FhoZlxAKBdZRW1+p895SqrqA3Y+rZ8dCxg65OkXYcqtGbTQe06VBGwdbYy9WOnK1as+M5lVVVVcjqdvtnpdKqystLMOAAstuXjMh04VqWJt7n0wLSoHq+vrKJWu0pOa1RYiKIiR/R4Xa9s/IeaWwwNHODQM3Nv79E6yypq9Z9vHVRLi0cDBgRp0ZwJPc6461CF/lT4uSSp9P+dkyTdNT6yR+u8nGXnIXg8HjkcDt9sGIbf3BlhYcO6/fpO5/BuP7e3kLHn7J5Psm/GNfn7deBolSaOdWnhQz/p8fre/HOpPvz7SUnSh38/qZCQQfpFcnS313e0vEavvn1Qzc0eDRwYpBULfqqxo6/r9vq2/J8Tam759nyB5hZDB8uqNXn8qG6vb1fJabW0eOQxJLV4dKq6vkfrk6TDJ2razA/MGNujdV7OskKIiIiQ2+32zWfPnr3qrqX2VFfXdeuED6dzuNzu811+Xm8iY8/ZPZ9k34x520u198i3W+y7/lGhxsZmzU/t/i9vSdq5/2Sb+d5JN3V7fXtLKtTU7JFhSM3NHu0tqVDY94K7vb6Ghktt5p58b0aFhWjAgCDJu4UwKiykx9/rmFuu08Fjbr+5K+sMCnK0+4e0ZYUQGRmpwYMH68CBA5o4caK2bdum+Ph4q+IAfV4gd8fsO1rVZu5pIThHDNW585f85p4Yc9NIDRwQ5NslM+amkT1a35SYG/TJ4dPytBgKGuDQlJgberS+qMgRWjRngk5V1wdkl5b0791DBz6v0sQxroDuLpIsKISsrCxlZ2crJiZGr776qpYuXaq6ujpFR0crMzOzt+MA14QtH5f57Y6R1KNSuHLLOxCXXkifFqVVGw7IMCSH49u5JwL9CzcqcoSenXu7Pj95TmNuGhmQX+BRkSM0efyogG4F3jU+MuBF0MphGEbfuMjGVbDLyFp2z2j3fFLgMj7+6i41NXt8c/DAIP3X/76r++v7z4/V1PLvn63gAQ7916JpPYko6dsDrYH8hSvZ//tsp3y23WUEIHAuL4OrzV01/Sc3+rY0WudAiIocEbAiQOBRCADaaN3ddOhf1Rr/g7CAfEQU9kchALiqB6ZF6Vf/McE2uztgPgoBsMiatw/qeEWtbo0coYUZE6yOA3C1U8AKa94+qNLyc7rU5FFp+Tmtefug1ZEACgGwwpEvzrU7d1X06JHtzkBnUAiABa78sHdPP/y9MGOCokePVPDAIEWPHskuKHQLxxCAawQlgJ5iCwEAIIlCAAB4UQgAAEkUAgDAi0IAAEiiEAAAXhQCAEAShQB0SllFrf5nT7nKKmqtjgKYhhPTgA6UVdT63enr+Ycnck1/XJPYQgA68N8Fpb5LSxjGtzNwLaIQgA5Ufd3Y7twdXIwOdkQhABZovRjdoGAuRgf74BgCYJGFGRNsdQN2gC0EAIAkCgEA4EUhAAAkUQgAAC8KAQAgiUIAAHhRCAAASRQCAMCLQgAASKIQAABephZCQUGBkpKSlJCQoPz8/DbLS0tLNXv2bKWmpurxxx/XN998Y2YcAEA7TCuEyspK5ebmauPGjdq6das2bdqksrIyv8esWLFC2dnZ2r59u77//e9r3bp1ZsUBAHTAtEIoLi5WXFycQkNDFRISosTERBUWFvo9xuPx6MKFC5KkhoYGDRkyxKw4AIAOmHa106qqKjmdTt/scrlUUlLi95jnnntOjzzyiFauXKmhQ4dq8+bNXXqNsLBh3c7ndA7v9nN7Cxl7zqx8gVxvf30PA8nuGe2er5VpheDxeORwOHyzYRh+c2Njo5YsWaI333xTsbGxWr9+vZ599lnl5eV1+jWqq+vk8RhdztYXLjlMxp4zM1+g1tuf38NAsXtGO+ULCnK0+4e0abuMIiIi5Ha7fbPb7ZbL5fLNx44d0+DBgxUbGytJevDBB7Vv3z6z4gAAOmBaIUyZMkV79uxRTU2NGhoaVFRUpPj4eN/ym2++WWfOnNGJEyckSR999JFiYmLMioN+pqyiVls+OqayilqrowB9hmm7jMLDw5WTk6PMzEw1NTUpPT1dsbGxysrKUnZ2tmJiYrRq1So9/fTTMgxDYWFhWrlypVlx0I+UVdRq1YYDMgzJ4ZCef3iioiJHWB0LsD1Tb6GZkpKilJQUv6+tXbvW9/+pU6dq6tSpZkZAP/TfBaUyvIeWDOPbefWCKd1eX5BDuvxQVZDjux8L9GWcqYxrTtXXje3OXfW/fhje7gxcK0zdQgCuBfNToyVJh09UK+aWMN8MXGsoBKATKAH0B+wyAgBIohAAAF4UAgBAEoUAAPCiEAAAkigEAIAXhQAAkEQhAAC8KAQAgCQKAQDgRSEAACRRCAAALwoBACCJQgAAeFEIAABJFAIAwItCAABIohAAAF4UAgBAEoUAAPCiEAAAkigEAIAXhQAAkEQhAAC8KAQAgCQKAQDgRSEAACSZXAgFBQVKSkpSQkKC8vPz2yw/ceKE5s2bp9TUVD366KOqra01Mw4AoB2mFUJlZaVyc3O1ceNGbd26VZs2bVJZWZlvuWEY+uUvf6msrCxt375dP/zhD5WXl2dWHABAB0wrhOLiYsXFxSk0NFQhISFKTExUYWGhb3lpaalCQkIUHx8vSVqwYIEeeughs+IAADow0KwVV1VVyel0+maXy6WSkhLffPLkSV1//fVavHixPvvsM91yyy164YUXzIoDm8vbXqrDJ6oVc0uY5qdGWx0H6JdMKwSPxyOHw+GbDcPwm5ubm7Vv3z5t2LBBMTEx+t3vfqfVq1dr9erVnX6NsLBh3c7ndA7v9nN7S3/JuCZ/v/YeqZQk7T1SqSFDBmrhQz/p8XovZ+f30s7ZJPvnk+yf0e75WplWCBEREdq/f79vdrvdcrlcvtnpdOrmm29WTEyMJCk5OVnZ2dldeo3q6jp5PEaXszmdw+V2n+/y83pTf8r4fw991WbOTOj+el2hQ1T1daPfbNf30u7fZ7vnk+yf0U75goIc7f4hbdoxhClTpmjPnj2qqalRQ0ODioqKfMcLJGnChAmqqanR0aNHJUk7d+5UdDS7CvqjK0u9OyV/ucdSotudAVydaVsI4eHhysnJUWZmppqampSenq7Y2FhlZWUpOztbMTEx+uMf/6ilS5eqoaFBEREReuWVV8yKg34kKnKEFs+bqFPV9RoVFqKoyBFWRwL6BNMKQZJSUlKUkpLi97W1a9f6/j9u3Di98847ZkZAPxUVOUKTx4+yzaY60BdwpjIAQBKFAADwohAAAJIoBACAF4UAAJBEIQAAvCgEAIAkCgEA4EUhAAAkUQgAAC8KAQAgiUIAAHhRCAAASZ0ohCeffFLFxcW9kQUAYKEOC2HGjBl64403lJiYqHXr1unrr7/ujVwAgF7WYSGkpqZqw4YNeuONN1RdXa309HQtWrRIJSUlvZEPANBLOnUMwePx6IsvvlB5eblaWloUFhamF198Ua+99prZ+QAAvaTDO6bl5ubqvffe04033qi5c+fq97//vYKDg1VfX69p06YpOzu7N3ICAEzWYSHU1NRo7dq1Gjt2rN/XQ0JCtGbNGtOCAQB6V4eF8NJLL33nsjvvvDOgYQAA1uE8BACAJAoBAOBFIQAAJHXiGAJwNXnbS1VaXqPo0ddpfmq01XEABACFgC7L216qvUcqJcn3L6UA9H3sMkKX7Tta1e7cVa7QIe3OAHoHhYAu83iMdueueiwlut0ZQO9glxEsFxU5QovnTdTnJ89pzE0jFRU5wupIQL9EIcAWoiJHUASAxdhlBACQRCEAALxMLYSCggIlJSUpISFB+fn53/m4Xbt26e677zYzCgCgA6YdQ6isrPRdOnvQoEHKyMjQpEmTFBUV5fe4s2fP6uWXXzYrBgCgk0zbQiguLlZcXJxCQ0MVEhKixMREFRYWtnnc0qVL9etf/9qsGACATjJtC6GqqkpOp9M3u1yuNrfd/NOf/qQf/ehHGjduXLdeIyxsWLfzOZ3Du/3c3tIXMraya1a75rqc3TPaPZ9k/4x2z9fKtELweDxyOBy+2TAMv/nYsWMqKirSm2++qTNnznTrNaqr67p1UpTTOVxu9/luvWZv6QsZL2fHrH3hPbR7Rrvnk+yf0U75goIc7f4hbdouo4iICLndbt/sdrvlcrl8c2Fhodxut2bPnq358+erqqpKc+fONSsOAKADphXClClTtGfPHtXU1KihoUFFRUWKj4/3Lc/OztaOHTu0bds25eXlyeVyaePGjWbFAQB0wLRCCA8PV05OjjIzM3XfffcpOTlZsbGxysrK0uHDh816WQBAN5l66YqUlBSlpKT4fW3t2rVtHjdq1Cjt3LnTzCgAgA5wpjIAQBKFAADwohAAAJIoBACAF4UAAJBEIQAAvCgEAIAkCgEA4EUhAAAkUQgAAC8KAQAgiUIAAHhRCAAASRQCAMCLQgAASKIQAABeFAIAQBKFAADwohAAAJJMvqcy7GHLx2U6cKxKE29z6YFpUVbHAWBTFMI1bsvHZfrw7yclyfcvpQDgathldI3764FT7c7dETJ4QLszgL6JQrjGNTV72p27I/2KLYwrZwB9E7uM0GV3jY+UJB0+UaOYW67zzQD6NgoB3XLX+Eg9MGOs3O7zVkcBECDsMgIASKIQAABeFAIAQBKFAADwohAAAJJMLoSCggIlJSUpISFB+fn5bZb/9a9/1axZs5Samqpf/epXqq2tNTMOAKAdphVCZWWlcnNztXHjRm3dulWbNm1SWVmZb3ldXZ1efPFF5eXlafv27RozZoxef/11s+IAADpgWiEUFxcrLi5OoaGhCgkJUWJiogoLC33Lm5qatHz5coWHh0uSxowZo9OnT5sVBwDQAdMKoaqqSk6n0ze7XC5VVlb65pEjR2rGjBmSpMbGRuXl5Wn69OlmxQEAdMC0M5U9Ho8cDodvNgzDb251/vx5PfHEExo7dqzuv//+Lr1GWNiwbudzOod3+7m9xayMgVyv3d9Hu+eT7J/R7vkk+2e0e75WphVCRESE9u/f75vdbrdcLpffY6qqqvToo48qLi5Oixcv7vJrVFfXyeMxuvw8p3O47S+5YGbGQK3X7u+j3fNJ9s9o93yS/TPaKV9QkKPdP6RN22U0ZcoU7dmzRzU1NWpoaFBRUZHi4+N9y1taWrRgwQLNnDlTS5YsuerWAwCg95i2hRAeHq6cnBxlZmaqqalJ6enpio2NVVZWlrKzs3XmzBkdOXJELS0t2rFjhyTpxz/+sVasWGFWJABAO0y92mlKSopSUlL8vrZ27VpJUkxMjI4ePWrmywMAuoAzlQEAkigEAIAXhQAAkEQhAAC8KAQAgCQKAQDgRSEAACRRCAAAL1NPTEP3rHn7oI5X1OrWyBFamDHB6jgA+gm2EGxmzdsHVVp+TpeaPCotP6c1bx+0OhKAfoJCsJkjX5xrd+6q6NEj250BoBWFYDOG0f7cVQszJih69EgFDwxS9OiR7IIC8J04htAPUAIAOoMtBACAJAoBAOBFIQAAJFEIAAAvCgEAIIlCAAB4UQgAAEkUAgDAi0IAAEiiEAAAXhQCAEAShQAA8OLidgGw5u2DOnaqVreN4oY2APouthB6qPWGNk3N3NAGQN9GIfRQoG9o4wod0u4MAGahEHoo0De0eSwlut0ZAMzCMQSbiYococXzJupUdb1GhYUoKnKE1ZEA9BMUgg1FRY7Q5PGj5HaftzoKgH7E1F1GBQUFSkpKUkJCgvLz89ss/+yzz5SWlqbExEQtWbJEzc3NZsYBALTDtEKorKxUbm6uNm7cqK1bt2rTpk0qKyvze8yiRYu0bNky7dixQ4ZhaPPmzWbF8cnbXqq5L3ygvO2lpr8WAPQlphVCcXGx4uLiFBoaqpCQECUmJqqwsNC3vKKiQo2NjRo/frwkKS0tzW+5GfK2l2rvkUqdr2/S3iOVlAIAXMa0QqiqqpLT6fTNLpdLlZWV37nc6XT6LTfD4RPV7c7dwcdEAVwrTDuo7PF45HA4fLNhGH5zR8s7IyxsWJcef8ePwrXrHxV+s9M5vEvruNKieXdo0euf+M09XWerQK3HTHbPaPd8kv0z2j2fZP+Mds/XyrRCiIiI0P79+32z2+2Wy+XyW+52u33z2bNn/ZZ3RnV1nTyezn/wPzNhjBobm1VaXqPo0dcpM2FMjz/JE/a9YC2eN1GfnzynMTeNVNj3ggPy6SCnc7jtP2Vk94x2zyfZP6Pd80n2z2infEFBjnb/kDatEKZMmaLXX39dNTU1Gjp0qIqKivTSSy/5lkdGRmrw4ME6cOCAJk6cqG3btik+Pt6sOD7zU6MD/g2KihzB+QIA+jzTjiGEh4crJydHmZmZuu+++5ScnKzY2FhlZWXp8OHDkqRXX31Vq1at0j333KP6+nplZmaaFQcA0AGHYfT0YgvW6eouo1Z22oT7LmTsObvnk+yf0e75JPtntFO+jnYZcS0jAIAkCgEA4EUhAAAk9fGL2wUFde28hUA9t7eQsefsnk+yf0a755Psn9Eu+TrK0acPKgMAAoddRgAASRQCAMCLQgAASKIQAABeFAIAQBKFAADwohAAAJIoBACAF4UAAJDUDwuhoKBASUlJSkhIUH5+vtVx2vjDH/6ge++9V/fee69eeeUVq+O06+WXX9Zzzz1ndYw2du7cqbS0NM2cOVO/+c1vrI5zVdu2bfN9n19++WWr4/jU1dUpOTlZp06dkiQVFxcrJSVFCQkJys3NtTjdt67MuGnTJiUnJyslJUXPP/+8Ll26ZKt8rTZs2KB58+ZZlKqTjH7kzJkzxrRp04xz584ZFy5cMFJSUozjx49bHctn9+7dxoMPPmhcvHjRuHTpkpGZmWkUFRVZHeuqiouLjUmTJhnPPvus1VH8nDx50rjzzjuN06dPG5cuXTLmzJlj7Nq1y+pYfurr64077rjDqK6uNpqamoz09HRj9+7dVscyDh06ZCQnJxvR0dHGl19+aTQ0NBhTp041Tp48aTQ1NRmPPPKI5e/llRlPnDhhzJgxwzh//rzh8XiMZ555xli/fr1t8rU6fvy48bOf/cx4+OGHLcvWGf1qC6G4uFhxcXEKDQ1VSEiIEhMTVVhYaHUsH6fTqeeee06DBg1ScHCwfvCDH+irr76yOlYbX3/9tXJzc7VgwQKro7Txl7/8RUlJSYqIiFBwcLByc3M1btw4q2P5aWlpkcfjUUNDg5qbm9Xc3KzBgwdbHUubN2/W8uXLffc2Lykp0c0336wbb7xRAwcOVEpKiuU/L1dmHDRokJYvX65hw4bJ4XDotttus/Rn5sp8knTp0iUtW7ZM2dnZluXqrD59tdOuqqqqktPp9M0ul0slJSUWJvJ36623+v5fXl6uDz/8UG+99ZaFia5u2bJlysnJ0enTp62O0sYXX3yh4OBgLViwQKdPn9Zdd92lp59+2upYfoYNG6annnpKM2fO1NChQ3XHHXfo9ttvtzqWVqxY4Tdf7eelsrKyt2P5uTJjZGSkIiMjJUk1NTXKz8/XqlWrrIgmqW0+SVqzZo1mz56tUaNGWZCoa/rVFoLH45HD8e/LvxqG4TfbxfHjx/XII4/omWee0ejRo62O42fLli264YYbNHnyZKujXFVLS4v27NmjlStXatOmTSopKdH7779vdSw/R48e1bvvvquPP/5Yn3zyiYKCgrRu3TqrY7XRV35eJKmyslI///nPNXv2bE2aNMnqOD67d+/W6dOnNXv2bKujdEq/KoSIiAi53W7f7Ha7/Tbt7ODAgQP6xS9+oYULF+r++++3Ok4bH3zwgXbv3q1Zs2bptdde086dO7Vy5UqrY/lcf/31mjx5sq677joNGTJE06dPt9VWoCT97W9/0+TJkxUWFqZBgwYpLS1N+/btszpWG33h50WS/vWvfykjI0P333+/nnjiCavj+Pnzn/+s48ePa9asWVq6dKn++c9/2m6L1Y/VBzF6U+tB5erqaqO+vt5ITU01Pv30U6tj+Xz11VfGpEmTjOLiYqujdMq7775ru4PKhw4dMhITE43a2lqjubnZePzxx43NmzdbHcvPJ598YqSmphoXLlwwPB6P8cILLxivvfaa1bF8pk2bZnz55ZdGY2OjER8fb5SXlxvNzc3Go48+anzwwQdWxzMM498Zz58/b0ydOtV4//33rY7kpzXf5fbu3Wv7g8r96hhCeHi4cnJylJmZqaamJqWnpys2NtbqWD7r1q3TxYsXtXr1at/XMjIyNGfOHAtT9S3jxo3TY489prlz56qpqUk//elPbbe5fuedd+rIkSNKS0tTcHCwYmJiNH/+fKtjtTF48GCtXr1aTz75pC5evKipU6fqnnvusTqWn3feeUdnz57V+vXrtX79eknS3XffraeeesriZH0Td0wDAEjqZ8cQAADfjUIAAEiiEAAAXhQCAEAShQAA8KIQAACSKAQAgBeFAATI+++/r+nTp+vChQuqr6/XzJkztXXrVqtjAZ3GiWlAAC1cuFDDhw/XpUuXNGDAAL300ktWRwI6jUIAAqiurk6zZs3SkCFD9N5779niPgdAZ7HLCAig6upqXbx4Ud98842qqqqsjgN0CVsIQIA0NTUpIyNDGRkZ8ng82rJli9566y0FBwdbHQ3oFLYQgAD57W9/q+uvv14PPPCAHnzwQY0cOdI2N6YHOoMtBACAJLYQAABeFAIAQBKFAADwohAAAJIoBACAF4UAAJBEIQAAvCgEAIAk6f8DKyzq/NQZSl8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute ECDF: x, y\n", "x, y = ecdf(n_defaults)\n", "\n", "# Plot the ECDF with labeled axes\n", "_ = plt.plot(x, y, marker='.', linestyle='none')\n", "_ = plt.xlabel('x')\n", "_ = plt.ylabel('y')\n", "\n", "# Compute the number of 100-loan simulations with 10 or more defaults: n_lose_money\n", "n_lose_money = np.sum(n_defaults >= 10)\n", "\n", "# Compute and print probability of losing money\n", "print('Probability of losing money =', n_lose_money / len(n_defaults))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probability distributions and stories: The Binomial distribution\n", "- Probability mass function (PMF)\n", " - The set of probabilities of discrete outcomes\n", "- Probability distribution\n", " - A mathmatical description of outcomes\n", "- Binomial distribution\n", " - The number r of successes in n Bernoulli trials with probability p of success, is Binomially distributed\n", " - The number r of heads in 4 coin flips with probability 0.5 of heads, is Binomially distributed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sampling out of the Binomial distribution\n", "Compute the probability mass function for the number of defaults we would expect for 100 loans as in the last section, but instead of simulating all of the Bernoulli trials, perform the sampling using ```np.random.binomial()```. This is identical to the calculation you did in the last set of exercises using your custom-written ```perform_bernoulli_trials()``` function, but far more computationally efficient. Given this extra efficiency, we will take 10,000 samples instead of 1000. After taking the samples, plot the CDF as last time. This CDF that you are plotting is that of the Binomial distribution.\n", "\n", "Note: For this exercise and all going forward, the random number generator is pre-seeded for you (with ```np.random.seed(42)```) to save you typing that each time." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYD0lEQVR4nO3dfXBU9aHG8WeDIZgLc6NhN+mgg21tpdIE++IQKE2EgURiNkiIw4s1dqABrBrNMGgqCE4tEKw0LbU6TQZwvAQnEeUlvTakbeSOJXEYmIuxvCiUiyjCZkkUCQmQZM/9o4etG/JC4p6cHPP9/OP+cnbPeWaH+Ozv/HLOugzDMAQAGPQi7A4AABgYKAQAgCQKAQBgohAAAJIoBACAiUIAAEiiEAAApuvsDvBlfPrpBQUCvb+MIjZ2uBoamixIZA0n5XVSVslZeZ2UVXJWXidllfqeNyLCpRtu+I8utzu6EAIBo0+FcOW1TuKkvE7KKjkrr5OySs7K66SskjV5OWUEAJBEIQAATBQCAEBSPxRCU1OTMjIy9PHHH1+17fDhw8rKylJaWpqWLVumtrY2q+MAALpgaSG8++67mjt3rk6cONHp9qVLl2rFihXatWuXDMNQeXm5lXEAAN2wtBDKy8u1cuVKeTyeq7adOnVKFy9e1B133CFJysrKUmVlpZVxADjc7gOntK7sf7X7wKmw7XN+YbW8S3ZofmF12PYpSYt+/ZbmF1Zr0a/fCut+i3ce1Lyn31TxzoNh3a9kcSGsWrVKP/zhDzvdVl9fL7fbHRy73W75fD4r4wBwsN0HTumVyvd18P8+1SuV74elFDqWQLhKYdGv31Jr+7/+LLS13QhbKRTvPKh3Dvl0vrlV7xzyhb0UbLsOIRAIyOVyBceGYYSMr0Vs7PA+H9/tHtHn19rBSXmdlFVyVl6nZPUu2RF8XLFuRlj2+Url+1eN75s2Jiz7/qJwvMdXyuCL43Ds9+CJxqvG4fw3YVshxMfHy+/3B8dnz57t9NRSdxoamvp0cYbbPUJ+//lev84uTsrrpKySs/I6JWvHT9neJTu0sWCKJcey4v0Ixz4jh7hCSiFyiCss+x17y41655AvZNyb/UZEuLr9IG3bn52OGjVKUVFR2r9/vyRpx44dSk5OtisOgEGoY1GFq7j+uHSyIof864xH5BCX/rh0clj2uzBzrJJuj9OI6Egl3R6nhZljw7LfK/p9hpCbm6u8vDwlJCTo+eef1/Lly9XU1KSxY8cqJyenv+MAg9oXP81b9Sk+XKKjhqj5UnvIOBw2FkyxZPYVrhLoaGHmWMtmi/1SCNXV//5HV1JSEnw8ZswYbd26tT8iAOigswXVgVwKL+Sn6JGi/1HzpXZFRw3RC/kpdkf6ynH0ze0ADC6UgLW4dQUAQBKFACDMrFqohfU4ZQQg7KxaqIW1mCEAACQxQwAcwUl/HgrnYoYADHBW3W8H6IhCAABIohCAQSvn7tu6HWPwYQ0BGKTuumOUJGn/+/X6wW2e4BiDF4UADGJ33TGKIkAQp4wAAJIoBACAiUIAAEiiEAAAJgoBACCJQgAAmCgEAIAkrkMAwo4b0cGpmCEAYcSN6OBkFAIAQBKFAAAwUQjAAMd3FKO/sKgMOADfUYz+wAwBACCJQgAAmCgEAIAkCgEAYKIQAACSKAQAgIlCAABIsrgQKioqlJ6ertTUVJWWll61/eDBg5o1a5YyMzO1aNEiff7551bGAQB0w7JC8Pl8Kioq0pYtW7R9+3aVlZXp2LFjIc9ZtWqV8vLytHPnTn3961/Xhg0brIoDAOiBZYVQU1OjpKQkxcTEKDo6WmlpaaqsrAx5TiAQ0IULFyRJLS0tGjZsmFVxAAA9sOzWFfX19XK73cGxx+NRXV1dyHMKCgo0f/58rV69Wtdff73Ky8t7dYzY2OF9zud2j+jza+3gpLxOyipZnzec++e9tY6TskrW5LWsEAKBgFwuV3BsGEbI+OLFi1q2bJlefvllJSYmatOmTXryySdVXFx8zcdoaGhSIGD0OpvT7gnjpLxOyir1T95w7Z/31jpOyir1PW9EhKvbD9KWnTKKj4+X3+8Pjv1+vzweT3D8wQcfKCoqSomJiZKk2bNna+/evVbFAQD0wLJCmDhxompra9XY2KiWlhZVVVUpOTk5uH306NE6c+aMjh8/Lkn629/+poSEBKviAAB6YNkpo7i4OOXn5ysnJ0etra3Kzs5WYmKicnNzlZeXp4SEBK1Zs0aPP/64DMNQbGysVq9ebVUcAEAPLP0+BK/XK6/XG/KzkpKS4OOUlBSlpKRYGQHoUvHOgzp4olFjb7lRCzPH2h0HsB1fkINBqXjnQb1zyCdJwf9SChjsuHUFBqUrJdDVGBiMKAQAgCQKAQirjQVTuh0DAxlrCECYUQJwKmYIAABJFAIAwEQhAAAkUQgAABOFAACQRCEAAEwUAgBAEoUAADBRCAAASRQCAMBEIQAAJFEIAAAThQAAkEQhAABMFAIAQBKFAAAwUQgAAEkUAgDARCEAACRRCAAAE4UAAJBEIQAATBQCAEAShQAHOHbqnP679oSOnTpndxTgK+06uwMA3Tl26pxW/9f+4PipB36gW0f9p42JgK8uS2cIFRUVSk9PV2pqqkpLS6/afvz4cT3wwAPKzMzUggULdO4cnwAR6otl0NkYQPhYVgg+n09FRUXasmWLtm/frrKyMh07diy43TAMPfTQQ8rNzdXOnTv1ne98R8XFxVbFAQD0wLJCqKmpUVJSkmJiYhQdHa20tDRVVlYGtx88eFDR0dFKTk6WJC1evFj333+/VXGAEBsLpnQ7BgYjy9YQ6uvr5Xa7g2OPx6O6urrg+OTJkxo5cqSeeuopHT58WN/4xjf09NNPWxUHuMrGgilyu0fI7z9vdxRgQLCsEAKBgFwuV3BsGEbIuK2tTXv37tXmzZuVkJCg3/72tyosLFRhYeE1HyM2dnif87ndI/r8Wjs4Ka/VWcO9f95b6zgpr5OyStbktawQ4uPjtW/fvuDY7/fL4/EEx263W6NHj1ZCQoIkKSMjQ3l5eb06RkNDkwIBo9fZnPap0El5+yNrOPfPe2sdJ+V1Ulap73kjIlzdfpC2bA1h4sSJqq2tVWNjo1paWlRVVRVcL5Ck733ve2psbNSRI0ckSdXV1Ro7dqxVcQAAPbBshhAXF6f8/Hzl5OSotbVV2dnZSkxMVG5urvLy8pSQkKA//OEPWr58uVpaWhQfH6/nnnvOqjgAgB5YemGa1+uV1+sN+VlJSUnw8bhx47R161YrIwAArhG3rgAASKIQAAAmCgEAIIlCAACYKAQAgCQKAQBgohAAAJIoBACAiUIAAEiiEAAAJgoBACCJQgAAmCgEAICkayiERx99VDU1Nf2RBQBgox4LYdq0aXrxxReVlpamDRs26LPPPuuPXACAftZjIWRmZmrz5s168cUX1dDQoOzsbC1dulR1dXX9kQ8A0E+uaQ0hEAjoww8/1IkTJ9Te3q7Y2Fg988wzWr9+vdX5AAD9pMdvTCsqKtIbb7yhm2++WfPmzdPvfvc7RUZGqrm5WZMnT1ZeXl5/5AQAWKzHQmhsbFRJSYnGjBkT8vPo6GitW7fOsmAAgP7VYyE8++yzXW6bNGlSWMMAAOzDdQgAAEkUAgDARCEAACRRCAAAE4UAAJBEIQAATD3+2SlwrR5at1uXWgOKiozQS0vusjsOgF5ihoCwuFIGknSpNaCH1u22NxCAXqMQEBZXyqCrMYCBj0IAAEiiEDDAbSyY0u0YQPhYuqhcUVGhl156SW1tbXrwwQd1//33d/q83bt365e//KWqq6utjAOHogSA/mFZIfh8vuCts4cOHao5c+Zo/PjxuvXWW0Oed/bsWa1du9aqGACAa2TZKaOamholJSUpJiZG0dHRSktLU2Vl5VXPW758uR555BGrYgAArpFlM4T6+nq53e7g2OPxXPW1m6+88opuv/12jRs3rk/HiI0d3ud8bveIPr/WDk7LKzkns1NySs7KKjkrr5OyStbktawQAoGAXC5XcGwYRsj4gw8+UFVVlV5++WWdOXOmT8doaGhSIGD0+nVu9wj5/ef7dEw7OC3vFU7I7KT31klZJWfldVJWqe95IyJc3X6QtuyUUXx8vPx+f3Ds9/vl8XiC48rKSvn9fs2aNUsLFy5UfX295s2bZ1UcAEAPLCuEiRMnqra2Vo2NjWppaVFVVZWSk5OD2/Py8rRr1y7t2LFDxcXF8ng82rJli1VxAAA9sKwQ4uLilJ+fr5ycHN17773KyMhQYmKicnNz9d5771l1WABAH1l6HYLX65XX6w35WUlJyVXPu+mmm7gGAQBsxpXKAABJFAIAwEQhAAAkUQgAABOFAACQRCEAAEwUAgBAEoUAADBRCAAASRQCAMBEIQAAJFEIAAAThQAAkEQhAABMFAIAQBKFAAAwUQgAAEkUAgDARCEAACRRCAAA03V2B4A95hdWBx9vLJhiYxIAAwUzhEHoi2XQ2RjA4EQhAAAkUQgAABOFgLDouA7BugTgPCwqI2w2FkyR2z1Cfv95u6MA6ANmCAAASRQCAMBEIQAAJFEIAACTpYVQUVGh9PR0paamqrS09Krtf/3rXzVjxgxlZmbq5z//uc6dO2dlHABANywrBJ/Pp6KiIm3ZskXbt29XWVmZjh07Ftze1NSkZ555RsXFxdq5c6duu+02/f73v7cqDgCgB5YVQk1NjZKSkhQTE6Po6GilpaWpsrIyuL21tVUrV65UXFycJOm2227T6dOnrYoDAOiBZYVQX18vt9sdHHs8Hvl8vuD4hhtu0LRp0yRJFy9eVHFxsaZOnWpVHABADyy7MC0QCMjlcgXHhmGEjK84f/68Hn74YY0ZM0YzZ87s1TFiY4f3OZ/bPaLPr7WD1XnDuX/eW+s4KavkrLxOyipZk9eyQoiPj9e+ffuCY7/fL4/HE/Kc+vp6LViwQElJSXrqqad6fYyGhiYFAkavX+e0q2n7I2+49s97ax0nZZWclddJWaW+542IcHX7QdqyU0YTJ05UbW2tGhsb1dLSoqqqKiUnJwe3t7e3a/HixZo+fbqWLVvW6ewBANB/LJshxMXFKT8/Xzk5OWptbVV2drYSExOVm5urvLw8nTlzRocOHVJ7e7t27dolSfrud7+rVatWWRUJANANS29u5/V65fV6Q35WUlIiSUpISNCRI0esPDwAoBe4UhkAIIlCAACYKAQAgCQKAQBgohAAAJIoBACAiUIAAEiiEAAAJgoBACDJ4iuV8eXNL6wOPt5YMMXGJAC+6pghDGBfLIPOxgAQThQCAEAShQAAMFEIg1DHtQjWJgBILCoPWpQAgI6YIQAAJFEIAAAThQAAkEQhAABMFAIAQBKFAAAwUQgAAEkUAgDARCEAACRRCAAAE7euCBO+twCA0zFDCAO+twDAVwGFAACQRCEAAEwUwgDG9xYA6E8sKg9wGwumyO0eIb//vN1RAHzFWTpDqKioUHp6ulJTU1VaWnrV9sOHDysrK0tpaWlatmyZ2trarIwjSVpQWC3vkh1awMIvAISwrBB8Pp+Kioq0ZcsWbd++XWVlZTp27FjIc5YuXaoVK1Zo165dMgxD5eXlVsWR9K8yMMzHhjkGAPyLZYVQU1OjpKQkxcTEKDo6WmlpaaqsrAxuP3XqlC5evKg77rhDkpSVlRWy3QpGD2MAGMwsW0Oor6+X2+0Ojj0ej+rq6rrc7na75fP5enWM2NjhXzqn2z3iS++jP/ZrVU4rOCmr5Ky8TsoqOSuvk7JK1uS1rBACgYBcLldwbBhGyLin7deioaFJgcCX+5wfjsXajQVTrrpSOZyLwE5aVHZSVslZeZ2UVXJWXidllfqeNyLC1e0HacsKIT4+Xvv27QuO/X6/PB5PyHa/3x8cnz17NmS7FTr7H3c49w0ATmbZGsLEiRNVW1urxsZGtbS0qKqqSsnJycHto0aNUlRUlPbv3y9J2rFjR8h2q2wsmKKKdTP4HzgAdGBZIcTFxSk/P185OTm69957lZGRocTEROXm5uq9996TJD3//PNas2aN7r77bjU3NysnJ8eqOACAHrgMw3DsH9v0dQ1hsJwvtIOTskrOyuukrJKz8jopq2TdGgK3rgAASKIQAAAmCgEAIMnhN7eLiOjddQvheq0dnJTXSVklZ+V1UlbJWXmdlFXqW96eXuPoRWUAQPhwyggAIIlCAACYKAQAgCQKAQBgohAAAJIoBACAiUIAAEiiEAAAJgoBACBpEBZCRUWF0tPTlZqaqtLSUrvjdOuFF17QPffco3vuuUfPPfec3XGuydq1a1VQUGB3jB5VV1crKytL06dP169+9Su74/Rox44dwX8La9eutTtOp5qampSRkaGPP/5YklRTUyOv16vU1FQVFRXZnC5Ux6xlZWXKyMiQ1+vVL37xC12+fNnmhKE65r1i8+bNeuCBB8J3IGMQOXPmjDF58mTj008/NS5cuGB4vV7j6NGjdsfq1J49e4zZs2cbly5dMi5fvmzk5OQYVVVVdsfqVk1NjTF+/HjjySeftDtKt06ePGlMmjTJOH36tHH58mVj7ty5xu7du+2O1aXm5mbjzjvvNBoaGozW1lYjOzvb2LNnj92xQhw4cMDIyMgwxo4da3z00UdGS0uLkZKSYpw8edJobW015s+fP2De445Zjx8/bkybNs04f/68EQgEjCeeeMLYtGmT3TGDOua94ujRo8aPf/xj4yc/+UnYjjWoZgg1NTVKSkpSTEyMoqOjlZaWpsrKSrtjdcrtdqugoEBDhw5VZGSkvvnNb+qTTz6xO1aXPvvsMxUVFWnx4sV2R+nRX/7yF6Wnpys+Pl6RkZEqKirSuHHj7I7Vpfb2dgUCAbW0tKitrU1tbW2KioqyO1aI8vJyrVy5Mvi96HV1dRo9erRuvvlmXXfddfJ6vQPmd61j1qFDh2rlypUaPny4XC6Xvv3tbw+o37WOeSXp8uXLWrFihfLy8sJ6LEff7bS36uvr5Xa7g2OPx6O6ujobE3XtW9/6VvDxiRMn9Oc//1mvvvqqjYm6t2LFCuXn5+v06dN2R+nRhx9+qMjISC1evFinT5/WXXfdpccff9zuWF0aPny4HnvsMU2fPl3XX3+97rzzTn3/+9+3O1aIVatWhYw7+13z+Xz9HatTHbOOGjVKo0aNkiQ1NjaqtLRUa9assSNapzrmlaR169Zp1qxZuummm8J6rEE1QwgEAnK5/n37V8MwQsYD0dGjRzV//nw98cQTuuWWW+yO06nXXntNX/va1zRhwgS7o1yT9vZ21dbWavXq1SorK1NdXZ22bdtmd6wuHTlyRK+//rreeustvf3224qIiNCGDRvsjtUtJ/6u+Xw+Pfjgg5o1a5bGjx9vd5wu7dmzR6dPn9asWbPCvu9BVQjx8fHy+/3Bsd/vD5mGDTT79+/XT3/6Uy1ZskQzZ860O06X3nzzTe3Zs0czZszQ+vXrVV1drdWrV9sdq0sjR47UhAkTdOONN2rYsGGaOnXqgJ0pStLf//53TZgwQbGxsRo6dKiysrK0d+9eu2N1y2m/a//85z81Z84czZw5Uw8//LDdcbr1pz/9SUePHtWMGTO0fPly/eMf/wjfDDdsqxEOcGVRuaGhwWhubjYyMzONd9991+5Ynfrkk0+M8ePHGzU1NXZH6ZXXX399wC8qHzhwwEhLSzPOnTtntLW1GYsWLTLKy8vtjtWlt99+28jMzDQuXLhgBAIB4+mnnzbWr19vd6xOTZ482fjoo4+MixcvGsnJycaJEyeMtrY2Y8GCBcabb75pd7wQV7KeP3/eSElJMbZt22Z3pG5dyftF77zzTlgXlQfVGkJcXJzy8/OVk5Oj1tZWZWdnKzEx0e5YndqwYYMuXbqkwsLC4M/mzJmjuXPn2pjqq2HcuHH62c9+pnnz5qm1tVU/+tGPLJl+h8ukSZN06NAhZWVlKTIyUgkJCVq4cKHdsboVFRWlwsJCPfroo7p06ZJSUlJ099132x2rU1u3btXZs2e1adMmbdq0SZI0ZcoUPfbYYzYn6398YxoAQNIgW0MAAHSNQgAASKIQAAAmCgEAIIlCAACYKAQAgCQKAQBgohCAMNm2bZumTp2qCxcuqLm5WdOnT9f27dvtjgVcMy5MA8JoyZIlGjFihC5fvqwhQ4bo2WeftTsScM0oBCCMmpqaNGPGDA0bNkxvvPHGgPveAqA7nDICwqihoUGXLl3S559/rvr6ervjAL3CDAEIk9bWVs2ZM0dz5sxRIBDQa6+9pldffVWRkZF2RwOuCTMEIEx+85vfaOTIkbrvvvs0e/Zs3XDDDQPuy+WB7jBDAABIYoYAADBRCAAASRQCAMBEIQAAJFEIAAAThQAAkEQhAABMFAIAQJL0/0guu2cvlDomAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Take 10,000 samples out of the binomial distribution: n_defaults\n", "n_defaults = np.random.binomial(n=100, p=0.05, size=10000)\n", "\n", "# Compute CDF: x, y\n", "x, y = ecdf(n_defaults)\n", "\n", "# Plot the CDF with axis labels\n", "_ = plt.plot(x, y, marker='.', linestyle='none')\n", "_ = plt.xlabel('x')\n", "_ = plt.ylabel('y')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the Binomial PMF\n", "Plotting a nice looking PMF requires a bit of matplotlib trickery that we will not go into here. Instead, we will plot the PMF of the Binomial distribution as a histogram with skills you have already learned. The trick is setting up the edges of the bins to pass to ```plt.hist()``` via the ```bins``` keyword argument. We want the bins centered on the integers. So, the edges of the bins should be -0.5, 0.5, 1.5, 2.5, ... up to ```max(n_defaults) + 1.5```. You can generate an array like this using ```np.arange()``` and then subtracting 0.5 from the array.\n", "\n", "You have already sampled out of the Binomial distribution during your exercises on loan defaults, and the resulting samples are in the NumPy array ```n_defaults```." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEJCAYAAAC61nFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df0yTd+IH8DdgQfulOwZrwXDLTE5P4w80CwuMc7idSgXpOhlmCBssHCjenbrGwdjpdDI2f0TEWya4sGWXzR+TibTXCwPcdsupcFG5qZi5Re/i5ikrFZwWBGnp8/3Dr893HQrlsQ9P0fcrMeHz/Oq7pfVNn/Z5ngBBEAQQERENU6DSAYiIaHRigRARkSQsECIikoQFQkREkrBAiIhIEhYIERFJwgIhIiJJxigdYCRdudINt1vew14iIkLR0dEl623cDX/PB/h/Rn/PBzCjL/h7PkD+jIGBAXjwwf+54/z7qkDcbkH2Arl1O/7M3/MB/p/R3/MBzOgL/p4PUDYjd2EREZEkLBAiIpKEBUJERJKwQIiISBIWCBERScICISIiSVggREQkyX11HAjdWzQPjMPYEN89hXtvuOC41uOz7RHd61ggNGqNDRkDw2qLz7ZnLTPC4bOtEd37uAuLiIgkYYEQEZEkLBAiIpKEBUJERJKwQIiISBIWCBERScICISIiSWQ9DsRqtaKyshIulws5OTnIysq67XJFRUWIj49HWloaOjo6kJubK85zOBy4cuUKvvrqKxw9ehQrVqxAVFQUAGDq1KnYuHGjnHeBiIjuQLYCsdlsKC8vx4EDBxAcHIyMjAzExcVh4sSJHsusX78ezc3NiI+PBwBERETAYrl5cJjb7UZOTg5MJhMA4PTp08jNzcWyZcvkik1ERF6SbRdWU1MT4uPjERYWBrVaDb1ej/r6eo9lrFYr5s6di+Tk5Ntuo6amBuPGjYPBYAAAtLa24vDhwzAYDCgoKEBbW5tc8YmIaAiyFUh7ezu0Wq041ul0sNlsHsvk5eVh8eLFt12/v78fO3fuxOrVq8VpGo0GL7zwAqxWK+bMmSO+MyEiopEn2y4st9uNgIAAcSwIgsd4KIcOHcKECRMwefJkcVpJSYn485IlS1BWVgaHwwGNRuPVNiMiQr2+/buh1XqXRyn+ng9QLqO3t8vH0Df8PaO/5wOUzShbgURFReH48ePi2G63Q6fTeb3+Z599hpSUFHHsdrvx7rvvYunSpQgKChKn//TnoXR0dMHtFrxeXgqtVgO73X9Pyefv+QDvM8rxwvH2du+Vx1BJ/p7R3/MB8mcMDAwY9A9v2XZhJSQkoLm5GZ2dnejp6UFjYyMSExO9Xv/EiROIjY0Vx4GBgTh48CAaGhoAAGazGTNnzoRarfZ5diIiGppsBRIZGQmTyYTs7Gw888wzSE1NRUxMDPLz89Ha2jrk+hcuXBC/rnvL5s2b8eGHH2LhwoWoqalBaWmpXPGJiGgIAYIgyLtPx49wF5ay+Xx9ASgAPr8eCHdhjRx/z+jv+QDld2HxglI0YuS4ABQRKYenMiEiIklYIEREJAkLhIiIJGGBEBGRJCwQIiKShAVCRESSsECIiEgSFggREUnCAiEiIklYIEREJAkLhIiIJGGBEBGRJCwQIiKShAVCRESSsECIiEgSFggREUnCAiEiIklYIEREJImsBWK1WpGSkoKkpCTs3r37jssVFRXhwIED4ri2thazZ8+G0WiE0WhEeXk5AODSpUvIysrCggULsHz5cnR3d8sZn4iIBiFbgdhsNpSXl2PPnj0wm83Yt28fzp07N2CZgoICNDQ0eEw/ffo0iouLYbFYYLFYYDKZAAAbNmxAZmYm6uvrMX36dFRUVMgVn4iIhiBbgTQ1NSE+Ph5hYWFQq9XQ6/Wor6/3WMZqtWLu3LlITk72mN7a2ora2loYDAa8/PLLuHr1KpxOJ44dOwa9Xg8ASEtLG7A9IiIaOWPk2nB7ezu0Wq041ul0OHXqlMcyeXl5AICWlhaP6VqtFrm5uXj00Uexbds2lJSU4JVXXkFoaCjGjBkjLmOz2YaVKSIiVMpdGTatVjMityOVv+dTkrePzWh4DJnx7vl7PkDZjLIViNvtRkBAgDgWBMFjPJgdO3aIP+fl5WH+/PkoKioasL6327ulo6MLbrcwrHWGS6vVwG53yHobd0PJfKPhxejNY+Pvv2OAGX3B3/MB8mcMDAwY9A9v2XZhRUVFwW63i2O73Q6dTjfkeg6HA3/5y1/EsSAICAoKQnh4OBwOB/r7+4e1PSIikodsBZKQkIDm5mZ0dnaip6cHjY2NSExMHHI9tVqN9957DydPngQA7Nq1C/Pnz4dKpUJsbCzq6uoAAGaz2avtERGRPGTbhRUZGQmTyYTs7Gw4nU6kp6cjJiYG+fn5WLlyJWbMmHHb9YKCgrB9+3a8/vrr6O3txYQJE7BlyxYAwPr161FcXIzKykqMHz8e27Ztkys+3Yf6nP0+/wyk94YLjms9dxOLyG/JViAAYDAYYDAYPKZVVVUNWG7Tpk0e49jYWNTW1g5YLjo6Gh999JFvQxL9n2BVEAyrLT7dprXMCP/ei04kHY9EJyIiSVggREQkCQuEiIgkYYEQEZEkLBAiIpJE1m9h0eileWAcxobw6UFEd8b/Iei2xoaMkeUrrUR07+AuLCIikoQFQkREkrBAiIhIEhYIERFJwgIhIiJJWCBERCQJC4SIiCRhgRARkSQsECIikoQFQkREkrBAiIhIEhYIERFJImuBWK1WpKSkICkpCbt3777jckVFRThw4IA4bmlpQXp6OoxGI3JycnDx4kUAwNGjRxEXFwej0Qij0YhXX31VzvhERDQI2c7Ga7PZUF5ejgMHDiA4OBgZGRmIi4vDxIkTPZZZv349mpubER8fL04vLCxERUUFpkyZgv3796O0tBSVlZU4ffo0cnNzsWzZMrliExGRl2R7B9LU1IT4+HiEhYVBrVZDr9ejvr7eYxmr1Yq5c+ciOTlZnNbX14dVq1ZhypQpAIDJkyejra0NANDa2orDhw/DYDCgoKBAnE5ERCNPtncg7e3t0Gq14lin0+HUqVMey+Tl5QG4ucvqluDgYBiNN68b4Xa78c4772DevHkAAI1Gg+TkZCQlJWHv3r0wmUz4+OOPvc4UEREq+f4Mh1arGZHbkcrf891rlHq8R8Pv2d8z+ns+QNmMshWI2+1GQECAOBYEwWM8lL6+PhQXF8Plcom7rEpKSsT5S5YsQVlZGRwOBzQa7x7Ajo4uuN2C1xmk0Go1sNsdst7G3fA232h44YwWSjwf/P15CPh/Rn/PB8ifMTAwYNA/vGXbhRUVFQW73S6O7XY7dDqdV+t2d3cjLy8PLpcLlZWVUKlUcLvdqKysRH9/v8eyQUFBPs1NRETeka1AEhIS0NzcjM7OTvT09KCxsRGJiYlerVtYWIhHHnkE27dvR3Bw8M2ggYE4ePAgGhoaAABmsxkzZ86EWq2W6y4QEdEgZNuFFRkZCZPJhOzsbDidTqSnpyMmJgb5+flYuXIlZsyYcdv1vv76a3z++eeYOHEiFi1aBODm5ydVVVXYvHkzXnvtNezYsQPh4eHYsmWLXPGJiGgIshUIABgMBhgMBo9pVVVVA5bbtGmT+PPUqVPx7bff3nZ7kyZNGtaH5kREJB8eiU5ERJKwQIiISBIWCBERScICISIiSVggREQkCQuEiIgkYYEQEZEkLBAiIpKEBUJERJKwQIiISBIWCBERScICISIiSWQ9mSLR/a7P2e/Ti3P13nDBca3HZ9sjuhssECIZBauCYFht8dn2rGVG+Pc18uh+wl1YREQkCQuEiIgkYYEQEZEkLBAiIpJE1gKxWq1ISUlBUlISdu/efcflioqKcODAAXF86dIlZGVlYcGCBVi+fDm6u7sBANeuXcPSpUuRnJyMrKws2O12OeMTEdEgZCsQm82G8vJy7NmzB2azGfv27cO5c+cGLFNQUICGhgaP6Rs2bEBmZibq6+sxffp0VFRUAAC2b9+O2NhYfPrpp1i8eDHefPNNueITEdEQZCuQpqYmxMfHIywsDGq1Gnq9HvX19R7LWK1WzJ07F8nJyeI0p9OJY8eOQa/XAwDS0tLE9b788ksYDAYAQGpqKv7xj3/A6XTKdReIiGgQshVIe3s7tFqtONbpdLDZbB7L5OXlYfHixR7Trly5gtDQUIwZc/MQFa1WK673022OGTMGoaGh6OzslOsuEBHRIGQ7kNDtdiMgIEAcC4LgMb6T2y13p/UEQUBgoPcdGBER6vWyd8OXRx7Lwd/z0eC8/f2Nht+zv2f093yAshmHLJAVK1ZgyZIlSEhIGNaGo6KicPz4cXFst9uh0+mGXC88PBwOhwP9/f0ICgryWE+n0+Hy5cuIioqCy+VCd3c3wsLCvM7U0dEFt1sY1v0YLq1WA7vdf48V9jbfaHjh3K+8/f358/MQ8P+M/p4PkD9jYGDAoH94D/nn+/z581FRUQG9Xo/3338fP/74o1c3nJCQgObmZnR2dqKnpweNjY1ITEwccj2VSoXY2FjU1dUBAMxms7jenDlzYDabAQB1dXWIjY2FSqXyKg8REfnWkAXy9NNPY9euXaioqEBHRwfS09NRWFiIU6dODbpeZGQkTCYTsrOz8cwzzyA1NRUxMTHIz89Ha2vroOuuX78e1dXVSElJwfHjx/HSSy8BAFatWoUTJ05g4cKF2LNnD9atWzeMu0pERL7k1Wcgbrcb3333Hc6fP4/+/n5ERETg9ddfx5NPPomVK1fecT2DwSB+a+qWqqqqActt2rTJYxwdHY2PPvpowHJhYWHYuXOnN5GJiEhmQxZIeXk5Dhw4gIcffhiZmZn485//DJVKhevXr+Opp54atECIiOjeNWSBdHZ2oqqqClOmTPGYrlarUVZWJlswIiLyb0MWyBtvvHHHebNnz/ZpGCIiGj14MkUiIpKEBUJERJKwQIiISBIWCBERScICISIiSVggREQkCQuEiIgkYYEQEZEkLBAiIpKEBUJERJKwQIiISBIWCBERSSLbNdFpZGkeGIexId79Onm5WiLyBRbIPWJsyBgYVlt8tj1rmdFn2yKiexN3YRERkSQsECIikkTWXVhWqxWVlZVwuVzIyclBVlaWx/wzZ85gzZo16O7uRmxsLDZs2ICrV68iNzdXXMbhcODKlSv46quvcPToUaxYsQJRUVEAgKlTp2Ljxo1y3gUiIroD2QrEZrOJ11MPDg5GRkYG4uLiMHHiRHGZwsJClJaWYtasWfjTn/6E6upqZGZmwmK5uS/f7XYjJycHJpMJAHD69Gnk5uZi2bJlcsUmIiIvybYLq6mpCfHx8QgLC4NarYZer0d9fb04/+LFi+jt7cWsWbMAAGlpaR7zAaCmpgbjxo2DwWAAALS2tuLw4cMwGAwoKChAW1ubXPGJiGgIshVIe3s7tFqtONbpdLDZbHecr9VqPeb39/dj586dWL16tThNo9HghRdegNVqxZw5c8R3JkRENPJk24XldrsREBAgjgVB8BgPNf/QoUOYMGECJk+eLE4rKSkRf16yZAnKysrgcDig0Xh3XENERKik+zJcPM6C5OTt82s0PA/9PaO/5wOUzShbgURFReH48ePi2G63Q6fTecy32+3i+PLlyx7zP/vsM6SkpIhjt9uNd999F0uXLkVQUJA4/ac/D6WjowtutzDs+zIcWq0GdrtD1tu40+3S/cGb55dSz8Ph8PeM/p4PkD9jYGDAoH94y7YLKyEhAc3Nzejs7ERPTw8aGxuRmJgozo+OjkZISAhaWloAABaLxWP+iRMnEBsb+/9BAwNx8OBBNDQ0AADMZjNmzpwJtVot110gIqJByFYgkZGRMJlMyM7OxjPPPIPU1FTExMQgPz8fra2tAICtW7di48aNWLBgAa5fv47s7Gxx/QsXLohf171l8+bN+PDDD7Fw4ULU1NSgtLRUrvhERDQEWY8DMRgM4jeobqmqqhJ/njJlCvbv33/bdU+ePDlg2qRJk/Dxxx/7NiQREUnCI9GJiEgSFggREUnCAiEiIkl4OneiUaTP2e/T40B6b7jguNZzt7HoPsUCIRpFglVBPr/ui38f6UD+jLuwiIhIEhYIERFJwgIhIiJJWCBERCQJC4SIiCRhgRARkSQsECIikoQFQkREkrBAiIhIEhYIERFJwgIhIiJJWCBERCQJC4SIiCRhgRARkSSyFojVakVKSgqSkpKwe/fuAfPPnDmDtLQ06PV6rFmzBi6XCwBQW1uL2bNnw2g0wmg0ory8HABw6dIlZGVlYcGCBVi+fDm6u7vljE9ERIOQrUBsNhvKy8uxZ88emM1m7Nu3D+fOnfNYprCwEOvWrUNDQwMEQUB1dTUA4PTp0yguLobFYoHFYoHJZAIAbNiwAZmZmaivr8f06dNRUVEhV3wiIhqCbAXS1NSE+Ph4hIWFQa1WQ6/Xo76+Xpx/8eJF9Pb2YtasWQCAtLQ0cX5raytqa2thMBjw8ssv4+rVq3A6nTh27Bj0ev2A5YmIaOTJViDt7e3QarXiWKfTwWaz3XG+VqsV52u1Wvz+97/HX//6V4wfPx4lJSW4cuUKQkNDMWbMmAHLExHRyJPtkrZutxsBAQHiWBAEj/Fg83fs2CFOz8vLw/z581FUVOSxPIAB46FERIQOa3mpvL1mNZE/UPL56u+vFX/PByibUbYCiYqKwvHjx8Wx3W6HTqfzmG+328Xx5cuXodPp4HA4UFNTgxdffBHAzWIJCgpCeHg4HA4H+vv7ERQUNGB73ujo6ILbLdzdHRuCVquB3T7yV5keDU908k9KPF8B5V4r3vL3fID8GQMDAwb9w1u2XVgJCQlobm5GZ2cnenp60NjYiMTERHF+dHQ0QkJC0NLSAgCwWCxITEyEWq3Ge++9h5MnTwIAdu3ahfnz50OlUiE2NhZ1dXUAALPZ7LE9IiIaWbK9A4mMjITJZEJ2djacTifS09MRExOD/Px8rFy5EjNmzMDWrVuxdu1adHV1Ydq0acjOzkZQUBC2b9+O119/Hb29vZgwYQK2bNkCAFi/fj2Ki4tRWVmJ8ePHY9u2bXLFJyKiIchWIABgMBhgMBg8plVVVYk/T5kyBfv37x+wXmxsLGprawdMj46OxkcffeT7oERENGw8Ep2IiCRhgRARkSQsECIikoQFQkREkrBAiIhIEhYIERFJwgIhIiJJWCBERCQJC4SIiCRhgRARkSQsECIikoQFQkREkrBAiIhIEhYIERFJIuvp3InIv/U5+316NcveGy44rvX4bHvk31ggRPexYFUQDKstPtuetcwI/74ILPkSC0QhmgfGYWwIH34iGr34P5hCxoaM8flffkREI4kfohMRkSSyFojVakVKSgqSkpKwe/fuAfPPnDmDtLQ06PV6rFmzBi6XCwDQ0tKC9PR0GI1G5OTk4OLFiwCAo0ePIi4uDkajEUajEa+++qqc8YmIaBCyFYjNZkN5eTn27NkDs9mMffv24dy5cx7LFBYWYt26dWhoaIAgCKiurhanl5aWwmKxwGAwoLS0FABw+vRp5ObmwmKxwGKxYOPGjXLFJyKiIchWIE1NTYiPj0dYWBjUajX0ej3q6+vF+RcvXkRvby9mzZoFAEhLS0N9fT36+vqwatUqTJkyBQAwefJktLW1AQBaW1tx+PBhGAwGFBQUiNOJiGjkyfYhent7O7RarTjW6XQ4derUHedrtVrYbDYEBwfDaLz5gbDb7cY777yDefPmAQA0Gg2Sk5ORlJSEvXv3wmQy4eOPP/Y6U0RE6N3eLa/48nv1RKPNcJ7//v5a8fd8gLIZZSsQt9uNgIAAcSwIgsd4qPl9fX0oLi6Gy+XCsmXLAAAlJSXi/CVLlqCsrAwOhwMajXcPYEdHF9xuQfJ98oZWq4HdPvQ34UfDE5NICm+e/4D3rxWl+Hs+QP6MgYEBg/7hLdsurKioKNjtdnFst9uh0+nuOP/y5cvi/O7ubuTl5cHlcqGyshIqlQputxuVlZXo7+/3uJ2goCC57gIREQ1CtgJJSEhAc3MzOjs70dPTg8bGRiQmJorzo6OjERISgpaWFgCAxWIR5xcWFuKRRx7B9u3bERwcfDNoYCAOHjyIhoYGAIDZbMbMmTOhVqvlugtERDQI2XZhRUZGwmQyITs7G06nE+np6YiJiUF+fj5WrlyJGTNmYOvWrVi7di26urowbdo0ZGdn4+uvv8bnn3+OiRMnYtGiRQBufn5SVVWFzZs347XXXsOOHTsQHh6OLVu2yBWfiIiGIOuR6AaDAQaDwWNaVVWV+POUKVOwf/9+j/lTp07Ft99+e9vtTZo0aVgfmhMRkXx4JDoREUnCAiEiIklYIEREJAkLhIiIJOHp3InIZ4Z7hUNvluVVDv0XC4SIfMbXVzgEeJVDf8ZdWEREJAkLhIiIJGGBEBGRJCwQIiKShAVCRESSsECIiEgSFggREUnCAiEiIklYIEREJAmPRPeS5oFxGBvi3cPF650T+c5wT48yFJ4axXdYIF4aGzLGp6dosJYZfbYtonuZr0+PwlOj+A53YRERkSQsECIikkTWArFarUhJSUFSUhJ27949YP6ZM2eQlpYGvV6PNWvWwOVyAQAuXbqErKwsLFiwAMuXL0d3dzcA4Nq1a1i6dCmSk5ORlZUFu90uZ3wiIhqEbAVis9lQXl6OPXv2wGw2Y9++fTh37pzHMoWFhVi3bh0aGhogCAKqq6sBABs2bEBmZibq6+sxffp0VFRUAAC2b9+O2NhYfPrpp1i8eDHefPNNueIT0T3q1ofyQ/0D4NVymgfGKXyPlCPbh+hNTU2Ij49HWFgYAECv16O+vh5//OMfAQAXL15Eb28vZs2aBQBIS0vD22+/jcWLF+PYsWPYsWOHOP35559HYWEhvvzyS/GdTGpqKkpKSuB0OqFSqbzKFBgYcFf3Sfegb58o99v25Nimv29Pjm3eb9vz9TaDVUH4XWmjz7ZX+cpcn3/z8sYNF7q6er1a9m7/X7urbQsy2blzp7Bt2zZxXF1dLaxdu1Yc/+tf/xIyMjLE8fnz54WkpCTBZrMJTzzxhDjd6XQK06ZNEwRBEKZNmyY4nU5x3hNPPCH88MMPct0FIiIahGy7sNxuNwIC/r+9BEHwGN9p/s+XAzBg/NN1AgP5PQAiIiXI9r9vVFSUx4fcdrsdOp3ujvMvX74MnU6H8PBwOBwO9Pf3D1hPp9Ph8uXLAACXy4Xu7m5xFxkREY0s2QokISEBzc3N6OzsRE9PDxobG5GYmCjOj46ORkhICFpaWgAAFosFiYmJUKlUiI2NRV1dHQDAbDaL682ZMwdmsxkAUFdXh9jYWK8//yAiIt8KEARBkGvjVqsV7777LpxOJ9LT05Gfn4/8/HysXLkSM2bMwDfffIO1a9eiq6sL06ZNw8aNGxEcHIyLFy+iuLgYHR0dGD9+PLZt24Zf/OIX+PHHH1FcXIwLFy5Ao9Fg69at+OUvfylXfCIiGoSsBUJERPcufgJNRESSsECIiEgSFggREUnCAiEiIklYID4y1Ikj/cE777yDhQsXYuHChdiyZYvSce5o8+bNKC4uVjrGbX3xxRdIS0tDcnIySktLlY4zgMViEX/HmzdvVjqOh66uLqSmpuK///0vgJunOzIYDEhKSkJ5ebnC6Qbm27dvH1JTU2EwGPDqq6+ir69P4YQDM96ya9cuvPDCCyMfSLmD4O8dP/zwg/DUU08JV65cEbq7uwWDwSCcPXtW6Vgejhw5Ijz33HPCjRs3hL6+PiE7O1tobGxUOtYATU1NQlxcnPDKK68oHWWA77//Xpg9e7bQ1tYm9PX1CUuWLBG+/PJLpWOJrl+/Ljz22GNCR0eH4HQ6hfT0dOHIkSNKxxIEQRBOnDghpKamCtOmTRMuXLgg9PT0CHPmzBG+//57wel0Crm5uYo+lj/P95///EeYP3++4HA4BLfbLRQVFQkffPCBYvlul/GWs2fPCk888YTw/PPPj3gmvgPxgZ+eOFKtVosnjvQnWq0WxcXFCA4Ohkqlwq9+9StcunRJ6VgefvzxR5SXl6OgoEDpKLd18OBBpKSkICoqCiqVCuXl5Zg5c6bSsUT9/f1wu93o6emBy+WCy+VCSEiI0rEAANXV1Vi/fr14VolTp07hkUcewcMPP4wxY8bAYDAo+pr5eb7g4GCsX78eoaGhCAgIwK9//WvFXy8/zwgAfX19WLduHVauXKlIJl7S1gfa29uh1WrFsU6nw6lTpxRMNNCkSZPEn8+fP49PP/0Ue/fuVTDRQOvWrYPJZEJbW5vSUW7ru+++g0qlQkFBAdra2vDkk0/ipZdeUjqWKDQ0FKtWrUJycjLGjRuHxx57DI8++qjSsQBgwKUXbveasdlsIx1L9PN80dHRiI6OBgB0dnZi9+7d2LhxoxLRRLe7fEVZWRmeffZZxQ6o5jsQHxjqxJH+5OzZs8jNzUVRUREmTJigdBzRJ598gvHjx+Pxxx9XOsod9ff3o7m5GW+99Rb27duHU6dOoba2VulYom+++QY1NTX4+9//jkOHDiEwMBDvv/++0rFua7S8Zmw2G3JycvDss88iLi5O6Tgejhw5gra2Njz77LOKZWCB+MBQJ470Fy0tLXjxxRexevVqLFq0SOk4Hurq6nDkyBEYjUa8/fbb+OKLL/DWW28pHcvDQw89hMcffxzh4eEYO3Ys5s2b51fvNA8fPozHH38cERERCA4ORlpaGo4ePap0rNsaDa+Zf//738jIyMCiRYvwhz/8Qek4A/ztb3/D2bNnYTQasXbtWpw+fXrk3xGP+Kcu96BbH6J3dHQI169fF55++mnh5MmTSsfycOnSJSEuLk5oampSOsqQampq/PJD9BMnTgh6vV64evWq4HK5hGXLlgnV1dVKxxIdOnRIePrpp4Xu7m7B7XYLr732mvD2228rHcvDU089JVy4cEHo7e0VEhMThfPnzwsul0v43e9+J9TV1SkdT8zncDiEOXPmCLW1tUpHGuBWxp/65z//qciH6PwMxAciIyNhMpmQnZ0tnjgyJiZG6Vge3n//fdy4cQObNm0Sp2VkZGDJkiUKphpdZs6ciby8PGRmZsLpdOI3v/mNorsPfm727Nn4+um6XS8AAAFcSURBVOuvkZaWBpVKhRkzZmDp0qVKx7qtkJAQbNq0CStWrMCNGzcwZ84cLFiwQOlYov379+Py5cv44IMP8MEHHwAAfvvb32LVqlUKJ/MvPJkiERFJws9AiIhIEhYIERFJwgIhIiJJWCBERCQJC4SIiCRhgRARkSQsECIikoQFQqSQ2tpazJs3D93d3bh+/TqSk5NhNpuVjkXkNR5ISKSg1atXQ6PRoK+vD0FBQXjjjTeUjkTkNRYIkYK6urpgNBoxduxYHDhwwG+u30HkDe7CIlJQR0cHbty4gWvXrqG9vV3pOETDwncgRApxOp3IyMhARkYG3G43PvnkE+zduxcqlUrpaERe4TsQIoVs27YNDz30EBYvXoznnnsODz74IMrLy5WOReQ1vgMhIiJJ+A6EiIgkYYEQEZEkLBAiIpKEBUJERJKwQIiISBIWCBERScICISIiSVggREQkyf8CYY53nwys7V0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute bin edges: bins\n", "bins = np.arange(0, max(n_defaults) + 2) - 0.5\n", "\n", "# Generate histogram\n", "_ = plt.hist(n_defaults, bins=bins, density=True)\n", "\n", "# Label axes\n", "_ = plt.xlabel('x')\n", "_ = plt.ylabel('y')\n", "plt.savefig('../images/bin-dist.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Poisson processes and the Poisson distribution\n", "- Poisson process\n", " - The timing of the next event is completely independent of when the previous event happened\n", "- Example\n", " - Natural births in a given hospital\n", " - Hit on a website during a given hour\n", " - Meteor strikes\n", " - Molecular collisions in a gas\n", " - Aviation incidents\n", " - Buses in Poissonville\n", "- Poisson Distribution\n", " - The number r of arrivals of a Poisson process in a given time interval with average rate of? arrivals per interval is Poisson distributed.\n", " - The number r of hits on a website in one hour with an average hit rate of 6 hits per hour is Poisson distributed.\n", " - Limit of the Binomial distribution for low probability of success and large number of trials\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Relationship between Binomial and Poisson distributions\n", "You just heard that the Poisson distribution is a limit of the Binomial distribution for rare events. This makes sense if you think about the stories. Say we do a Bernoulli trial every minute for an hour, each with a success probability of 0.1. We would do 60 trials, and the number of successes is Binomially distributed, and we would expect to get about 6 successes. This is just like the Poisson story we discussed in the video, where we get on average 6 hits on a website per hour. So, the Poisson distribution with arrival rate equal to np approximates a Binomial distribution for n Bernoulli trials with probability p of success (with n large and p small). Importantly, the Poisson distribution is often simpler to work with because it has only one parameter instead of two for the Binomial distribution.\n", "\n", "Let's explore these two distributions computationally. You will compute the mean and standard deviation of samples from a Poisson distribution with an arrival rate of 10. Then, you will compute the mean and standard deviation of samples from a Binomial distribution with parameters n and p such that np=10." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Poisson: 9.9895 3.176159591393354\n", "n = 20 Binom: 10.0621 2.240277569856021\n", "n = 100 Binom: 10.0432 2.988567175085747\n", "n = 1000 Binom: 10.0106 3.1406508306400442\n" ] } ], "source": [ "# Draw 10,000 samples out of Poisson distribution: samples_poisson\n", "samples_poisson = np.random.poisson(10, 10000)\n", "\n", "# Print the mean and standard deviation\n", "print('Poisson: ', np.mean(samples_poisson), np.std(samples_poisson))\n", "\n", "# Specify values of n and p to consider for Binomial: n, p\n", "n = [20, 100, 1000]\n", "p = [0.5, 0.1, 0.01]\n", "\n", "# Draw 10,000 samples for each n,p pair: samples_binomial\n", "for i in range(3):\n", " samples_binomial = np.random.binomial(n[i], p[i], size=10000)\n", " \n", " # Print results\n", " print('n = ', n[i], 'Binom:', np.mean(samples_binomial), np.std(samples_binomial))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Was 2015 anomalous?\n", "1990 and 2015 featured the most no-hitters of any season of baseball (there were seven). Given that there are on average 251/115 no-hitters per season, what is the probability of having seven or more in a season?" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Probability of seven or more no-hitters: 0.0071\n" ] } ], "source": [ "# Draw 10,000 samples out of Poisson distribution: n_nohitters\n", "n_nohitters = np.random.poisson(251/115, 10000)\n", "\n", "# Compute number of samples that are seven or greater: n_large\n", "n_large = np.sum(n_nohitters >= 7)\n", "\n", "# compute probability of getting seven or more: p_large\n", "p_large = n_large / 10000\n", "\n", "# Print the result\n", "print('Probability of seven or more no-hitters:', p_large)" ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }