{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Class 1: Probability & Statistics\n", "\n", "In this class we will review how statistics are used to summarize data, special probability distributions, their use in simple applications using Frequentist and Bayesian methods, and Monte Carlo techniques. At the end of this class you should be able to:\n", "\n", "- determine summary statistics for datasets and their errors\n", "- optimally combine data\n", "- apply probability distributions for Gaussian, Binomial and Poisson statistics\n", "- compare the Frequentist and Bayesian frameworks for statistical analysis\n", "- solve statistical problems using Monte Carlo techniques\n", "\n", "## Summary statistics and their errors\n", "\n", "A **statistic** is a quantity which summarizes our data. Example statistics for a sample of $N$ independent estimates of a quantity are the **mean** (average value), **median** (middle value) and **standard deviation** $\\sigma$ (spread).\n", "\n", "- Mean $= \\overline{x} = \\frac{1}{N} \\sum_i x_i$\n", "- Median $=$ middle value when ranked\n", "- Variance $= \\sigma^2 = \\frac{1}{N-1} \\sum_i \\left( x_i - \\overline{x} \\right)^2$\n", "\n", "Here is an example of how to calculate these in python using the scipy.stats library:\n", "\n", "_We have 10 measurements of a variable, $x_i = (7.6, 5.8, 8.0, 6.9, 7.2, 7.5, 6.4, 8.1, 6.3, 7.0)$. Estimate the mean, variance and median of this dataset. What are the errors in your estimates?_" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mean = 7.08\n", "Variance = 0.5662222222222222\n", "Standard deviation = 0.7524773898411979\n", "Median = 7.1\n" ] } ], "source": [ "import numpy as np\n", "from scipy import stats\n", "x = [7.6, 5.8, 8.0, 6.9, 7.2, 7.5, 6.4, 8.1, 6.3, 7.0] # 10 measurements of a variable\n", "s = stats.describe(x)\n", "mu = s.mean # mean\n", "var = s.variance # variance\n", "sig = np.sqrt(var) # standard deviation\n", "median = np.median(x) # median\n", "print('Mean =',mu)\n", "print('Variance =',var)\n", "print('Standard deviation =',sig)\n", "print('Median =',median)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can quote an error in each of these statistics:\n", "\n", "- Error in the mean = $\\sigma/\\sqrt{N}$\n", "- Error in the variance = $\\sigma^2 \\sqrt{2/(N-1)}$\n", "- Error in the median = $1.25 \\sigma/\\sqrt{N}$\n", "\n", "The error in the mean holds for all distributions, the other relations assume a Gaussian distribution." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error in mean = 0.23795424396766326\n", "Error in variance = 0.2669197153278997\n", "Error in median = 0.2974428049595791\n" ] } ], "source": [ "n = len(x)\n", "print('Error in mean =',sig/np.sqrt(n))\n", "print('Error in variance =',var*np.sqrt(2./(n-1)))\n", "print('Error in median =',1.25*sig/np.sqrt(n))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Optimal combination of data\n", "\n", "Suppose we have $N$ independent estimates $x_i$ of some quantity $y$, which have varying errors $\\sigma_i$. What is our best combined estimate of $y$?\n", "\n", "A simple average $\\hat{y} = \\sum_i x_i/N$ is not the optimal combination, because we want to give **more weight to the more precise estimates**. Let's weight each estimate by $w_i$: $\\hat{y} = \\sum_i w_i x_i / \\sum_i w_i$.\n", "\n", "The weights which minimize the combined error are the **inverse-variance weights**, $w_i = 1/\\sigma_i^2$. In this case, the variance in the combined estimate is given by, $1/{\\rm Var}(\\hat{y}) = \\sum_i 1/\\sigma_i^2$.\n", "\n", "Note that this approach is only useful if the errors in the data are dominated by statistical, not systematic, errors.\n", "\n", "Here is an example:\n", "\n", "_We have 5 measurements of a quantity: $(7.4 \\pm 2.0, 6.5 \\pm 1.1, 4.3 \\pm 1.7, 5.5 \\pm 0.8, 6.0 \\pm 2.5)$. What is the optimal estimate of this quantity, and the error in that estimate?_" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Weights = [0.25 0.82644628 0.34602076 1.5625 0.16 ]\n", "Combined estimate = 5.807227819726063\n", "Error in combined estimate = 0.5638868290446053\n" ] } ], "source": [ "x = np.array([7.4, 6.5, 4.3, 5.5, 6.0]) # 5 measurements of a quantity\n", "sig = np.array([2.0, 1.1, 1.7, 0.8, 2.5]) # errors in these measurements\n", "wei = 1./(sig**2) # weights\n", "print('Weights =',wei)\n", "print('Combined estimate =',np.sum(wei*x)/np.sum(wei))\n", "print('Error in combined estimate =',1./np.sqrt(np.sum(wei)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probability distributions\n", "\n", "A **probability distribution**, $P(x)$, is a function which assigns a probability for each particular value (or range of values) of a discrete or continuous variable $x$.\n", "\n", "- A probability distribution must be normalized, $\\int_{-\\infty}^\\infty P(x) \\, dx = 1$ (or $\\sum_i P(x_i) = 1$ for a discrete variable)\n", "- Probability in the range $[x_1, x_2]$ is $\\int_{x_1}^{x_2} P(x) \\, dx$\n", "\n", "A probability distribution may be quantified by its...\n", "\n", "- Mean, $\\mu = \\overline{x} = \\langle x \\rangle = \\int_{-\\infty}^\\infty x \\, P(x) \\, dx$\n", "- Variance $\\sigma^2 = \\int_{-\\infty}^\\infty (x - \\mu)^2 \\, P(x) \\, dx = \\langle x^2 \\rangle - \\langle x \\rangle^2$\n", "\n", "Certain important types of variables have known probability distributions:\n", "\n", "- **Binomial** distribution\n", "- **Poisson** distribution\n", "- **Gaussian** or **Normal** distribution\n", "\n", "## Binomial distribution\n", "\n", "The binomial distribution applies where there is a random process with **two possible outcomes** with probabilities $p$ and $1-p$ (for example, tossing a coin).\n", "\n", "If we have $N$ trials, and the probability of \"success\" in each trial is $p$, then the probability of $n$ successes is:\n", "\n", "$P_{Binomial}(n) = \\frac{N!}{n! (N-n)!} p^n (1-p)^{N-n}$\n", "\n", "The mean and variance of this distribution are $\\overline{n} = p N$, ${\\rm Var}(n) = N p (1-p)$.\n", "\n", "Here is an example Binomial distribution for $N = 10$ and $p = 0.2$ (you can change these choices in the notebook of course!):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 0.4)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZyNdf7H8dfHjJ2YiXXXHSp+LItEaFYq3dhlWlKSoVWhomLRpqKU3OS+tjs3KUOrm5lCucmYH1psKzKaqZDJmPULm5WbpWEwxuf3xzmjYxrm7pxznbnO5/l4zMOc61zXud7mUfN2Xd/r+l6iqhhjjDGlVcHpAMYYY8o3KxJjjDFlYkVijDGmTKxIjDHGlIkViTHGmDKxIjHGGFMmQS0SEekkIhkikikiw8+z3l0ioiLS2mfZCO92GSLSMTiJjTHGFCUyWDsSkQhgGvB7YDewUUQWq+rWAutdCAwBNvgsawL0BJoClwErReQ3qpoXrPzGGGMKF8wjkmuBTFXNUtWTQCLQtZD1xgKTgOM+y7oCiap6QlX/BWR6P88YY4zDgnZEAtQGdvm83g3E+q4gItcAdVX1ExF5osC26wtsW7vgDkSkP9AfIDo6ulXjxo39FN0YY8LDpk2b9qtqrZJsE8wiOS8RqQC8BPQp7Weo6ixgFkDr1q01NTXVP+HKsV27PN1dt25dh5MYY8oDEfm/km4TzCLZA/j+NqvjXZbvQqAZsFpEAC4BFovI7cXY1pzDvffeC8Dq1audDWKMca1gFslGoKGI1MNTAj2Be/LfVNXDQM381yKyGhimqqkikgO8JyIv4Rlsbwh8EcTs5dbIkSOdjmCMcbmgFYmqnhKRQUAKEAEkqOoWERkDpKrq4vNsu0VEPgC2AqeAgXbFVvF06NDB6QjGGJcTt04jb2MkHllZWQDUr1/f4STGmPJARDapauui1/xZyAy2m8Do168fYGMkxpjAsSJxudGjRzsdwRjjclYkLte+fXunIxhjXM4mbXS5jIwMMjIynI5hjHExOyJxuQEDBgA2RmKMCRwrEpcbP3680xGMMS5nReJy1113ndMRjDEuZ2MkLrd582Y2b97sdAxjjIvZEYnLDRo0CLAxEmNM4FiRuNyUKVOcjmCMcTkrEpdr06aN0xGMMS5nYyQul56eTnp6utMxjDEuZkckLjd06FDAxkiMMYFjReJyL7/8stMRjDEuZ0Xici1atHA6gjHG5WyMxOU2btzIxo0bnY5hjHExOyJxuSeeeAKwMRJjTOBYkbjc66+/7nQEY4zLBfXUloh0EpEMEckUkeGFvP+wiHwjIuki8pmINPEuv1JEcrzL00VkZjBzl2fNmjWjWbNmTscwxrhY0I5IRCQCmAb8HtgNbBSRxaq61We191R1pnf924GXgE7e93aoqo0cl9C6desAm7zRGBM4wTy1dS2QqapZACKSCHQFzhSJqh7xWT8a0CDmc6Wnn34asDESY0zgBLNIagO7fF7vBmILriQiA4G/AL8CbvF5q56IpAFHgJGq+o8AZnWNN954w+kIxhiXC7nBdlWdBkwTkXuAkcD9wA/A5ap6QERaAR+LSNMCRzCISH+gP8Dll18e5OShqVGjRk5HMMa4XDAH2/cAdX1e1/EuO5dE4A4AVT2hqge8328CdgC/KbiBqs5S1daq2rpWrVp+C16erVmzhjVr1jgdwxjjYsE8ItkINBSRengKpCdwj+8KItJQVbd7X/4R2O5dXgs4qKp5IlIfaAhkBS15OTZq1CjAxkiMMYETtCJR1VMiMghIASKABFXdIiJjgFRVXQwMEpEOQC5wCM9pLYAbgTEikgucBh5W1YPByl6eJSQkOB3BGONyourOC6Nat26tqampTscwxphyRUQ2qWrrkmxjc2253MqVK1m5cqXTMYwxLhZyV20Z/xo3bhwAHTp0cDiJMcatrEhcbt68eU5HMMa4nBWJy9WtW7folYwxpgxsjMTlli9fzvLly52OYYxxMTsicbmJEycC0KlTpyLWNMaY0rEicbnExESnIxhjXM6KxOUuueQSpyMYY1zOxkhcbsmSJSxZssTpGMYYF7MjEpd78cUXAejSpYvDSYwxbmVF4nLz5893OoIxxuWsSFyuZs2aTkcwxricjZG43MKFC1m4cKHTMYwxLmZHJC736quvAtCtWzeHkxhj3MqKxOUWLVrkdARjjMtZkbhc1apVnY5gjHE5GyNxuaSkJJKSkpyOYYxxMTsicbkZM2YAEB8f73ASY4xb2RGJi+Xl5TFo0CBuuOEGli5dSl5entORjDEuFNQjEhHpBLwCRABvqerEAu8/DAwE8oBsoL+qbvW+NwJ4wPveYFVNCWb28iYvL4+OHTuyYcMGjh49SnR0NLGxsaSkpBAREeF0PGOMiwTtiEREIoBpQBzQBOglIk0KrPaeql6lqi2AycBL3m2bAD2BpkAnYLr388w5LF26lM8++4zs7GxUlezsbDZs2EBycrLT0YwxLhPMU1vXApmqmqWqJ4FEoKvvCqp6xOdlNKDe77sCiap6QlX/BWR6P88UoKosWrSIBx98kBMnTpz13tGjR0lPT3comTHGrYJZJLWBXT6vd3uXnUVEBorIDjxHJINLuG1/EUkVkdQff/zRb8HLi/Xr13PjjTdyxx13EBUVxQUXXHDW+5UqVaJFixYOpTPGuFXIDbar6jRV/R/gKWBkCbedpaqtVbV1rVq1AhMwBH333Xd0796dtm3bkpmZycyZM9mxYwft2rUjJiYGEQGgSpUqxMXFOZzWGOM2wSySPUBdn9d1vMvOJRG4o5TbhoX//Oc/DBw4kKZNm5KSksLo0aPZvn07AwYMICoqipSUFB544AHuuOMO4uPj2bt3L2vXrnU6tjHGZYJZJBuBhiJST0R+hWfwfLHvCiLS0OflH4Ht3u8XAz1FJEpE6gENgS+CkDkkZWdnM2bMGBo0aMCsWbPo378/mZmZPPfcc8TExJxZLyIigvT0dA4ePMicOXOoX78+jz76KCdPnnQwvTHGdVQ1aF/AbcB3wA7gGe+yMcDt3u9fAbYA6cDfgaY+2z7j3S4DiCtqX61atVK3yc3N1ZkzZ+rFF1+sgN51112akZFR7O2XLVumgI4fPz6AKY0x5RmQqiX83S6e7dyndevWmpqa6nQMv1DvlVjDhw8nIyOD66+/nsmTJ9O2bdsSf1b37t1ZtmwZW7ZsoV69egFIa4wpz0Rkk6q2Lsk2ITfYbs62bt06rr/+eu68805EhEWLFrF27dpil8ibb77Jm2++eeb1X//6VypUqMDgwYNx6z8ijDHBZUUSojIyMujWrRvt2rUjKyuLWbNm8c0333D77befuQqrOApO2li3bl1Gjx7N0qVLWbx48Xm2NMaY4rFTWyFm7969jB49mjfffJNKlSrx1FNP8dhjjxEdHe23feTm5nLNNddw5MgRtm7d6tfPNsaUb3Zqqxz76aefeP7552nQoAFvvfUWjzzyCDt27GDkyJF+/0VfsWJFZs6cyffff8+YMWP8+tnGmPBjReKw3NxcZsyYQYMGDRg9ejS33XYbW7du5bXXXuOiiy4q8+dPnz6d6dOn/2J5u3bt6NevHy+99BJbtmwp836MMeHLisQhqsrChQtp1qwZjz76KI0bN2b9+vV88MEHNGzYsOgPKKYlS5awZMmSQt+bNGkSVapU4dFHH7WBd2NMqVmROOCzzz6jXbt23HXXXURERLB48WJWr15NbGys3/eVnJx8zhl/a9asyaRJk1i7di3z5s3z+76NMeHBiiSIvv32W+644w5uuOEGdu7cyZtvvsnXX39Nly5dSnQllj/169ePtm3bMmzYMA4ePOhIBmNM+WaP2g2QvLw8kpOTSUtL44orruAf//gHCQkJREdHM27cOIYOHRqUq6VeeeUVAIYMGVLo+xUqVGDGjBm0atWKp59+mpkzZwY8kzHGXaxIAiD/6YTr16/n6NGjZ5YPHDiQUaNGEcyZiVetWgWcu0gArr76agYPHszLL79M3759A3KKzRjjXnYfSQAsXbqU+Ph4jh07dmZZ5cqVSUpKonPnzo5kKspPP/1E48aNufjii/niiy+IjLR/YxgTjuw+khCRlpZ2VokA5OTkhPTTCS+88EJefvll0tLSCr1c2BhjzsWKJACaN2/+i8Hz6OhoR55OOHXqVKZOnVqsdbt3707Hjh0ZOXIkP/zwQ4CTGWPcwookACIjI1FVoqKiEBFiYmKIjY115OmEn3/+OZ9//nmx1hURXn/9dU6ePMlf/vKXACczxriFnQgPgL/97W9Uq1aN2bNns2XLFlq0aEFcXBwRERFBz7JgwYISrd+gQQNGjBjB888/T79+/fj9738foGTGGLewwXY/O3DgAJdddhkPP/zwmUtvy5vjx49z1VVXISJ88803REVFOR3JGBMkNtgeAt577z1OnjxJv379nI4CwMSJE5k4cWKJtrnggguYNm0a27dvZ/LkyQFKZoxxCysSP0tISOCaa67h6quvdjoKAOnp6aW6WuwPf/gDPXr04IUXXmDHjh0BSGaMcYugFYmIdBKRDBHJFJHhhbz/FxHZKiJfi8gqEbnC5708EUn3foXs05jS0tJIT08PmaMRgMTERBITE0u17UsvvUTFihX585//bJM6GmPOKShFIiIRwDQgDmgC9BKRJgVWSwNaq2pzYD7ge04lR1VbeL9uD0bm0khISCAqKopevXo5HcUvateuzdixY0lOTmbhwoVOxzHGhKhgHZFcC2SqapaqngQSga6+K6jq31U1/y6+9UCdIGXzi+PHj/Puu+9y5513Ur16dafjnDF27FjGjh1b6u0HDRrE1VdfzZAhQ/jpp5/8mMwY4xbBKpLawC6f17u9y87lAcB37vMLRCRVRNaLyB3n2khE+nvXS/3xxx/LlriEFi1axKFDh0LqtBZ4nv2ekZFR6u0jIyOZMWMGe/bsYfTo0X5MZoxxi6Bc/isi3YFOqvqg9/W9QKyqDipk3d7AIKC9qp7wLqutqntEpD7wKXCrqp53BDjYl/927NiRb7/9ln/961+O3C8SaP379ychIYEvv/yS5s2bOx3HGBMgoXz57x6grs/rOt5lZxGRDsAzwO35JQKgqnu8f2YBq4GWgQxbUrt27WLFihX06dPHlSUCMGHCBKpVq8YjjzzC6dOnnY5jjAkhJS4SEYn2Dp6XxEagoYjUE5FfAT2Bs66+EpGWwBt4SmSfz/JqIhLl/b4m0A7YWtLcgfT222+jqvTp08fpKL/w3HPP8dxzz5X5c2rUqMHkyZNZt24dc+fOLXswY4xrFFkkIlJBRO4RkU9EZB+wDfjBe6nuFBFpUNRnqOopPKerUoBvgQ9UdYuIjBGR/KuwpgAxwIcFLvP9LZAqIl8BfwcmqmrIFMnp06eZM2cON998M/Xr13c6zi/s2rWLXbt2Fb1iMdx///1cf/31PPnkkxw4cMAvn2mMKf+KHCMRkTXASmARsFlVT3uXVwduBu4BPlLVdwKctUSCNUayevVqbr75ZubNm0fv3r0Dvj+nffPNN7Rs2ZK+ffvy5ptvOh3HGONngRoj6aCqY1X16/wSAVDVg6q6QFXvApJKGtYtEhISqFKlCt26dXM6SlBcddVVPPbYY7z11lusW7fO6TjGmBBQZJGoaq4/1nGjw4cPM3/+fHr16kXlypWdjlOoESNGMGLECL9+5qhRo6hTpw6PPPIIp06d8utnG2PKn2IPtovILSIyW0ReFJG+ItIqfxA8XCUlJZGTkxNy9474OnDggN/HM2JiYnjllVf4+uuvee211/z62caY8qfY95GIyE5gKFARaO79aqqqRQ62OyEYYyS/+93vyM7O5ptvvvnFExHdTlXp3Lkza9eu5dtvv6VOnXI1EYEx5hwCfR/J/6nqx6r6oao+q6pdQ7VEgmHr1q1s2LCBfv36hV2JgOdpiq+99hqnTp2ypykaE+ZKUiRrReQxCcffmoWYM2cOkZGRIX+l1rBhwxg2bFhAPrt+/fo888wzfPjhh6SkpARkH8aY0FeSImkCPILnHpJPROQFEbk7QLlCWm5uLn/729/o3LkzF110kdNxzisnJ4ecnJyAff4TTzzBb37zGwYOHBjQ/RhjQlexi0RV71LV3wD1gOeA7cDvAhUslC1btox9+/aF9CB7vmnTpjFt2rSAfX5UVBTTp09nx44dTJo0KWD7McaEruLckChaxErFWSfYAjnY3rVrV7744gt27dpFZGRkQPZR3txzzz0sWLCAzZs307BhQ6fjGGNKKVCD7X8XkT+LyOUFdvYr7yXBbwP3l2Sn5dnevXv55JNPuO+++8pFiQwdOpShQ4cGfD8vvvgiF1xwAQMHDrSnKRoTZopTJJ2APOB9Efm3d46tf+E5tdULeFlV5wYwY0iZN28eeXl59O3b1+koIeXSSy9l3LhxrFixgg8//NDpOMaYICrR80hEpCJQE8+jb/8bsFR+EIhTW6pKkyZNqF69Ov/85z/9+tlukJeXR5s2bdi7dy/btm2jSpUqTkcyxpRQwJ9H4p0KpSVwnYjclv9Vks8oz9avX8+2bdvKxSC7EyIiIpg5cyZ79+5l1KhRTscxxgRJaR5sVcv7VdPnKywkJCRQuXJlevTo4XSUYhs4cCADBw4M2v6uvfZaBgwYwKuvvkp6enrQ9muMcU6Ji0RV3wYygFg808jf5OdMIeno0aMkJSVx9913c+GFFzodp9gqVapEpUqVgrrP8ePHU6NGDXuaojFhorSP2n0QOASMArL8Fyd0LViwgJ9++qncndaaOnUqU6dODeo+q1WrxtSpU1m/fj2zZ88O6r6NMcFX2iL5D3ABcBq4xH9xQldCQgINGjTghhtucDpKuXDvvffSvn17nnrqKX788Uen4xhjAqi0RfIuMBN4Es/TE10tMzOTNWvW0Ldv33I3QWP//v3p379/0PcrIkyfPp0jR47Qq1cvxo4dy9KlS8nLywt6FmNMYJX4jjoRaQdU9b5MAYp9/bCIdAJeASKAt1R1YoH3/4LntNkp4Eegn6r+n/e9+4GR3lXHecdqgmLu3LlUqFCB++67L1i79JsaNWo4tu9GjRpRp04dVq1axaeffkp0dDSxsbGkpKQQERHhWC5jjH+V5tbs6t6vfMUqEhGJAKYBvwd2AxtFZLGqbvVZLQ1orarHROQRYDIQ730+/CigtXd/m7zbHipF/hLJy8tj7ty5dOzYsVw+c2PChAmO7Ts5OZn9+/cDnntwsrOz2bBhA8nJyXTu3NmxXMYY/yrNVVtL8Jza+gr4FthWzE2vBTJVNUtVTwKJQNcCn/13VT3mfbkeyP/N3RFY4X1O/CFgBZ477gNuxYoV7Nmzp9wNsoeCtLQ0jh07dtayo0eP2mXBxrhMacdIPsDzi/wG4PpiblMb2OXzerd32bk8ACSXZFsR6S8iqSKS6q8B3oSEBGrUqEGXLl388nnB1rdvX8emc2nZsiXR0dFnLatcuTItWrRwJI8xJjBKO+vgloLjG/4kIr3xnMZqX5LtVHUWMAs8U6SUNceBAwdYtGgRjzzyCFFR5fPx9HXr1nVs33FxccTGxrJhwwaOHj2KqlKrVi3i4uIcy2SM8b/SFkmuiKzAMyCOqt5TjG32AL6/1ep4l51FRDoAzwDtVfWEz7Y3Fdh2dYlTl9B7773HyZMny/UEjWPGjHFs3xEREaSkpJCcnEx6ejopKSls2rSJffv2cemllzqWyxjjXyWatPHMRiLTVfXREm4TCXwH3IqnGDYC96jqFp91WgLzgU6qut1neXVgE3CNd9GXQCtVPXiu/flj0saWLVtSoUIFNm3aVKbPMR6ZmZn89re/pV+/frzxxhtOxzHGFCLgkzb6qCwiPUsyaaOqngIG4blk+FvgA1XdIiJjROR272pTgBjgQxFJF5HF3m0PAmPxlM9GYMz5SsQf0tLSSE9PL/eD7L179w6Z58o3aNCARx55hNmzZ7NtW3Gv0TDGhLrSntr6O/ArSjhho6ouA5YVWPacz/cdzrNtApBQspill5CQQFRUFL169QrWLgOiUaNGTkc4y7PPPsvcuXMZMWIEH330kdNxjDF+UNoiOQT0xnNPx3v+ixMajh8/zrvvvsudd95J9erVi94ghD377LNORzhLrVq1ePLJJ3n22Wf55z//Sbt27ZyOZIwpo9Ke2uqsqj1UNZ4g3c8RTIsWLeLQoUPl/rRWqHrssce49NJLeeKJJ+yxvMa4QGmLpJKIXO59jnt0kWuXMwkJCVx++eXccsstTkcps549e9KzZ0+nY5wlOjqa0aNH8/nnn/Pxxx87HccYU0alLZLngT97v8b6LU0I+P7771mxYgV9+vRxxXxQLVq0CMkbAPv27Uvjxo0ZMWIEp06dcjqOMaYMSjVGoqo7gCf8nCUkvP3226gqffr0cTqKXwwfPtzpCIWKjIxk4sSJ3HHHHcyePZsBAwY4HckYU0qlvY/kaaApniMaLeYNiUFVmvtITp8+TcOGDbnyyitZtWpVgJKZfKrKDTfcwI4dO9i+fTsxMTFORzIm7AXzPhJR1T+paq9QLJHSWrt2LVlZWeX6TvaC7rrrLu666y6nYxRKRJgyZQp79+7lpZdecjqOMaaUSlsk/yMid5fkhsTyICEhgSpVqtCtWzeno/hN27Ztadu2rdMxzqlt27Z069aNKVOmsG/fPqfjGGNKobRFsgaoDNSihDclhqrDhw8zf/58evXqReXKlZ2O4zfDhg1j2LBhTsc4r/Hjx5OTk+PovGDGmNIrcZGIyPV4Jmvc7/NV7iUlJZGTk2P3jjigUaNGPPTQQ7zxxhts37696A2MMSGlNEck1fAchdQA6gNX+jOQUxISEmjatClt2rRxOopf3X777dx+++1Fr+iwUaNGERUVxdNPP+10FGNMCZWmSG7F87jcVnhuRrzOr4kcsGXLFjZs2EC/fv0QEafj+NWtt97Krbfe6nSMIl1yySUMGzaM+fPns2HDBqfjGGNKoDRFEqmq9wKVvA+3CugsvMEwZ84cIiMjQ2aWXH8aMmQIQ4YMcTpGsTz++ONcdNFFPPnkkzZ1ijHlSGmK5ArvlVoXef+80r+Rgis3N5d58+bRpUsXLrroIqfjhLULL7yQUaNGsXbtWj755BOn4xhjiqk0RTIfz9VaH3n/XODXREG2bNky9u3b59pB9ri4uHL1aNuHHnqIhg0b8tRTT9nUKcaUE0VOkSIi9wMv4imdpcBAVf0p0MGCJSEhgUsuuYROnVw3iTEAXbp0cTpCiVSsWJEJEybQvXt33n77bR544AGnIxljilDkFCkikgncjefxuH8Gfq2qfw5CtjIpzhQpe/fupU6dOjz++ONMmjQpSMlMUVSV6667ju+//57t27e76r4eY0JdoKZIOaKqaaq6T1WfBa4tXbzQM2/ePPLy8lw1JYobiAiTJ0/m3//+N6+88orTcYwxRShOkVwqIv1F5EYRqQVULM2ORKSTiGSISKaI/GJKWu/nfykip0Ske4H38rzPcD/zHPeyUlUSEhK47rrraNy4sT8+MiR16NCBDh3O+QTjkHXDDTfQpUsXJk6cyP79rrjn1RjXKk6RjAKuwvPckQygmYgsE5EJIlKsB5qLSAQwDYgDmgC9RKRJgdW+B/pQ+KN7c1S1hffLL3fXrV+/nm3btrl2kD1ffHw88fHxTscolYkTJ5Kdnc24ceOcjmKMOY8iB9tVdZbvaxGpg6dYmgO3Ae8XYz/XApmqmuX9jESgK7DVZz87ve+dLmb2MklISKBy5cr06NEjGLtzzEMPPeR0hFJr0qQJ/fr1Y/r06QwePJj69es7HckYU4gSX/6rqrtVNVlVJ3lvTCyO2sAun9e7vcuK6wIRSRWR9SJyx7lW8p6CSxWR1B9//PGcH3b06FESExPp0aMHF154YQlimGAbPXo0kZGRjBw50ukoxphzKO3sv8F2hfcqgnuAl0XkfwpbSVVnqWprVW1dq1atc37Y/Pnzyc7Odv1pLYCbbrqJm266yekYpXbZZZfx2GOP8f7777Np0yan4xhjChGsItkD1PV5Xce7rFhUdY/3zyxgNdCyLGHmzJlDgwYNuP7668vyMeVCnz59yv1jg5988klq1qxpU6cYE6KCVSQbgYYiUk9EfgX0BIp19ZWIVBORKO/3NYF2+IytlFRmZiZr1qxx5QSNhXFDkVStWpVnn32WTz/9lJSUFKfjGGMKCEqRqOopYBCQAnwLfKCqW0RkjIjcDiAibURkN56bH98QkS3ezX8LpIrIV8DfgYmqWuoimTt3LhUqVOC+++4ry1+p3MjNzSU3N9fpGGX28MMPU79+fZ566iny8vKcjmOM8VHkne3lVWF3tufl5XHFFVfQvHlzli1b5lCy4MofH1m9erWjOfwhMTGRXr168fbbb4fNPwSMCbZA3dnuGitWrGDPnj1hMcie78EHH+TBBx90OoZf9OjRg9atWzNy5EiOHz/udBxjjFdYFUlCQgI1atQodxMZlkXv3r1d85yVChUqMHnyZHbt2sVrr73mdBxjjFfYFMn+/fv5+OOP6d27N1FRUU7HCZpjx45x7Ngxp2P4zc0330xcXBzjx4/n4MFy/0w1Y1whbIrkvffeIzc3N+wmaLztttu47bbbnI7hVxMnTuTw4cNMmDDB6SjGGMJosL1FixZERkZS1NTybpOUlARQbufbOpc+ffqQmJhIRkYGV1xxhdNxjHENG2w/h7S0NL766quwGmTPV54nbTyfsWPHAvDss886nMQYExZFkpCQQFRUFL16FWuyYlc5fPgwhw8fdjqG39WtW5chQ4bwzjvv8NVXXzkdx5iw5voiOX78OO+++y7dunWjWrVqTscJuq5du9K1a1enYwTE8OHD+fWvf81TTz3ldBRjwprri2TRokUcOnQo7AbZ8w0ePJjBgwc7HSMgqlWrxjPPPENKSgorV650Oo4xYcv1g+0dO3Zk27ZtZGVlERER4XQs42fHjx+ncePGVK9endTUVCpUcP2/jYwJKBtsL+D7779nxYoV9OnTJ2xLZP/+/a5+VO0FF1zAuHHjSEtLIzEx0ek4xoQlVx+RdO3aleeee46srCzq1avndCRHuGmurXM5ffo0rVq14hpenvIAABFUSURBVL///S/btm0LqxtOjfE3OyIpYM6cOdxyyy1hWyIAjz/+OI8//rjTMQKqQoUKTJo0iZ07dzJjxgyn4xgTdlx7RNKoUSP97rvveOedd/jTn/7kdBwTBH/4wx/YtGkTO3bs4Ne//rXTcYwpl+yIxMeuXbuoXLmyay99La69e/eyd+9ep2MExaRJkzh48CCTJk1yOooxYcW1RyQiopGRkbRv356UlJSwHWwPhzESX71792bBggVs376dOnXqOB3HmHKnNEckri4SgJiYGN5//306d+7sdCRHLF++HIBOnTo5nCQ4du7cSaNGjejduzezZ892Oo4x5Y6d2irE0aNHSU9PdzqGYzp16hQ2JQJw5ZVXMmjQIObOncvmzZudjmNMWAhqkYhIJxHJEJFMERleyPs3isiXInJKRLoXeO9+Ednu/bq/uPuMjo6mRYsW/ohfLu3atYtdu3Y5HSOonn76aS688EKGD//Ff2LGmAAIWpGISAQwDYgDmgC9RKRJgdW+B/oA7xXYtjowCogFrgVGiUiRE2fFxMQQGxtLXFxc2f8C5dS9997Lvffe63SMoKpRowYjRozgk08+Yc2aNU7HMcb1gnlEci2QqapZqnoSSATOuqRKVXeq6tfA6QLbdgRWqOpBVT0ErADOe77msssu4/333w/rgXaAkSNHMnLkSKdjBN3gwYOpU6cOTz75JG4dBzQmVASzSGoDvudYdnuX+W1bEekvIqkikhoZGUnnzp3DukQAOnToQIcOHZyOEXSVKlVizJgxfPHFF8yfP9/pOMa4mqsG21V1lqq2VtXWtWrVcjpOSMjKyiIrK8vpGI647777aNasGSNGjODkyZNOxzHGtYJZJHuAuj6v63iXBXrbsNavX7+wfDIkQEREBJMmTWLHjh3MmjXL6TjGuFbQ7iMRkUjgO+BWPCWwEbhHVbcUsu5cYKmqzve+rg5sAq7xrvIl0EpVD55rfwWf2R6u8geb27dv73ASZ6gqt9xyC1u2bCEzM5MqVao4HcmYkBbyNySKyG3Ay0AEkKCqL4jIGCBVVReLSBvgI6AacBzYq6pNvdv2A572ftQLqjrnfPuyIjH5Nm7cyLXXXkt8fDxNmzalZcuWxMXFhf34mTGFCfkiCSYrEo+MjAwAGjVq5HAS5+Tl5XHZZZexb98+RITo6GhiY2PD/oo+Ywpjd7abXxgwYAADBgxwOoajkpOTyc7OBjynurKzs9mwYQPJyckOJzPGHSKdDmACa/z48U5HcFxaWho5OTlnLcufOidc52Azxp+sSFzuuuuuczqC41q2bEl0dPSZoxLwXNHVvHlzB1MZ4x52asvlNm/eHPaTF8bFxREbG0tMTAwiQsWKFTl16hSLFi3i9OmCkygYY0rKjkhcbtCgQUD4PI+kMBEREaSkpJCcnEx6ejpXX301n3/+ORMmTCAiIoKZM2dSoYL9m8qY0rIicbkpU6Y4HSEkRERE0Llz5zNjIp07d0ZEzowhWZkYU3pWJC7Xpk0bpyOEJBFh3LhxAFYmxpSRFYnL5T/UK5yfyXIuVibG+IcVicsNHToUCO8xkvPJLxNVZcKECYgIM2bMsDIxpgSsSFzu5ZdfdjpCyBMRXnjhBQAmTJgAYGViTAlYkbicndIqHisTY0rPisTlNm7cCNige3EULBMRYfr06VYmxhTBisTlnnjiCcDGSIorv0xUlYkTJwJYmRhTBCsSl3v99dedjlDu+N5fYmViTNGsSFyuWbNmTkcolwqWiYgwbdo0KxNjCmFF4nLr1q0DbPLG0ijsyMTKxJhfsiJxuaef9jxU0sZISie/TFSVSZMmAVYmxhRkReJyb7zxhtMRyj0ROXNJsJWJMb8U1CIRkU7AK3ie2f6Wqk4s8H4U8DegFXAAiFfVnSJyJfAtkOFddb2qPhys3OVZOD9i158KlomI8Prrr1uZGEMQi0REIoBpwO+B3cBGEVmsqlt9VnsAOKSqDUSkJzAJiPe+t0NV7e66ElqzZg0A7du3dzhJ+VfYkYmViTHBPSK5FshU1SwAEUkEugK+RdIVeN77/XzgdRGRIGZ0nVGjRgE2RuIv+WWiqkyePBnwnOay/0xNOAtmkdQGdvm83g3EnmsdVT0lIoeBGt736olIGnAEGKmq/whwXldISEhwOoLriMiZq7isTIwpP4PtPwCXq+oBEWkFfCwiTVX1iO9KItIf6A9w+eWXOxAz9NSvX9/pCK5kZWLMz4JZJHuAuj6v63iXFbbObhGJBKoCB1RVgRMAqrpJRHYAvwFSfTdW1VnALIDWrVtrIP4S5c3KlSsB6NChg8NJ3Ce/TFT1zJMorUxMOApmkWwEGopIPTyF0RO4p8A6i4H7gc+B7sCnqqoiUgs4qKp5IlIfaAhkBS96+ZX/4CYrksAQkTMD71OmTDlzNZeViQknQSsS75jHICAFz+W/Caq6RUTGAKmquhiYDcwTkUzgIJ6yAbgRGCMiucBp4GFVPRis7OXZvHnznI7gegXLBLAyMWElqGMkqroMWFZg2XM+3x8H7i5kuwXAgoAHdKG6desWvZIpMysTE87Ky2C7KaXly5cD0KlTJ4eTuF9+magqU6dOBaxMTHiwInG5/CuLrEiCQ0TOXMU1depURITXXnvNysS4mhWJyyUmJjodIewULBPAysS4mhWJy11yySVORwhLViYmnFiRuNySJUsA6NKli8NJwk9+magqL774IqpKp06dSE9Pp2XLlsTFxREREeF0TGPKTDz3+rlP69atNTU1tegVXe6mm24CbK4tJ6kqjz/+OH/961+pWLEip06dIjo6mtjYWFJSUqxMTEgRkU2q2rok29i0pS43f/585s+f73SMsCYi3HzzzVSsWJHc3FxUlezsbDZs2EBycrLT8YwpMysSl6tZsyY1a9Z0OkbYS09P59SpU2cty87O5t133yUnJ8ehVMb4hxWJyy1cuJCFCxc6HSPstWzZkujo6F8sT0xM5KKLLuJPf/oTixcv5sSJEw6kM6ZsrEhc7tVXX+XVV191OkbYi4uLIzY2lpiYGESEmJgYbr75ZpYvX06vXr1Yvnw5Xbt25eKLL6ZPnz4kJyeTm5vrdGxjisUG213u8OHDAFStWtXhJCYvL4/k5GTS09Np0aLFWVdt5ebmsmrVKpKSkvjoo484fPgw1atXp1u3bsTHx3PTTTcRGWkXWZrAK81guxWJMSHmxIkT/O///i9JSUksWrSI7OxsatWqRffu3YmPj+f666+3K71MwFiR+LAi8UhKSgIgPj7e4SSmNHJyckhOTiYpKYklS5aQk5PDpZdeeqZU2rZta8+MN35lReLDisTD7iNxj6NHj7J06VKSkpJYtmwZJ06coG7dutx9993Ex8fTpk0bu3PelJkViQ8rEo9jx44BULlyZYeTGH86cuQIixcvJikpiZSUFHJzc6lXrx49evQgPj6eFi1aWKmYUrEi8WFFYsLFoUOH+Pjjj0lKSmLlypXk5eXRsGFD4uPjiY+Pp1mzZk5HNOWIFYkPKxKPd955B4DevXs7nMQEw/79+1m4cCFJSUmsXr2a06dP06RJkzOl0qhRozNXj6WlpdmcX+YXrEh8WJF42BhJ+Nq7dy8LFiwgKSmJzz77DFWlefPmHD16lB9++IGcnByb88v8ghWJDysSj/yb2ipWrOhwEuOkPXv28OGHHzJz5kwyMjLOeq9ChQrExsZy1VVXUbNmTWrUqHHWn/nfV61a1W/jLnZUFLpKUyRBvcNJRDoBrwARwFuqOrHA+1HA34BWwAEgXlV3et8bATwA5AGDVTUliNHLLSsQA1C7dm2GDh3KTz/9xKhRo/D9B+Tp06fJzMxkx44dHDhwgLy8vEI/IzIykurVq59VLoUVju+fVatW/cXlyXl5eXTs2JENGzZw9OhRR4+KQqHQQiGDbw7g0pJuG7QiEZEIYBrwe2A3sFFEFqvqVp/VHgAOqWoDEekJTALiRaQJ0BNoClwGrBSR36hq4f/FmzPmzp0LQJ8+fRzNYUJD/pxf2dnZZ5bFxMSQkJBA586dUVUOHz7MgQMH2L9/P/v37z/zfcFl3333HevWrePAgQO/mJAyX0RExJnyyS+XY8eOsXbt2jNHy9nZ2Xz22WcMHz6cNm3aEBEREZQvgD/+8Y988cUXjhVaqJSqbw48v2NLJGintkSkLfC8qnb0vh4BoKoTfNZJ8a7zuYhEAnuBWsBw33V91zvX/uzUloeNkRhfgfjFpaocOXKkWOWzf/9+du7ceVaRhRoR+cVRVMFTer6vz/deUeuePn260Ik6o6KigjolzqlTp87KoaolOocZzFNbtYFdPq93A7HnWkdVT4nIYaCGd/n6AtvWLrgDEekP9Pe+PCEim/0TvdyrKSL7nQ4RImoC9rOAqkDN7Ozs/atWrToc5Hm8qgL1OXvS2NNAFnA4iDkupZB/favqv/Py8n5wMsOJEyf+feLEiWBlOGeO4nLVLHCqOguYBSAiqSUdMHIr+1n8zH4WP7Ofxc/sZ/EzESnxqZxgTtKzB6jr87qOd1mh63hPbVXFM+henG2NMcY4IJhFshFoKCL1RORXeAbPFxdYZzFwv/f77sCn6hnEWQz0FJEoEakHNAS+CFJuY4wx5xG0U1veMY9BQAqey38TVHWLiIwBUlV1MTAbmCcimcBBPGWDd70PgK3AKWBgMa7YmhWov0s5ZD+Ln9nP4mf2s/iZ/Sx+VuKfhWtvSDTGGBMc9iADY4wxZWJFYowxpkxcWSQi0klEMkQkU0SGO53HKSJSV0T+LiJbRWSLiAxxOpOTRCRCRNJEZKnTWZwkIr8Wkfkisk1EvvXeLByWROQx7/8bm0XkfRG5wOlMwSQiCSKyz/eeOxGpLiIrRGS7989qRX2O64rEZyqWOKAJ0Ms7xUo4OgU8rqpNgN8BA8P4ZwEwBPjW6RAh4BVguao2Bq4mTH8mIlIbGAy0VtVmeC4C6ulsqqCbC3QqsGw4sEpVGwKrvK/Py3VFAlwLZKpqlqqeBBKBrg5ncoSq/qCqX3q//wnPL4xfzAgQDkSkDvBH4C2nszhJRKoCN+K5QhJVPamq/3U2laMigUre+9YqA/92OE9QqepaPFfI+uoKvO39/m3gjqI+x41FUthULGH5y9OXiFwJtAQ2OJvEMS8DT+KZiiOc1QN+BOZ4T/O9JSLRTodygqruAaYC3wM/AIdV9X+dTRUSLlbV/OlZ9gIXF7WBG4vEFCAiMcACYKiqHnE6T7CJSGdgn6pucjpLCIgErgFmqGpL4CjFOHXhRt5z/13xlOtlQLSI2KNEfXhvCC/yHhE3FolNp+JDRCriKZF3VXWh03kc0g64XUR24jnVeYuIvONsJMfsBnarav6R6Xw8xRKOOgD/UtUfVTUXWAhc53CmUPAfEbkUwPvnvqI2cGORFGcqlrAgnjmrZwPfqupLTudxiqqOUNU6qnolnv8ePlXVsPyXp6ruBXaJSCPvolvxzBgRjr4Hficilb3/r9xKmF54UIDvVFX3A4uK2sBVs//CuadicTiWU9oB9wLfiEi6d9nTqrrMwUzGeX8G3vX+QysL6OtwHkeo6gYRmQ98iecKxzTCbKoUEXkfuAnPoyZ2A6OAicAHIvIA8H9AjyI/x6ZIMcYYUxZuPLVljDEmiKxIjDHGlIkViTHGmDKxIjHGGFMmViTGGGPKxIrEGGNMmViRGGOMKRMrEmOCSEQWisg4EVkrIt+LSAenMxlTVlYkxgTXVcB/VfVGPM9H+ZPDeYwpMysSY4JERCoDVYG/ehdVBML5WSDGJaxIjAmeJsAmVc3zvm4ObD7P+saUC1YkxgTPVUC6z+vmwNcOZTHGb6xIjAmegkXSDDsiMS5gs/8aY4wpEzsiMcYYUyZWJMYYY8rEisQYY0yZWJEYY4wpEysSY4wxZWJFYowxpkysSIwxxpTJ/wPqdXdrYK+LTgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "n = 10 # number of trials\n", "p = 0.2 # probability of success in each trial\n", "x = np.arange(n+1) # number of successes\n", "y = stats.binom.pmf(x,n,p) # generate Binomial probability function\n", "mu = n*p # mean of distribution\n", "xmin,xmax = 0.,float(n)\n", "ymin,ymax = 0.,0.4\n", "plt.plot(x,y,marker='o',markersize=5,color='black')\n", "plt.plot([mu,mu],[ymin,ymax],color='black',linestyle='dotted')\n", "plt.xlabel(r'$n$')\n", "plt.ylabel(r'$P_{\\rm Binomial}(n)$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Poisson distribution\n", "\n", "The Poisson distribution applies to a **discrete random process where we are counting** a quantity in a fixed interval (of space or time). Example: radioactive decay, photons arriving at a CCD. If the mean number of events expected in some interval is $\\mu$, the probability of observing $n$ events is:\n", "\n", "$P_{\\rm Poisson}(n) = \\frac{\\mu^n \\, e^{-\\mu}}{n!}$\n", "\n", "The mean and variance of this distribution are equal, $\\overline{n} = {\\rm Var}(n) = \\mu$.\n", "\n", "The variance of the distribution allows us to set the **Poisson error** for a counting experiment: if a bin contains $N$ events, by assuming the mean count is the observed count we can place an error:\n", "\n", "Count = $N \\pm \\sqrt{N}$\n", "\n", "Here is an example Poisson distribution for $\\mu = 5$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 0.25)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3iUVfr/8fed0EyCCmIDgoIFFZHkJ0sUXUWKEqUouAIr2FDABQVFXGFVliYoIKBShUhRygookDUEUNoKBnATkLpSpCmuCLISWsr9+yMD3xBCmCSTOTPJ/bquuTLlKZ91Q+45z3nOOaKqGGOMMRcS4jqAMcaY4GAFwxhjjFesYBhjjPGKFQxjjDFesYJhjDHGK1YwjDHGeMVvBUNEmorINhHZLiKv5fL5yyKyWUQ2iMiXInJNts8yRCTF85jvr8zGGGP+j/hjHIaIhAL/AZoA+4C1QDtV3Zxtm/uAJFU9JiLPAw1UtY3ns6OqGlHkQY0xxpyXv1oY9YDtqrpTVU8BM4GW2TdQ1aWqeszz8hugqp+yGWOM8UIpP52nCrA32+t9QEwe23cEErK9Lici64B0YIiqfp7bTiLSCegEEB4efvtNN91UqNDGGFPSfPvttwdV9fLcPvNXwfCaiLQH6gL3Znv7GlXdLyI1gK9E5DtV3ZFzX1WdAEwAqFu3rq5bt84vmY1v7d2b9d0iMjLScRJjSh4R2X2+z/xVMPYD2f/1V/W8dxYRaQz8DbhXVU+efl9V93t+7hSRZUA0cE7BMMVDhw4dAFi2bJnbIMaYs/irYKwFbhCR6mQVirbAn7NvICLRwHigqar+N9v7FYBjqnpSRCoBdwHv+Cm3ceD11193HcEYkwu/FAxVTReRbkAiEArEqeomEekPrFPV+cBQIAL4VEQA9qhqC+BmYLyIZJLVST8k+91Vpvhp3Lix6wjGmFz45bZaF6wPI3jt3LkTgBo1ajhOYkzJIyLfqmrd3D4LuE5vY5555hnA+jCMCTRWMEzA6devn+sIxphcWMEwAefee++98EbGGL+zyQdNwNm2bRvbtm1zHcMYk4O1MEzA6dy5M2B9GMYEGisYJuC89dZbriMYY3JhBcMEnPr167uOYIzJhfVhmICzceNGNm7c6DqGMSYHa2GYgNOtWzfA+jCMCTRWMEzAGTp0qOsIxphcWMEwAecPf/iD6wjGmFxYH4YJOCkpKaSkpLiOYYzJwVoYJuD06NEDsD4MYwKNFQwTcEaOHOk6gjEmF1YwTMCJiopyHcEYkwvrwzABZ+3ataxdu9Z1DGNMDtbCMAGnV69egPVhGBNorGCYgPPBBx+4jmCMyYUVDBNwbr31VtcRjDG5sD4ME3BWrVrFqlWrXMcwxuRgLQwTcPr06QNYH4YxgcYKhgk448ePdx3BGJMLKxgm4NSsWdN1BGNMLqwPwwSc5cuXs3z5ctcxjDE5WAvDBJy+ffsC1odhTKCxgmECTlxcnOsIxphcWMEwAadGjRquIxhjcmF9GCbgLFmyhCVLlriOYYzJwVoYJuAMHDgQgMaNGztOYozJzgqGCTjTpk1zHcEYkwsrGCbgREZGuo5gjMmF9WGYgLNw4UIWLlzoOoYxJgdrYZiAM2TIEACaNm3qOIkxJjsrGCbgzJw503UEY0wurGCYgHPVVVe5jmCMyYVf+zBEpKmIbBOR7SLyWi6fvywim0Vkg4h8KSLXZPvsSRH53vN40p+5jX8tWLCABQsWuI5hjMnBby0MEQkFRgNNgH3AWhGZr6qbs22WDNRV1WMi8jzwDtBGRCoCfYG6gALfevY97K/8xn+GDx8OQPPmzR0nMcZk589LUvWA7aq6E0BEZgItgTMFQ1WXZtv+G6C95/kDwGJVPeTZdzHQFJjhh9zGz2bPnu06gjEmF/4sGFWAvdle7wNi8ti+I5CQx75Vcu4gIp2ATgDVqlUrTFbjUKVKlVxHMMbkIiDHYYhIe7IuPw3Nz36qOkFV66pq3csvv7xowpkiN3fuXObOnes6hjEmB3+2MPYD2YfwVvW8dxYRaQz8DbhXVU9m27dBjn2XFUlK49x7770HQKtWrRwnMcZk58+CsRa4QUSqk1UA2gJ/zr6BiEQD44GmqvrfbB8lAm+JSAXP6/uB3kUf2bgwb9481xGMMbnwW8FQ1XQR6UbWH/9QIE5VN4lIf2Cdqs4n6xJUBPCpiADsUdUWqnpIRAaQVXQA+p/uADfFzyWXXOI6gjEmF6KqrjMUibp16+q6detcxzAFMGvWLADatGnjOIkxJY+IfKuqdXP7zEZ6m4AzduxYwAqGMYHGCoYJOF988YXrCMaYXFjBMAEnLCzMdQRjTC4CchyGKdk+/vhjPv74Y9cxjDE5WAvDBJyJEycC0L59+wtsaYzxJysYJuAsXrzYdQRjTC6sYJiAU7p0adcRjDG5sD4ME3AmT57M5MmTXccwxuRgBcMEHCsYxgQmuyRlAs6yZctcRzDG5MJaGMYYY7xiBcMEnA8//JAPP/zQdQxjTA52ScrkW0ZGBgkJCSQnJxMdHU1sbCyhoaE+O/bo0aP5/fffufrqq316bGNM4VjBMPmSkZHBAw88QFJSEqmpqYSHhxMTE0NiYmKh/7CfPvaOHTtITU2lXbt2Pju2MabwrGCYfElISCApKYmjR48CcPToUVasWMEjjzxC9erVC3XsXbt2sWLFCtLS0s4cOykpiYSEBJo1a1bo7MaYwrGCYfIlOTmZ1NTUs95LS0tj8eLFlCtXrlDHPnHixJlicVpqaiopKSlWMIwJAFYwTL5ER0dTqlSps/6wR0REMGPGjEL/UY+Pj6ddu3ZnWi+QNXNtVFRUoY5rjPENu0vK5Et6ejppaWmULl0aESEiIoKYmBhiY2MLfezY2FhiYmKIiIjAs0QvYWFh3H///YU+tjGm8GyJVuO1/fv3U6dOHSIjI3nzzTfZtGkTUVFRPr9LKiEhgZSUFH744QcmTZrESy+9xLvvvuuT4xtj8mZLtJpCy8jI4PHHH+fEiRPMmjWLG2+8kUceecTn5wkNDWXHjh2UL1+eDz/8kPDwcEaMGME111xD9+7dfX4+Y4z37JKU8cqgQYNYvnw5o0eP5sYbbyzSc3355Zd8+eWXiAjvvvsujzzyCC+99BJz5swp0vMaY/Jml6TMBa1cuZIGDRrQrl07pk2bdqZ/wV+OHz9Ow4YNSUlJ4csvv6R+/fp+Pb8xJUlel6SsYJg8HTp0iKioKMqUKUNycjLly5d3kuOXX36hfv36HD58mFWrVhV5K8eYkiqvgmGXpMx5qSodO3bkwIEDzJw502/FYtiwYQwbNuys9y6//HISEhIQEWJjY/nvf//rlyzGmP9jBcOc19ixY/n8888ZMmQIdevm+oWjSKxevZrVq1ef8/71119PfHw8P/30E82bN+fYsWN+y2SMsUtS5jw2bNhAvXr1aNiwIfHx8YSEBM53i88//5xWrVrRokUL5syZY/NMGeNDdknK5Etqaipt2rShQoUKTJ48OaCKBcDDDz/Me++9x7x58+jevTvF9UuPMYHGxmGYc3Tv3p1t27axePFirrjiCr+ff8iQIQC89tpr592mW7du7N69m2HDhnHttdfyyiuv+CueMSWWFQxzllmzZjFp0iR69+5No0aNnGRISUnxaru3336bPXv20KtXLyIjI2nTpk0RJzOmZLM+DHPGrl27iIqKolatWixfvpzSpUu7jnRBJ06coEmTJqxZs4bFixdzzz33uI5kTFCzPgxzQWlpabRr1w4RYfr06UFRLADKlSvHvHnzqF69Oi1btmTLli2uIxlTbFnBMAC88cYbJCUlMXHiRK699lqnWQYMGMCAAQO83r5ixYokJCRQtmxZYmNjOXDgQBGmM6bksoJhWLRoEW+//TadOnXi0UcfdR2Hbdu2sW3btnztU716deLj4/nll1946KGHzlpTwxjjG9aHUcL9/PPP1KlTh0qVKrFmzRrCwsJcRyqU+Ph4WrZsSdOmTZk3bx6lStl9HcbkR0D0YYhIUxHZJiLbReSc+yVF5B4R+beIpIvIozk+yxCRFM9jvr8yF3eZmZk88cQTHDlyhJkzZwZ9sQBo1qwZY8aM4YsvvqBr1642RsMYH/LL1y8RCQVGA02AfcBaEZmvqpuzbbYHeArI7Yb646pq63T62PDhw1m0aBHjxo3j1ltvdR3njDfffBOA/v37F2j/zp07s3v3bgYPHsw111xDnz59fBnPmBLLX+31esB2Vd0JICIzgZbAmYKhqj94Psv0U6YSbc2aNfTp04fWrVvTqVMn13HOsnfv3kIfY+DAgezevZu//e1vVKtWjfbt2/sgmTElW74LhoiEAydUNSMfu1UBsv8V2AfE5GP/ciKyDkgHhqjq5+fJ1gnoBFCtWrV8HL5kOXLkCO3ataNy5cp8+OGHfl/f4kI++uijQh8jJCSEuLg4fvzxR5555hkqV65Mw4YNfZDOmJLrggVDREKAtsDjwB+Ak0BZETkI/BMYr6rbizQlXKOq+0WkBvCViHynqjtybqSqE4AJkNXpXcSZgpKq0qVLF3bv3s2KFSuoUKGC60hFpmzZsnz22WfcddddPPzwwwwePJhDhw4RHR3t03XIjSkpvGlhLAWWAL2BjaqaCSAiFYH7gLdF5DNV/TiPY+wHIrO9rup5zyuqut/zc6eILAOigXMKhrmwjz76iJkzZzJo0KCAXbmud+/eAAwePLjQx7r00kuJj4/npptu4oUXXgAgPDycmJgYEhMTrWgYkw/e3CXVWFUHqOqG08UCQFUPqeocVW0NzLrAMdYCN4hIdREpQ1aLxau7nUSkgoiU9TyvBNxFtr4P470tW7bwwgsv0LBhQ/7617+6jnNev/76K7/++qvPjrdp0yZCQ0NRVVSVo0ePkpSUREJCgs/OYUxJcMEWhqqmFXYbVU0XkW5AIhAKxKnqJhHpD6xT1fki8gfgM6AC0FxE+qlqLeBmYLynMzyErD4MKxj5dOLECdq2bUtYWBjTpk0L6G/WEyZM8OnxkpOTOXHixFnvpaamkpKSQrNmzXx6LmOKM687vUWkIVn9GL8BG4ENZF2iOunN/qr6BfBFjvfezPZ8LVmXqnLutwqo7W1Ok7tevXqxYcMG/vnPf1K5cmXXcfwqOjqa8PDws0Z/h4WFERVld2obkx/5GbgXBywAvgFqAG8Cm4oilPGNjIwM4uPj+fOf/8wHH3xA9+7defDBB13HuqBXXnnFp+tbxMbGEhMTQ0RExJk7wi6++GJiY2N9dg5jSoL83Fa7O9vtrJ8WRRjjOxkZGTzwwAOsXr2aY8eOERISwoYNG8jIyAjoy1EAx48f9+nxQkNDSUxMJCEhgZSUFFJSUpgzZw6rVq3ij3/8o0/PZUxx5vVcUiIyADgEjNQgmG+hpM8lFR8fT7t27c66DBMREcGMGTNK/HX71NRUatWqRUREBMnJyUEzlbsx/uCruaRuAZ4HfhKRf4rIIBH5k08SGp9LTk4mNTX1rPdOd/SWdOHh4bz//vts2rSJESNGuI5jTNDwumCoamtVvRGoTlb/xffAHUUVzBROdHT0OZeewsPDg6Kjt0ePHvTo0aNIz9G8eXNatmxJv3792L17d5Gey5ji4oIFQ3LMG6Gqx1X1W1WdrKo9c9vGuHfrrbeSnp5O6dKlEREiIiKIiYmxjt5s3nvvPQC6d+/uOIkxwcGrkd4iMgeYp6p7Tr/pGYB3N/AkWaPBJxdJQlMg77//PiEhIYwbN44ff/yRqKiooJkOY+TIkX45T7Vq1ejbty9//etfmT9/Pi1atPDLeY0JVhfs9BaRcsAzZI3BqE7WOIyLyGqdLALGqGpyEefMt5Lc6f3bb78RGRlJixYt+OSTT1zHCWhpaWlER0fz+++/s3nzZsLDw11HMsapQnV6q+oJVR2jqncB1wCNgGhVvUZVnwvEYlHSjR8/nqNHj9KrVy/XUQqka9eudO3a1S/nKl26NOPGjWPPnj35WkfcmJIoXyvuqWqaqv4EPC8in4jIdBGZXkTZTAGcPHmSUaNG0bhx46Do4M7NRRddxEUXXeS389199908/fTTDB8+nE2bbCyqMedT0AWUQlT1cZ8mMT4xffp0fvrpJyZPnuw6SoENGzbM7+d85513mDdvHs8//zzLly8PuDVCjAkEBV3T+zoR+ZOIPCgigT/XRAmRmZnJsGHDqFOnDk2aNHEdJ6hUqlSJd955h5UrVzJlyhTXcYwJSAUtGMuBMOByz8MEgISEBDZv3swrr7wS1N+QO3Xq5GTZ2Keffpr69evTq1cvn06vbkxxUaCCoapTsj98HcoUzNChQ4mMjKRNmzauoxTKZZddxmWXXeb384aEhDB27FgOHz58ZhEnY8z/KVAfhog0BzoAmcAMVZ3n01Qm39auXcvy5csZPnx40M+N5IuV9grqtttuo0ePHgwfPpynnnoqYFclNMYFrycfPGsnkQmq2snzfKyqPu/zZIVU0sZhPPbYYyxatIi9e/dSvnx513GC2tGjR7n55pupWLEi3377LaVKFfTeEGOCj68mH8zuIhGpJiLVABvp5NjOnTuZM2cOXbp0KRbF4umnn+bpp592dv6IiAjee+89NmzYcGb6EGNMwQvGAOAFYDhg/6Ice/fddwkNDeXFF190HcUnIiMjiYyMdJrh4Ycf5qGHHuLNN99k7969TrMYEygKeklqKDAKGApkqGp7XwcrrJJySergwYNUq1aNtm3bEhcX5zpOsbJr1y5q1apFbGwsc+bMcR3HGL8oiktSFwMtgMHAjwUNZgpvzJgxHD9+3KdLmpos1atX54033mDu3Ll88cUXF97BmGKuoAVjGVBBVTeQtS6GceD48eN88MEHPPTQQ9xyyy2u4/hM+/btad8+MBqtPXv25Oabb6Zbt24cO3bMdRxjnCpowZitqoNE5DpsfW9npkyZwi+//BK0kwyeT82aNalZs6brGACUKVOGMWPGsGvXLt566y3XcYxxqjB9GCOBYVgfhhMZGRncdNNNVKhQgaSkpKAe2R0MnnjiCWbOnMn69eu5+eabXccxpsgUVR9GS6wPw5l58+axfft2evXqZcXCD4YNG0Z4eDh/+ctfKMiXLGOKA+vDCEKqytChQ6lRowatWrVyHcfn2rZtS9u2bV3HOMsVV1zBkCFDWLZsmS1KZUqsfBcMEbkLOAKkeGaq3e/zVCZPX3/9Nd988w0vv/xyUCy5ml9RUVEBuZbHc889R0xMDD179uTw4cOu4xjjd/nuw/DMI1UROLOjqk71ca5CK859GC1btuTrr79mz549hIWFuY5ToiQnJ1O3bl06derE2LFjXccxxud82oehqguAw0Az4CGyWhvGT7Zu3cr8+fPp2rWrFQsHoqOjefHFFxk/fjxr1qxxHccYvypoH0YzVX1MVdsATX0ZyORt+PDhlCtXjm7durmOUmRat25N69atXcc4r/79+3P11VfTpUsX0tPTXccxxm8KM/ngNTb5oH8dOHCAqVOn8tRTT3H55cV33ao777yTO++803WM8ypfvjwjR44kOTmZMWPGuI5jjN8UdBzGdUAXz8sJqhpwd0oVxz6Mv/3tbwwePJht27Zxww03uI5ToqkqsbGxrFq1iq1bt1K5cmXXkYzxibz6MArS6d0LqAdMUdV4H+QrEsWtYBw9epRq1apx33332UR4AWL79u3UqlWLevXqcf/99xMdHU1sbGyxvHPNlBx5FYyCrAxzi6r+SUTGAwFbMIqbSZMmcfjw4WI3DUhuWrRoAcD8+fMdJ8lb9erVqVq1Kv/617/4+uuvCQ8PJyYmhsTERCsaplgqSB9GJc/4iytF5EHPc1OE0tPTGTFiBHfffTd33HGH6zhFrlGjRjRq1Mh1jAtKSEjg559/BrIuUR09epSkpCQSEhIcJzOmaBSkhTEbuBz4zPPT5kkoYp9++im7d+8uMau/de/e3XUEryQnJ58zg21qaiopKSk0a9bMUSpjis4FWxgi8qSIHBSRQyIyFZirqlOyPbwetCciTUVkm4hsF5HXcvn8HhH5t4iki8ijueT43vN40ttzBrvT04DcdNNN9kcowERHRxMefvZNguXKlQvIUerG+II3l6TeAJoANwG7gQLN8SwiocBoIBa4BWgnIjkXcdgDPAVMz7FvRaAvEENWh3tfEalQkBzB5quvviI5OZmePXsSElLQu6CDS2xsLLGxsa5jXFBsbCwxMTFEREQgIogIISEhNGzY0HU0Y4qEN5ek/qeqyZ7nb4hIUgHPVQ/Yrqo7AURkJlkz3m4+vYGq/uD5LDPHvg8Ai1X1kOfzxWQNGJxRwCxBY+jQoVx55ZUBs6CQPzRv3tx1BK+EhoaSmJhIQkICKSkppKWl0b9/f95++2369evnOp4xPudNwbhaRDoBW4EtQOkCnqsKsDfb631ktRgKum+VnBt5cnYCqFatWsFSBpANGzaQmJjIoEGDKFeunOs4fvOXv/zFdQSvhYaG0qxZszOXC08vtNSqVSvq1KnjOJ0xvuXNNY6+QG1gALANuFVEvhCRwSLSrkjT5ZOqTlDVuqpatziMhD69BsPzzz/vOorx0siRI7nssst4+umnSUtLcx3HGJ+6YMHw/BF+QVXvVdWKQA3gfeA3ID+31O4HIrO9ror3U6MXZt+gtHfvXmbMmMGzzz5LhQolorvmjMaNG9O4cWPXMQqkYsWKjB49muTkZIYNG+Y6jjE+le/balV1H1mXhPJ7s/la4AYRqU7WH/u2wJ+93DcReCtbR/f9QO98nj+ojBo1ClXlpZdech3F79q0aeM6QqG0bt2aRx99lH79+vHwww/bkq6m2CjQXFIFPlnWIL+RQCgQp6qDRKQ/sE5V54vIH8ga31EBOAEcUNVann2fAfp4DjVIVT/K61zBPDXIkSNHiIyMpFmzZkyfPv3CO5iA8/PPP3PLLbdQs2ZNVq5caSO/TdAoijW9C0RVv1DVG1X1OlUd5HnvTVWd73m+VlWrqmq4ql52ulh4PotT1es9jzyLRbAbP348v//+e4mYBqS4uvLKK3nvvfdYvXo177//vus4xviEX1sY/hSsLYxTp05RvXp1br75ZpYsWeI6jhMNGjQAYNmyZU5zFJaq0rx5c7766iu+++47rrvuOteRjLkgX08+aIrQ9OnT+fHHH4mLi3MdxZmnnnrKdQSfEBHGjRtHrVq1eO6551iyZEmJGXxpiidrYQQQVaV27dqEhoaSkpKCiLiOZHxg4sSJPPfcc4wbN47OnTu7jmNMngKmD8PkLSEhgU2bNvHKK6+U6GKRlpZWrMYwdOzYkUaNGtGrVy/27NnjOo4xBWYFI4AMHTqUqlWr0rZtW9dRnGrSpAlNmjRxHcNnRIQPP/yQjIwMOnfuTHFt1ZvizwpGAMjIyGDEiBEsW7aM+++/v8Rf53722Wd59tlnXcfwqerVqzNkyBAWLlzItGnTXMcxpkCsD8OxjIwMHnjgAZYvX056ejrh4eHccccdtmpbMZSZmck999zDpk2b2Lx5M1dffbXrSMacw/owAlhCQgKrV68mPT0dyFqAp6Sv2nbs2LFzFiYqDkJCQoiLi+PEiRN07drVLk2ZoGMFw7G8Vm0rqR588EEefLB4rvx744030q9fPz777DNmz57tOo4x+WLjMByrVKnSOe+Fh4eX6FXbivvsvC+//DKffvopXbt25b777sv1d8CYQGQtDMeWLVtGSEgI4eHhiAgRERHExMQExYpzRaVNmzZBPwFhXkqVKkVcXBy//fZb0KxfbgxYC8OpjRs38umnn/Lqq69y9913k5KSQlRUFLGxsSW6w/vIkSMAXHLJJY6TFJ3atWvz+uuv07dvX9q2bRs0qwyaks3uknLoT3/6E4mJiezatYvLLrvMdZyAUVzmkrqQU6dOUbduXX799Vc2bdrEpZde6jqSMXaXVCBav349s2fPpkePHlYscnjxxRd58cUXXccocmXKlOGjjz7i559/5pVXXnEdx5gLshaGI4888ghLly5l165dJW5FPXO23r17M2TIEBYtWlSsRrib4GQtjADz7bff8vnnn/Pyyy9bscjFwYMHOXjwoOsYftO3b19q1qzJc889x9GjR13HMea8rGA48Pe//50KFSrYHTLn8eijj/Loo4+6juE35cqVIy4ujj179tC7d7FeedgEObtLys/WrFlDfHw8gwYNKtZ3ARVGz549XUfwu/r16/Piiy8yatQoHnvsMf74xz+6jmTMOawPw89iY2NZu3Ytu3btonz58q7jmACSmpp6Zj2U9evXExYW5jqSKYGsDyNArFq1ioULF/Lqq69ascjDgQMHOHDggOsYfhceHs7EiRPZvn07ffv2dR3HmHNYC8OPmjRpwvr169m1axfh4eGu4wSskjIO43y6dOnChx9+yOrVq6lXr57rOKaEsTW9A8CKFStYsmQJw4cPt2JxAa+99prrCE698847/POf/+Tpp59m4MCBbNy4kejo6BI/A4Bxz1oYfnLfffexdetWduzYYdemzQUtWLCAFi1aULp06TPrpMTExNg6KabIWR+GY0uXLmXZsmX07t3bioUX9u7dy969e13HcEpEKFWqFGlpaagqR48eLfHrpBj3rGAUMVXlzTffpHLlynTq1Ml1nKDQoUMHOnTo4DqGU8nJyWRkZJz1XklfJ8W4Z30YRWzJkiX861//YvTo0ZQrV851nKDw+uuvu47gXHR0NOHh4WeN/A4LCyvR66QY96wPowipKvXr12f//v18//33lC1b1mkeEzxOr/WelJREamoqqsoll1zCvn37iIiIcB3PFGN2l5QjCxcu5JtvvmHChAlWLPJh586dANSoUcNxEndCQ0NJTEwkISGBlJQUUlNTGTJkCE888QSffvqpdXwbJ6yFUURUlXr16vHrr7+ybds2Spcu7SxLsCnp4zDOZ9SoUfTo0YMXXniBUaNGISKuI5liyFoYDsTHx7Nu3Tri4uKsWORTv379XEcISN27d2fPnj28++67VKtWzdbQMH5nLYwikJmZye23387vv//O1q1bKVXK6rLxjczMTNq1a8c//vEPpk+fTrt27VxHMsWMtTD87PPPPyclJYWpU6dasSiAbdu2AVCzZk3HSYqcAbIAABMiSURBVAJPSEgIU6ZM4cCBAzz55JNcddVV3Hfffa5jmRLCWhg+lpmZSZ06dUhLS2Pjxo1WMArA+jAu7PDhw9x9993s37+flStXUrt2bdeRTDFhLQw/mj17Nhs3bmT69OlWLArorbfech0h4FWoUIGEhATuuOMOHnzwQVavXk3VqlVdxzLFnF9bGCLSFBgFhAITVXVIjs/LAlOB24FfgTaq+oOIXAtsAbZ5Nv1GVbvkdS4XLYyMjAxq166NiLBhwwa79dEUufXr1/PHP/6Ra6+9lpUrV9qiXKbQAmIuKREJBUYDscAtQDsRuSXHZh2Bw6p6PTACeDvbZztUNcrzyLNYuDJr1iy2bNnC3//+dysWhbBx40Y2btzoOkZQqFOnDnPnzmXLli20atWKU6dOuY5kijF/ziVVD9iuqjtV9RQwE2iZY5uWwBTP89lAIwmSm83T09Pp168ftWvXpnXr1q7jBLVu3brRrVs31zGCRuPGjYmLi+Orr77imWeeITMz03UkU0z58yJ7FSD7FKT7gJjzbaOq6SJyBLjM81l1EUkG/ge8rqorc55ARDoBnQCqVavm2/QXMH36dP7zn/8wd+5cQkJsTsfCGDp0qOsIQadDhw7s27ePPn36ULVqVYYMGXLhnYzJp2Dplf0JqKaqv4rI7cDnIlJLVf+XfSNVnQBMgKw+DH+FS0tLo3///kRHR/Pwww/767TF1h/+8AfXEYLSa6+9xp49e3j77beJjIyka9euriOZYsafBWM/EJntdVXPe7lts09ESgGXAL9qVs/8SQBV/VZEdgA3AgGxQtK0adPYsWMH8+fPt+kafOD0FN42M2v+iAjvv/8++/fv54UXXqBKlSr2Bcb4lN/ukvIUgP8AjcgqDGuBP6vqpmzbdAVqq2oXEWkLtFLVx0TkcuCQqmaISA1gpWe7Q+c7n7/ukjp16hQ1a9akUqVKrFmzxgqGD9g4jMI5duwYDRs2ZP369Xz11VfceeedriOZIBIQ4zA8fRLdgESybquNU9VNItIfWKeq84FJwDQR2Q4cAtp6dr8H6C8iaUAm0CWvYuFPkydP5ocffmDMmDFWLHxk5MiRriMEtbCwMBYsWED9+vVp3rw5q1at4sYbb3QdyxQDNtK7EE6ePMkNN9xAlSpVWLVqlRUME1B27NjBnXfeSUREBKtXr+bKK690HckEgYAYh1EcTZo0ib1799K/f38rFj60du1a1q5d6zpG0LvuuuuIj4/n559/5qGHHjpr9T5jCsJaGAV04sQJrrvuOmrUqMGKFSusYPiQ9WH4Vnx8PC1btqRp06bMmzfPpqwxeQqIPoziZsKECfz44498/PHHVix87IMPPnAdoVhp1qwZY8aMoUuXLjz//PNMmDDBfmdNgVjBKIBjx44xePBgGjRoYFNLF4Fbb73VdYRip3Pnzuzdu5dBgwZRpUoV6tatS3JyMtHR0cTGxtpUNsYrVjAKYNy4cRw4cIBZs2a5jlIsrVq1CoD69es7TlK8DBgwgN27d9OvXz/Kli3LqVOnCA8PJyYmhsTERCsa5oKs09tLGRkZxMfH88Ybb9CvXz8aNWrEPffc4zpWsdSnTx/69OnjOkaxIyK0atWK0NBQTp48iapy9OhRkpKSSEhIcB3PBAFrYXghIyODBx54gKSkpDN3mvz2229kZGTYt7IiMH78eNcRiq2NGzeeMzlhamoqKSkpNGvWzFEqEyysheGFhISEs4oFZC0jat/KikbNmjVtedYiEh0dTXh4+FnviQgVKlRwlMgEEysYXkhOTiY1NfWs905/KzO+t3z5cpYvX+46RrEUGxtLTEwMERERiAjlypUjNDSUl19+mVGjRtnU6CZPdknKC3Xq1CE0NJT09PQz74WHh9vkeEWkb9++gI3DKAqhoaEkJiaSkJBASkoKUVFR3H777Tz33HP06NGDL774gsmTJ3P11Ve7jmoCkA3cuwBVpWvXrowdO5YyZcqQlpZmd5YUsZ07dwJQo0YNx0lKDlVl/PjxvPzyy4SFhTFx4kSb6baEymvgnhWMPKgqvXv35u233+all17ivvvuY/369URFRdm966ZY2rp1K48//jj//ve/efbZZxkxYgQRERGuYxk/yqtgoKrF8nH77bdrYQ0cOFAB7dKli2ZmZhb6eMY7ixcv1sWLF7uOUWKdPHlSX3vtNRURvf766zUpKcl1JONHZM0enuvfVev0Po+RI0fy+uuv06FDB0aPHm1TKfjRwIEDGThwoOsYJVaZMmUYPHgwS5cu5dSpU9SvX58BAwac1YdnSia7JJWLiRMn8txzz9GqVStmzZplk7X52d69WUu/R0ZGXmBLU9R+++03unbtyvTp07nrrruYNm0a1atXdx3LFCGb3jwfZsyYQadOnWjatCnTp0+3YuFAZGSkFYsAcemll/LJJ5/w8ccf891331GnTh2mTp1Kcf2iafJmBSObefPm0aFDB+655x7mzJlD2bJlXUcqkRYuXMjChQtdxzDZPP7442zYsIGoqCiefPJJ2rZty+HDh13HMn5mBcNj8eLFPPbYY9StW5cFCxYQFhbmOlKJNWTIEIYMGeI6hsnhmmuuYenSpbz11lvMnTuX2267jaVLl7qOZfzICgawcuVKWrZsyc0330xCQgLly5d3HalEmzlzJjNnznQdw+QiNDSU3r17s3r1asLCwmjUqBGvvvoqJ0+edB3N+EGJ7/Ret24dDRs2pHLlyqxYsYIrrrjCD+mMCX6pqam88sorjBs3jqioKKZOncru3bttnY0gZwP3zmPjxo3ce++9XHzxxaxcuZKqVav6KZ3Jy4IFCwBo3ry54yTGGwsWLOCZZ57h119/pUyZMrbORpCzu6Ry8f3339O4cWPKlSvHl19+acUigAwfPpzhw4e7jmG81Lx5c959911CQkLOWmfjm2++sRmdi5kSec/o7t27adSoEZmZmSxbtszmLAows2fPdh3B5NMPP/yQ6zobQ4cOpWbNmtxwww2OkhlfKnEtjJ9++olGjRrx+++/s2jRIm666SbXkUwOlSpVolKlSq5jmHzIbZ2N0NBQVq5cyY033si9997L1KlTOXbsmKOExhdKVME4ePAgTZo04cCBAyQkJNj05AFq7ty5zJ0713UMkw8519mIiIigQYMG7N69m7feeosff/yRJ598kquuuorOnTuzZs0aG/wXhEpMp/eRI0do2LAhmzdvJiEhgQYNGrgLZ/J0+v8bWw8juGRkZJy1zkb2u6RUlZUrVzJp0iQ+/fRTjh8/Tq1atejYsSPt27fn8ssvd5zenFbi75JKTU3l/vvvZ+3atXz++ec8+OCDjtOZvBw5cgSASy65xHESUxSOHDnCrFmzmDRpEmvWrKF06dK0aNGCjh07cv/999tdVY6V6IJx4sQJmjVrxtKlS/nHP/5B69atXUczxnhs3LiRuLg4pk2bxsGDB6lSpQpPPfUUzzzzjN2M4kiJLRirV6+mdevWLFiwgClTpvDEE0+4jmW8MGvWLADatGnjOInxl1OnTjF//nzi4uJITEwkMzOTBg0a0LFjR1q3bk2ZMmVISEiwQYF+UCILRuXKlfX6669n5cqVjBkzhueff951JOMl68Mo2fbt28eUKVOIi4tj586dXHzxxZQvX55Dhw5x4sQJGxRYxEpkwRARBbj++uvZunWr/WIFkdO3XtoEkCVbZmYmy5cvp3///ud8eShVqhSdO3fmiSee4NZbb7XfFR8q0QUjIiKCGTNm0KxZM9eRjDEFMGDAAPr27Xve23BDQkKoWbMmUVFRZz1sXriCyatgFPuR3qmpqaSkpFjBCCIff/wxAO3bt3ecxASC04MCjx49eua9iIgIRowYwWWXXUZKSgopKSl8/fXXzJgx48w2V1999TlF5Prrryck5P+Gn52+Fdj6RrxjLQwTcKwPw2SXkZHBAw88QFJSEqmpqXn2YRw6dIj169efKSIpKSls3rz5zHrk4eHh3HbbbURFRXHbbbcxadIktmzZwrFjx6xvxKPEXpKKiIiwX4AglJaWBkDp0qUdJzGBIq9BgRdy8uRJNm/efFYRSUlJ4X//+98524aGhnL33XcTFRXF5ZdfnuujQoUKiIhXeYui5VLUraKAKRgi0hQYBYQCE1V1SI7PywJTgduBX4E2qvqD57PeQEcgA3hRVRPzOleVKlV0/Pjx1sQ0xpxDVenZsycjR448p2+kYsWKpKen51pQIKvDvVKlSuctKBUrVuSdd95h69atHD9+nLCwMOrVq8eiRYsoVapwvQD5aW0V5NgJCQk0b978R1Wtkts2fisYIhIK/AdoAuwD1gLtVHVztm3+Atymql1EpC3wiKq2EZFbgBlAPaAysAS4UVUzznc+bxdQMoFn8uTJADz11FNOc5jiLT4+nnbt2p3TN3L6EvbJkyf55ZdfvH54s8b5RRddxEUXXUS5cuW8/pn9+Y4dO5g6depZKxyWLVuWnj17cscddxAaGprno1SpUrm+D/DEE0+QkpJCamoqqpprE8qfBeNO4O+q+oDndW8AVR2cbZtEzzarRaQUcAC4HHgt+7bZtzvf+axgBC/rwzD+4Otv62lpaRw8eJCBAwcyduzYc1ouDRo0oG7dupw4cYLjx4/n+TPne6dOnfLV/2yvnK9g+PMuqSrA3myv9wEx59tGVdNF5Ahwmef9b3Lse06TSUQ6AZ08L0+KyEbfRPebSsBB1yHyoUjzXug6cQEF239jCL7MwZb3EqDS0aNHD3755ZdHCnvZyHO8Gpw9G3jmsmXLdi5btuxIURwb2AkU5thXk3X1Jk/F6rZaVZ0ATAAQkXXn67gJVMGWOdjygmX2h2DLC5bZW/5cD2M/EJntdVXPe7lu47kkdQlZnd/e7GuMMaYI+bNgrAVuEJHqIlIGaAvMz7HNfOBJz/NHga8060LgfKCtiJQVkerADcAaP+U2xhiDHy9JefokugGJZN1WG6eqm0SkP7BOVecDk4BpIrIdOERWUcGz3T+AzUA60DWvO6Q8JhTV/5YiFGyZgy0vWGZ/CLa8YJm9UmwH7hljjPGtErWmtzHGmIKzgmGMMcYrxbJgiEhTEdkmIttF5DXXefIiIpEislRENovIJhHp7jqTt0QkVESSRSTedRZviMilIjJbRLaKyBbPYNKAJSIveX4nNorIDBEp5zpTTiISJyL/zT7mSUQqishiEfne87OCy4w5nSfzUM/vxQYR+UxELnWZMbvc8mb7rKeIqIhU8keWYlcwPFOQjAZigVuAdp6pRQJVOtBTVW8B7gC6Bnje7LoDW1yHyIdRwEJVvQmoQwBnF5EqwItAXVW9lawbRdq6TZWryUDTHO+9BnypqjcAX3peB5LJnJt5MXCrqt5G1hRGvf0dKg+TOTcvIhIJ3A/s8VeQYlcwyJpvaruq7lTVU8BMoKXjTOelqj+p6r89z38n649YrhN/BRIRqQo8BEx0ncUbInIJcA9Zd+KhqqdU9Te3qS6oFHCRZ0xSGPCj4zznUNUVZN3RmF1LYIrn+RTgYb+GuoDcMqvqIlVN97z8hqyxXgHhPP+NAUYArwJ+u3OpOBaM3KYgCfg/wAAici0QDSS5TeKVkWT9sma6DuKl6sAvwEeey2gTRSTcdajzUdX9wDCyvj3+BBxR1UVuU3ntSlX9yfP8AHClyzAF8AyQ4DpEXkSkJbBfVdf787zFsWAEJRGJAOYAPVQ193mVA4SINAP+q6rfus6SD6WA/weMVdVoIJXAu1Ryhue6f0uyCl1lIFxEgm4JQs/A26C5d19E/kbWZeJPXGc5HxEJA/oAb/r73MWxYATdNCIiUpqsYvGJqs51nccLdwEtROQHsi75NRSRj91GuqB9wD5VPd16m01WAQlUjYFdqvqLqqYBc4H6jjN562cRuRrA8/O/jvN4RUSeApoBj2tgD1C7jqwvEus9/warAv8WkauK+sTFsWB4MwVJwJCsKVknAVtU9V3Xebyhqr1VtaqqXkvWf9+vVDWgv/2q6gFgr4jU9LzViKyZAwLVHuAOEQnz/I40IoA76XPIPsXPk8A8h1m84lnc7VWghaoec50nL6r6napeoarXev4N7gP+n+d3vEgVu4Lh6bg6PQXJFuAfqrrJbao83QV0IOtbeorn8aDrUMXUC8AnIrIBiALecpznvDwtodnAv4HvyPq3GnDTV4jIDGA1UFNE9olIR2AI0EREvierpTQkr2P423kyfwCUBxZ7/g2Ocxoym/PkdZMlsFtexhhjAkWxa2EYY4wpGlYwjDHGeMUKhjHGGK9YwTDGGOMVKxjGGGO8YgXDGGOMV6xgGGOM8YoVDGP8SETmishAEVkhIntEpLHrTMZ4ywqGMf5VG/hNVe8haz2Rxx3nMcZrVjCM8RPPLKOXkLWOAUBpINDX5DDmDCsYxvjPLcC3qprheX0bcM6ym8YEKisYxvhPbSAl2+vbgA2OshiTb1YwjPGfnAXjVqyFYYKIzVZrjDHGK9bCMMYY4xUrGMYYY7xiBcMYY4xXrGAYY4zxihUMY4wxXrGCYYwxxitWMIwxxnjl/wNA2k+1HlutzgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mu = 5. # mean number of events expected\n", "nmax = 15 # maximum observed number of events to consider\n", "x = np.arange(nmax+1) # variable for number of events\n", "y = stats.poisson.pmf(x,mu) # generate Poisson probability function\n", "xmin,xmax = 0.,float(nmax)\n", "ymin,ymax = 0.,0.25\n", "plt.plot(x,y,marker='o',markersize=5,color='black')\n", "plt.plot([mu,mu],[ymin,ymax],color='black',linestyle='dotted')\n", "plt.xlabel(r'$n$')\n", "plt.ylabel(r'$P_{\\rm Poisson}(n)$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gaussian distribution\n", "\n", "The **Gaussian (or \"normal\") probability distribution** for a variable $x$, with mean $\\mu$ and standard deviation $\\sigma$, is:\n", "\n", "$P_{\\rm Gaussian}(x) = \\frac{1}{\\sigma \\sqrt{2\\pi}} e^{-\\frac{1}{2} \\left( \\frac{x-\\mu}{\\sigma} \\right)^2}$\n", "\n", "The Gaussian distribution is an ubiquitous and important probability distribution because:\n", "\n", "- It is the **high-$N$ limit** for the Binomial and Poisson distributions\n", "- The **central limit theorem** says that if we average together variables drawn many times from any probability distribution, the resulting average will follow a Gaussian distribution (see example below!)\n", "\n", "Here is an example Gaussian distribution for $\\mu = 2$ and $\\sigma = 0.5$, marking in intervals of standard deviations:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 1.0)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hUVd7A8e8hCSgoXTpIiyhKWUGwB1GkKEUWBVSsEIMgooArawNBVmUXcaWFIiAo8MqyUpTQS1AWk9UICRLAQCAiIkKooYSc948Z2BiSyWTmztx7Jr/P8+QxIXduvpwYTu49M/cqrTVCCCGEN0rYHSCEEMIcMmkIIYTwmkwaQgghvCaThhBCCK/JpCGEEMJrMmkIIYTwWtAmDaXUx0qpQ0qp5AI+r5RS/1RK7VZKbVVK3RysNiGEEN4J5pHGLKCDh893BCLdb9HA5CA0CSGEKIKgTRpa643AEQ+bdAU+0S7/AcorpaoHp04IIYQ3wu0OyKUmsD/XxxnuP/sl74ZKqWhcRyOUKVOmxfXXXx+UQCGECAX//e9/D2utr/HlsU6aNLymtZ4KTAVo2bKlTkxMtLnIs/37XXNh7dq1bS4pmAmNIJ1Wk05rmdKplEr39bFOmjR+BnKPdC33nxmvT58+AKxfv97eEA9MaATptJp0WsuUTn84adJYAgxUSs0HWgPHtNaXnZoy0euvv253QqFMaATptJp0WsuUTn+oYF3lVik1D2gDVAZ+Bd4CIgC01lOUUgqYgOsZVqeBp7XWhZ53MuH0lBBCOIlS6r9a65a+PDZoRxpa696FfF4DA4KUE1RpaWkA1K9f3+aSgpnQCNJpNem0limd/gjakUagmHCk0aZNG8DZ5zlNaATptJp0WsuUTiOONIqzkSNH2p1QKBMaQTqtJp3WMqXTH3KkIYQQxYw/RxpywcIgSE1NJTU11e4Mj0xoBOm0mnRay5ROf8iRRhCYcJ7ThEaQTqtJp7VM6ZQ1DYcbM2aM3QmFMqERpNNq0mktUzr9IUcaQghRzMiahsMlJyeTnJzvbUQcw4RGkE6rSae1TOn0hxxpBIEJ5zlNaATptJp0WsuUTlnTcLixY8fanVAoExpBOq0mndYypdMfcqQhhBDFjKxpOFxSUhJJSUl2Z3hkQiNIp9Wk01qmdPpDjjSCwITznCY0gnRaTTqtZUqnrGk43Pjx4+1OKJQJjSCdVpNOa5nS6Q850hBCiGJG1jQcLiEhgYSEBLszPDKhEaTTatJpLVM6/SFHGkFgwnlOExpBOq0mndYypVPWNBxuwoQJdicUyoRGkE6rSae1TOn0hxxpCCFEMSNrGg73zTff8M0339id4ZEJjSCdVpNOa5nS6Q850ggCE85zmtAI0mk16bSWKZ3+HGnIpBEEF+/k1ahRI5tLCmZCI0in1aTTWqZ0yqTh8ElDCCGcRNY0HG7Dhg1s2LDB7gyPTGgE6bSadFrLlE5/yJFGEJhwntOERpBOq0mntUzplNNTDp800tLSAKhfv77NJQUzoRGk02rSaS1TOmXScPikIYQQTiJrGg63evVqVq9ebXeGRyY0gnRaTTqtZUqnP+RIIwhMOM9pQiNIp9Wk01qmdMrpKYdPGvv37wegdu3aNpcUzIRGkE6rSae1TOmUScPhk4YQQjiJrGk4XFxcHHFxcXZneGRCI0in1aTTWqZ0+kOONILAhPOcJjSCdFpNOq1lSqecnnL4pHHw4EEAqlWrZnNJwUxoBOm0mnRay5ROmTQcPmkIIYSTGLOmoZTqoJRKVUrtVkq9ms/n6yil1imlvldKbVVKdQpmX6AsXbqUpUuX2p3hkQmNIJ1Wk05rmdLpj6AdaSilwoCdQDsgA0gAemutt+faZirwvdZ6slKqMfCV1rqup/2acKRhwnlOExpBOq0mndYypdOUe4S3AnZrrdMAlFLzga7A9lzbaKCs+/1ywIEg9gXMwoUL7U4olAmN4PzO06dPs2PHDnr27Mm5c+f4z3/+ww033EC5cuXsTsuX08fzIul0jmBOGjWB/bk+zgBa59lmBLBSKfUCUAa4L78dKaWigWiAOnXqWB5qtcqVK9udUCgTGsGZndnZ2Xz66af84x//IDk5mfyO3uvVq8cLL7zAc889R+nSpW2ozJ8TxzM/0ukcwTw91QPooLXu6/64D9Baaz0w1zYvu5v+oZS6DZgB3KS1zilovyacnlq0aBEA3bt3t7mkYCY0grM6s7Oz+eSTT3jnnXdIS0ujWbNmPPTQQ9x4440cOHCAsLAwateuzfbt21mxYgXr16+nSpUqDB06lAEDBjhi8nDSeHoindby5/QUWuugvAG3AStyfTwcGJ5nmxSgdq6P04AqnvbbokUL7XRRUVE6KirK7gyPTGjU2jmdv//+u27Xrp0GdMuWLfWSJUt0Tk7Opc/n1xkfH6/vv/9+DeimTZvqtLS0IFdfzinjWRjptBaQqH39t9zXBxb5C7lOhaUB9YCSwA/AjXm2WQ485X7/BlxrGsrTfk2YNDIzM3VmZqbdGR6Z0Ki1MzqTk5N1/fr1dcmSJfX06dP/MFlc5Knzq6++0uXLl9eVKlXSa9euDXSuR04YT29Ip7WMmDRcnXTC9Qyqn4DX3H/2NtDF/X5j4Gv3hJIE3F/YPk2YNEToWLZsmb7qqqt0tWrV9ObNm33ez65du/QNN9ygw8LC9OTJky0sFKJw/kwa8uK+IFiwYAEAPXv2tLmkYCY0gr2d69evp3379jRp0oTFixdTs2bNArf1pvP48eM89thjLFu2jFmzZvHkk09a3lwY+b5by5ROeUW4wycNE567bUIj2Ne5detW7rrrLmrWrMmmTZuoWLGix+297Tx37hwPPPAA69atY+nSpXTs2NGiYu/I991apnTKpOHwSeP06dMAjni2TEFMaAR7OtPT07n99ttRSvHNN9949TTvonQeP36cNm3akJqayrp162jVqpXfzd6S77u1TOmUScPhk4Yw14kTJ2jdujUHDhxg06ZN3HTTTQH5OgcPHuT222/nxIkTJCYmcu211wbk6wgBBl17qriaO3cuc+fOtTvDIxMaIfidgwYNIjU1lX//+99FmjCK2lmtWjXi4uI4e/Ysffr04cKFC77kFpl8361lSqdffF1Bd8qbCc+eMuG52yY0ah3czvnz52tAv/HGG0V+rK+dc+bM0YAeNWpUkR/rC/m+W8uUTuTZU84+PXX+/HkAIiIibC4pmAmNELzO9PR0mjVrxg033EB8fDzh4UW74o4/nY899hgLFixg06ZN3HrrrUV+fFHI991apnTKmobDJw1hlgsXLtCmTRt++OEHkpKSqF+/flC//rFjx2jevDlKKZKSkihbtmzhDxKiCGRNw+FmzZrFrFmz7M7wyIRGCE7nBx98wKZNm5g0aZLPE4Y/neXKlWPu3Lmkp6czdOhQn/bhLfm+W8uUTn/IkUYQmPDcbRMaIfCdGRkZXH/99bRt25YlS5b4vB8rOocMGcIHH3zA5s2bad067wWhrSHfd2uZ0imnpxw+aQhzPPLIIyxdupTt27dTr149W1tOnDjB9ddfT7Vq1fj2228JCwuztUeEDjk9JYQFVq5cyeeff85rr71m+4QBcPXVVzNu3Di+++47YmNj7c4RApAjjaCYNm0aAP369bO5pGAmNELgOs+ePUvTpk3Jyclh27ZtXHHFFX7tz6pOrTX3338/iYmJpKamUqVKFb/2l1dx/75bzZROOdJwuAULFly6kJlTmdAIgescN24cO3fuZMKECX5PGGBdp1KKCRMmcOrUKV555RW/95dXcf++W82UTn/IkYYo9g4dOkSDBg1o167dpTuvOc2rr77Ke++9x/fff0/z5s3tzhGGkyMNIfwwZswYsrKy+Nvf/mZ3SoFeffVVKlSowGuvvWZ3iijmZNIIgkmTJjFp0iS7MzwyoRGs79y7dy+TJ0/m6aefplGjRpbt1+rO8uXLM3z4cL766is2btxo2X6L6/c9UEzp9IecngqCi/dIWL58uc0lBTOhEazvfOqpp5g/fz67d++mVq1aluwTAjOeWVlZNGzYkLp167Jp0yaUUn7vs7h+3wPFlE55nYbDJw3hTMnJyTRt2pQhQ4YwduxYu3O8Mm3aNKKjo1myZAmdO3e2O0cYSiYNmTSED7p168a6detIS0ujUqVKdud4JTs7m8aNG1OqVCmSkpLkBX/CJ7IQ7nAffvghH374od0ZHpnQCNZ1JiQksHjxYoYNGxaQCSNQ4xkeHs7o0aNJTk625Kmdxe37HmimdPpDJo0gWLNmDWvWrLE7wyMTGsG6znfeeYfy5cszaNAgC6ouF8jx7NGjB40bN2bMmDHk5OT4ta/i9n0PNFM6/SGnp0Sxs23bNpo2bcpbb73FiBEj7M7xyaeffsrjjz/Ov//9b7p162Z3jjCMrGnIpCGK4NFHH2Xp0qWkp6dTsWJFu3N8kp2dTaNGjahYsSLffvutJc+kEsWHrGk43N///nf+/ve/253hkQmN4H/nrl27WLBgAf379w/ohBHo8QwPD2f48OEkJiaycuVKn/dTXL7vwWJKpz+Kdg9L4ZPNmzfbnVAoExrB/853332XiIgIXn75ZYuK8heM8XziiScYOXIk77zzDu3bt/dpH8Xl+x4spnT6Q05PiWJj3759NGjQgJiYGD766CO7cyzxz3/+kxdffJENGzZw9913250jDCGnp4Twwrhx4wAYNmyYzSXW6du3L1WqVOHdd9+1O0UUEzJpBMG7777r+B9qExrB987MzExmzJhBr169qFOnTgDK/ihY41m6dGkGDhzI8uXL+fHHH4v8+FD/vgebKZ3+kDWNIEhKSrI7oVAmNILvndOmTePkyZMBX8u4KJjjGRMTw5gxY/jggw+YOnVqkR4b6t/3YDOl0x+ypiFC3vnz56lfvz6RkZGsXbvW7pyAiImJYdasWezbt8/yu/uJ0CNrGkJ48Pnnn5ORkcGQIUPsTgmYwYMHc/bsWSZPnmx3ighxcqQRBKNGjQLgjTfesLmkYCY0QtE7tda0bNmSU6dOsX37dkqUCM7vSXaMZ+fOndmyZQvp6elceeWVXj0mVL/vdjGl058jDVnTCILU1FS7EwplQiMUvXPjxo189913xMbGBm3CAHvG8+WXX6Zt27Z8+umn9O3b16vHhOr33S6mdPpDjjRESOvatSvffPMN+/bt8/q3b1NprWnRogVnzpwhJSVFLi0iCmTEmoZSqoNSKlUptVsp9WoB2zyilNqulEpRSn0WrDYRmtLS0li6dCkxMTEhP2EAKKUYPHgwP/74Y8hfaVXYp8iThlKqjFKqSHd+cW8/EegINAZ6K6Ua59kmEhgO3KG1vhEYXNQ2p3rzzTd588037c7wyIRGKFrnpEmTKFGiBDExMQGuupxd49mzZ0+uueYar1/xHorfdzuZ0umPQtc0lFIlgF7AY8AtwFmglFLqMPAlEKu13l3IbloBu7XWae59zge6AttzbdMPmKi1PgqgtT5UxL+LY+3fv9/uhEKZ0Ajed54+fZoZM2bQvXt3atasGeCqy9k1nqVKlSI6OpoxY8awd+9e6tat63H7UPu+282UTn8UuqahlNoArAYWA8la6xz3n1cE7gEeBf6ttZ7rYR89gA5a677uj/sArbXWA3Nt8wWwE7gDCANGaK3jCthfNBANUKdOnRbp6ene/W1FsTF9+nT69evHxo0bueuuu+zOCaqMjAzq1q3LkCFDeO+99+zOEQ4U0PtpKKUitNbn/dnGy0ljGXAeeASoBWwEmmitMz19bVkIF3lprWnevDngeoVucVwQfvjhh1m7di0ZGRnFYj1HFE1AF8IvTgZKqQ9VAT99hU0qwM9A7Vwf13L/WW4ZwBKt9Xmt9R5cRx2RhfWZYPjw4QwfPtzuDI9MaATvOjdt2sTWrVt54YUXbJsw7B7PgQMHcuTIEebNm+dxO7s7vSWdzlGUhfATwBKlVBkApVR7pdTXXj42AYhUStVTSpXEtUayJM82XwBt3PuuDFwHpBWhz7F+//13fv/9d7szPDKhEbzr/Oijj6hQoQKPPvpokKouZ/d43n333TRp0oSPPvoIT2cT7O70lnQ6iNba6zdc6xcJwNfACuCuIjy2E66jh5+A19x/9jbQxf2+AsbhWhzfBvTyZr8tWrTQQlyUkZGhw8LC9NChQ+1OsV1sbKwGdHx8vN0pwmGARF2Ef/tzv3n94j6l1L3A6+5/3Ku7/7G3/eWPsqYhchs5ciQjRoxg9+7dNGjQwO4cW506dYqaNWvy4IMPMndugc9TEcVQsF7c9xrwhta6DdADWKCUauvLFy1uhg4dytChQ+3O8MiERvDcmZ2dzbRp02jfvr3tE4YTxrNMmTI88cQTfP755/z222/5buOETm9Ip3N4PWlordtqrTe539+G64V6owMVFkqysrLIysqyO8MjExrBc+fSpUv5+eef6d+/f5CrLueU8YyJieHcuXPMnDkz3887pbMw0ukc3jzlVukCNlJKXam1zvK0TaDJ6SlxUfv27dm+fTt79uwhPFyuxXlRVFQUGRkZ7Nq1K6gXbRTOFejTU+uUUi8opf5wj0z3s6BuU0rNBp705YsLYZXdu3ezcuVK+vXrJxNGHv379yctLY1Vq1bZnSJCgDeTRgfgAjBPKXXAfUHBNGAX0BsYr7WeFcBG4w0ePJjBg519KS0TGqHgztjYWMLCwry+JHigOWk8H3roIa655pp8b9DkpE5PpNM5Cv2VTGt9BpgETFJKRQCVgSxdyCu1hQiWM2fOMHPmTLp160aNGjXsznGcUqVK8eyzz/L++++TkZFBrVq17E4SBpP7aQjjzZ07lz59+rB69Wruvfdeu3Mcac+ePTRo0IA33niDkSNH2p0jbBbQa08V8AX/CtyE6zUbWmtt20tvZdIQd911FwcPHiQ1NVUWej3o1KkTP/zwA+np6bLuU8zZcRMmpbV+VGvd284JwxQDBgxgwIABdmd4ZEIjXN6ZkpLCpk2biI6OdtSE4cTxjI6O5sCBA3z55ZeX/syJnfmRTufw9deNBkqph4FTAFrrr6xLCj0mXGXUhEa4vHPatGmULFmSp556yp6gAjhxPB988EFq1KjB1KlT6dq1K+DMzvxIp3P4enrqSeDSA7XWn1gZVRRyeqr4ysrKokaNGnTo0KHQq7kKlzfffJPRo0ezZ88err32WrtzhE3sOD2VCrTGdROmNj7uQwi/LFy4kMzMTKKjo+1OMcazzz4LuG5SJYQvfJ00+gJHgbcIkcuXB1J0dLTj/2EzoRH+2BkbG0tkZCRt2rSxNyofTh3Pa6+9lo4dOzJjxgyys7Md25mXdDqHr2savwJXADlANetyQlOlSpXsTiiUCY3wv86UlBS+/vprxo4d68g78zl5PKOjo+nWrRvLli1zdGdu0ukcvq5pNAbOAYOAtVrrL6wO85asaRRPL774IlOmTOHnn3+mcuXKducYJTs7m2uvvZamTZuyfPlyu3OEDYK6pqGUugOoi+vOeitwTR5CBE1WVhaffPIJ3bt3lwnDB+Hh4Tz77LOsWLGCvXv32p0jDOPLmkZF4Br3W2X3m/Dg6aef5umnn7Y7wyMTGsHV2b59ezIzM3nuuefszimQ08fz4jW6unbt6ujOi5w+nheZ0umPIq9paK2XKqUitNbnlVINgBC/Ia7/ateubXdCoUxoBFfnmjVriIyMJCoqyu6cAjl9POvUqUOHDh2Ij4+nc+fOducUyunjeZEpnf7wdU1jLDAe+DtwQWv9uNVh3pI1jeLlxx9/pHHjxrz//vsMGzbM7hyjffHFFzz00EMsXryYLl262J0jgsiO12mUBboCfwMO+LgPIYps2rRpRERE8OSTcgsXfz3wwANUr16dqVOn2p0iDOLrpPE1rrWNU8Ah63JC0+OPP87jj9t2MOYVExrPnDnDxIkTqV69OlWqVLE7xyMTxjMiIoJKlSrx5Zdfsm/fPrtzPDJhPMGcTn/4+jqNJsCHwFhcN2gSHjRq1MjuhEKZ0Lho0SLOnTtH27Zt7U4plAnjCdCuXTuSk5P5+OOPGTFihN05BTJlPE3p9IevaxqxwA/AJuBxrfUrVod5S9Y0io82bdqQkZHBzp07HXVFW9N16NCBlJQUubd6MWLHmsZ6oKLWeiuu274KEVCpqals2LCBfv36yYRhsejoaDIyMoiLi7M7RRjAp58+rfU8rfVo9/vTrE0KPb169aJXr152Z3jk9MZp06YRHh7O5s2bHd15kdPH86JevXoxf/58qlat6ugFcZPG04ROf/h0LKqUmofr0uhXAeW01s59wrwDNG/e3O6EQjm58cyZM8yaNYtu3brRokULu3O84uTxzO1iZ8OGDXnvvffYv3+/I19rYNp4hjK/7xGulBqstR5vUU+RyZpG6Js3bx6PPvooK1eupF27dnbnhKS0tDQaNGjAiBEjeOutt+zOEQEW9DUNpVQn91tX4GZf9iGEt2JjY6lfvz733nuv3Skhq379+tx///1Mnz6dCxfkCZGiYL6uKF687tQVgG3PnDLFn//8Z/785z/bneGRUxvzLoA7tTMvEzsvLog78cq3Jo5nqPL1+XWpQB+gNNAeeMayohB022232Z1QKKc2Tp06lfDw8EsXgXNqZ14mdnbp0uXSgviDDz5oY9XlTBzPUOXr6zSmAweBqcCTWutRVod5S9Y0QteZM2eoVasW99xzD59//rndOcXCX//6V9577z327t3ryAVxYQ07XqeR+859VX3chxAeLVq0iN9//93Rl0APNf369SMnJ4cZM2bYnSIcyt87970ArJM793l28QqiS5YssbmkYE5svPvuuzlw4MAfXgHuxM78mNzZoUMHkpOT2bt3r2NeIW7yeDqRP0caRf4/Qin1FDAHmIBrMfw6X75wcWLCs36c1piSkkJ8fDzvv//+H14B7rTOgpjcGRMTw0MPPcSXX35J165dbai6nMnjGWqKfKShlPpIa/2CUuoDrfVLSqnxWuvBAeorlAlHGqLoBg0aRGxsLBkZGVxzzTV25xQr2dnZ1K1blyZNmjjymVTCf8Fe07g4y7zs/m+Ytw9USnVQSqUqpXYrpV71sN2flVJaKeXTX0qY7dSpU8yePZuHH35YJgwbhIeH069fP1asWEFaWprdOcJhfJo0lFLtgQpKqQ64FsMLpZQKAyYCHYHGQG/32kje7a4GXgS2+NDmSB07dqRjx452Z3jkpMb58+dz/PhxYmJiLvuckzo9Mb2zb9++lChRwjHXozJ9PEOJL6tcQ4F+QBcg2f2xN1oBu7XWaQBKqfm47v63Pc92o4D3gJC5l6cJ92B2UuPkyZO58cYbueOOOy77nJM6PTG9s2bNmnTu3JkZM2YwcuRISpUqFeSyPzJ9PEOK1trjG/AkcBg4AnwCXF3YYwrYTw9geq6P+wAT8mxzM/Av9/vrgZYF7CsaSAQS69Spo0XoSEhI0ICeMGGC3SnF3ooVKzSgP/vsM7tThMWARO3Dv+Naa69OT70BtAOuB9KBMRbMVZdRSpUAxgFDCttWaz1Va91Sa91SznmHlsmTJ1O6dOmQv2WmCe677z4aNGjA5MmT7U4RDuLNpHFca/291vqQ1voNXKeZfPEzkPslprXcf3bR1cBNwHql1F7gVmBJKCyG33fffdx33312Z3jkhMajR49euqJtuXLl8t3GCZ3eCIXOEiVK8NxzzxEfH09ycnKQy/4oFMYzVHizplFdKRUN7AB+BCJ8/FoJQKRSqh6uyaIX8OjFT2qtj+F63QcASqn1wFCttfHPp+3Zs6fdCYVyQuPMmTPJyspiwIABBW7jhE5vhErnM888w5tvvsnEiRNtPeIIlfEMBYW+TsM9YTQBmrr/exWwGtc9wrdqred5/cWU6gSMx/U03Y+11u8opd7GdX5tSZ5t1+PFpCGv0wgNOTk5XHfddVSrVo1NmzbZnSNyefrpp/n888/5+eefCzwCFGYJ6Os03OsHL2ito7TWFYH6wEdAJtCpKF9Ma/2V1vo6rXUDrfU77j97M++E4f7zNqFwlCG8s2LFCn766ScGDhxod4rIY+DAgZdeOyOE33fus5sJRxpt2rQBYP369bZ2eGJ344MPPkhiYiL79u2jZMmSBW5nd6e3Qq3z1ltv5ejRo/z4449/uKxLsITaeNotqNeeEkX31FNP2Z1QKDsb09LS+Oqrr3jjjTc8ThhgxlhC6HUOHDiQPn36sGbNGltuuRtq42kyOdIQths2bBgffPAB6enp1KxZ0+4ckY+zZ89Su3Ztbr/9dr74wraLWguL2HE/DVEE58+f5/z583ZneGRX4+nTp5kxYwbdu3f3asIwYSwh9DpLlSpFv379WLp0Kenp6UEo+6NQG0+TyaQRBO3atbPlkL4o7Gr87LPPOHr0qMen2eZmwlhCaHbGxMSglGLixIkBrrpcKI6nqWRNIwj69u1rd0Kh7GjUWjN+/HiaN2/O3Xff7dVjTBhLCM3O2rVr0717d6ZNm8Zbb71FmTJlAlj2R6E4nqaSNQ1hm9WrV9OuXTtmzpxZLBYQQ8HXX3/NnXfeyaRJk+jfv7/dOcJH/qxpyKQRBKdPnwagdOnSNpcUzI7GBx98kISEBNLT07niiiu8eowJYwmh26m1plWrVpw4cYLt27cH7em3oTqedpGFcIfr1KkTnToV6XWQQRfsxl27dvHll1/Sv39/rycMMGMsIXQ7lVIMHjyY1NRUVq5cGcCyPwrV8TSRrGkEgQmH8cFu/Oc//0nJkiXzvdGSJyaMJYR258MPP8ywYcMYP348HTp0CEDV5UJ5PE0jp6dE0GVmZlKrVi169OjBrFmz7M4RPnjnnXd4/fXXSUlJoXHjy27AKRxOTk853LFjxzh27JjdGR4Fs3H69OmcOnWKF198sciPNWEsIfQ7o6OjKVWqFB9++GEAqi4X6uNpEjnSCAITrkcTrMbz589Tv359GjRo4NPXMmEsoXh09uvXjzlz5rBv3z6qVKlibVgexWE8g0muPeVwgwYNsjuhUMFqnD9/PhkZGcTGxvr0eBPGEopH59ChQ5kxYwYTJkzg7bfftrDqcsVhPE0hRxoiaLTWNGvWjJycHLZt24ZSyu4k4adu3boRHx/Pvh1qtYkAABecSURBVH37gvpiP+EfWdNwuMOHD3P48GG7MzwKRuOKFSvYtm0bw4YN83nCMGEsofh0vvLKKxw5coSPP/7YwqrLFZfxNIEcaQSBCec5g9HYtm1bdu7cSVpaWqGXQC+ICWMJxavzzjvv5Oeff2bXrl2EhwfmjHdxGs9gkDUNhxsyZIjdCYUKdGNiYiLr1q1j7NixPk8YYMZYQvHqHDZsGN26dWPhwoX06tXLgqrLFafxdDo50hBB0bNnT+Li4ti/fz9ly5a1O0dYKCcnh8aNG1O6dGn++9//ylqVAWRNw+EOHjzIwYMH7c7wKJCNu3btYuHChcTExPg9YZgwllC8OkuUKMGwYcP4/vvvA3ZpkeI0nk4nRxpBYMJ5zkA2Pv300yxYsIA9e/ZQtWpVv/ZlwlhC8es8d+4cDRs2pE6dOsTHx1t+tFHcxjPQZE3D4V599VW7EwoVqMY9e/YwZ84cBg4c6PeEAWaMJRS/zpIlS/KXv/yFgQMHsn79eu655x5L9ntRcRtPJ5MjDRFQzz33HLNmzSItLU3u/x3izpw5Q7169bjhhhtYu3at3TnCA1nTcLj9+/ezf/9+uzM8CkRjRkYGM2fO5JlnnrFswjBhLKF4dl5xxRUMGzaMdevW8fXXX1uyz4uK43g6lRxpBIEJ5zkD0Tho0CAmT57M7t27ufbaay3ZpwljCcW389SpU9StW5eWLVuyfPlyS/YJxXc8A0XWNBzu9ddftzuhUFY3Hjx4kGnTpvHEE09YNmGAGWMJxbezTJkyDBkyhOHDh5OQkMAtt9xiyX6L63g6kRxpiIAYNGgQkyZNYseOHTRs2NDuHBFEJ06coG7durRq1crSow1hHVnTcLi0tDTS0tLszvDIysb09HRiY2N55plnLJ8wTBhLKN6dV199Na+++ipxcXHEx8dbss/iPJ5OI0caQWDCeU4rG5999lnmzp3L7t27qV27tt/7y82EsQTpPH36NA0aNKBhw4Zs3LjR79dtFPfxtJqsaTjcyJEj7U4olFWNqampzJ49mxdeeMHyCQPMGEuQztKlS/PGG28wYMAAVqxY4fe9xIv7eDqJHGkIS/Xq1Ytly5aRlpYW8Lu5CWc7d+4cjRo1okKFCiQmJlKihJwNdwpZ03C41NRUUlNT7c7wyIrGpKQkFixYwODBgwM2YZgwliCd4HqV+IgRI/j+++9ZtGiRX/uS8XQOOdIIAhPOc1rR2LFjR/7zn/+wZ88eypcvb01YHiaMJUjnRRcuXKBJkyZcuHCB5ORkIiIifNqPjKe1ZE3D4caMGWN3QqH8bYyLiyMuLo5//OMfAZswwIyxBOm8KCwsjPfee48uXbowZcoUXnjhBZ/2I+PpHEE70lBKdQA+BMKA6Vrrd/N8/mWgL5AN/AY8o7VOL2y/JhxphLrs7GyaNWvG2bNnSUlJoVSpUnYnCQfRWnPfffeRlJTE7t27qVChgt1JxZ7j1zSUUmHARKAj0BjorZRqnGez74GWWuumwELg/WC0BUNycjLJycl2Z3jkT+P06dPZvn07Y8eODfiEYcJYgnTmppRi3LhxHD16lNGjR/u0DxlP5wjKkYZS6jZghNa6vfvj4QBa678VsP2fgAla6zsK27cJRxomnOf0tfHYsWNERkbSuHFj1q1bF/C7tpkwliCd+enXrx+zZ88mJSWFyMjIIj1WxtNaJqxp1ARyX/oxA2jtYftngQKvP6CUigaiAerUqWNFX0CNHTvW7oRC+do4ZswYDh8+zLhx44Jym08TxhKkMz+jRo1i/vz5/OUvfynys6lkPJ0jWEcaPYAOWuu+7o/7AK211gPz2fZxYCAQpbU+W9i+TTjSCFU7d+6kSZMm9O7dm1mzZtmdIwwwZswYXnvtNVatWsV9991nd06x5fg1DeBnIPfLg2u5/+wPlFL3Aa8BXbyZMEyRlJREUlKS3RkeFbVRa83zzz/PlVdeybvvvlv4AyxiwliCdBbk5ZdfJjIykueff54zZ854/TgZTwfRWgf8DddpsDSgHlAS+AG4Mc82fwJ+AiKLsu8WLVpop4uKitJRUVF2Z3hU1MZPP/1UA3rixImBi8qHCWOptXR6smrVKg3oESNGeP0YGU9rAYnax3/Pg/mU207AeFxPuf1Ya/2OUuptd/wSpdRqoAnwi/sh+7TWXQrbrwmnpy7+5tG8eXObSwpWlMbMzEyuv/566tSpw+bNmwkLCwt03iUmjCVIZ2EeffRRFi1axLZt27xaFJfxtJY/p6fkFeGiyAYMGMCUKVNISEjg5ptvtjtHGOjgwYM0atSI1q1bs2LFiqA8iUL8jwlrGsVaQkICCQkJdmd45G3jli1bmDx5MgMHDrRlwjBhLEE6C1OtWjXGjBnDqlWrmDdvXqHby3g6hxxpBIEJz932pjErK4ubb76ZkydPkpKSQtmyZYMTl4sJYwnS6Y0LFy5wxx13sGvXLlJSUqhWrVqB28p4WsuE12kUaxMmTLA7oVDeNL7++uvs2LGDlStX2jJhgBljCdLpjbCwMGbPnk3z5s3p168fS5YsKfA0lYync8iRhvBKfHw8UVFRxMTEMGnSJLtzRAgZP348L730EjNnzuSpp56yO6dYkIVwh08a33zzDQC33367zSUF89R48uRJmjVrBsAPP/zAVVddFdS23EwYS5DOosjJyeGee+4hKSmJ5OTkfO/46IROb5jSKZOGwycNE85zemqMiYlh6tSprF+/nrvvvju4YXmYMJYgnUWVlpZG06ZNufXWW1mxYsVlT+N2SmdhTOn0Z9IIyov7Avlmwov7duzYoXfs2GF3hkcFNc6fP18DetiwYTZUXc6EsdRaOn0xbdo0DehRo0Zd9jkndXpiSicmvLgvUEw40jDVrl27aNGiBTfddBMbNmzw+a5rQnhDa02fPn2YN28ea9asufRbu7CevE7D4TZs2MCGDRvszvAob+OZM2d45JFHiIiIYMGCBY6ZMEwYS5BOXyilmDJlCpGRkfTu3Ztff/310uec1OmJKZ3+kCONIDDhPGfexv79+zNlyhSWLVvGAw88YF9YHiaMJUinP7Zu3Urr1q254447Lq1vOLEzP6Z0ypqGw/3000/6p59+sjvDo9yNsbGxGtCvvPKKzVWXM2EstZZOf02fPl0DesiQIVpr53bmZUonsqbh7CMNk6xdu5b27dvTrl07li5dGtSLEQqR28CBA5k4cSLTp0/n2WeftTsnpMgrwh1u9erVAI6+6czq1avJyMjgpZdeolGjRsyfP9+RE4YJYwnSaYXx48eza9cuYmJiyMzMpFmzZo7szM3J42kVOdIIAhPOc95555189913lClThm+//ZZ69erZnZQvE8YSpNMqmZmZ3HbbbezevZubb76ZLVu22J3kkdPH8yI50nC4OXPm2J3g0YkTJzh16hTZ2dl88cUXjp0wwPljeZF0WqN8+fIsW7aMli1bcuDAATIyMqhVq5bdWQVy+nhaQY40irmsrCw6depEfHw8//rXv+jatavdSUJcJjExkbZt21KjRg02btxIlSpV7E4ymrxOw+Hi4uKIi4uzO+My586d4+GHH2bDhg0MGTKEUqVK2Z1UKKeOZV7Saa3Dhw/z5ptvsm/fPu6//36OHj1qd1K+TBlPf8iRRhA48Tzn2bNnL91yMzY2ls8++wxwVmN+nDiW+ZFOa13sHD58OJ07d+ZPf/oTy5cvp2LFivaG5WHKeMoFCx0+aRw8eBDA401mgunkyZN0796dVatWMX78eF588UXHNRZEOq1lYufixYt55JFHuO6661i5ciXVq1e3ue5/TBlPmTQcPmk4yZEjR3jggQf49ttvmTFjhty/QBhpzZo1dO3alapVq7Jq1Srq169vd5JRZE3D4ZYuXcrSpUvtzmDv3r1ERUXx3Xff8a9//esPE4ZTGgsjndYytfPee+9lzZo1ZGZmcuedd+KUXxxNGU9/yJFGEDjhPOf69evp0aMH2dnZLFq0iLZt2/7h805o9IZ0Wsv0zpSUFDp16sShQ4eYPn06jz32WPDjcjFlPOX0lMMnjcOHDwNQuXLloH9trTUTJ05k8ODBXHfddSxevJjIyMjLtrOzsSik01qh0Pnbb7/Ro0cPNm7cyLBhw/jb3/5m29UMTBlPmTQcPmnY5ejRozz//PPMnz+fzp07M3fuXMqWLWt3lhCWO3/+PC+99BITJ06kTZs2zJ49mzp16tid5ViypuFwixYtYtGiRUH9muvWraNZs2YsXLiQ0aNH88UXX3icMOxo9IV0WitUOiMiIpgwYQKzZs0iMTGRpk2bMm/evCAWupgynn7x9fK4Tnkz4dLoUVFROioqKihf6/jx43rw4MFaKaWvu+46nZCQ4NXjgtnoD+m0Vih2/vTTT/r222/XgO7Zs6f+5ZdfAhuXiynjiR+XRrf9H31/30yYNDIzM3VmZmZAv0ZOTo6eP3++rlGjhlZK6eeff16fPHnS68cHo9EK0mmtUO08f/68Hj16tC5ZsqQuW7as/vDDD/X58+cDWOhiynjKpFHMJSQk6LZt22pAt2jRQm/ZssXuJCEcYefOnbp9+/Ya0M2aNdNxcXE6JyfH7izb+TNpyJpGECxYsIAFCxZYvt9t27bx0EMPccstt/DDDz8wadIktmzZQqtWrRzTaDXptFaod0ZGRrJ8+XIWLlxIZmYmHTp0ICoqio0bNwag0pzx9Iuvs41T3kw40rDyPGdOTo5eu3at7tKli1ZK6XLlyulRo0bp48ePO6YxkKTTWsWp8+zZs3rSpEm6Ro0aGtB33XWXXrRokc7OzrYmUpsznsjtXp39lNvTp08DULp0aZ/3cezYMf7v//6PCRMmsHXrVipXrkz//v0ZPHiwJRdts6IxGKTTWsWxMysri9jYWMaPH096ejp169ZlwIABPP74435fM8qU8ZTXaTh80vDV2bNnWb16NXPmzGHx4sWcOXOGJk2aMHjwYHr37s2VV15pd6IQxsrOzmbJkiWMHz+e+Ph4wsLCuP/++3niiSd44IEHuPrqq+1ODBiZNBw+acydOxeAxx9/vNBtDx06xMqVK1m8eDFxcXGcPHmSSpUq0bt3b5544glatmyJUsrWRjtJp7Wk02XHjh3MmTOHOXPmsH//fkqWLMm9995Lly5d6NSpk9cvFDRlPGXScPikUdD1aLTW7Nu3j2+//ZaNGzeybt06UlJSANellbt06UKXLl1o164dJUuWtKXRaaTTWtL5Rzk5OcTHx7N48WIWL15MWloaAPXq1eOee+4hKiqK1q1bExkZSYkSlz+PyJTxlEnD4ZPG+fPnOX36NHv27CElJYXk5GS2bt1KYmIihw4dAlznQO+8807uuece2rZtS8uWLfP9nzKQjeB6Za2TSae1pLNgWmu2b9/O6tWrWbduHRs2bCAzMxOAcuXK0bJlS5o2bcpNN93ETTfdxPXXX3/plLHTx9OYSUMp1QH4EAgDpmut383z+VLAJ0AL4Hegp9Z6r6d9OmHSOHv2LIcOHeLXX3/l119/JSMjg/3797N//3727NnD7t27+eWXXy5tHxERQaNGjWjRogWtWrXilltuoVmzZgE/mhBC+O7ChQts376dhISES2/bt28nKyvr0jaVK1emYcOG1K9fn9q1a1OrVi1q165NtWrVqFq1KlWqVHHEIrkRk4ZSKgzYCbQDMoAEoLfWenuubZ4HmmqtY5RSvYCHtNY9Pe23RYsWevPmzeTk5Fz2duHChUtv2dnZZGdnc+HCBc6fP3/p7dy5c5w9e/bSf8+cOUNWVhZZWVmcPn2aU6dOXXo7fvw4x44d4/jx4xw9epQjR45w5MgRTpw4cVlXWFgY1atXp169euTk5FClShV69uxJkyZNiIyMdNxvIrNmzQJw/E2ZpNNa0umfCxcusGfPHpKTk9m5cydfffUVhw4d4syZM2RkZFw6QsqtdOnSVKxYkYoVK1KhQgXKlStHuXLlKFu2LFdddRVlypShTJkylC5dmiuvvJIrrriCK6+8klKlSlGyZMlL/42IiLj0Fh4eTnh4OGFhYZe9lShR4rK3iIgIIyaN24ARWuv27o+HA2it/5ZrmxXubTYrpcKBg8A12kOkUirgf4Hw8HDKlClz6RtbtmxZKlSocOkbX6lSJapWrXrprXbt2lStWpXw8HDAjPOcJjSCdFpNOq2VuzMnJ4fffvuN/fv3XzoL8euvv3L48OFLv3AeOXKE48ePX/qF9NSpU5w7dy4YqT5PGuFWl3hQE9if6+MMoHVB22its5VSx4BKwOHcGymlooFo94dngeRABF+UnZ3NsWPHOHbsmD+7qayUOlz4ZrYyoRGk02rSaS0TOhv5+sBgThqW0VpPBaYCKKUSfZ0xg8mEThMaQTqtJp3WMqFTKeXzQnAwrz31M1A718e13H+W7zbu01PlcC2ICyGEcIBgThoJQKRSqp5SqiTQC1iSZ5slwJPu93sAaz2tZwghhAiuoJ2ecq9RDARW4HrK7cda6xSl1Nu4Lp61BJgBzFFK7QaO4JpYCjM1YNHWMqHThEaQTqtJp7VM6PS50fgX9wkhhAgeuZ+GEEIIr8mkIYQQwmtGTBpKqQ5KqVSl1G6l1Kv5fL6UUmqB+/NblFJ1g1/pVedTSqnflFJJ7re+NnV+rJQ6pJTK9/UtyuWf7r/HVqXUzQ5sbKOUOpZrLN8MdqO7o7ZSap1SartSKkUp9WI+2zhhPL3ptH1MlVJXKKW+VUr94O4cmc82tv68e9noiJ91d0uYUup7pdSyfD5X9LH09e5NwXrDtWj+E1AfKAn8ADTOs83zwBT3+72ABQ7tfAqY4IAxvRu4GUgu4POdgOWAAm4FtjiwsQ2wzAFjWR242f3+1bgulZP3++6E8fSm0/YxdY/RVe73I4AtwK15trH1593LRkf8rLtbXgY+y+9768tYmnCk0QrYrbVO01qfA+YDXfNs0xWY7X5/IXCvUgG46YRn3nQ6gtZ6I65npxWkK/CJdvkPUF4pVT04dS5eNDqC1voXrfV37vdPAD/iurJBbk4YT286beceo5PuDyPcb3mfrWPrz7uXjY6glKoFPABML2CTIo+lCZNGfpcfyfs/+x8uPwJcvPxIMHnTCfBn9ymKhUqp2vl83gm8/bvY7Tb3KYLlSqkb7Y5xH9r/Cddvnrk5ajw9dIIDxtR9OiUJOASs0loXOJ52/bx70QjO+FkfD7wC5BTw+SKPpQmTRihZCtTVWjcFVvG/GV4U3XfAtVrrZsBHwBd2xiilrgL+BQzWWh+3s8WTQjodMaZa6wta6+a4rhrRSil1kx0dnnjRaPvPulLqQeCQ1vq/Vu7XhEnDlMuPFNqptf5da33W/eF0XPcNcSJvxtxWWuvjF08RaK2/AiKUUpXtaFFKReD6h/hTrfWifDZxxHgW1umkMXU3ZALrgA55PuWEn3eg4EaH/KzfAXRRSu3Fdbq8rVJqbp5tijyWJkwaplx+pNDOPOexu+A6r+xES4An3M/6uRU4prX+pbAHBZNSqtrFc69KqVa4/l8O+j8c7oYZwI9a63EFbGb7eHrT6YQxVUpdo5Qq737/Slz339mRZzNbf969aXTCz7rWerjWupbWui6uf4/Waq3z3ry8yGPp+Kvc6sBdfsSOzkFKqS5AtrvzqWB3Aiil5uF6pkxlpVQG8BauxTy01lOAr3A942c3cBp42oGNPYD+SqlsIAvoZcMvCuD6ba4PsM19jhvgr0CdXK22j6eXnU4Y0+rAbOW6aVsJ4P+01ssc9vPuTaMjftbz4+9YymVEhBBCeM2E01NCCCEcQiYNIYQQXpNJQwghhNdk0hBCCOE1mTSEEEJ4TSYNIYQQXpNJQwghhNdk0hDCYsp134p27vdHK6U+srtJCKs4/hXhQhjoLeBtpVQVXFeT7WJzjxCWkVeECxEASqkNwFVAG/f9K4QICXJ6SgiLKaWa4Lo+0TmZMESokUlDCAu5r276Ka47op1USuW9rLcQRpNJQwiLKKVKA4uAIVrrH4FRuNY3hAgZsqYhhBDCa3KkIYQQwmsyaQghhPCaTBpCCCG8JpOGEEIIr8mkIYQQwmsyaQghhPCaTBpCCCG89v8NpKe5TVESCAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mu = 2. # mean\n", "sig = 0.5 # standard deviation\n", "xmin,xmax,nx = 0.,4.,101 # range of variable\n", "x = np.linspace(xmin,xmax,nx) # values of variable\n", "y = stats.norm.pdf(x,mu,sig) # generate Gaussian probability function\n", "ymin,ymax = 0.,1.\n", "plt.plot(x,y,color='black')\n", "for i in range(-3,4):\n", " x1 = mu + i*sig\n", " plt.plot([x1,x1],[ymin,ymax],color='black',linestyle='dotted')\n", "plt.xlabel(r'$x$')\n", "plt.ylabel(r'$P_{\\rm Gaussian}(x)$')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(ymin,ymax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Gaussian distribution is often used as shorthand for the **\"confidence\" of a statement**. For example, the probability contained within $\\pm 1,2,3$ standard deviations is $(68.27, 95.45, 99.73)\\%$. If a statement is true with \"3$\\sigma$ confidence\", it implies it is true with a probability of $99.73\\%$. The following snippet of code generates this for a general number of standard deviations." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Probability = 0.9721931049730028\n" ] } ], "source": [ "nsigma = 2.2 # number of standard deviations\n", "prob = stats.norm.cdf(nsigma) # 1-tailed probability integrated up to this location\n", "prob = 1. - 2.*(1.-prob) # add in other tail!\n", "print('Probability =',prob)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Frequentist and Bayesian frameworks\n", "\n", "It is important to understand different methods of using probability in science. We describe them here in a practical way.\n", "\n", "First let's introduce the concept of **conditional probability**. $P(A|B)$ means \"the probability of $A$ on the condition that $B$ has occurred\".\n", "\n", "**Frequentist statistics** assign probabilities to a measurement, i.e. they determine $P(data|model)$.\n", "\n", "- We are defining probability by imagining a series of hypothetical experiments, repeatedly sampling the population, which have not actually taken place\n", "- Our philosophy of science is to attempt to \"rule out\" or falsify models, if $P(data|model)$ is too small\n", "\n", "**Bayesian statistics** assign probabilities to a model, i.e. they give us tools for calculating $P(model|data)$.\n", "\n", "- We update the model probabilities in the light of each new dataset\n", "- Our philosophy of science is we do not \"rule out\" models, just determine their relative probabilities\n", "\n", "Bayesian statistics uses Bayes' theorem to derive $P(model|data) = P(data|model) \\, P(model) / P(data)$.\n", "\n", "## Worked example\n", "\n", "The following is a worked example solved by both frequentist and Bayesian methods.\n", "\n", "_I observe 100 galaxies, 30 of which are AGN. What is the best estimate of the AGN fraction and its error?_\n", "\n", "### Solution by frequentist approach\n", "\n", "- Estimate the AGN fraction as $p = N_{AGN}/N_{total} = 30/100 = 0.3$\n", "- There are 2 possible outcomes (\"AGN\" or \"not an AGN\") so the binomial distribution applies\n", "- Estimate the error in $N_{AGN}$ as the standard deviation of the binomial distribution, $\\sqrt{p(1-p)N_{total}} = 4.6$. Hence, the error in $p$ is $4.6/100 = 0.046$\n", "- Answer: $p = 0.3 \\pm 0.046$\n", "\n", "### Solution by Bayesian approach\n", "\n", "- Use Bayes' theorem, $P(p|D) \\propto P(D|p) \\, P(p)$\n", "- $P(p|D)$ is the probability distribution of $p$ given the data $D$, which we are aiming to determine\n", "- $P(D|p)$ is the probability of the data for a given value of $p$, which is given by the binomial distribution as $P_{\\rm Binomial}(n=30|N=100,p)$\n", "- $P(p)$ is the prior in $p$, which we take as a uniform distribution between $p=0$ and $p=1$.\n", "\n", "Here is the resulting calculation of $P(p|D)$:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$P(p|D)$')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXiU5d328e8vCQFCAUGQRUEUEQQEl7TFBaEiCKgsgooJCEpZWmx5bXFBacEXC/rW2qK4IZuyyxZQQIQiKCogCMrqgyxGkCVssieQXO8fCT6RRSYkM9dk5vwcxxyZfU5vmZy5rnsz5xwiIiKBiPEdQERECg+VhoiIBEylISIiAVNpiIhIwFQaIiISMJWGiIgELGSlYWYjzWy3ma3JdV9ZM5tnZhtzfpYJVR4REcm7UI40RgPNT7vvKeC/zrkawH9zbouISJiyUO7cZ2bVgPedc3Vzbn8DNHbO7TCzSsBC51zNkAUSEZE8ifP8+RWccztyru8EKpzriWbWHegOUKJEiRtr1aoVgngiIpFjxYoVe5xz5fPzHr5L4yfOOWdm5xz2OOeGAcMAEhMT3fLly0OWTUQkEpjZd/l9D99bT+3KmZYi5+duz3lEROQX+C6NmUDnnOudgRkes4iIyHmEcpPbCcDnQE0z22ZmXYHngaZmthG4I+e2iIiEqZCt03DOPXiOh5qEKoOIiOSP7+kpEREpRFQaIiISMJWGiIgETKUhIiIBU2mIiEjAVBoiIhIwlYaIiARMpSEiIgFTaYiISMBUGiIiEjCVhoiIBEylISIiAVNpiIhIwFQaIiISMJWGiIgETKUhIiIBU2mIiEjAVBoiIhIwlYaIiARMpSEiIgGL8x1ApLDbtWsXH3/8MRs2bCAzM5OyZcty6623ct111xETo7/LJLKoNEQu0KpVq+jfvz+zZs0iMzPzjMerVKlCnz596NatG8WLF/eQUKTg6c8gkTxKT0+nd+/e3HjjjXz22Wf06dOHL774gqNHj3Ly5Em2b9/OO++8Q7Vq1ejduzfXXXcdS5cu9R1bpECoNETyYOfOndx+++28/PLL9OzZk40bN/L888+TmJhI8eLFiY2NpXLlynTq1ImPP/6YuXPncuzYMW655RbefPNN3/FF8k3TUyIB2rFjB40aNWL79u28++673Hfffed9TbNmzVi9ejVJSUn07NmTHTt2MGDAgOCHFQkSlYZIAPbs2UOTJk3YsWMH8+bN4+abbw74taVLl2bGjBl069aNZ599loSEBJ544okgphUJHpWGyHmcOHGC9u3bs2XLFubOnZunwjglLi6OESNGcPz4cZ588kkqVKhA586dg5BWJLhUGiLn8de//pVFixYxduxYbrvttgt+n5iYGN5++23S0tLo0aMHderUITExsQCTigSfVoSL/IIZM2bwyiuv8Nhjj5GcnJzv94uPj2fixIlUqFCBdu3asXfv3gJIKRI6Kg2Rc9i1axfdunXjuuuu4/nnny+w9y1XrhxTp05lx44d/OEPf8A5V2DvLRJsKg2Rc/jjH//IwYMHGTt2LPHx8QX63omJiQwYMIDJkyczceLEAn1vkWBSaYicxaxZs5g2bRr9+/enTp06QfmMJ554ggYNGtCrVy/S0tKC8hkiBU2lIXKao0eP8uijj3LNNdfw17/+NWifExcXx/Dhwzl06BBPP/100D5HpCCpNERO8+9//5utW7fy2muvFfi01Onq1KlD7969GTFihA41IoWCFcaVcImJiW758uW+Y0gE2rt3L1deeSWNGzdmxowZIfnMQ4cOUbNmTSpXrszSpUuJjY0NyedK9DGzFc65fG3nHRYjDTN7zMzWmtkaM5tgZsV8Z5LoNGjQIA4fPsygQYNC9pklS5bkxRdfZMWKFYwaNSpknytyIbyPNMzsUmAxUNs5d8zM3gVmO+dGn+s1GmlIMKSmplKjRg2SkpJC/svbOcctt9xCamoq3377LcWK6e8mKXgRM9Ige8/04mYWByQAP3jOI1FowIABmBnPPvtsyD/bzBg0aBDbt2/ntddeC/nniwTKe2k457YDLwKpwA7gR+fch6c/z8y6m9lyM1uuzROloG3dupV33nmHP/zhD1StWtVLhsaNG9O0aVMGDRrEwYMHvWQQOR/vpWFmZYDWwBVAZaCEmXU8/XnOuWHOuUTnXGL58uVDHVMi3L/+9S9iYmKCuoltIAYNGsTevXv5z3/+4zWHyLl4Lw3gDmCLcy7NOXcCmAbk/TCiIhdo9+7dDB8+nE6dOnHZZZd5zZKYmMi9997Liy++yIEDB7xmETmbcCiNVKCBmSWYmQFNgPWeM0kUeeWVV0hPT+fxxx/3HQWAv/3tbxw6dIhXX33VdxSRM3gvDefcUmAK8CWwmuxMw7yGkqhx6NAhhg4dStu2balVq5bvOABcd911tGjRgiFDhnD06FHfcUR+xntpADjn+jvnajnn6jrnOjnn0n1nkugwcuRIDhw4wJNPPuk7ys/07duXtLQ0Ro4c6TuKyM9430/jQmg/DSkIWVlZ1KpVi3LlyvHZZ5/5jnOGW2+9le+//55vv/2WIkWK+I4jESCS9tMQCbl58+axceNGHn30Ud9Rzuqpp54iNTWVCRMm+I4i8hONNCRq3XPPPXzxxRekpqYG/cCEF8I5R926dSlatCgrVqwgezsRkQunkYbIBdqyZQuzZs2ie/fuYVkYkL2XeO/evVm5ciWLFy/2HUcEUGlIlHr99deJiYmhR48evqP8oo4dO1K2bFmGDBniO4oIoNKQKJSens7IkSNp27Ytl156qe84vyghIYFu3boxffp0tm7d6juOiEpDos+MGTPYu3cv3bp18x0lIL169cLMtLOfhAWVhkSd4cOHU7VqVe644w7fUQJSpUoV2rVrx/Dhw7Wzn3in0pCo8t133zF//nwefvhhYmIKzz//Xr16ceDAAd59913fUSTKFZ5vjUgBOHVypYcffthzkrxp2LAhtWrV4s033/QdRaKcSkOiRmZmJqNGjeKOO+7g8ssv9x0nT8yM7t27s2TJEr7++mvfcSSKqTQkaixYsIDU1FS6du3qO8oF6dy5M0WLFtVoQ7xSaUjUGDNmDBdddBGtW7f2HeWClC1blvvuu4+xY8dy5MgR33EkSqk0JCocOXKEadOm0b59e4oVK+Y7zgXr0aMHBw8eZOLEib6jSJRSaUhUmDlzJkeOHKFjxzPOJFyo3HLLLdSuXVtTVOKNSkOiwtixY6lSpQoNGzb0HSVfzIwePXrwxRdfsHLlSt9xJAqpNCTipaWlMXfuXJKSkgrVvhnn0qlTJ4oWLfrT5sMioVT4v0Ei5zFp0iQyMzML/dTUKWXKlKFNmzaMGzeO9HSd5FJCS6UhEW/s2LHUq1ePunXr+o5SYLp06cK+ffuYNWuW7ygSZVQaEtG+/fZbli5dGjGjjFOaNm1KpUqVGD16tO8oEmVUGhLRxo0bh5nx4IMP+o5SoGJjY+nUqROzZ89m165dvuNIFFFpSMRyzjFx4kQaNWrEZZdd5jtOgevcuTOZmZmMGzfOdxSJIioNiVhr1qxhw4YN3H///b6jBEXt2rX5zW9+w+jRo3HO+Y4jUUKlIRFr8uTJxMTEcO+99/qOEjRdunRh9erVrFq1yncUiRIqDYlIzjkmT57MbbfdRoUKFXzHCZoOHToQHx+vFeISMioNiUhr165lw4YN3Hfffb6jBFXufTYyMjJ8x5EooNKQiDR58mTMLKKnpk7p0qULe/fu5f333/cdRaKASkMi0qmpqYoVK/qOEnRNmzalYsWKjBkzxncUiQIqDYk4a9euZf369RE/NXVKXFwcHTp0YPbs2ezfv993HIlwKg2JOKemptq1a+c7SsgkJyeTkZHBlClTfEeRCKfSkIgTTVNTp9x4443UrFlTO/pJ0Kk0JKKsW7eOdevWRc3U1ClmRnJyMosWLSI1NdV3HIlgKg2JKNE4NXVKUlISABMmTPCcRCKZSkMiypQpU2jYsGFUTU2dUr16dRo0aKApKgmqsCgNM7vIzKaY2QYzW29mN/nOJIXPxo0bWbNmTVSOMk5JTk5m9erVrF692ncUiVBhURrAEOAD51wtoD6w3nMeKYRSUlIAaNOmjeck/tx///3ExsZqtCFB4700zKw0cBswAsA5l+GcO+A3lRRGKSkp3HDDDVStWtV3FG8uueQS7rzzTsaPH09WVpbvOBKBvJcGcAWQBowys5VmNtzMSpz+JDPrbmbLzWx5Wlpa6FNKWNu5cyeff/55VI8yTklOTub7779n8eLFvqNIBAqH0ogDbgBed85dDxwBnjr9Sc65Yc65ROdcYvny5UOdUcLczJkzcc6pNIDWrVtTokQJxo4d6zuKRKBwKI1twDbn3NKc21PILhGRgKWkpFC9enXq1q3rO4p3JUqUoE2bNkyePJn09HTfcSTCeC8N59xO4Hszq5lzVxNgncdIUsgcPHiQ//73v7Rp0wYz8x0nLCQnJ3PgwAHmzJnjO4pEGO+lkeNPwDgz+xq4DhjkOY8UInPmzCEjI0NTU7k0bdqU8uXLaysqKXBxvgMAOOdWAYm+c0jhlJKSQvny5bnpJu3ec0pcXBwPPPAAw4cP5+DBg5QqVcp3JIkQ4TLSELkg6enpzJo1i1atWhEbG+s7TlhJSkri+PHjTJ8+3XcUiSAqDSnUFi5cyKFDh2jbtq3vKGGnQYMGXHHFFYwfP953FIkgeS4NMythZvqTTsLC9OnTKVGiBE2aNPEdJeyYGUlJScyfP59du3b5jiMR4rylYWYxZpZkZrPMbDewAdhhZuvM7J9mdlXwY4qcKSsrixkzZtCiRQuKFSvmO05YSkpKIisri0mTJvmOIhEikJHGR0B1oC9Q0TlXxTl3CXArsAR4wcw6BjGjyFktW7aMnTt3aqupX1C7dm3q16+vKSopMIGUxh3OuYHOua+dcz8dzMY5t885N9U51w7QnzEScikpKcTFxXHXXXf5jhLWkpOTWbp0KZs2bfIdRSLAeUvDOXcCwMwuM7Nrz3ZcqFPPEQkV5xzTp0/nd7/7HRdddJHvOGGtQ4cOgE7OJAUjkHUa1czsS2ApkALsNrP3zOzqoKcTOYcNGzbwP//zP5qaCkCVKlW47bbbGDduHM4533GkkAtkeuoF4E3n3KXOuepAaeA9YI6Z1QhqOpFzOHXujFatWnlOUjgkJyezYcMGVq1a5TuKFHKBlMbVzrk3T91wzp10zg0D/gD8PWjJRH5BSkoKv/71r7nssst8RykU2rVrR5EiRbRCXPItkNI463jWOfchcE3BxhE5v+3bt7Ns2TLt0JcHF198Mc2bN2fChAk6OZPkSyClUdHMuprZb83sV6c9pglSCbkZM2YA0X1a1wuRlJTE9u3b+fjjj31HkUIskNIYQPaRZ/8f8J2ZbTWz983seaBiMMOJnE1KSgpXX301tWrV8h2lUGnVqhUlSpTQFJXkSyCb3A5zzv3JOdfIOXcx2Tv1vQrsBxYFO6BIbvv37+ejjz7SuTMuQEJCAm3btmXKlCk6OZNcsEA2ua2a+5LzmrXABODpXI/p2MsSdLNnz+bkyZNan3GBkpKS2L9/Px988IHvKFJIBXI+jbfPcf+p9RmWc3008E4BZBI5p5SUFCpWrMhvfvMb31EKpTvuuINy5coxfvx4Wrdu7TuOFELnLQ3n3O9CEUTkfI4fP86cOXPo2LEjMTE6qv+FKFKkCA888AAjRozg0KFDlCxZ0nckKWT0zZNCY/78+Rw5ckRbTeXTqZMzndpBUiQvAi4NM2thZkvN7Bsze9fMdG5NCamUlBRKlSrF7bff7jtKoXbTTTdRrVo1nT9cLkheRhqvAX8BGgDDgH+a2YNBSSVymszMTGbOnEnLli2Jj4/3HadQMzMefPBBnZxJLkheSmO3c+5T59x+59x84E7gmSDlEvmZzz77jLS0NE1NFZDk5GQyMzOZPHmy7yhSyOSlNLaY2XNmdurPvBPAySBkEjnD9OnTiY+Pp0WLFr6jRIQ6depQr1497egneZaX0sgC2gLfm9li4FtgoY50K8HmnCMlJYUmTZpQqpR2ByooSUlJfP7552zevNl3FClEAtm5zwCcc0nOuTpAVaA32YcXMWCYmW0PZkiJbqtXr2bLli3aoa+APfhg9ipJnZxJ8iKgc4Sb2Z9y9gbHOZfunFsBjAVmAt8D/YKYUaLc9OnTMTOdO6OAVa1alYYNG+rkTJIngZRGcyATmGBmO8xsnZltBjYCHYB/O+dGBTOkRLeUlBRuvvlmKlSo4DtKxElKSmL9+vV89dVXvqNIIRHIAQuPO+dec87dQvbUVBPgBufc5c65bs65lUFPKVFry5YtrFq1SltNBUn79u2Ji4vTCnEJWCDrNDqb2R4z2wcMBw475w4EP5qIzp0RbOXKldPJmSRPApme+hvQFKgFpAKDgppIJJfp06dTt25drrrqKt9RIlZSUhLbtm1j8eLFvqNIIRBIaRx0zq10zu12zv0N0OFFJSTS0tJYvHixtpoKslatWpGQkKDDikhAAimNSmbW3cxuM7PyQJFghxIBeO+998jKytLUVJCVKFGCNm3aMHnyZDIyMnzHkTAXSGn0B64FBgLfAHXNbLaZDdaxpySYUlJSqFq1Ktdff73vKBEvOTmZ/fv3M2fOHN9RJMzl9XSvZYErgVeAA0DLYAeU6HT48GE+/PBDndY1RJo2bUqFChV4++1znXNNJFsgZ+77GefcNmAboD9JJGjmzp1Lenq6pqZCpEiRIiQnJ/PKK6+wZ88eypUr5zuShCmdhEnCUkpKCmXLlqVhw4a+o0SNzp07c+LECR1WRH5R2JSGmcWa2Uoze993FvHrxIkTvP/++9xzzz3ExeV5MCwXqF69elx//fWMHj3adxQJY2FTGmQfBHG97xDi30cffcSBAwe0qa0HXbp04csvv2T16tW+o0iYCovSMLPLgLvI3uNcotyUKVP41a9+RbNmzXxHiToPPvggcXFxWiEu5xQWpQH8B3iC7HN2nFXOviLLzWx5Wlpa6JJJSJ08eZLp06dzzz33ULx4cd9xok758uW56667GDt2LCdP6hxrcibvpWFmd5N9KtkVv/S8nE1/E51zieXLlw9ROgm1RYsWsWfPHtq3b+87StTq0qULu3btYu7cub6jSBjyXhrALUArM9sKTARuN7OxfiOJL1OmTCEhIYHmzZv7jhK1WrZsycUXX6wpKjkr76XhnOvrnLvMOVeN7PNzLHDOdfQcSzzIzMxk2rRp3H333SQkJPiOE7Xi4+NJTk5mxowZ7Nu3z3ccCTPeS0PklMWLF7N7925NTYWBzp07k5GRoX025AxhVRrOuYXOubt95xA/Jk+eTPHixWnZUken8e3666+nfv36DB8+XKeClZ8Jq9KQ6JWVlcXUqVNp2bIlJUqU8B0n6pkZ3bt3Z9WqVaxY8YvbqEiUUWlIWPjss8/YuXOnpqbCSHJyMsWLF+ett97yHUXCiEpDwsKUKVMoWrQod911l+8okqN06dLcf//9jB8/nsOHD/uOI2FCpSHeZWVlMWXKFFq0aEHJkiV9x5FcunXrxuHDh5k0aZLvKBImVBri3dKlS9m+fbumpsLQzTffzDXXXKMpKvmJSkO8mzhxIkWLFuWee+7xHUVOc2qF+NKlS3UQQwFUGuLZyZMnmTRpEnfffTelSpXyHUfOolOnTsTHx2u0IYBKQzxbuHAhu3bt4sEHdbr5cHXxxRfTrl07xowZw7Fjx3zHEc9UGuLVhAkTKFmypHboC3PdunXjwIEDTJ061XcU8UylId6kp6czdepU2rZtq8Ogh7nGjRtTo0YNXn/9dd9RxDOVhnjzwQcf8OOPP2pqqhAwM/74xz/y2WefsXLlSt9xxCOVhngzYcIEypUrR5MmTXxHkQB06dKFhIQEXn31Vd9RxCOVhnhx+PBhZs6cyf3330+RIkV8x5EAXHTRRXTq1Ilx48bpkOlRTKUhXsyYMYNjx45paqqQ6dWrF8ePH2fkyJG+o4gnKg3xYty4cVSpUoWbb77ZdxTJg2uvvZZGjRrx2muvkZmZ6TuOeKDSkJDbsWMHc+fOpWPHjsTE6J9gYfPoo4+yZcsW5syZ4zuKeKBvrITcuHHjyMrKonPnzr6jyAVo3bo1l156Ka+88orvKOKBSkNCyjnH22+/TYMGDahZs6bvOHIBihQpQs+ePfnwww9Zt26d7zgSYioNCamVK1eyZs0ajTIKuZ49e1K8eHFeeukl31EkxFQaElKjR4+maNGiPPDAA76jSD6UK1eOLl26MGbMGHbu3Ok7joSQSkNCJiMjg/Hjx9O6dWvKlCnjO47k02OPPcaJEycYOnSo7ygSQioNCZnZs2ezd+9eTU1FiBo1atC6dWtef/11jhw54juOhIhKQ0Jm9OjRVKxYkWbNmvmOIgWkT58+7Nu3j9GjR/uOIiGi0pCQ2LVrF7NmzSI5OZm4uDjfcaSA3Hzzzfz2t7/l3//+t3b2ixIqDQmJUaNGcfLkSX7/+9/7jiIFyMzo06cPmzZtYtq0ab7jSAiYc853hjxLTEx0y5cv9x1DApSVlUWNGjWoUqUKCxcu9B1HClhmZia1a9emePHirFy5EjPzHUnOwcxWOOcS8/MeGmlI0M2fP5/NmzfTs2dP31EkCGJjY3n66af56quveP/9933HkSDTSEOC7t577+WTTz5h27ZtFC1a1HccCYITJ05Qs2ZNypUrx9KlSzXaCFMaaUjY++GHH5g5cyYPP/ywCiOCFSlShL59+/LFF1/w4Ycf+o4jQaTSkKAaOXIkmZmZdO/e3XcUCbLOnTtTpUoVBg4cSGGcwZDAqDQkaDIzM3nrrbdo0qQJV111le84EmTx8fE8+eSTfPrpp9rgIYKpNCRoPvjgA1JTU7UCPIp07dqVSpUq0b9/f402IpRKQ4Lm5ZdfpnLlyrRq1cp3FAmRYsWK0a9fPz755BOdpClCqTQkKNauXcuHH35Ir169iI+P9x1HQqhbt25Ur16dvn37kpWV5TuOFDDvpWFmVczsIzNbZ2Zrzay370ySfy+//DLFihXTCvAoVKRIEQYOHMjXX3/NxIkTfceRAuZ9Pw0zqwRUcs59aWYlgRVAG+fcOU8Jpv00wtvevXu57LLL6NSpE8OGDfMdRzzIysrihhtu4NChQ6xfv16jzTAREftpOOd2OOe+zLl+CFgPXOo3leTHsGHDOH78OL17a9AYrWJiYhg8eDCbN2/mrbfe8h1HCpD3kUZuZlYN+Bio65w7eNpj3YHuAFWrVr3xu+++C3k+Ob/09HSuvPJKateuzbx583zHEY+cczRu3JgNGzawceNGSpUq5TtS1IuIkcYpZvYrYCrwf04vDADn3DDnXKJzLrF8+fKhDygBGTNmDD/88AOPP/647yjimZnxz3/+k927d/Pcc8/5jiMFJCxGGmZWBHgfmOucO++Z6rVOIzxlZmZyzTXXULJkSZYvX67jDwkAjzzyCGPHjmXNmjVcffXVvuNEtYgYaVj2b5YRwPpACkPC19SpU9m4cSN9+/ZVYchPBg8eTPHixXnsscd8R5EC4L00gFuATsDtZrYq59LSdyjJG+ccgwcP5uqrr6Zt27a+40gYqVChAv3792f27NnMmjXLdxzJp7CYnsorTU+Fnzlz5tCyZUtGjBjBI4884juOhJmMjAzq16/PyZMnWbNmjY547ElETE9J4eec4+9//zuXX345HTt29B1HwlB8fDxDhgzh22+/5YUXXvAdR/JBpSH5NnPmTJYvX07//v21E5ecU7NmzejQoQPPPfcc69adc99dCXOanpJ8ycrK4rrrruP48eOsW7eOuLg435EkjO3evZvatWtTo0YNFi9eTGxsrO9IUUXTU+Ld5MmTWb16NQMGDFBhyHldcsklDBkyhCVLljB06FDfceQCaKQhF+zEiRNce+21xMXF8dVXX+mvRgmIc467776bhQsXsmbNGq644grfkaKGRhri1RtvvME333zD4MGDVRgSMDPjjTfeIDY2lkceeYTMzEzfkSQPVBpyQfbv38+AAQNo0qQJd999t+84UshUqVKF//znPyxcuJAXX3zRdxzJA5WGXJCBAweyf/9+XnrpJe39LRfk4Ycfpn379vTr1w9NNxceKg3Js40bNzJ06FC6du1KvXr1fMeRQsrMePPNN6lYsSJJSUkcPnzYdyQJgEpD8sQ5R69evShWrBgDBw70HUcKubJlyzJmzBg2bdrE73//ewrjhjnRRqUheTJ+/HjmzZvH4MGDqVixou84EgEaN27MP/7xDyZNmsSQIUN8x5Hz0Ca3ErB9+/ZRq1YtrrzySj799FNtMSUFxjnHvffey3vvvceCBQu47bbbfEeKSNrkVkLqiSeeYN++fQwbNkyFIQXKzBg9ejTVq1enffv2bN682XckOQeVhgRk9uzZjBgxgj59+mjltwRF6dKlmTlzJidPnuSuu+7iwIEDviPJWag05LzS0tJ45JFHuPbaa3n22Wd9x5EIVrNmTaZPn86mTZto164dGRkZviPJaVQa8oucc3Tv3p39+/czduxYnQdBgq5Ro0YMHz6cBQsW0KlTJ+0xHmZ0hDn5RW+99RYpKSn885//1LSUhMxDDz3E7t27efzxxylZsiRvvfWWdiINEyoNOadly5bxpz/9iWbNmvGXv/zFdxyJMn369OHgwYMMHDiQhIQEhgwZouIIAyoNOavdu3fTrl07KlWqxPjx44mJ0UymhN6zzz7LkSNHeOmll0hPT+f111/Xv0XPVBpyhoyMDDp06MCePXv49NNPufjii31HkihlZrz44osUK1aMQYMGcfToUUaNGqVzt3ikJS8/45yja9eufPTRR7zzzjvccMMNviNJlDMz/vGPf5CQkEC/fv3Ys2cPkyZNolSpUr6jRSWN8+Rnnn76acaOHctzzz1Hp06dfMcR+ckzzzzDsGHDmDdvHg0bNuT777/3HSkqqTTkJ//61794/vnn6dmzJ08//bTvOCJn6NatG3PmzGHr1q0kJiaycOFC35GijkpDAHj++efp06cP999/P0OHDtVWKhK2mjZtyueff07ZsmVp0qQJL7zwAllZWb5jRQ2VRpRzzjFw4ED69t/r4iwAAAo/SURBVO1LUlIS48aN03GlJOzVrl2bZcuWce+99/LUU0/RrFkzUlNTfceKCiqNKHbixAm6d+/O3//+dx566CHeeecdbZUihUbJkiV59913eeONN1iyZAl169Zl5MiROidHkKk0otT+/ftp0aIFw4cPp1+/fowaNUojDCl0zIwePXqwevVqbrjhBrp27UqzZs1Yu3at72gRS6URhT7//HOuv/56Pv74Y95++20GDhyoHaakULviiitYsGABQ4cOZfny5dSvX58///nP7Nu3z3e0iKPfFFHkxIkTDBo0iIYNGxITE8Mnn3zCQw895DuWSIGIiYmhV69ebNy4ke7du/Pqq69SvXp1nn32Wfbv3+87XsRQaUSJpUuX8utf/5pnnnmG9u3bs3LlSn7729/6jiVS4MqVK8drr73GypUradSoEQMGDKBatWr069ePnTt3+o5X6Kk0IlxqaiqPPPIIN910E3v27GHatGlMmDCB0qVL+44mElT16tUjJSWFVatW0axZMwYNGkSVKlV44IEHWLRokVaYXyCVRoT67rvv6N27NzVq1GD8+PE89thjrFu3jrZt22ofDIkq9evXZ/LkyXzzzTf8+c9/Zt68eTRu3Jirr76av/3tb1ppnkdWGNs2MTHRLV++3HeMsJOVlcWiRYt4/fXXmTp1KmZG586dGTBgAFWqVPEdTyQsHDt2jEmTJjFu3DgWLFhAVlYWtWvXpmXLljRv3pxbb701Yk82ZmYrnHOJ+XoPlUbhlpmZybJly5g5cybjx48nNTWV0qVL06NHDx599FGVhcgv2LlzJ1OmTGH69OksXryYjIwMEhISaNy4MbfeeisNGjQgMTGRkiVL+o5aIFQaUSgrK4v169fz+eefs2jRIj744AP27NlDbGwsTZs2pXPnzrRq1YqEhATfUUUKlcOHD7Nw4ULmzp3L/Pnz2bBhA5C9VVadOnWoX78+tWvXpk6dOtSpU4dq1aoVun2bIqY0zKw5MASIBYY7557/pedHQ2lkZGSwY8cONm7cyPr163+6rFy5kh9//BHI3kqkefPm3HXXXdx5552UKVPGc2qRyLF//36WLVvGkiVLWLJkCWvWrGHbtm0/PV6sWDGqVavG5Zdf/tOlatWqVKpUifLly3PJJZdQrly5sDrKQkSUhpnFAv8DNAW2AV8ADzrn1p3rNQVZGs65ny65b+fnekZGBunp6T+75L7v8OHD/Pjjj/z4448cOHDgp+s7d+7khx9+YPv27aSlpf0s50UXXcQ111xDvXr1uOmmm7jpppuoUaOGVmqLhNDBgwdZv349a9euZd26dWzdupXvvvuO77777ozv7CllypShfPnylC5dmlKlSlGyZElKlSr1s+slSpSgaNGiZ1zi4+PPuB0bG0tMTAyxsbF5vh4TExMRpXETMMA5d2fO7b4AzrnB53pNbGysK1asWL5+uYeTEiVKULp0aSpUqEDlypW59NJLqVy5MpUrV+aqq67immuuoUKFCioIkTB29OhRUlNT2b17N7t37yYtLe2nn2lpaRw8eJCDBw9y6NChn/08efJkSPIVL16cY8eO5bs0wmHcdCmQ+2wq24Az9jozs+5A95yb6UePHl0Tgmz5VQ7Yc74nHTlyhCNHjvDDDz+wcuXKEMQ6Q0A5PSsMGUE5C5pyFpBjx44B1Mzv+4RDaQTEOTcMGAZgZsvz25ahoJwFpzBkBOUsaMpZsMws3/P64bBz33Yg93ahl+XcJyIiYSYcSuMLoIaZXWFm8UAHYKbnTCIichbep6eccyfN7FFgLtmb3I50zp1vv/5hwU9WIJSz4BSGjKCcBU05C1a+c3rfekpERAqPcJieEhGRQkKlISIiAQur0jCz5mb2jZl9a2ZPneXx28zsSzM7aWbtT3uss5ltzLl0DuOcmWa2KucS1BX+AeT8i5mtM7Ovzey/ZnZ5rsfCaXn+Us5wWp49zWx1TpbFZlY712N9c173jZndGY45zayamR3LtTzf8Jkz1/PamZkzs8Rc94VkeV5oxnBblmbWxczScuX5fa7H8vZdz30YDZ8XsleCbwKuBOKBr4Dapz2nGlAPeAdon+v+ssDmnJ9lcq6XCbecOY8dDqPl+TsgIef6H4BJYbo8z5ozDJdnqVzXWwEf5FyvnfP8osAVOe8TG4Y5qwFrwmV55jyvJPAxsARIDOXyzGfGsFqWQBdg6Flem+fvejiNNH4DfOuc2+ycywAmAq1zP8E5t9U59zWQddpr7wTmOef2Oef2A/OA5mGYM5QCyfmRc+5ozs0lZO8jA+G3PM+VM5QCyXkw180SwKmtTFoDE51z6c65LcC3Oe8XbjlD6bw5cwwEXgCO57ovVMszPxlDKdCcZ5Pn73o4lcbZDidyaQhem1f5/axiZrbczJaYWZuCjfYzec3ZFZhzga/Nj/zkhDBbnmbWy8w2Af8P+HNeXhsGOQGuMLOVZrbIzBoGKWNAOc3sBqCKc25WXl8bBhkhjJZljnY5U7xTzOzUDtV5Xpbe99OIQpc757ab2ZXAAjNb7Zzb5DOQmXUEEoFGPnOczzlyhtXydM69CrxqZklAPyCo64Mu1Dly7gCqOuf2mtmNQIqZ1TltZBISZhYDvET2tEpYOk/GsFmWOd4DJjjn0s2sB/A2cPuFvFE4jTTycziRUB6KJF+f5ZzbnvNzM7AQuL4gw+USUE4zuwN4BmjlnEvPy2vDIGfYLc9cJgKnRj5htzxz+SlnznTP3pzrK8ieJ7/aU86SQF1goZltBRoAM3NWNIdqeV5wxjBbljjn9ub63gwHbgz0tWcIxYqaAFfmxJG9EuYK/ndlTp1zPHc0Z64I30L2ipwyOdfLhmHOMkDRnOvlgI2cZcVaqHKS/Qt2E1DjtPvDann+Qs5wW541cl2/B1iec70OP19xu5ngrQjPT87yp3KRvVJ1ezh8j3Kev5D/XckckuWZz4xhtSyBSrmutwWW5FzP83e9wP8D8vkf35LsEzJtAp7Jue//kv3XJcCvyZ5zOwLsBdbmeu0jZK8Q+xZ4OBxzAjcDq3P+p64GunrOOR/YBazKucwM0+V51pxhuDyHAGtzMn6U+4tL9ihpE/AN0CIccwLtct3/JXCPz5ynPXchOb+QQ7k8LzRjuC1LYHBOnq9y/p/XyvXaPH3XdRgREREJWDit0xARkTCn0hARkYCpNEREJGAqDRERCZhKQ0REAqbSEBGRgKk0REQkYDr2lEgBMLMJZP8RdgVQAfijO/tB7EQKNY00RApGfWCzc+43QDLQ33MekaDQHuEi+WRmxcg+vHQV59xxMysLLHXO1fAcTaTAaaQhkn91gY3OuVMn4bmB7GP8iEQcrdMQyb/6QNWcEUcs8CzwhN9IIsGh0hDJv/rANGApUAQY5Jz71G8kkeDQOg2RfDKzRUB359w3vrOIBJtKQySfzGwb2af2zPKdRSTYVBoiIhIwbT0lIiIBU2mIiEjAVBoiIhIwlYaIiARMpSEiIgFTaYiISMBUGiIiErD/D+PeQG0nMp+rAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Observed number of AGN\n", "n = 30\n", "# Total number of objects\n", "nmax = 100\n", "# Range of values for the AGN fraction to evaluate\n", "# This is the same as the width of the uniform prior in P(frac)\n", "xmin,xmax,nx = 0.1,0.5,400\n", "# Array of values for the AGN fraction\n", "dx = (xmax-xmin)/nx\n", "x = np.linspace(xmin+0.5*dx,xmax-0.5*dx,nx)\n", "# Compute the (unnormalised) probability of each AGN fraction value\n", "# This calculation uses Bayes' theorem: P(frac|n) \\propto P(n|frac) P(frac)\n", "# P(n|frac) is determined using the Binomial probability distribution\n", "# P(frac) is the prior, which we take as uniform and therefore do not include in the calculation\n", "prob = np.empty(nx)\n", "for i in range(nx):\n", " prob[i] = stats.binom.pmf(n,nmax,x[i])\n", "# Normalise the probability distribution. This ensures that \\int P(x) dx = 1\n", "# Expressing the integral as a sum, we are finding N such that N \\sum P(x) dx = 1, or N = 1/\\sum P(x) dx\n", "prob /= np.sum(prob)*dx\n", "# Plot the normalised posterior probability distribution for P(frac|n)\n", "ymax = 10.\n", "plt.plot(x,prob,color='black')\n", "plt.xlim(xmin,xmax)\n", "plt.ylim(0.,ymax)\n", "plt.xlabel(r'$p$')\n", "plt.ylabel(r'$P(p|D)$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could then quote our answer as the mean and standard deviation of this distribution:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mean = 0.3039158356556867\n", "Standard deviation = 0.04530760743511314\n" ] } ], "source": [ "# These calculations work because the probability distribution is already normalised such that \\int P(x) dx = 1\n", "x1 = np.sum(x*prob*dx) # determine by evaluating \\int x P(x) dx = \\sum x P(x) dx\n", "x2 = np.sum((x**2)*prob*dx) # determine by evaluating \\int x^2 P(x) dx = \\sum x^2 P(x) dx\n", "mu = x1 # mean = \n", "sig = np.sqrt(x2 - x1**2) # standard deviation = \\sqrt( - ^2)\n", "print('Mean =',mu)\n", "print('Standard deviation =',sig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "Consider the following question:\n", "\n", "_A survey of area $A = 1$ deg$^2$ finds $N=20$ quasars. What is the number of quasars per sq degree, $\\sigma$?_\n", "\n", "Solve this question using both a Frequentist and a Bayesian approach.\n", "\n", "## Monte Carlo simulations\n", "\n", "A **Monte Carlo simulation** is a computer model of an experiment, in which many random realizations of the results are created and analysed like the real data.\n", "\n", "This powerful technique allows determination of errors from the **distribution of results over the realizations**.\n", "\n", "## Activity: Monte Carlo methods\n", "\n", "Solve the following problem by Monte Carlo methods: I'm dealt 5 playing cards from a normal deck (i.e. 13 different values in 4 suits). What is the probability of obtaining \"3 of a kind\" (i.e. 3 of my 5 cards having the same value?)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity: central limit theorem\n", "\n", "Write a code that draws $n$ values of $x$ from an exponential distribution $P(x) \\propto e^{-x}$ (where $0 < x < \\infty$) and computes their arithmetic mean $\\mu$. Repeat this process $m$ times, and plot the probability distribution of $\\mu$ across the $m$ realisations. Run this experiment for values $n = 1,2,5,10,20,50$.\n", "\n", "Hint: to do a single draw, select a uniform random number $y$ in the range $0 < y < 1$, then $x = -\\ln{y}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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": 2 }