day-3-data-mining/01 - Bootstrapping MC chains.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:33:01.579122Z",
     "start_time": "2019-07-03T11:33:01.562072Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunss.otf');\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        font-weight: bold;\n",
       "        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunsx.otf');\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        font-style: oblique;\n",
       "        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunsi.otf');\n",
       "    }\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        font-weight: bold;\n",
       "        font-style: oblique;\n",
       "        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunso.otf');\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: Helvetica, serif;\n",
       "    }\n",
       "    h4{\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "    div.text_cell_render{\n",
       "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 145%;\n",
       "        font-size: 130%;\n",
       "        margin-left:auto;\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\", source-code-pro,Consolas, monospace;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 22pt;\n",
       "        color: #4057A1;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"../Data/www/styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synopsis\n",
    "\n",
    "In this unit we will discuss how to randomly generate possible data sets from a set of sample data (known as bootstrapping). To do that we will cover:\n",
    "\n",
    "* Why/when you would need to bootstrap data\n",
    "* How to bootstrap data\n",
    "* How to statistically evaluate bootstrapped data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:33:10.485263Z",
     "start_time": "2019-07-03T11:33:08.523523Z"
    }
   },
   "outputs": [],
   "source": [
    "import json\n",
    "import scipy.stats as stat\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.special import erf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:33:12.380707Z",
     "start_time": "2019-07-03T11:33:12.372082Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:33:13.701800Z",
     "start_time": "2019-07-03T11:33:13.692915Z"
    }
   },
   "outputs": [],
   "source": [
    "def half_frame(sub, xaxis_label, yaxis_label, font_size = 15, padding = -0.02):\n",
    "    \"\"\"Formats frame, axes, and ticks for matplotlib made graphic with half frame.\"\"\"\n",
    "\n",
    "    # Format graph frame and tick marks\n",
    "    sub.yaxis.set_ticks_position('left')\n",
    "    sub.xaxis.set_ticks_position('bottom')\n",
    "    sub.tick_params(axis = 'both', which = 'major', length = 7, width = 2, direction = 'out', pad = 10,\n",
    "                    labelsize = font_size)\n",
    "    sub.tick_params(axis = 'both', which = 'minor', length = 5, width = 2, direction = 'out', labelsize = 10)\n",
    "    for axis in ['bottom','left']:\n",
    "        sub.spines[axis].set_linewidth(2)\n",
    "        sub.spines[axis].set_position((\"axes\", padding))\n",
    "    for axis in ['top','right']:\n",
    "        sub.spines[axis].set_visible(False)\n",
    "\n",
    "    # Format axes\n",
    "    sub.set_xlabel(xaxis_label, fontsize = 1.6 * font_size)\n",
    "    sub.set_ylabel(yaxis_label, fontsize = 1.6 * font_size)\n",
    "    \n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap methods\n",
    "\n",
    "When doing science, ideally, one would have many replicates of an experiment.  The replicates would enable us to calculate both an expected outcome and the variability in outcomes.  Those quantities would then enable us, under certain assumptions, to estimate the chance that some observed difference from a null model could be explained by chance or not.\n",
    "\n",
    "But what can we do when we only have one replicate?  This is frequently the case with observational studies, but it can occur in other contexts too.\n",
    "\n",
    "In order to learn more about these situations, let us consider the data we obtained when analyzing the sentiment score of the lines of Othello and Iago.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:34:15.172243Z",
     "start_time": "2019-07-03T11:34:15.162801Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "537 470\n",
      "[1, 1, -3, 3, 1, 2, -3, 2, 3, 2]\n"
     ]
    }
   ],
   "source": [
    "with open('../Data/Day7-Structured-Data-Analysis/iago.json', 'r') as file1:\n",
    "    iago_sentiment = json.load(file1)\n",
    "    \n",
    "\n",
    "with open('../Data/Day7-Structured-Data-Analysis/othello.json', 'r') as file1:\n",
    "    othello_sentiment = json.load(file1)\n",
    "    \n",
    "print(len(iago_sentiment), len(othello_sentiment))\n",
    "\n",
    "print(iago_sentiment[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:34:23.819010Z",
     "start_time": "2019-07-03T11:34:23.574377Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/lgaalves/miniconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6521: MatplotlibDeprecationWarning: \n",
      "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n",
      "  alternative=\"'density'\", removal=\"3.1\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAFKCAYAAADc7h9KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcVNWd9/HPVwREVARB0LiAmmhQJ6MxImoCBnEhMe7bJGM00VFjopMY48QYbZ2YjCsx8RHjiibuG+4yoiJqcMvyPFFClLjvoCDIooK/549zC4qiqrtvdXVXNf19v173VV3nnnPvr6ob6lfnnnOuIgIzMzOzRrVKvQMwMzMza46TFTMzM2toTlbMzMysoTlZMTMzs4bmZMXMzMwampMVMzMza2hOVszMzKyhOVkxMzOzhuZkxczMzBqakxUzMzNraE5WzMzMrKE5WWkHkkKSb7pkZmZWA6vWO4AuzgmNmZl1JaqmkXtWzMzMrKE5WTEzM7OG5mTFzMzMGpqTFTMzM2toTlbMzMysoTlZMTMzs4bmZMXMzMwampMVMzMza2hOVszMzKyhOVkxMzOzhuZkxczMzBqa7w1kZmat19S0cp7LGlqunhVJf5C0S3sFY2ZmZlYq72WgfwMmSXpR0s8lbdgeQZmZmZkV5L0M9Htgf2Aw0AScLmkScCUwISI+rml0tqJ6dYu28bzjx49n3LhxPPvss3Tr1o1tttmGE088kW984xtL61x66aWsu+667LPPPsu1HTx4MAcccADnnXdem2IAkMRvf/tbvv/97wMwcuRI+vfvzy233NLmY5uZWfvI1bMSEd8GBgH/ATyZtd8NuB54S9JvJG1T8yitUzv22GM58sgjGTZsGBMmTODGG29k8ODB7L333px99tlL61166aVMmDChjpGamVkjyj3ANiI+BC4HLpf0OeC7wLeA9YDjgOMk/V/gCuC6iJhdw3itk5kwYQKXXHIJ48aN45hjjllavueeezJo0CBOOeUURo8ezbbbblvHKM3MrJG1aepyRDwfEScDGwLfAO4AFgP/CvwGeFPS9ZJ2a3Ok1ildeOGFbLbZZhx11FEr7DvllFNYc801ueiiixg5ciR/+tOfuPrqq5GEJMaPH79c/bFjx7LBBhvQt29fDjnkEObMmbPc/vfff5+jjz6agQMHstpqq7Hjjjvy5JNP5o75oYceYtiwYay22moMHDiQ733ve3z44Ye5j2NmZrVRk3VWIuLTiLg7IvYDNgUeBwT0BA4C7ssG5Z4gqUctzmmNb/HixUydOpW99tqLbt26rbC/T58+7LLLLkyZMoWLL76YLbbYgjFjxjB16lSmTp3K1772taV1b7rpJh588EEuvfRSzj77bO6++25OOeWUpfs/+ugjdt11Vx544AHOPfdcJkyYwIABA9h11115++23Wx3ztGnT2GOPPejfvz+33norZ5xxBtdddx0HHHBA294MMzOrWs3WWZG0LXAEcCjQNyv+CHgE2IE0KPcC4GhJu0fEa7U6tzWmWbNm8dFHH7HxxhtXrLPxxhtz//33M3ToUHr37s2AAQPYYYcdVqjXvXt3JkyYwKqrpj/ZadOmccMNN3DxxRcD8Ic//IFnn32W5557js9+9rMA7Lrrrmy++eacf/75nHvuua2K+cwzz2TjjTfmzjvvXJpg9evXj4MPPpipU6cyfPjwXO+BmZm1XZt6ViT1l/Sf2RiVp0ljVvoBzwH/CawfEXuQxrMcBbwBbA60fVqHdSm77LLL0kQFYOjQobz77rt8/HGagDZp0iS++MUvMmTIEBYvXszixYsBGDFiBM8880yrz/PUU0+x7777LtcTtP/++7Pqqqvy2GOP1ejVmJlZHrl7ViStAowh9aJ8DehOuuTzIXAjcHlELDdQICIWAldIehCYAXy1jXFbJ9C/f3969uzJK6+8UrHOK6+8wmc+85kWj7X22msv97xHjx5EBB9//DE9evRg1qxZPPHEE3Tv3n2FtptuummrY37rrbcYOHDgcmXdunVjnXXW4f3332/1ccys7ZqaHu+g8+zUIeex6uVKViSdQ5r5M5CUoAA8RZoddH1EzG+ufUS8LOktYP0qYrVOZtVVV2X48OHcc889nHfeeayyyvIdeXPnzmXy5Mnsu+++bT5Xv3792G677Rg3btwK+3r27Nnq46y33nq8++67y5UtWbKE9957j379+rU5TjMzyy/vZaAfk9ZZmU2a7fMvEbFDRFzeUqJS5HFgSs7zWid1wgkn8Pzzz3P55ZevsO9//ud/mDt37tIF2nr06MGiRYuqOs+oUaOYMWMGG220Edttt91y29Zbb93q4wwbNozbb7+dJUuWLC277bbbWLx4MTvvvHNVsZmZWdvkvQw0GbgMuLXa1Woj4pCW6kgaCvwWGA7MIfXcnBERS5pp8yXge8CXST03rwHXAWdHxKKSujuRBvv+C/A2MDYiflPN67Hm7bPPPhxzzDEcd9xxTJs2ja9//essXryYG2+8kfHjx/OrX/1q6RorW2yxBRMnTmTixImss846DBkyhHXWWadV5znssMO45JJLGDlyJD/+8Y/ZZJNNeO+993jqqacYNGgQP/zhD1t1nFNPPZVtttmGffbZh2OPPZbXX3+dk08+md13392Da83M6iRXshIR7T7WRFJfYBIwDdibNBX6fFIv0KnNND04q3s28AIpEfnv7HH/ouNvBkwE7gZ+CmwPXCBpQUSs+PW/0XTCu5BefPHFDBs2jHHjxnHZZZexyiqrsO2223LHHXcst9z+qaeeyquvvspBBx3E3Llzueqqqzj88MNbdY7VVluNhx9+mNNOO43TTz+dd955h3XXXZftt99+uXO0ZMstt+S+++7jlFNOYb/99mOttdbi0EMP5Zxzzsn7ss3MrEYUEa2vLD0EvBcRB7ay/vXAuhExKsc5fgr8BNg4IuZmZT8h3YtoUKGsTLsBETGzpOw/gN8BgyPilazsd8AuwNCIWJyVXQzsBWwUed6Qyq8hACJCLVRt87nMzDpUB35hamJ0x5zHA2w7Ukufi2XlHbMyEsjzW90ha5PHnsDEkqTkBqAXMKJSo9JEJfOX7HHdkuPfVkhUio6/AbBVzljNzMysndVkBdtmdCN/78EWwPTigoh4FViQ7ctjR+BT4B8AknqTbg0wvaTe34vO3WqSotyWM0YzMzNrRrslK5J6kno0yl62aUZf0qDaUrNZtjJua84/CPgZ8PuiXprCYh2lxy/cbLHVxzczM7OO0ewAW0kbkZbJL9ZD0pepfN1JpKTgUKAH8Mcq4irXO6EK5StWTPcfuom0UF25aSCVjpOrV6TSmBT3rpiZmdVOS7OBjgBOKynrS5rC3JLCB/mvc8Y0m2U9IMX6UL7HZfmTSgKuAbYEdoqI2UW7C+1Lj9+3ZL+ZmZk1iJaSlTnAq0XPNyaNAXm9mTafki79PAdcEREP54xpOiVjRyRtCPRmxbEm5YwlTXkeHRGlY1/mS3qt9PhFz1tzfDMzM+tAzSYrEXEhcGHhuaRPgZkRMaQdY7oPOEnSmhExLys7GFhIuoNzRdm05x8AB0VEpbvO3QfsK+nUokXmDiYtIvdsm6M3MzOzmsq7gu0ZpHEg7ekS4HjgNklnA5uQ1li5oHg6s6QZwCMR8d3s+b8BvwTGA29I2qHomP8smtp8LvBN4PeSLgO+BBwNHFuLNVbMzMystvKuYHtGewVSdI7ZkkYBFwF3kS5FjSUlLMVWJU2NLtgtezw824odQUpiiIgZkvYgLbd/H2m5/RM7xeq1ZmZmXVDenpUOERHTgGaX9o+IwSXPD2fFJKVS28dIy+xbO2tqauKiiy5i1qxZ9Q7FzMw6qYrJiqQrsx/fioiflZTlEYVLNdZ2TU2P1+m8Xo7azMzqo7melcOzx+mkxdUKZUG+tf0DcLJiZmZmVWkuWSmMT5lVpswst/nz53PyySfzwAMP8NprrzFw4EDGjBnDr371K9Zaa62l9WbPns2xxx7LXXfdRZ8+fTjhhBOYOXMmt9xyCy+//PLSen/961858cQTmTp1Kj179mTMmDFccMEFDBw4sA6vzszM2kvFZKXcYNqOGGBrK68FCxawZMkSzjrrLAYMGMBrr73GWWedxYEHHsjEiROX1jv88MN57LHHuPDCCxk0aBBjx47l+eefp1u3ZeOpZ86cyciRI/n85z/Pddddx4cffsh//dd/MXr0aJ555hl69OhRj5doZmbtoCEH2NrKacCAAYwbN27p88WLFzNkyBB23nlnXn31VTbaaCOeffZZ7rzzTm666SYOPPBAAEaNGsWGG27IGmussbTt+eefD8DEiROX9sp87nOfY9iwYdx6660ceuihHfjKzMysPdX8RoaSeknqU+vj2srh97//Pdtssw1rrLEG3bt3Z+eddwbg+eefB+CZZ54BYK+99lraplevXuy6667LHeepp55it912W+7y0fbbb8/gwYN57LFK6wGamVlnlCtZkbShpP+Q9I0y+7aW9CQwD3hf0lRJW9YqUOv8br/9dg477DCGDx/OzTffzBNPPMHtt98OwKJFiwB4++23WXPNNVlttdWWaztgwIDlnr/11ltlx6YMHDiQ999/v51egZmZ1UPenpUjgXHAF4sLs56UScB22TEFDAMelNS/BnHaSuDmm29m2LBhXHzxxey5554MGzaMvn37Lldn0KBBzJs3b2nyUjBz5szlnq+33nq8++67K5zjnXfeoV+/frUP3szM6iZvslLoi7+xpPwoYADppod7ACOAv2Vl/9mWAG3lsXDhQnr27Llc2bXXXrvc8+222w6AO++8c7l2DzzwwHL1hg0bxsSJE5k3b97SsqeffpqXX3556aUlMzNbOeRNVjYkrZvyQkn5vln5yRHxvxHxKCmBEfC1NkdpK4XRo0czZcoUzjrrLCZNmsSPfvQjHnzwweXqbLXVVuy1114ce+yxXHHFFdxzzz18/etfZ/XVV2eVVZb9uf7oRz8CYPfdd+eOO+7g2muvZb/99mPrrbdm//3379DXZWZm7SvvbKABwJyI+KRQIGk10s0APyHdyweAiHhK0ifAprUI1JLOvJLs0UcfzYsvvsiFF17IokWLGD16NNdddx077LDDcvXGjx/Psccey/HHH88aa6zBcccdxyabbMLTTz+9tM6AAQN4+OGHOfHEEzn00EPp0aMHY8aMYezYsZ62bGa2klGeGw1L+hDoERE9ispGAg8BUyNip5L6s4DeEdGrNuF2DpICICJaWunXd3luhcWLF7PVVlsxbNgwrr766nqHY9a1NTV13KkY3THn6cRfAjuhPCvgL5W3Z+UlYKikHSPij1nZAaQP3SnLRSN1B/oAb1QTmHVdN998M2+++SZbb701c+fO5bLLLuOFF17gmmuuqXdoZmZWB3mTlfuBLYGrJJ0KrEeaIQRwe0ndLwDdSINuzVqtd+/eXHXVVcyYMYMlS5aw9dZbc9ddd7H99r5RtplZV5Q3WTkH+CbwWeCGrEzAHRHxVEndwqDbKZjlMGbMGMaMGVPvMMzMrEHkSlYiYqakHYAm0joqc4F7gbOL62WXgA7M9k/EzMzMrEq57w0UEa8C32mhzifA56oNyszMzKyg5vcGMjMzM6slJytmZmbW0HJfBgKQ9Hlgf2AroC/QvZnqERGjqjmPmZmZWe5kRdIFwPGkWUCtWdzFC5+ZmZlZ1XIlK5KOY9mNCf8G3EFa9G1RxUZmZmZmbZC3Z+UoUk/JbyPCd1M2MzOzdpd3gG1hOvJptQ7EzMzMrJy8PSvzgUURMbc9gjEzMzMrlbdn5UlgLUkD2iMYMzMzs1J5k5Vfkcas/KwdYjEzMzNbQa5kJSIeJ91l+WhJl0ga3B5BmZmZmRXknbr8YvbjEtLMoKMkvQ/Ma6ZZRMSmVcZnZmZmXVzeAbaDy5Stk22VeFE4MzMzq1reZGWXdonCzMzMrIJcyUpEPNJegZitDJqaHu+g8+zUIecxM2sEvuuymZmZNbSq7roMIGlV4IvAhsDqEXFNzaIyMzMzy1SVrEg6GTgJ6FtUfE3R/rWBx4GewA4RMastQZqZmVnXlfsykKRrgV+SEpUXgcWldSJiDjAZGALs27YQzczMrCvLlaxIOgQ4FHgLGB4RnwXer1D9OkDA3m2K0MzMzLq0vD0r3yWtm3JCRDzVQt1ngE+Bf6kmMDMzMzPIn6xsQ0pA7mqpYkR8BHwA+KaHZmZmVrW8ycoawPyI+LiV9XuSluY3MzMzq0reZGUmsKaktVqqKGlLYHXg9WoCMzMzM4P8yUphec5DWlH3NNL4lodzngNJQyU9KGmBpDclnSmpWwttekg6V9KjkhZKKntPIknjJUWZbYu8cZqZmVn7y5us/JY0w+dMSV8sV0FSX0mXAweSkpWL8pxAUl9gUtZ2b+BM4ETgjBaarg4cCSwA/thC3enA8JLt5TxxmpmZWcfIe2+gxyWdS1oQ7o+SHgPWApB0HjAUGAGsljU5LSKeyxnTMUAvYL+ImAs8kF12apJ0TlZWLrY5kvpFREj6PvDVZs4xPyKeyBmXmZmZ1UHuReEi4mTgh8BHpLsw9yL1tvwQ2CN7vgA4PiJ+WUVMewITS5KSG7LjjmghtrKXfszMzKzzqmq5/Yi4UNJ4YH9gR2A9UuLzDjAVuDkiKi0W15ItgIdKzveqpAXZvhanTbfCUElzSbOVngZ+Vs0dpSuNizEzM7PaqfpGhhHxAXBlttVSX2BOmfLZLH8vomr9BXgSmEZaA+ZE0qWmnVux0J2ZmZl1sKqTlXZWrsdCFcrzHTjiwuUOKt1DSlxOAfbJeSyVK3ePi5mZWe3kHrPSAWYDa5cp70P5Hpc2iYiFwL3AtrU+tpmZmbVdxZ4VSQ9V2pdTRMSoHPWnk8amFMeyIdA729de3BtiZmbWgJq7DDSyhbZBujRTaR9Ud+nmPuAkSWtGxLys7GBgIZB7EGxLJPUizUD6U62PbWZmZm3XXLJyRIXyfqTVafsAU0gJxBukxGQ90vTir5BuYngm6bJOHpcAxwO3STob2ARoAi4ons4saQbwSER8t6hsT1IPzL9mzw/Idj0dEa9I6gPcDfwBmAH0J025/gxwUM44zczMrANUTFYi4urSsuzD/mnSGitfiYjHyrWVtCNwK2mBt+3zBBQRsyWNIq18exdpnMpYUsJSGnvpEvzjgI2Lnt+cPR4BjM/ingmcCqwLLCJNtR4REc/kidPMzMw6Rt7ZQKcBmwLfqJSoAETEHyUdSUo2fk5a8bbVImIaza9AS0QMbk1Zyf5FwH55YjEzM7P6yjsbaB9gYUTc04q695LGmeybOyozMzOzTN5kZX3g09ZUzJa+X5K1MTMzM6tK3stA7wHrSdopIh5vrqKknYA1gDerDc6sJpqaOvBkozvwXGbWUDrq/5oO/T+tMeTtWbmXNOvnKkmbVaokaVPgKtK05dZcMjIzMzMrK2/PyumkcSubAn+TdBtp6nKh92R90rTl/YDVgHezNmZmZmZVyZWsRMRbkkYAtwCfBw7JtlIi3W/nwIh4u81RmpmZWZeV+0aGEfF3SV8gJSkHkO6pMyDbPRP4M2l9kxsjYnGtAjUzM7Ouqaq7LmdJyB+yzczMzKzdNOJdl83MzMyWcrJiZmZmDa2qy0BmZmZWH01NzS5zVuNz7dRh52qOe1bMzMysoTlZMTMzs4bmZMXMzMwampMVMzMza2hOVszMzKyh5UpWJK3ZXoGYmZmZlZO3Z+UtSVdLGtkewZiZmZmVypusrA58C3hQ0gxJP5O0QTvEZWZmZgbkT1a+ClwHLAQ2Ac4EXpJ0r6QDJHWvdYBmZmbWteVKViJickT8O7AecAzwFNAN2AO4EXhT0q+zuzKbmZmZtVlVs4EiYl5EXBoRw4GhwPnAO8A6wA+AP0t6RtKxkvrULlwzMzPrato8dTkipkfEScCGwD7AHcBiYBvgItKg3Gs9KNfMzMyqUbN1ViJiCXAfcAPwl6xYwGrAoaRBuX+WtEutzmlmZmYrv5okK5K2kfQb4C3gemB74BPgFuDfgCuA+cC/Ag9I2qsW5zUzM7OVX9XJiqR1JJ0g6a/AM8BxQD/gH8CPgQ0i4qCIuCEijgI2AK7Kznla20M3MzOzrmDVPJUlrQLsCRwBfB3oTrrUswC4Gbg8Ih4v1zYi5ko6BjgI2LItQZuZmVnXkStZAV4DBpESFIA/A5cD10bEvJYaR8Qnkt4jDcY1MzMza1HeZGU94APSwnCXRcRfqzjnicAaVbQzMzOzLihvsvJt4OaIWFTtCSPi1mrbmpmZWdeTd4DtI6SF31pF0vqSNsp5DjMzM7Ol8vasvEyanvyZVtZ/nDQ+Je95zMzMzIDqpi6r5Sptqm9mZma2VM1WsK1gddLS+2ZmZmZVabdkRdJmQH/g7fY6h5mZma38mh1LImlvYO+S4j6SrmyuGbA2sHP2/OHqwzMzM7OurqWBr/8KHF5S1qtMWSX/BH6eLyQzMzOzZVpKViaXPD8d+BA4v5k2nwJzgeeAyRHhMStmZmZWtWaTlYh4hLS2CgCSTgc+jIgz2jswMzMzM8i//skQYEl7BGJmZmZWTq7ZQBHxSkS83l7BFEgaKulBSQskvSnpTEndWmjTQ9K5kh6VtFBSNFN3b0l/k7RI0jRJB9f+VZiZmVkttPc6K7lJ6gtMAoI0E+lM0s0PW7r0tDpwJLAA+GMzx98ZuJU0S2lP4B7gekm7tTl4MzMzq7mKl4EkFS73TI+ILUvK8oiIyHO56RjSjKP9ImIu8ICktYAmSedkZeVOMkdSv4gISd8Hvlrh+D8HpkTE8dnzhyVtCZwG/G+OOM3MzKwDNNezoqKtXFlrt7y9N3sCE0uSkhtICcyI5hpGRMVLPwCSegK7ADeV7LoBGC6pT85YzczMrJ011+MxJHv8pExZe9oCeKi4ICJelbQg23dXG469KdAdmF5S/ndSUvU54Ok2HN/MzMxqrGKyEhGvtKasHfQF5pQpn53ta+uxKXP82SX7W6W5QbxmZmZWGw03wDZTLglQhfJaHF8Vys3MzKzO8q6z0hFmk+4tVKoP5Xtc8h6bMscvPM91/IhQuXL3uJiZmdVOc7OBDqvVSSLimhzVp5PGphTHsiHQmxXHmuT1T9IYnC0oWpk3e/4p8Hwbj29mK4umppXrPGadWHM9K+Op3WWRPMnKfcBJktaMiHlZ2cHAQpZPMHKLiI8kPQwcCPyuaNfBwNSI+KAtxzczM7Paay5ZmUJ9xnBcAhwP3CbpbGAToAm4oHg6s6QZwCMR8d2isj1JPTD/mj0/INv1dNHg4P8GJkv6NTABGJNte7Tni+oQ/iZoZmYroeZmA43swDiKzztb0ijgItI05TnAWFLCUmxVoHQJ/nHAxkXPb84ejyD1FBERj2VJzC+AY4GXgH+LCC8IZ2Zm1oAacYAtETGNyivQFuoMbk1ZhbYTSL0qZmZ11dT0eAedZ6cOOY9Ze2jUqctmZmZmgJMVMzMza3DNTV0uLHn/SkQcUVKWR0TEqGqCMzMzM2tuzMrI7HF6mbI8vECamZmZVa25ZOWI7PGDMmVmZmZmHaK5qctXt6bMzMzMrD15gK2ZmZk1NCcrZmZm1tCqWhROUjfS/XQOALYFBmS7ZgJ/Bm4Cbo6IJbUI0hpLRy1ilc7lhazMzLq63MmKpM1Jy9hvCahk90bZtjfwU0kHRcQ/2hylmZmZdVm5khVJg0g3OBwAfAzcQroT8hukxGU9YASpx2Vr0g0Dt4mIt2sZtJmZmXUdeXtWziAlKi8CYyLi+TJ1Lpd0JnAv6Y7Jp5NuGGhmZmaWW94BtmNIi7wdUSFRASAiXgC+Q+pt+Xr14ZmZmVlXlzdZ6Q/Mj4hHW6qY1fkwa2NmZmZWlbzJyps523TL2piZmZlVJW+ycifQS9KeLVXM6vQCJlQTmJmZmRnkT1bOAF4CrpS0Y6VKknYArgRmAGdWH56ZmZl1dRVnA0k6rMKui4GfA1MkPQpMJk1dBlifNHV5BDAXOIe05so1NYrXzMzMupjmpi6PJ838KaewGNwI4CsV9vUBzst+drJiZmZmVWkuWZlC5WTFzMzMrENUTFYiYmQHxmFmZmZWlu+6bGZmZg3NyYqZmZk1NCcrZmZm1tDy3sgQAElfAo4BdiJNV+7dTPWIiKrOY2ZmZpY7iZD0X8AvaH2vjFquYmZmZlZerstAknYBfkma0nwasG22ayawGamn5XRgVrbtDQypVbBmZmbW9eQds/IDUqJyekT8IiL+mpUviYgXI2JqRPw38AVgNnAFsLh24ZqZmVlXk/cy0LDs8dKS8uWSnoh4S9L3gAeAU4DjqgvPzFYGTU2Pd+C5duqwc5lZx8jbs9IfmB8Rs4rKFgOrl6n7ELAQaPEOzWZmZmaV5E1WZrNib8xsoLekPsWFERHAp8B61YdnZmZmXV3eZOV1oKekAUVl07LHkcUVJX2BNKV5ftXRmZmZWZeXN1kpXHjerqjsTtL05PMkfUlSd0nbAleTBuM+0vYwzczMrKvKm6zcTkpMvl1UNg54AdgUeAJYBDwN/AtpzEpTm6M0MzOzLitvsjIF2Br4eaEgIhYBI4CbgY9ZtgjcVOCrEfG3GsRpZmZmXVSuqcsR8SnwXJnyt4GDJXUnzRiaGxEeq2JmZmZtVtN79kTEJ8BbtTymmZmZdW2+waBZV9bU1EEnGt1B5zGzlVHeMSsASOom6d8k3SbpZUnzs+3lrOwQSd1qHayZmZl1PdXcdXlz0mDaLVnxjsobZdvewE8lHRQR/2hzlGZmZtZl5UpWJA0izQgaQJr5cwtpHZU3SInLeqSZQQeQZg1NlrRNNgDXzMzMLLe8l4HOICUqLwL/EhHfiojLIuLeiLgnIi6PiH8nrbHyT2Bd4PS8QUkaKulBSQskvSnpzNZcVpLUR9JVkmZL+kDStZLWKakzXlKU2bbIG6eZmZm1v7yXgcaQVqU9IiKer1QpIl6Q9B1Sr8vXgWNbewJJfYFJpGX89yYtNnc+KbE6tYXmNwKbA0eS7kt0NjAB+HJJvenAESVlL7c2RjMzM+s4eZOVwl2XH22pYkQ8KunDrE0exwC9gP0iYi7wgKS1gCZJ52RlK5A0HNgdGBERU7KyN4AnJe3ZuMDHAAAZSUlEQVQaEZOKqs+PiCdyxmVmZmZ1kPcy0Js523TL2uSxJzCxJCm5gZTAjGih3TuFRAUgIp4CXsr2mZmZWSeUN1m5E+glqcUP/6xOL9JlmDy2IF2mWSoiXgUWZPta3S7z9zLthkqaK+kjSY9Jai4JqqjC2Jeo5lhmZmZWXjUDbF8CrpS0Y6VKknYArgRmAGfmPEdfYE6Z8tnZvra2+wtwIrAX8E1S788DkrbPGaeZmZl1gIpjViQdVmHXxaQbGU6R9CgwmTR1GWB90qWaEcBc4BzSINlrcsZVrndCFcpztYuIC5fbKd1DGsx7CrBPriAjSteZKRzTvStmZmY10twA2/FUTg4KH9IjgK9U2NcHOC/7OU+yMhtYu0x5H8r3nBS3G1CmfO3m2kXEQkn3knpazMzMrME0l6xMoeWejPYwnZIxJpI2BHpTfkxKcbvSKcpkx2rNuBn3hpiZmTWgislKRIzswDiK3QecJGnNiJiXlR0MLCSt29Jcu59L2jkiHgOQtB2wSbavLEm9SLOF/lSL4M3MzKy2qrqRYTu7BPgIuE3SrpL+A2gCLiiezixphqQrCs8jYiowEbhG0n6S9gGuBR4rrLGSrXD7qKSjJY2SdDDwMPAZ4Jcd9QLNzMys9XLfyLC9RcRsSaOAi4C7SONNxpISlmKrkmbyFDskq3slKRG7Gzi+aP9HwEzSSrjrAouAqaSF5J6p6QsxMzOzmqg6WZG0CemGhduybGDrTODPwC0R8WK1x46IacBXW6gzuEzZHNIy+qVL6Rf2LwL2qzYuMzMz63i5k5VsjMeFwHdIM39Kp+8eCPxS0uXADyNiYZujNDMzsy4rV7IiaRXgDmAUKUl5g7TOyutZlQ2AkaQxIEcBQyTtERGeaWNmZmZVyduzcgSwK2msxwnA5aWJiCSREpULs7pHkMaQmJmZmeWWdzbQYaT1SI6PiMvK9ZhEcilpYKuAb7c9TDMzM+uq8iYrWwOfAFe3ou7VWd2t8wZlZmZmVpA3WekFLIiIT1qqGBEfA/OzNmZmZmZVyZusvAn0kbRZSxUlfY50X543qwnMzMzMDPInK5NI41B+J2m1SpWyfZeQxrc8UH14ZmZm1tXlTVbOJs0EGgn8P0nHSNpC0pqS+kv6oqQfAy+Q7si8CDinphGbmZlZl5Jr6nJEvCjpIOB6YDPg/1SoKtJ4lUPbspKtmZmZWe4bGUbE3cAXgKuAuSxbxbawfUBaV+ULWV0zMzOzqlV1b6Cst+S7wHezewQtvTeQe1LMzMyslvIut/+N7Mc/RsQsWJq4OEExMzOzdpG3Z2UCsBjo1w6xmJmZma0gb7LyPkBEfNgOsZiZmZmtIO8A2+dIi8Kt1R7BmJmZmZXKm6xcCnQDftAOsZiZmZmtIO86K9dK2h44I1uldmxEvN8+oZmZmZnlnw30UPbjAuAU4GRJM4CZwJIKzSIiRlUfopmZmXVleQfYjizTfotsqyRynsPMzMxsqbzJyhHtEoWZmZlZBXnHrFzdXoGYmZmZlZP73kBmZmZmHalVPSuSegL7AF8E1gLmAE8Cd0XE4vYLz8zMzLq6FpMVSTsCNwODyux+WdI+EfG3mkdmZmZmRguXgSR9BriblKiINLNnZmE3MAS4V1Kf9gzSzMzMuq6WxqycAKxNuuxzGLB6RAwCegPHAwuB9YHvtmeQZmZm1nW1lKyMJvWmHB8Rf4iIjwEiYlFEXAScTuph2a19wzQzM7OuqqVkZRNSsnJrhf03F9UzMzMzq7mWkpU1gZkRsajczoh4Jfuxd02jMjMzM8u0Zp2V1iyXr7YGYmZmZlaOF4UzMzOzhtaaReH6Fd1tuZo6vuuymZmZVa01yUoPVrzbcp46vuuymZmZVa2lZMU3LjQzM7O6ajZZiYgjOioQMzMzs3I8wNbMzMwampMVMzMza2hOVszMzKyhOVkxMzOzhuZkxczMzBpaQyYrkoZKelDSAklvSjpTUrdWtOsj6SpJsyV9IOlaSeuUqbe3pL9JWiRpmqSD2+eVmJmZWVs1XLIiqS8wibSY3N7AmcCJwBmtaH4jaXG6I4HDgS8BE0qOvzPpLtIPA3sC9wDXS9qtJi/AzMzMaqo1K9h2tGOAXsB+ETEXeEDSWkCTpHOyshVIGg7sDoyIiClZ2RvAk5J2jYhJWdWfA1Mi4vjs+cOStgROA/63/V6WmZmZVaPhelZIvR0TS5KSG0gJzIgW2r1TSFQAIuIp4KVsH5J6ArsAN5W0vQEYLqlP28M3MzOzWmrEZGULYHpxQUS8CizI9rW6XebvRe02BbqXqfd30nvxuTyBSopyW55jmJmZWfMU0VifrZI+AU6KiF+XlL8OXBMRp1Ro9wAwPyL2KSn/A7BJROwoaSfgMWCbiPhrUZ3NgBeA3SOi1ZeCWkpMIkKtPZaZmZmV14hjVqD8nZpVobyadqXPVaG8+ZM5GTEzM2t3jXgZaDawdpnyPsCcKtqtXdRudlFZaR1aOL6ZmZnVQSMmK9MpGZsiaUOgN+XHpFRslykey/JP4JMy9bYAPgWeryJeMzMza0eNeBnoPuAkSWtGxLys7GBgIfBIC+1+LmnniHgMQNJ2wCbZPiLiI0kPAwcCvytqezAwNSI+qO1LqcwDcc3MrCuqZghFIw6w7QtMA54FziYlGxcAv46IU4vqzQAeiYjvFpXdT5rR82NST8nZwLsR8eWiOjsDk4GLSAvGjcnq75FncG1bOVkxM7OuaKVIViAtt09KJoaTxpFcDjRFxJKiOi8DkyPi8KKytYGxwL6kS1x3A8dHxKyS4+8D/AL4LGkdlqaIuKEdX1KbFBKblWlAr19T5+DX1Dn4NXUOK9tr6sjX05DJii1vZfsDB7+mzsKvqXPwa+ocVrbX1JGvpxEH2JqZmZkt5WTFzMzMGpqTFTMzM2toTlbMzMysoXmArZmZmTU096yYmZlZQ3OyYmZmZg3NyYqZmZk1NCcrZmZm1tCcrJiZmVlDc7JiZmZmDc3JipmZmTU0JytmZmbW0JysmJmZWUNzsmJmZmYNzclKJyfpPyWFpFvqHUu1JF0iabqkDyXNljRF0q71jqtaktaSdIakpyR9IOltSbdL+ly9Y2sLSQdLuk3SW9nf3OH1jikPSUMlPShpgaQ3JZ0pqVu942oLSZtJ+p2k/ytpiaTJ9Y6pLSQdKOlOSW9k/x/8SdKh9Y6rLSQdIOmPkt6TtEjSPySdKqlHvWOrBUmfyX5XIWmN9jqPk5VOTNK6wGnAzHrH0ka9gIuAfYFvAbOA+yTtUNeoqrcRcBQwETgAOBpYD3hS0ob1DKyNDgAGA3fXOY7cJPUFJgEB7A2cCZwInFHPuGpgS2AM8Hy2dXY/Aj4Efgh8A3gYuE7SD+oaVdusQ3odRwJ7AlcCPwMuqGdQNXQu6XfWrnwjw05M0hVAD2BDYFZEHFDnkGoi+7b7EjAhIo6vdzx5SeoNfBoRC4vK+gGvAudGRKf8gJS0SkR8mn17mgccERHj6xxWq0j6KfATYOOImJuV/QRoAgYVyjqbwu8k+/kWoH9EjKxvVNWT1D8iZpWUXQcMj4ghdQqr5iSdBRwH9I1O/CEs6cvAHcAvSUnLmhHRLomLe1Y6KUlfAg4C/qvesdRaRCwB5pASsU4nIuYXJypZ2fvAK8C69Ymq7Qofip3UnsDEkqTkBlKv3oj6hNR2nfx3soLSRCXzFzrxv5sK3qOT/v9WkH2p/C2pl7Lc762mnKx0QpJEumxyTkS8Ue94akHJqpLWkfRD4LOk7tKVgqQBwGbAtHrH0kVtAUwvLoiIV4EF2T5rXDuyEvy7kdRN0uqSdgaOB8Z15l4V4BhgNeD/dMTJVu2Ik1jNHQEMAs6rdyA1dDBwffbzfODgiHiqjvHU2vmk67o31DuQLqovqbeu1OxsnzUgSaNIY4y+U+9YamA+0DP7+RrgpDrG0iaS1gH+G/hWRHySvj+3LycrDUBSH9IAzGZFxPSs7i+B40svNTSKPK+n6OlE4EtAf+CbwA2SxkTE5HYJMqcqX1Oh7bGkgcP7R8R77RBeVdrymjqpct9iVaHc6kzSYOA64I7OMjaqBTsCqwPbkyZGXAR8r64RVe8s4MmIuLejTuhkpTEcCFzWinoCTgFeA/5X0tpZ+apA9+z5vGzMRz3leT0ARMRs4Jns6f2S1iddC/1K7cOrSu7XBCDpG6TruidHxO3tEVgbVPWaOqnZwNplyvtQvsfF6igbkH4faVD6t+ocTk1ExJ+zHx+TNAu4WtL5EfHPesaVl6QtST1dXyn6DFo9e+wjaUl7fJH2mJUGEBGXR4Ra2rLqmwPbkf7zLWw7kab5zQaG1+VFFMn5eir5C7BJR8TbGtW8Jkk7ki77XBIR59Yn8spq9HvqLKZTMjYlm0bem5KxLFZfklYnTY/vAXwtIubXOaT2UEhcOuMMp88C3YGpLPsMKoxbeZ305azm3LPS+ZwK/Lqk7NfAB8DpwN86PKIaywYQDydNX+6Usm8fdwP3kwbTWX3dB5wkac2ImJeVHQwsBB6pX1hWTNKqwM2kD8SdIuLdOofUXnbKHjvj/3GPAbuUlO0BnExa8+fF9jipk5VOJiKeLS2TNIe0zsrkjo+obbJ5+icCt5G6fNcBvg3sAOxVx9Cqli3Wdz9pQO1vgO2LBqDNjYhOObNB0lBgKGkGAMB2kj4EZkZEo3/gX0JKGm+TdDap164JuKCzrrECS3shxmRPPwOsJamw3tK9EbGgPpFV7WLS6zkB6FeyMORfIuKj+oRVPUn3kxYkfA5YQkpUTgRu7GyXgGDp9PLJxWXZ+CKAR9trnRUnK1ZvrwGLSYOGB5BW4/0rsHNETK1nYG0wFNgg+/nhkn2PACM7NJraOYjUe1dwXLY1/GuKiNnZzJKLgLtI41TGkhKWzmxdUk9EscLzIcDLHRpN2+2WPV5YZl9nfD0ATwOHk1Z/XkzqefgpKYG2VvIKtmZmZtbQPMDWzMzMGpqTFTMzM2toTlbMzMysoTlZMTMzs4bmZMXMzMwampMVMzMza2hOVsysZiRFtg2udyxmtvLwonBmdZAtK/4t4BDgC6SVe+cDb5MWjZoCPBQRT9ctyCKSRpIWfvtrREyobzSNIUvIDgfmRETpLTDMrIa8KJxZB5M0ALiXdEPKgkXAR8BaLLvL8QcRUe5OwR1OUhNp9dqrI+LwZuoVbgo4KiLe6IDQ6iZL4B4GXomIwfWNxmzl5stAZh3vD6REZR7wE2C9iOiVJSZ9gNGke6TMqV+I1YmILbJtpU5UzKxj+TKQWQeStAXL7n/ynYi4pXh/dkfgScAkST/u6PjMzBqRe1bMOtbWRT/f3VzFiFhYaZ+kNSSdIulpSR9IWiTpBUm/kbRhhTaTs8Gvh0vqJalJ0j8kLZT0rqQbJH22pM1gScGyGxh+u2gQ7QqDaSsNsM3OFZLGKzlO0l8kzZf0lqSrJW1QVP+zWdnr2Wt7VtJRzb1fklaR9O+SHpA0U9LHkt6UdKOkYRXaLI0re/5tSU9KmidprqSHJY0u0+5llt2kcuMy78nhzcVacqwhksZJej77XSyQ9Er2+/qppP4V2m0o6fzsvZmXbdMkXSFplwptNpX0O0kvZu/rbElTJB0pqVuFNsV/N2tLOlvS9CzOFXr/JG0l6UpJL2XnmCPpcUnHSOre2vfFbDkR4c2btw7agAOByLZNqzzG50l3ny0c5xPgw6Ln7wM7lWk3Odt/PPDn7OdFwIKitu8VxwVsSBr0Wzj+wux58bZhUf3CcQaXnLspKx8P3JD9/FFJ3C+S7ry9AzA7K5sDfFpU56QK78mawANF9T4FPih6vgT4fpl2xXFdnv28uEzb/UvaPZ29z4X9pe/Jwa38XW4LzC0618dFr72w7VGm3f4lv7eFpMuKhecvl2nz9axeoc6c7HyF5w8AvZv5uzkJ+GfR381c0uDi4rrfz96PwjE/zN7PwvOHgdXr/e/QW+fb6h6AN29daQM2KfqP+35gQM72fYCXsva3A9sAq2b7BgPXZPveBtYuaVv40JmdHWN3oBuph/XLwGvZ/pvKnHfph3oL8bWUrMzJPlS/CfQgDSb+MvBWtv8SUiJ2F7BJ1nYtYFzRh/I6Zc57e7b//wJjgF5Z+drAT0mJ0RJKkriiuGZnxz6m8GEKDAEeyfa/WXifi9qOrJQY5Ph9PpQd4wlgm6Ly1UnjmsYCw0vaDCclqJG1/xLLJksMAPYBrixpsynLEsPJwOZZeU/gP0jJRwCXl4mx8HczD3gV2ANYJdu3WVG9vVmWoPwUWDcr704ahzU92/+7ev879Nb5troH4M1bV9uAq4s+1D8ijVH5RfaffbPJS1YvgAmFD6gyde7J6vy4pLzwobOg+EOmaP/+LPvW3KNkX+FDfXwL8bWUrATw7TLt/r1o//NlEoNVgBey/YeV7Ns1K38J6Fchrp9kde5uJq5vlmm3XvY7CuArJftG0vZkpdA7MixHmyezNo8A3VvZ5oqszQzK9GxkCUuhR2qzkn2Fv5uPga0qHL8by3r79q1QZ0iWyHxCGlRe93+L3jrP5jErZh3vKOAC0n/+PYBRwM9ICci7kp6S9E1JKtP229nj2IiICse/PntcYaxF5paImFGm/E7Sh01PYLOWX0ZVXgd+X6Z8UtHP50bE4uKdEfEpy8aIbFXStvCejI+I9yuc97rscZcKYzNeLapTfN63gKcqnLcW5maP67WmcjZAe/vs6U8i4pNWtBEpEYX0d7OgTLXLgTdIPV0HVDjUfRHxbIV9I4GNSYnb7eUqRMRLpB6kVbP6Zq3m2UBmHSwiPgZOlHQ2sC8wgtTlvxnpw+JLpOnNe0s6JPugJhs4WxiEerOkTyucokf2WHagLWm8Rbm4PpH0LjAQ6JvvVbXatMLrKfFu0c+VPhDfyR5LY9sxe/yhpGNbOP/qpAX43i0pf6aZ5K8wDbs93pN7gSOAayRdTEpY/9RMErJD9vh+RDzZynNsQrp8CMsSvuVExKeSJpMuz21b4ThTmzlH4XewvqS3m6lXiKPS36ZZWU5WzOokIt4FfpdtSBoI7AWcRvrP/EDgceDCrEnxt+8BrTjF6hXK5zXTZlH22F6zNt4qVxgRS4o6ksrWIY05gRVjK7wvfVj2Ydiccu9Lvd6Tk4DNSR/2J2fbIklTgZtJvUXFs8IGZo+v5jhH8d9Kc+vfvF6mfrGZzbQt/A56sCzG5lT62zQry5eBzBpERLwTEZeTvtkWehG+U1Sl+N9rn4hQC9vgjoq9zgrvy96teE8UES/XM9hiEfEesDPpkt1vgL+QPvB3IS0M+GzxlG6WrW5crZ5taLukmX2F38HtrfwdNLUhDuuCnKyYNZiImAXckT39XNGud4p+HtpxETW8wvvSKd+TSCZFxAkRsS3QHziaNDV6E9KMoILCJZaNcpyiuEdk42bqFZKi5npQKunUvwNrfE5WzBrT/Ozx40JBNkCx8KGwXwfHUxhn0tZv9u2hMJZi/2Zr1V67vCcRMTsiLgVOyYpGFO1+InvsJ2kHWudFlt26odJicauwbNDrn1sf7VKF38Hmkrasor1Zs5ysmHWgbLXSTVuoszpprQyAv5bsHp89fk/S55s5hiS1ZvxGaxVmrTTEjRVLjM8et5N0WHMVJdVykGzhPanqfc5W3G1u3GBhrMrSSzcRMZ1ls5POac2KsNnA4duypydkf1+ljgQ+Q5oNdkuZ/S15kGXjaMZWWg0Xav47sC7CyYpZx9oS+Iek2yQdJGnpoFlJvSXtBTxKWpMClg2uLfgf0jfl3sAj2fLwaxQdY8NsWfo/kWYa1cpz2ePOKlmSv94i4n6WfRhfKemMkve1r6S9Jd1BmjJeKy+Q1gzpI6maXp21gBmSfiZp68IHfJbEjALOyupNLGn3I9KqsF8G7pe09O7dkvpLOkTStSVtfknqrVsfuEfS5ln9ntnfy2+yeldUmNberGz20g9Iyc5o4H8lDStMv5e0qqQvSir8/Zrl4tlAZh3rE9ICWvtmG5IWki73FH9DXwKcFhG3FTeOiDmSdietifJ5Uq/Cldk9Wnpl29LqNYx7Mmmp9U1JydYs0oJmADtHxOuVGnaQw0hfvvYhzaY6TdIHpEs0axXVG1+rE0bEfEnXZ+e+JTtf4XLLj6PkJpUVbExa6O8XwCeS5pH+Dgo9Ey+SkpPi8z4u6d9Jr+WrwNPZ39ASoJC4vlLS5p+SDgVuIl3umZ79zfRm2SynB4H/bEXMZUXEnZK+S1qF+KukS1aLJM0n9chV7G0xa4mTFbMOFBETs2+1e5FmgWxF6n5fg/RB9yIwhbTs+XMVjjFD0jakmUIHkm6OuDbpssH/I61segvwWA3j/iT7tv/fpHEPg1g2xbXu/49ExHxgX0lfI70vw0jxfUpatfUpUu/LvTU+9TGk6cD7kRKPwgDWNSq2WGYu6X49u5KmLm9Aink+8A/Smiu/jXQn7uVExA2SniQlMruRprp/CvydNN19hYX3IuIuSVuTVvMdTeplWQD8jXSbhisjorkZPy2KiKskPQyckJ1jY1Ly9R6pd+4eqrvMZF1c4X4SZmZmZg3JY1bMzMysoTlZMTMzs4bmZMXMzMwampMVMzMza2hOVszMzKyhOVkxMzOzhuZkxczMzBqakxUzMzNraE5WzMzMrKE5WTEzM7OG9v8BaXi3SMj3CDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x324 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure( figsize = (8, 4.5) )\n",
    "sub1 = fig.add_subplot(1,1,1)\n",
    "my_font_size = 15\n",
    "\n",
    "half_frame(sub1, \"Sentiment score\", \"Probability density\", font_size = my_font_size)\n",
    "\n",
    "# Calculate and plot histogram\n",
    "#\n",
    "sub1.hist([othello_sentiment, iago_sentiment], 6, normed = 1, color = ['r', 'darkblue'], rwidth = 0.75, alpha = 0.5, \n",
    "          histtype = \"bar\", label = ['Othello', 'Iago'], cumulative = False)\n",
    "\n",
    "# Format legend\n",
    "sub1.legend(loc = \"best\", frameon = False, markerscale = 1.8, fontsize = my_font_size)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, are these distributions different? Let's look at some statistical metrics that might tell us something about this..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:34:51.410695Z",
     "start_time": "2019-07-03T11:34:51.401100Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DescribeResult(nobs=470, minmax=(-4, 4), mean=0.4085106382978723, variance=5.120609717370594, skewness=-0.41543139442792626, kurtosis=-1.3060163694576852)\n",
      "\n",
      "DescribeResult(nobs=537, minmax=(-4, 4), mean=0.32774674115456237, variance=4.907307040217906, skewness=-0.29978164467913276, kurtosis=-1.4307035152922654)\n"
     ]
    }
   ],
   "source": [
    "results_othello = stat.describe(othello_sentiment)\n",
    "print(results_othello)\n",
    "print()\n",
    "results_iago = stat.describe(iago_sentiment)\n",
    "print(results_iago)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:35:25.099636Z",
     "start_time": "2019-07-03T11:35:25.093633Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ttest_ind: \t t =  0.5706  \t p = 0.5684\n"
     ]
    }
   ],
   "source": [
    "# Using t-test\n",
    "\n",
    "t, p = stat.ttest_ind(othello_sentiment, iago_sentiment, equal_var = False)\n",
    "print('ttest_ind: \\t t =  {0:6.4f}  \\t p = {1:6.4f}'.format(t, p) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:35:41.691341Z",
     "start_time": "2019-07-03T11:35:41.685128Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KS 2 samp: \t  D =  0.0272  \t p = 0.9916\n"
     ]
    }
   ],
   "source": [
    "# Using Kolmogorov-Smirnov test\n",
    "\n",
    "D, p = stat.ks_2samp(othello_sentiment, iago_sentiment)\n",
    "print('KS 2 samp: \\t  D =  {0:6.4f}  \\t p = {1:6.4f}'.format(D, p) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:35:43.542459Z",
     "start_time": "2019-07-03T11:35:43.531193Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 9.57446809 14.04255319  5.90425532  9.73404255 21.70212766 14.04255319]\n",
      "[ 8.51955307 16.34078212  6.70391061  9.3575419  20.81005587 13.26815642]\n",
      "[0.13062254 0.32323154 0.09538441 0.0151485  0.03824075 0.04519771]\n",
      "[0.13062254 0.45385409 0.5492385  0.564387   0.60262775 0.64782546]\n",
      "0.0 nan\n"
     ]
    }
   ],
   "source": [
    "# Using chi_square test\n",
    "\n",
    "hist_othello, edges = np.histogram(othello_sentiment, bins = 6, density = True)\n",
    "hist_iago, edges = np.histogram(iago_sentiment, bins = 6, density = True)\n",
    "\n",
    "hist_othello = 100 * np.array(hist_othello)\n",
    "hist_iago = 100 * np.array(hist_iago)\n",
    "print(hist_othello)\n",
    "print(hist_iago)\n",
    "\n",
    "b = (hist_othello - hist_iago)**2 / hist_iago\n",
    "print(b)\n",
    "print(np.cumsum( b ))\n",
    "\n",
    "chi2 , p = stat.chisquare(hist_othello, hist_othello, ddof = 5)\n",
    "print(chi2, p)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-03T11:36:00.988213Z",
     "start_time": "2019-07-03T11:36:00.980971Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 range(1, 7)\n",
      "\n",
      "52 ['2 Clubs', '2 Diamonds', '2 Hearts', '2 Spades', '3 Clubs', '3 Diamonds', '3 Hearts', '3 Spades', '4 Clubs', '4 Diamonds', '4 Hearts', '4 Spades', '5 Clubs', '5 Diamonds', '5 Hearts', '5 Spades', '6 Clubs', '6 Diamonds', '6 Hearts', '6 Spades', '7 Clubs', '7 Diamonds', '7 Hearts', '7 Spades', '8 Clubs', '8 Diamonds', '8 Hearts', '8 Spades', '9 Clubs', '9 Diamonds', '9 Hearts', '9 Spades', '10 Clubs', '10 Diamonds', '10 Hearts', '10 Spades', 'Jack Clubs', 'Jack Diamonds', 'Jack Hearts', 'Jack Spades', 'Queen Clubs', 'Queen Diamonds', 'Queen Hearts', 'Queen Spades', 'King Clubs', 'King Diamonds', 'King Hearts', 'King Spades', 'Ace Clubs', 'Ace Diamonds', 'Ace Hearts', 'Ace Spades']\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "# Playing games of chance\n",
    "\n",
    "# Define outcomes of play\n",
    "one_die = range(1, 7)\n",
    "\n",
    "cards = []\n",
    "card_values = ['2', '3', '4', '5', '6', '7', '8', '9', '10', \n",
    "               'Jack', 'Queen', 'King', 'Ace']\n",
    "suites = ['Clubs', 'Diamonds', 'Hearts', 'Spades']\n",
    "for value in card_values:\n",
    "    for suit in suites:\n",
    "        cards.append(value + \" \" + suit)\n",
    "\n",
    "        \n",
    "no_one_die = len(one_die)\n",
    "no_cards = len(cards)\n",
    "\n",
    "print( no_one_die, one_die )\n",
    "print()\n",
    "print( no_cards, cards )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-06-23T21:58:56.601711Z",
     "start_time": "2019-06-23T21:58:56.594525Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16\n",
      "4481381406320.0\n"
     ]
    }
   ],
   "source": [
    "#How many possibilities\n",
    "def number_states_dice(n, no_one_die = 6):\n",
    "    \"\"\" Returns number of possible sums for n dice\"\"\" \n",
    "    return 6*n - (n-1)\n",
    "\n",
    "def number_states_cards(n, no_cards = 52):\n",
    "    \"\"\" Returns number of possible hands of n cards from a single deck\"\"\" \n",
    "    number = 1\n",
    "    for i in range(no_cards, no_cards-n, -1):\n",
    "        number *= i\n",
    "    return number / math.factorial(n)\n",
    "\n",
    "\n",
    "print( number_states_dice(3) )\n",
    "print( number_states_cards(15) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-06-23T21:59:00.106775Z",
     "start_time": "2019-06-23T21:58:57.160023Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1\n",
      "0.07\n",
      "0.057\n",
      "0.055088\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def prob_of_outcome_tossing_n_dice(sum_dots, n, n_attempts = 1000):\n",
    "    \"\"\" Probability of the sum of the dots in n dice being sum_dots\"\"\"\n",
    "    count_successes = 0\n",
    "    for a in range(n_attempts):\n",
    "        # Toss dice\n",
    "        sum_temp = 0\n",
    "        for i in range(n):\n",
    "            sum_temp += np.random.randint(6) + 1\n",
    "        if sum_temp == sum_dots:\n",
    "            count_successes += 1\n",
    "    \n",
    "    return float(count_successes) / n_attempts\n",
    "\n",
    "print( prob_of_outcome_tossing_n_dice(3, 2, 10) )\n",
    "print( prob_of_outcome_tossing_n_dice(3, 2, 100) )\n",
    "print( prob_of_outcome_tossing_n_dice(3, 2) )\n",
    "print( prob_of_outcome_tossing_n_dice(3, 2, 1000000) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-06-23T21:59:17.262657Z",
     "start_time": "2019-06-23T21:59:00.110075Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.02334\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEKCAYAAADaa8itAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHodJREFUeJzt3XuUHWWZ7/Hvz4SAqJBgGsQk2sGJOMHjJbYYx8sgOCGAh6ALNIxzyNGMGTVe57A0jC7ggKwD3nA4KkyUHAKDhIgXciQYIkTRdSRJc0+4mDZkSJtI2km4KBoMPOePejcWO3t3qrtr78qG32etvXbVU29VPVW9u5+uy35LEYGZmVkZnld1AmZm9uzhomJmZqVxUTEzs9K4qJiZWWlcVMzMrDQuKmZmVhoXFTMzK42LipmZlcZFxczMSjO66gTabfz48dHd3V11GmZmHeXWW2/9XUR07andc66odHd309vbW3UaZmYdRdJ/FGnn019mZlYaFxUzMyuNi4qZmZXGRcXMzErjomJmZqVxUTEzs9K4qJiZWWlcVMzMrDQuKmZmVprn3DfqbWi6F1xXyXo3nX9CJes1s5Fp2ZGKpEWStklaVxf/uKT7Ja2X9MVc/AxJfWnasbn4zBTrk7QgF58sabWkDZKuljSmVdtiZmbFtPL012XAzHxA0juAWcBrIuII4MspPhWYDRyR5vmmpFGSRgHfAI4DpgKnprYAFwAXRsQUYAcwt4XbYmZmBbSsqETEzcD2uvBHgPMjYmdqsy3FZwFLImJnRDwA9AFHpldfRGyMiCeAJcAsSQKOBq5J8y8GTmrVtpiZWTHtvlD/SuBt6bTVzyS9McUnAJtz7fpTrFn8xcDDEbGrLm5mZhVq94X60cA4YDrwRmCppMMANWgbNC56MUj7hiTNA+YBvOxlLxtiymZmVlS7j1T6ge9HZg3wFDA+xSfl2k0EtgwS/x0wVtLounhDEbEwInoioqera4/PmDEzs2Fqd1H5Idm1ECS9EhhDViCWAbMl7StpMjAFWAOsBaakO73GkF3MXxYRAawCTk7LnQNc29YtMTOz3bTs9Jekq4CjgPGS+oGzgEXAonSb8RPAnFQg1ktaCtwD7ALmR8STaTkfA1YAo4BFEbE+reKzwBJJXwBuBy5t1baYmVkxLSsqEXFqk0n/0KT9ecB5DeLLgeUN4hvJ7g4zM7O9hLtpMTOz0riomJlZaVxUzMysNC4qZmZWGhcVMzMrjYuKmZmVxkXFzMxK46JiZmalcVExM7PSuKiYmVlpXFTMzKw0LipmZlYaFxUzMyuNi4qZmZXGRcXMzErT7mfUmxXSveC6yta96fwTKlu3WafzkYqZmZWmZUVF0iJJ29Kjg+unnS4pJI1P45J0kaQ+SXdJmpZrO0fShvSak4u/QdLdaZ6LJKlV22JmZsW08kjlMmBmfVDSJODvgAdz4eOAKek1D7g4tT2I7Nn2byJ7dPBZksaleS5ObWvz7bYuMzNrr5YVlYi4GdjeYNKFwGeAyMVmAZdH5hZgrKRDgWOBlRGxPSJ2ACuBmWnaARHxy4gI4HLgpFZti5mZFdPWayqSTgR+ExF31k2aAGzOjfen2GDx/gbxZuudJ6lXUu/AwMAItsDMzAbTtqIiaX/gc8CZjSY3iMUw4g1FxMKI6ImInq6uriLpmpnZMLTzSOUVwGTgTkmbgInAbZJeQnakMSnXdiKwZQ/xiQ3iZmZWobYVlYi4OyIOjojuiOgmKwzTIuK3wDLgtHQX2HTgkYjYCqwAZkgaly7QzwBWpGmPSZqe7vo6Dbi2XdtiZmaNtfKW4quAXwKHS+qXNHeQ5suBjUAf8C3gowARsR04F1ibXuekGMBHgG+neX4NXN+K7TAzs+Ja9o36iDh1D9O7c8MBzG/SbhGwqEG8F3j1yLI0M7My+Rv1ZmZWGhcVMzMrjYuKmZmVxkXFzMxK46JiZmalcVExM7PSuKiYmVlpXFTMzKw0LipmZlYaFxUzMyuNi4qZmZXGRcXMzErjomJmZqVxUTEzs9K4qJiZWWlcVMzMrDStfPLjIknbJK3Lxb4k6T5Jd0n6gaSxuWlnSOqTdL+kY3PxmSnWJ2lBLj5Z0mpJGyRdLWlMq7bFzMyKaeWRymXAzLrYSuDVEfEa4FfAGQCSpgKzgSPSPN+UNErSKOAbwHHAVODU1BbgAuDCiJgC7AAGe1yxmZm1QcuKSkTcDGyvi90QEbvS6C3AxDQ8C1gSETsj4gGy584fmV59EbExIp4AlgCzJAk4Grgmzb8YOKlV22JmZsVUeU3lg8D1aXgCsDk3rT/FmsVfDDycK1C1eEOS5knqldQ7MDBQUvpmZlavkqIi6XPALuDKWqhBsxhGvKGIWBgRPRHR09XVNdR0zcysoNHtXqGkOcC7gGMiolYI+oFJuWYTgS1puFH8d8BYSaPT0Uq+vZmZVaStRyqSZgKfBU6MiMdzk5YBsyXtK2kyMAVYA6wFpqQ7vcaQXcxflorRKuDkNP8c4Np2bYeZmTXWyluKrwJ+CRwuqV/SXODrwIuAlZLukHQJQESsB5YC9wA/BuZHxJPpKORjwArgXmBpagtZcfpnSX1k11gubdW2mJlZMS07/RURpzYIN/3DHxHnAec1iC8HljeIbyS7O8zMzPYS/ka9mZmVxkXFzMxK46JiZmalcVExM7PSuKiYmVlpXFTMzKw0LipmZlYaFxUzMyuNi4qZmZXGRcXMzErjomJmZqVxUTEzs9K4qJiZWWkKFRVJr251ImZm1vmKHqlcImmNpI9KGtvSjMzMrGMVKioR8Vbg/WSP9u2V9B1Jf9fSzMzMrOMUvqYSERuAz5M9cfFvgYsk3SfpPY3aS1okaZukdbnYQZJWStqQ3seluCRdJKlP0l2SpuXmmZPab0jPt6/F3yDp7jTPRZI09M03M7MyFb2m8hpJF5I90vdo4L9GxF+n4QubzHYZMLMutgC4MSKmADemcYDjyJ5LPwWYB1yc1nsQcBbwJrKnPJ5VK0SpzbzcfPXrMjOzNit6pPJ14DbgtRExPyJuA4iILWRHL7uJiJuB7XXhWcDiNLwYOCkXvzwytwBjJR0KHAusjIjtEbEDWAnMTNMOiIhfRkQAl+eWZWZmFSn6jPrjgT9GxJMAkp4H7BcRj0fEFUNY3yERsRUgIrZKOjjFJwCbc+36U2yweH+DuJmZVajokcpPgOfnxvdPsbI0uh4Sw4g3Xrg0T1KvpN6BgYFhpmhmZntStKjsFxG/r42k4f2Hsb6H0qkr0vu2FO8nu7OsZiKwZQ/xiQ3iDUXEwojoiYierq6uYaRtZmZFFC0qf6i7I+sNwB+Hsb5lQO0OrjnAtbn4aekusOnAI+k02QpghqRx6QL9DGBFmvaYpOnprq/TcssyM7OKFL2m8ingu5JqRwOHAu8bbAZJVwFHAeMl9ZPdxXU+sFTSXOBB4JTUfDnZdZs+4HHgAwARsV3SucDa1O6ciKhd/P8I2R1mzweuTy8zM6tQoaISEWslvQo4nOx6xn0R8ec9zHNqk0nHNGgbwPwmy1kELGoQ7wXcfYyZ2V6k6JEKwBuB7jTP6yUREZe3JCszM+tIhYqKpCuAVwB3AE+mcO37IWZmZkDxI5UeYGo6TWVmZtZQ0bu/1gEvaWUiZmbW+YoeqYwH7pG0BthZC0bEiS3JyszMOlLRonJ2K5MwM7Nnh6K3FP9M0suBKRHxE0n7A6Nam5qZmXWaol3ffwi4Bvi3FJoA/LBVSZmZWWcqeqF+PvAW4FF4+oFdBw86h5mZPecULSo7I+KJ2oik0QzSK7CZmT03FS0qP5P0L8Dz07Ppvwv839alZWZmnahoUVkADAB3A/9E1gFkwyc+mpnZc1fRu7+eAr6VXmZmZg0V7fvrARpcQ4mIw0rPyMzMOtZQ+v6q2Y/sOSgHlZ+OmZl1sqKnv/6zLvQ1Sb8Aziw/JavXveC6qlMwMyuk6Jcfp+VePZI+DLxouCuV9GlJ6yWtk3SVpP0kTZa0WtIGSVdLGpPa7pvG+9L07txyzkjx+yUdO9x8zMysHEVPf30lN7wL2AS8dzgrlDQB+ARZV/p/lLQUmE32OOELI2KJpEuAucDF6X1HRPyVpNnABcD7JE1N8x0BvBT4iaRXRsSTDVZrZmZtUPT01ztasN7nS/ozsD+wFTga+Ps0fTFZJ5YXA7P4S4eW1wBfl6QUXxIRO4EHJPUBRwK/LDlXMzMrqOjdX/882PSI+GrRFUbEbyR9GXgQ+CNwA3Ar8HBE7ErN+sn6FyO9b07z7pL0CPDiFL8lt+j8PGZmVoGiX37sAT5C9kd7AvBhYCrZdZUhXVuRNI7sKGMy2WmrFwDHNWhau4VZTaY1izda5zxJvZJ6BwYGhpKumZkNwVAe0jUtIh4DkHQ28N2I+MdhrPOdwAMRMZCW9X3gb4Cxkkano5WJwJbUvh+YBPSnPscOBLbn4jX5eZ4hIhYCCwF6enrcZ5mZWYsUPVJ5GfBEbvwJoHuY63wQmC5p/3Rt5BjgHmAVcHJqMwe4Ng0vS+Ok6TdFRKT47HR32GRgCrBmmDmZmVkJih6pXAGskfQDslNM7wYuH84KI2K1pGuA28juJLud7CjiOmCJpC+k2KVplkuBK9KF+O1kd3wREevTnWP3pOXM951fZmbVKnr313mSrgfelkIfiIjbh7vSiDgLOKsuvJHs7q36tn8i+wZ/w7yA84abh5mZlavo6S/Ibv19NCL+lez6xuQW5WRmZh2q6DfqzwI+C5yRQvsA/96qpMzMrDMVPVJ5N3Ai8AeAiNjCCLppMTOzZ6eiReWJdMdVAEh6QetSMjOzTlW0qCyV9G9k3yX5EPAT/MAuMzOrU/Tury+nZ9M/ChwOnBkRK1uamZmZdZw9FhVJo4AVEfFOwIXEzMya2mNRiYgnJT0u6cCIeKQdSZlVqaqHom06/4RK1mtWpqLfqP8TcLeklaQ7wAAi4hMtycrMzDpS0aJyXXqZmZk1NWhRkfSyiHgwIha3KyEzM+tce7ql+Ie1AUnfa3EuZmbW4fZUVPIPwjqslYmYmVnn21NRiSbDZmZmu9nThfrXSnqU7Ijl+WmYNB4RcUBLszMzs44yaFGJiFHtSsTMzDrfUJ6nUhpJYyVdI+k+SfdKerOkgyStlLQhvY9LbSXpIkl9ku6SNC23nDmp/QZJc5qv0czM2qGSogL8K/DjiHgV8FrgXmABcGNETAFuTOMAx5E9f34KMA+4GEDSQWRPj3wT2RMjz6oVIjMzq0bbi4qkA4C3k55BHxFPRMTDwCyg9n2YxcBJaXgWcHlkbiHrKflQ4FhgZURsj4gdZP2SzWzjppiZWZ0qjlQOAwaA/yPpdknfTs9nOSQitgKk94NT+wnA5tz8/SnWLG5mZhWpoqiMBqYBF0fE68n6ElswSHs1iMUg8d0XIM2T1Cupd2BgYKj5mplZQVUUlX6gPyJWp/FryIrMQ+m0Ful9W679pNz8E4Etg8R3ExELI6InInq6urpK2xAzM3umtheViPgtsFnS4Sl0DHAPsAyo3cE1B7g2DS8DTkt3gU0HHkmnx1YAMySNSxfoZ6SYmZlVpGgvxWX7OHClpDHARuADZAVuqaS5wIPAKantcuB4oA94PLUlIrZLOhdYm9qdExHb27cJZmZWr5KiEhF3AD0NJh3ToG0A85ssZxGwqNzszMxsuKr6noqZmT0LuaiYmVlpXFTMzKw0LipmZlYaFxUzMyuNi4qZmZXGRcXMzErjomJmZqVxUTEzs9K4qJiZWWlcVMzMrDQuKmZmVhoXFTMzK42LipmZlcZFxczMSuOiYmZmpXFRMTOz0lRWVCSNknS7pB+l8cmSVkvaIOnq9KhhJO2bxvvS9O7cMs5I8fslHVvNlpiZWU2VRyqfBO7NjV8AXBgRU4AdwNwUnwvsiIi/Ai5M7ZA0FZgNHAHMBL4paVSbcjczswYqKSqSJgInAN9O4wKOBq5JTRYDJ6XhWWmcNP2Y1H4WsCQidkbEA0AfcGR7tsDMzBqp6kjla8BngKfS+IuBhyNiVxrvByak4QnAZoA0/ZHU/ul4g3meQdI8Sb2SegcGBsrcDjMzy2l7UZH0LmBbRNyaDzdoGnuYNtg8zwxGLIyInojo6erqGlK+ZmZW3OgK1vkW4ERJxwP7AQeQHbmMlTQ6HY1MBLak9v3AJKBf0mjgQGB7Ll6Tn8fMzCrQ9iOViDgjIiZGRDfZhfabIuL9wCrg5NRsDnBtGl6WxknTb4qISPHZ6e6wycAUYE2bNsPMzBqo4kilmc8CSyR9AbgduDTFLwWukNRHdoQyGyAi1ktaCtwD7ALmR8ST7U/bzMxqKi0qEfFT4KdpeCMN7t6KiD8BpzSZ/zzgvNZlaGZmQ+Fv1JuZWWlcVMzMrDQuKmZmVhoXFTMzK42LipmZlcZFxczMSuOiYmZmpXFRMTOz0riomJlZaVxUzMysNC4qZmZWGhcVMzMrjYuKmZmVxkXFzMxK46JiZmalcVExM7PStL2oSJokaZWkeyWtl/TJFD9I0kpJG9L7uBSXpIsk9Um6S9K03LLmpPYbJM1ptk4zM2uPKo5UdgH/IyL+GpgOzJc0FVgA3BgRU4Ab0zjAcWTPn58CzAMuhqwIAWcBbyJ7YuRZtUJkZmbVaHtRiYitEXFbGn4MuBeYAMwCFqdmi4GT0vAs4PLI3AKMlXQocCywMiK2R8QOYCUws42bYmZmdSp9Rr2kbuD1wGrgkIjYClnhkXRwajYB2JybrT/FmsXNOlL3gusqW/em80+obN327FLZhXpJLwS+B3wqIh4drGmDWAwSb7SueZJ6JfUODAwMPVkzMyukkqIiaR+ygnJlRHw/hR9Kp7VI79tSvB+YlJt9IrBlkPhuImJhRPRERE9XV1d5G2JmZs9Qxd1fAi4F7o2Ir+YmLQNqd3DNAa7NxU9Ld4FNBx5Jp8lWADMkjUsX6GekmJmZVaSKaypvAf4bcLekO1LsX4DzgaWS5gIPAqekacuB44E+4HHgAwARsV3SucDa1O6ciNjenk0wM7NG2l5UIuIXNL4eAnBMg/YBzG+yrEXAovKyMzOzkfA36s3MrDQuKmZmVhoXFTMzK42LipmZlcZFxczMSlNpNy2dpspuNMzMOoGPVMzMrDQuKmZmVhoXFTMzK42LipmZlcZFxczMSuOiYmZmpXFRMTOz0riomJlZaVxUzMysNP5GvZlV1lvEpvNPqGS91jodf6Qiaaak+yX1SVpQdT5mZs9lHV1UJI0CvgEcB0wFTpU0tdqszMyeuzq6qABHAn0RsTEingCWALMqzsnM7Dmr06+pTAA258b7gTdVlIuZDZGv5Tz7dHpRUYNY7NZImgfMS6O/l3T/MNc3HvjdMOdtB+c3Ms5vZDomP11QcSaN7e377+VFGnV6UekHJuXGJwJb6htFxEJg4UhXJqk3InpGupxWcX4j4/xGxvmNzN6eX1Gdfk1lLTBF0mRJY4DZwLKKczIze87q6COViNgl6WPACmAUsCgi1leclpnZc1ZHFxWAiFgOLG/T6kZ8Cq3FnN/IOL+RcX4js7fnV4gidruubWZmNiydfk3FzMz2Ii4qDeyp6xdJ+0q6Ok1fLam7jblNkrRK0r2S1kv6ZIM2R0l6RNId6XVmu/JL698k6e607t4G0yXporT/7pI0rY25HZ7bL3dIelTSp+ratHX/SVokaZukdbnYQZJWStqQ3sc1mXdOarNB0pw25vclSfeln98PJI1tMu+gn4UW5ne2pN/kfobHN5m35d08Ncnv6lxumyTd0WTelu+/0kWEX7kX2QX/XwOHAWOAO4GpdW0+ClyShmcDV7cxv0OBaWn4RcCvGuR3FPCjCvfhJmD8INOPB64n+57RdGB1hT/r3wIvr3L/AW8HpgHrcrEvAgvS8ALgggbzHQRsTO/j0vC4NuU3Axidhi9olF+Rz0IL8zsbOL3Az3/Q3/VW5Vc3/SvAmVXtv7JfPlLZXZGuX2YBi9PwNcAxkhp9EbN0EbE1Im5Lw48B95L1LNBJZgGXR+YWYKykQyvI4xjg1xHxHxWs+2kRcTOwvS6c/4wtBk5qMOuxwMqI2B4RO4CVwMx25BcRN0TErjR6C9l3xCrRZP8V0ZZungbLL/3deC9wVdnrrYqLyu4adf1S/0f76TbpF+sR4MVtyS4nnXZ7PbC6weQ3S7pT0vWSjmhrYlmvBjdIujX1ZlCvyD5uh9k0/2Wucv8BHBIRWyH7RwI4uEGbvWU/fpDsyLORPX0WWulj6fTcoianD/eG/fc24KGI2NBkepX7b1hcVHZXpOuXQt3DtJKkFwLfAz4VEY/WTb6N7JTOa4H/DfywnbkBb4mIaWS9R8+X9Pa66XvD/hsDnAh8t8HkqvdfUXvDfvwcsAu4skmTPX0WWuVi4BXA64CtZKeY6lW+/4BTGfwopar9N2wuKrsr0vXL020kjQYOZHiH38MiaR+ygnJlRHy/fnpEPBoRv0/Dy4F9JI1vV34RsSW9bwN+QHaaIa9Q9zotdhxwW0Q8VD+h6v2XPFQ7JZjetzVoU+l+TDcGvAt4f6QLAPUKfBZaIiIeiognI+Ip4FtN1lv1/hsNvAe4ulmbqvbfSLio7K5I1y/LgNqdNicDNzX7pSpbOgd7KXBvRHy1SZuX1K7xSDqS7Of8n23K7wWSXlQbJrugu66u2TLgtHQX2HTgkdqpnjZq+h9ilfsvJ/8ZmwNc26DNCmCGpHHp9M6MFGs5STOBzwInRsTjTdoU+Sy0Kr/8Nbp3N1lv1d08vRO4LyL6G02scv+NSNV3CuyNL7K7k35FdmfI51LsHLJfIID9yE6b9AFrgMPamNtbyQ7R7wLuSK/jgQ8DH05tPgasJ7ub5Rbgb9qY32FpvXemHGr7L5+fyB6u9mvgbqCnzT/f/cmKxIG5WGX7j6y4bQX+TPbf81yya3Q3AhvS+0GpbQ/w7dy8H0yfwz7gA23Mr4/sekTtM1i7G/KlwPLBPgttyu+K9Nm6i6xQHFqfXxrf7Xe9Hfml+GW1z1yubdv3X9kvf6PezMxK49NfZmZWGhcVMzMrjYuKmZmVxkXFzMxK46JiZmalcVGxvZakkPSV3Pjpks4uadmXSTq5jGXtYT2nKOtRelVdvFvS3+fGeyRd1Op8ipD0UknXpOHX5Xv4Tb3/nl5gGT9O3dysl3SJpFGtzNn2Hi4qtjfbCbyngm+zD2qIfyDnAh+NiHfUxbuBp4tKRPRGxCdKSG/EImJLRNQK7uvIvssxVO+NrJubVwNdwCll5Wd7NxcV25vtInvE6qfrJ9QfaUj6fXo/StLPJC2V9CtJ50t6v6Q16bkUr8gt5p2Sfp7avSvNP0rZs0LWps4I/ym33FWSvkP2pbr6fE5Ny18n6YIUO5Psy6qXSPpS3SznA29Lz8n4dFr+j9J8Z0taLOkGZc/TeI+kL6bl/zh104OkN6RtvVXSily3Lp+QdE/Kf0mDXJdLek0avj3liaRzJf1jOopal75lfg7wvpTn+9Iipkr6qaSNkhoWwvhLf3SjybqVjyK52bNA1d++9MuvZi/g98ABZM+UOBA4HTg7TbsMODnfNr0fBTxM9tyZfYHfAP8zTfsk8LXc/D8m+8dqCtk3nfcD5gGfT232BXqByWm5fwAmN8jzpcCDZP+RjwZuAk5K035Kgx4DqHtmS36c7FkgvwD2AV4LPA4cl6b9gKwb/H2A/wd0pfj7gEVpeAuwbxoe22DdC4D5ad+uBVak+CrgcLKjqHUp9t+Br+fmPTutd19gPFnPBPs0+fmtAHYA3wFGFcnNr85/+UjF9mqR/cd7OTCUU0NrI3vuzE6y7jduSPG7yf5g1iyNiKci63Z8I/Aqsv6VTlP2JL7VZN2lTEnt10TEAw3W90bgpxExENmjEK4kezDTSFwfEX9OOY8iK4D5bTic7NTSypTr5/nLM03uAq6U9A9kR3v1fp7yeytwHfBCSfsD3RFxf4HcrouInRHxO7KOLg9p1CgijuUvxf3ogrlZh3NRsU7wNbJrEy/IxXaRPr+p88cxuWk7c8NP5cafIjuSqKnvoyjI+iX7eES8Lr0mR0StKP2hSX6teEDbToDIetn9c0TUcq1tg4D1uTz/S0TMSG1OIOtb7Q3Arcp6w81bS9aH2NuAm4HbgQ8Btw4lt+RJnrlPnyEi/kTW91bt4Vd7ys06nIuK7fUiYjuwlKyw1Gwi+8ME2R+sfYax6FMkPS9dZzkMuJ/slM1HctctXpl6iB3MauBvJY1PF/FPBX62h3keI3sc9HDdD3RJenPKcx9JR0h6HjApIlYBnwHGAi/MzxjZUw43kz1x8BayI5fT0/uI85T0wtz1ndFkF/rvK5KbdT7/l2Cd4itkvQfXfAu4VtIasl58mx1FDOZ+sj/+h5D1FvsnSd8mO710WzoCGqDxo3yfFhFbJZ1Bdk1CZL3MNuqqPu8uYJekO8mu79w+lMQj4ol0o8JFkg4k+13+GlmPu/+eYgIujIiHGyzi58AxEfG4pJ+TnTprVFRWAQvSKbb/VTC9FwDLJO1LduruJuCSNFwkN+tg7qXYzMxK49NfZmZWGhcVMzMrjYuKmZmVxkXFzMxK46JiZmalcVExM7PSuKiYmVlpXFTMzKw0/x+N4KVD3onIrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#0.055917 * 100\n",
    "Number_of_times_I_observed_3 = 11\n",
    "\n",
    "count_3s = []\n",
    "extreme = 0\n",
    "for i in range(50000):\n",
    "    outcome = 100 * prob_of_outcome_tossing_n_dice(3, 2, 100)\n",
    "    count_3s.append( outcome )\n",
    "    if outcome >= 11:\n",
    "        extreme += 1\n",
    "\n",
    "print( extreme / float(50000) )\n",
    "plt.hist(count_3s)\n",
    "plt.xlabel(\"Number of times with 3s\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-06-23T21:59:17.276532Z",
     "start_time": "2019-06-23T21:59:17.265853Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7853981633974483\n",
      "0.8\n"
     ]
    }
   ],
   "source": [
    "# Calculate pi by calculating area of quarter circle\n",
    "\n",
    "def estimate_pi(n_attempts = 10000):\n",
    "    \"\"\" Estimate pi from area of quarter circle \"\"\" \n",
    "    count_successes = 0\n",
    "    for i in range(n_attempts):\n",
    "        x_rand = np.random.random()\n",
    "        y_rand = np.random.random()\n",
    "        dist = math.sqrt( x_rand **2 + y_rand**2 )\n",
    "        if dist <= 1.:\n",
    "            count_successes += 1\n",
    "    return float(count_successes) / n_attempts\n",
    "\n",
    "\n",
    "print( math.pi / 4. )\n",
    "print( estimate_pi(10) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-06-23T21:59:17.819967Z",
     "start_time": "2019-06-23T21:59:17.279677Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGQ9JREFUeJzt3XuwZWV95vHvYyOo8cKtUdIXG017IZYgOSIRTSF4ATGANTJiObFLSTpOwMuYjIJOiWPGKpwxos4kTlBawHFAxAs9SoItosSq4dIIclXpYA8cIXYrNxGFAX7zx37PuDl9unuvPmefvZv+fqpO7bXe/e69fr2qznn6XZd3paqQJGlQjxt1AZKk7YvBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKmToQVHklVJNiS5flr7O5L8KMkNSf5zX/vJSda1917T1354a1uX5KRh1StJGkyGded4kj8C7gPOrqoXtLZXAB8AjqyqB5LsVVUbkuwLnAMcCPwu8C3gOe2rfgy8CpgErgTeVFU3bmnbe+65Zy1btmwI/ypJeuy66qqrfl5VC7fWb6dhFVBVlyZZNq353wKnVtUDrc+G1n40cG5r/0mSdfRCBGBdVd0CkOTc1neLwbFs2TLWrl07J/8OSdpRJPk/g/Sb73MczwFenuTyJN9N8uLWvgi4ra/fZGvbXLskaUSGNuLYwvZ2Aw4CXgycl+RZQGboW8wcbDMeW0uyElgJsHTp0jkpVpK0qfkecUwCX6meK4BHgD1b+5K+fouB27fQvomqOr2qJqpqYuHCrR6ikyRto/kOjq8BhwIkeQ6wM/BzYDVwXJJdkuwDLAeuoHcyfHmSfZLsDBzX+kqSRmRoh6qSnAMcAuyZZBI4BVgFrGqX6D4IrKjeZV03JDmP3knvh4ATqurh9j0nAhcBC4BVVXXDsGqWJG3d0C7HHaWJiYnyqipJ6ibJVVU1sbV+3jkuSerE4JAkdWJwSJI6me/7OKShWnbSNwbqt/7UI4dcifTY5YhDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqROhhYcSVYl2dCeLz79vb9KUkn2bOtJ8qkk65Jcm+SAvr4rktzcflYMq15J0mCGOeI4Ezh8emOSJcCrgFv7mo8AlreflcCnW9/dgVOAlwAHAqck2W2INUuStmJowVFVlwJ3zvDWacB7geprOxo4u3ouA3ZNsjfwGmBNVd1ZVXcBa5ghjCRJ82dez3EkOQr4aVX9YNpbi4Db+tYnW9vm2mf67pVJ1iZZu3HjxjmsWpLUb96CI8mTgA8AH5zp7RnaagvtmzZWnV5VE1U1sXDhwm0vVJK0RfM54ng2sA/wgyTrgcXA95M8g95IYklf38XA7VtolySNyLwFR1VdV1V7VdWyqlpGLxQOqKp/AVYDb2lXVx0E3FNVdwAXAa9Osls7Kf7q1iZJGpFhXo57DvC/gecmmUxy/Ba6XwjcAqwDPgP8BUBV3Qn8NXBl+/lwa5MkjchOw/riqnrTVt5f1rdcwAmb6bcKWDWnxUmStpl3jkuSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOjE4JEmdDPOZ46uSbEhyfV/bf0nywyTXJvlqkl373js5ybokP0rymr72w1vbuiQnDateSdJghjniOBM4fFrbGuAFVfVC4MfAyQBJ9gWOA36/febvkixIsgD4W+AIYF/gTa2vJGlEhhYcVXUpcOe0tm9W1UNt9TJgcVs+Gji3qh6oqp8A64AD28+6qrqlqh4Ezm19JUkjMspzHG8D/qEtLwJu63tvsrVtrn0TSVYmWZtk7caNG4dQriQJRhQcST4APAR8Yapphm61hfZNG6tOr6qJqppYuHDh3BQqSdrETvO9wSQrgNcBh1XVVAhMAkv6ui0Gbm/Lm2uXJI3AvI44khwOvA84qqru73trNXBckl2S7AMsB64ArgSWJ9knyc70TqCvns+aJUmPNrQRR5JzgEOAPZNMAqfQu4pqF2BNEoDLqurtVXVDkvOAG+kdwjqhqh5u33MicBGwAFhVVTcMq2ZJ0tYNLTiq6k0zNJ+xhf4fAT4yQ/uFwIVzWJokaRa8c1yS1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1MnQngAoaWbLTvrGQP3Wn3rkkCuRts3QRhxJViXZkOT6vrbdk6xJcnN73a21J8mnkqxLcm2SA/o+s6L1vznJimHVK0kazDAPVZ0JHD6t7STg4qpaDlzc1gGOAJa3n5XAp6EXNMApwEuAA4FTpsJGkjQaQwuOqroUuHNa89HAWW35LOCYvvazq+cyYNckewOvAdZU1Z1VdRewhk3DSJI0j+b75PjTq+oOgPa6V2tfBNzW12+ytW2ufRNJViZZm2Ttxo0b57xwSVLPuFxVlRnaagvtmzZWnV5VE1U1sXDhwjktTpL0WwMFR5IXzNH2ftYOQdFeN7T2SWBJX7/FwO1baJckjcigI47/nuSKJH+RZNdZbG81MHVl1Arggr72t7Srqw4C7mmHsi4CXp1kt3ZS/NWtTZI0IgPdx1FVL0uyHHgbsDbJFcDnqmrN5j6T5BzgEGDPJJP0ro46FTgvyfHArcCxrfuFwGuBdcD9wFvbdu9M8tfAla3fh6tq+gl3SdI8GvgGwKq6Ocl/ANYCnwJelCTA+6vqKzP0f9NmvuqwGfoWcMJmtrsKWDVonZKk4Rr0HMcLk5wG3AQcCvxxVT2/LZ82xPokSWNm0BHHfwM+Q2908eupxqq6vY1CJEk7iEGD47XAr6vqYYAkjwOeUFX3V9Xnh1adJGnsDHpV1beAJ/atP6m1SZJ2MIMGxxOq6r6plbb8pOGUJEkaZ4MGx6+mzVj7B8Cvt9BfkvQYNeg5jncDX0oyddf23sAbh1OSJGmcDXoD4JVJngc8l978UT+sqv871MokSWOpyxMAXwwsa595URKq6uyhVCVJGlsDBUeSzwPPBq4BHm7NBRgckrSDGXTEMQHs26YGkSTtwAa9qup64BnDLESStH0YdMSxJ3BjmxX3ganGqjpqKFVJksbWoMHxoWEWIUnafgx6Oe53kzwTWF5V30ryJGDBcEuTJI2jQadV/zPgfODvW9Mi4GvDKkqSNL4GPTl+AnAwcC/0HuoE7DWsoiRJ42vQcxwPVNWDvQf+QZKd6N3HsU2S/DvgT9t3XEfvUbF7A+cCuwPfB/6kbXMXeveL/AHwC+CNVbV+W7et2Vl20jcG6rf+1COHXImkURl0xPHdJO8HnpjkVcCXgP+1LRtMsgh4JzBRVS+gd67kOOCjwGlVtRy4Czi+feR44K6q+j16Txv86LZsV5I0NwYNjpOAjfRGB38OXAjM5sl/O9ELoZ3oTc9+B73H0J7f3j8LOKYtH93Wae8flqmhjyRp3g16VdUj9B4d+5nZbrCqfprkY8Ct9KZm/yZwFXB3VT3Uuk3SOwFPe72tffahJPcAewA/n20tkqTuBp2r6ifMcE6jqp7VdYNJdqM3itgHuJveYa8jZug6tb2ZRheb1JJkJbASYOnSpV3LkiQNqMtcVVOeABxL7yT2tngl8JOq2giQ5CvAS4Fdk+zURh2Lgalnf0wCS4DJdmjracCd07+0qk4HTgeYmJhwTi1JGpKBznFU1S/6fn5aVZ+gd05iW9wKHJTkSe1cxWHAjcAlwBtanxXABW15dVunvf9tJ1uUpNEZ9FDVAX2rj6M3AnnKtmywqi5Pcj69S24fAq6mN1L4BnBukv/U2s5oHzkD+HySdfRGGsdty3YlSXNj0ENVf9O3/BCwHvjX27rRqjoFOGVa8y3AgTP0/Q29Q2OSpDEw6FVVrxh2IZKk7cOgh6res6X3q+rjc1OOJGncdbmq6sX0TlQD/DFwKe3+CknSjqPLg5wOqKpfAiT5EPClqvrTYRUmSRpPg045shR4sG/9QWDZnFcjSRp7g444Pg9ckeSr9O7afj29GWslSTuYQa+q+kiSfwBe3preWlVXD68s6dEGnc59GN/nFPHSow16qAp6s9jeW1WfpDf9xz5DqkmSNMYGfXTsKcD7gJNb0+OB/zGsoiRJ42vQEcfrgaOAXwFU1e1s45QjkqTt26DB8WCbWLAAkvzO8EqSJI2zQa+qOi/J39Ob+vzPgLcxBw91kub6pLek4Rv0qqqPtWeN3ws8F/hgVa0ZamWSpLG01eBIsgC4qKpeCRgW2uE4KpIebavnOKrqYeD+JE+bh3okSWNu0HMcvwGuS7KGdmUVQFW9cyhVSZLG1qDB8Y32I0nawW0xOJIsrapbq+qs+SpIkjTetnaO42tTC0m+PFcbTbJrkvOT/DDJTUn+MMnuSdYkubm97tb6JsmnkqxLcu20559LkubZ1oIjfcvPmsPtfhL4x6p6HrAfcBNwEnBxVS0HLm7rAEcAy9vPSuDTc1iHJKmjrQVHbWZ5myV5KvBHwBkAVfVgVd0NHA1MHRI7CzimLR8NnF09l9G7CXHvuahFktTd1k6O75fkXnojjye2Zdp6VdVTt2GbzwI2Ap9Lsh9wFfAu4OlVdQe9L74jyV6t/yIe/YjaydZ2xzZsW5I0S1sMjqpaMKRtHgC8o6ouT/JJfntYaiaZoW2T0U+SlfQOZbF06dK5qFOz4E1z0mNXl+dxzJVJYLKqLm/r59MLkp9NHYJqrxv6+i/p+/xi4PbpX1pVp1fVRFVNLFy4cGjFS9KObt6Do6r+BbgtyXNb02HAjcBqYEVrWwFc0JZXA29pV1cdBNwzdUhLkjT/Br0BcK69A/hCkp2BW4C30gux85IcD9wKHNv6Xgi8FlgH3N/6SpJGZCTBUVXXABMzvHXYDH0LOGHoRUmSBjKKcxySpO2YwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUiejugFQ0jwbdP6w9aceOeRKtL1zxCFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1MrLgSLIgydVJvt7W90lyeZKbk3yxPY+cJLu09XXt/WWjqlmSNNoRx7uAm/rWPwqcVlXLgbuA41v78cBdVfV7wGmtnyRpREYSHEkWA0cCn23rAQ4Fzm9dzgKOactHt3Xa+4e1/pKkERjViOMTwHuBR9r6HsDdVfVQW58EFrXlRcBtAO39e1p/SdIIzPu06kleB2yoqquSHDLVPEPXGuC9/u9dCawEWLp06RxUKo2W06BrXI1ixHEwcFSS9cC59A5RfQLYNclUkC0Gbm/Lk8ASgPb+04A7p39pVZ1eVRNVNbFw4cLh/gskaQc278FRVSdX1eKqWgYcB3y7qt4MXAK8oXVbAVzQlle3ddr7366qTUYckqT5MU73cbwPeE+SdfTOYZzR2s8A9mjt7wFOGlF9kiRG/OjYqvoO8J22fAtw4Ax9fgMcO6+FSZI2a5xGHJKk7YDBIUnqxOCQJHUy0nMckmZv0Ps9pLniiEOS1InBIUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTbwAU4E1kkgbniEOS1InBIUnqxOCQJHVicEiSOjE4JEmdzHtwJFmS5JIkNyW5Icm7WvvuSdYkubm97tbak+RTSdYluTbJAfNdsyTpt0Yx4ngI+Muqej5wEHBCkn2Bk4CLq2o5cHFbBzgCWN5+VgKfnv+SJUlT5j04quqOqvp+W/4lcBOwCDgaOKt1Ows4pi0fDZxdPZcBuybZe57LliQ1Iz3HkWQZ8CLgcuDpVXUH9MIF2Kt1WwTc1vexydYmSRqBkQVHkicDXwbeXVX3bqnrDG01w/etTLI2ydqNGzfOVZmSpGlGEhxJHk8vNL5QVV9pzT+bOgTVXje09klgSd/HFwO3T//Oqjq9qiaqamLhwoXDK16SdnDzPldVkgBnADdV1cf73loNrABOba8X9LWfmORc4CXAPVOHtLRlzj8laRhGMcnhwcCfANcluaa1vZ9eYJyX5HjgVuDY9t6FwGuBdcD9wFvnt1xJUr95D46q+h4zn7cAOGyG/gWcMNSiJP1/XUaq6089coiVaFx557gkqRODQ5LUicEhSerE4JAkdeKjY7dDXmYraZQccUiSOnHEMUYcSUjaHhgcksbGoP958v6R0TI4JG0z/9DvmDzHIUnqxOCQJHXioSpJQ+eFH48tjjgkSZ0YHJKkTgwOSVInBockqRODQ5LUicEhSepkuwmOJIcn+VGSdUlOGnU9krSj2i7u40iyAPhb4FXAJHBlktVVdeNoKxuM17BLc2uuf6ecEqWb7SI4gAOBdVV1C0CSc4GjgZEGh4EgPTYM43f5sRxG20twLAJu61ufBF4yrI0ZCJJm67E8AeT2EhyZoa0e1SFZCaxsq/cl+dEstrcn8PNZfH7YrG92rG92rG92HlVfPjrCSjb1zEE6bS/BMQks6VtfDNze36GqTgdOn4uNJVlbVRNz8V3DYH2zY32zY32zM+71DWJ7uarqSmB5kn2S7AwcB6wecU2StEPaLkYcVfVQkhOBi4AFwKqqumHEZUnSDmm7CA6AqroQuHCeNjcnh7yGyPpmx/pmx/pmZ9zr26pU1dZ7SZLUbC/nOCRJY8Lg6JNkfZLrklyTZO2o6wFIsirJhiTX97XtnmRNkpvb625jVt+Hkvy07cdrkrx2RLUtSXJJkpuS3JDkXa19LPbfFuobl/33hCRXJPlBq+8/tvZ9klze9t8X2wUr41TfmUl+0rf/9h9FfX11LkhydZKvt/Wx2H+zYXBs6hVVtf8YXS53JnD4tLaTgIurajlwcVsflTPZtD6A09p+3L+dnxqFh4C/rKrnAwcBJyTZl/HZf5urD8Zj/z0AHFpV+wH7A4cnOQj4aKtvOXAXcPyY1Qfw7/v23zUjqm/Ku4Cb+tbHZf9tM4NjzFXVpcCd05qPBs5qy2cBx8xrUX02U99YqKo7qur7bfmX9H55FzEm+28L9Y2F6rmvrT6+/RRwKHB+ax/l/ttcfWMjyWLgSOCzbT2Myf6bDYPj0Qr4ZpKr2p3o4+rpVXUH9P74AHuNuJ6ZnJjk2nYoa2SH0qYkWQa8CLicMdx/0+qDMdl/7TDLNcAGYA3wz8DdVfVQ6zLJCMNuen1VNbX/PtL232lJdhlVfcAngPcCj7T1PRij/betDI5HO7iqDgCOoHfY4I9GXdB26tPAs+kdPrgD+JtRFpPkycCXgXdX1b2jrGUmM9Q3Nvuvqh6uqv3pzdZwIPD8mbrNb1V9G55WX5IXACcDzwNeDOwOvG8UtSV5HbChqq7qb56h61iNkgZhcPSpqtvb6wbgq/R+UcbRz5LsDdBeN4y4nkepqp+1X+hHgM8wwv2Y5PH0/ih/oaq+0prHZv/NVN847b8pVXU38B1652J2TTJ1D9gm0/+MQl99h7dDgFVVDwCfY3T772DgqCTrgXPpHaL6BGO4/7oyOJokv5PkKVPLwKuB67f8qZFZDaxoyyuAC0ZYyyam/ig3r2dE+7EdTz4DuKmqPt731ljsv83VN0b7b2GSXdvyE4FX0jsPcwnwhtZtlPtvpvp+2PefgtA7fzCS/VdVJ1fV4qpaRm+apG9X1ZsZk/03G94A2CR5Fr1RBvTuqP+fVfWREZYEQJJzgEPozaj5M+AU4GvAecBS4Fbg2KoayQnqzdR3CL3DLAWsB/586pzCPNf2MuCfgOv47THm99M7jzDy/beF+t7EeOy/F9I7ebuA3n8yz6uqD7fflXPpHQa6Gvg37X/341Lft4GF9A4LXQO8ve8k+kgkOQT4q6p63bjsv9kwOCRJnXioSpLUicEhSerE4JAkdWJwSJI6MTgkSZ1sNw9ykoYhyR70JjoEeAbwMLCxrd9fVS+d4+3tD/xu14kL22XPv0/vhrbdgEur6ltzWZs0KINDO7Sq+gW9eyZI8iHgvqr62BA3uT8wQYenWSZ5BvDSqnrmVvotqKqHZ1mftFUeqpI2I8l97fWQJN9Ncl6SHyc5Ncmb27Mgrkvy7NZvYZIvJ7my/Rw87ft2Bj4MvLE9J+KN6T0b5GttQr7L2k1t030T2Kt95uXteRNvaN+5PskHk3wPODbJs5P8Y5uo85+SPG+oO0k7JEcc0mD2ozfB353ALcBnq+rA9B6+9A7g3cAn6T1n4XtJlgIX0TcpYFU9mOSDwERVnQiQ5L8CV1fVMUkOBc6mjYD6HAV8vU3mR5Lpz2/4TVW9rL13Mb07pW9O8hLg7+jNkSTNGYNDGsyVU9N+JPlneqMA6E0X8oq2/Epg394USQA8NclT2rM2NudlwL8CqKpvJ9kjydOq6p4OtX2x1fVk4KXAl/pqGOWU4nqMMjikwfTPJfRI3/oj/Pb36HHAH1bVrzt871xMs/2rvu3fPTUykYbFcxzS3PkmcOLUSmZ+1vUvgaf0rV8KvLn1PwT4+bY+M6R97idJjm3flyT7bct3SVticEhz553ARDvRfSPw9hn6XELvcNY1Sd4IfGjqM8Cp/Ha69231ZuD4JD8AbqD3mFxpTjk7riSpE0cckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnfw/VM/u8dCeyYkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate distribution of firing times of a model neuron\n",
    "\n",
    "def time_to_fire( threshold , lamb, scale_size):\n",
    "    \"\"\" Simulate distribution of times for neuron that sums\n",
    "    inputs and fires when sum exceeds threshold \"\"\"\n",
    "    state = 0.\n",
    "    time = 0.\n",
    "    while state < threshold: \n",
    "        # Update time\n",
    "        delay_between_inputs = np.random.poisson(lamb)\n",
    "        time += delay_between_inputs\n",
    "        \n",
    "        # Update state\n",
    "        size_of_input = scale_size * np.random.random()\n",
    "        state += size_of_input\n",
    "        \n",
    "    return time    \n",
    "\n",
    "# Initialize variables\n",
    "N_neurons = 10000\n",
    "threshold = 1.\n",
    "scale_size = 0.1\n",
    "lamb = 1.\n",
    "\n",
    "list_times = []\n",
    "for neuron in range(N_neurons):\n",
    "    list_times.append( time_to_fire( threshold , lamb, scale_size) )\n",
    "\n",
    "    \n",
    "#print list_times\n",
    "\n",
    "plt.hist(list_times, 30)\n",
    "plt.xlabel(\"Time to fire\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}