Python_and_distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Python to visualize distributions\n",
    "\n",
    "Using Python and some of its libraries can be a great way to visualize, compare and understand probability distributions.  Combining statistics calculation functionality with plotting tools allows us to easily see the shape of probability density functions and to answer thought experiments.  This blog posts looks at some of the common ways this is done.\n",
    "\n",
    "Let's start by importing the numpy, pandas and matplotlib libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete distributions\n",
    "\n",
    "The first example looks a rolling a dice.  Here we can simulate the rolling of a dice by using the numpy random integer generator, randint. In order to be able to replicate the results in the future, we'll pick a seed for the random number generator.\n",
    "\n",
    "If we run this for 10 samples, we get the following output: array([4, 2, 6, 4, 2, 6, 1, 5, 1, 6]) where these numbers represent the results of a six-sided dice.  Now we can run it again with a sample size of 1000 and plot the results in a histogram. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4, 2, 6, 4, 2, 6, 1, 5, 1, 6])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Pick a seed for the random number generator\n",
    "np.random.seed(101)\n",
    "\n",
    "# Generate random integers between 1 and 6 for a sample size of 10\n",
    "np.random.randint(1, 7, 10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEICAYAAAC55kg0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGOBJREFUeJzt3X+cXHV97/HX2wQQIeFHs6UhiW6EQAs82gBrxOuP0qISIhVsLSZt+WHRyBWs1HoteG8Ltebe/gBteWhjA6QQlWAEEWpjJVJbLq0ImxghCSBLSB7ZbUwWIoYIRQKf/nG+qyfLzuzszOyeXb7v5+Mxjz3nO99zzmfOJPOe8z1nZhQRmJlZnl5RdQFmZlYdh4CZWcYcAmZmGXMImJllzCFgZpYxh4CZWcYcAi9DkjZKOrXqOqok6V2StknaI+nEESz3OUl/Mko1bZH01tFY9xDbukHSJxvs2ykpJE1O81+XdP7oVjhsTRdIuqc0H5KOrrKml6vJVRdgIyNpC/C+iPhmqe2C1PYmgIg4voH1dAKPA/tFxN7RqLViVwGXRMTtI1koIi4apXpGRFIAcyKiZ6y3HRFnjPU2rTo+ErBRMfCuskKvATZWXIPVIGlS1TVYwSHwMlQedpA0T1K3pN2Sdkj6VOp2d/r7VBoyeYOkV0j6P5K2StopaYWkQ0rrPS/d96SkPxm0nSsl3SLpC5J2AxekbX9b0lOStkv6jKT9S+sLSR+U9KikpyX9uaSjJP1HqndVuf+gxzhkrZIOkLQHmAR8T9JjQywrSZ9Oy+2W9KCkE9J9Px1GkXSqpF5JH0t9t0s6W9ICSd+XtEvSx0vr3WcIZmD5GvXX3DeSBp6b76Xn5j2p/UxJ69My/yHpl0vrO1HSurQfvwS8cqjtpr6TJF0l6QlJm4F3DLr/XyW9rzT/fkkPpXVvknRSaj9S0q2S+iU9LukP6mzzBklLJa2W9GPg19LztSItvzU9n8O+JqX9vynV0yfpo8MtY3VEhG8T6AZsAd46qO0C4J6h+gDfBs5N0wcDp6TpTiCAyaXlfh/oAV6b+n4F+Hy67zhgD/AmYH+K4ZbnS9u5Ms2fTfHm4kDgZOAUimHHTuAh4NLS9gK4HZgKHA88B9yVtn8IsAk4v8Z+qFlrad1H11j2dGAtcCgg4JeA6em+G4BPpulTgb3AnwL7Ae8H+oGbgCmp5meB2YOXLS3fW+N5aWTfHF2aPxHYCbyeIuDOT+s7ID0fW4E/THW+Oz0Xn6zx+C8CHgZmAYcD3yr/WwD+lWJ4EeC3gT7gdWlfHU1xlPWKtA//NG3/tcBm4PQa27wB+BHwxrTsK4EV6fmfkvbB94ELa/yb/un+ALYDb07ThwEnVf3/ciLffCQwMX01vRt8StJTwN/V6fs8cLSkaRGxJyLurdP3d4FPRcTmiNgDXA4sVDG0827gHyPinoj4CcV//sFfPPXtiPhqRLwYEc9GxNqIuDci9kbEFuDvgV8dtMxfRcTuiNgIbADuTNv/EfB1ihe/kdY6nOcpXnh+EVBEPBQR2+v0XRIRzwM3A9OAv42Ip1PNm4BfaWCb+2hw35QtBv4+Ir4TES9ExI0UoXlKuu0H/E1EPB8RtwD311nXOanvtojYBfy/On3fR/Ec3R+FnojYShEKHRHxiYj4SURsBq4FFtZZ1+0R8e8R8SLFfl0IXJ725RbgauDcOssPeB44TtLUiPhhRKxrYBmrwSEwMZ0dEYcO3IAP1ul7IXAM8LCk+yWdWafvkRTvKAdspXinekS6b9vAHRHxDPDkoOW3lWckHSPpa5J+kIaI/i/Fi2jZjtL0s0PMH9xErXVFxL8AnwE+C+yUtEzS1Brdn4yIF0r1DFVzrRpranDflL0G+KNB4T+LYj8cCfRFRDmUtw61kmSf53KYvrOAlwyppXqOHFTPx6m//8vbnEYRXIOfwxl1lh/wW8ACYKukf5P0hgaWsRocAi9zEfFoRCwCfh74S+AWSQfx0nfxAP9J8Z97wKsphkN2UByCzxy4Q9KBwM8N3tyg+aUUww5zImIqxYuEmn80Ddc6rIi4JiJOphjmOgb4X22o6cfAq0rzv1Cn70j3zTaKI5JDS7dXRcRKiudmhqTy8q+us67tFC/ujfTdBhxVo/3xQfVMiYgFddZV/vfxBMU7+sHPYV+d5YuVFEclZ1H8m/4qsGq4Zaw2h8DLnKTfk9SRDsGfSs0vUoxtv0gxljtgJfCHkmZLOpji3emXoriE9BbgNyT9j3QC80qGf0GfAuwG9kj6ReB/tutxDVNrXZJeJ+n1kvajeOH+L4p90ar1wAJJh0v6BeDSOn2H2zc72Pe5uRa4KNUtSQdJeoekKRTnffYCfyBpP0m/Ccyrs+1Vqe9MSYcBl9Xpex3wUUknp+0eLek1wH3A05L+WNKB6WTzCZJeV2ddP5WOrlYBSyRNSev8CPCFestJ2l/S70o6JA3R7aY9z122HAIvf/OBjSqumPlbYGEar38GWAL8ezqcPwVYDnye4sqhxyleHD8EkMa/P0QxLr6d4iTxTopx6Vo+CvwO8DTFi9iX2vi4atbagKmpnh9SDEE8Cfx1G2r6PPA9ihO2d1L/8Q63b64EbkzPzTkR0U1xYvozqe4eipOnpHM0v5nmdwHvoThRXsu1wDdSrevq9Y2IL1P8O7kp1fpV4PD0In4mMJdi/z9BERiH1FjVUD5EEcKbgXvSNpY3sNy5wJY0jHYRxfkha5L2HUY0a0x69/0UxXDG41XXY2bN8ZGANUzSb0h6VTqncBXwIMW7XjOboBwCNhJnUZyQ/U9gDsXQkg8lzSYwDweZmWXMRwJmZhmr+ku+hjVt2rTo7Oysugwzswlj7dq1T0RERyN9x30IdHZ20t3dXXUZZmYThqR6nwLfh4eDzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGON/KjzchU/sr2h1PYlFT94vV7Fj42vT+2dkp4t3fe50jInq/hB7x5J1wz6AQwzM6tAIx8Wu4HiO8xXDDRExHsGpiVdTfED0gMei4i5Q6xnKcX3oX8HWE3xPfdfH3nJZmbWLsOGQETcLalzqPvSu/lzgF+vtw5J04GpAz9yLmkFcDYOATNrUudl/zSq69/yF+8Y1fWPF62eE3gzsCMiHi21zZb03fQD0G9ObTOA3lKfXur8oLSkxZK6JXX39/e3WKKZmdXSaggsovit1wHbgVdHxIkUvxd6k6SpI11pRCyLiK6I6OroaOg7kMzMrAlNf4GcpMkUv2t68kBbRDxH+s3ZiFgr6THgGKAPmFlafGZqMzOzCrVyJPBW4OGI+Okwj6QOSZPS9Gspfn1qc0RsB3ZLOiWdRzgPuL2FbZuZWRs0conoSuDbwLGSeiVdmO5ayL5DQQBvAR5Il4zeAlwUEbvSfR8ErgN6gMfwSWEzs8o1cnXQohrtFwzRditwa43+3cAJI6zPzMxGkT8xbGaWMYeAmVnGxv3PS5o1yx8mMhueQ8BsnHKI2VhwCJjZqBjtELP28DkBM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OM+eqgUTTRL/Gb6PWb2fB8JGBmljGHgJlZxhwCZmYZcwiYmWXMIWBmljGHgJlZxhwCZmYZcwiYmWXMIWBmlrFhQ0DSckk7JW0otV0pqU/S+nRbULrvckk9kh6RdHqpfX5q65F0WfsfipmZjVQjRwI3APOHaP90RMxNt9UAko4DFgLHp2X+TtIkSZOAzwJnAMcBi1JfMzOr0LDfHRQRd0vqbHB9ZwE3R8RzwOOSeoB56b6eiNgMIOnm1HfTiCs2M7O2aeWcwCWSHkjDRYelthnAtlKf3tRWq31IkhZL6pbU3d/f30KJZmZWT7MhsBQ4CpgLbAeubltFQEQsi4iuiOjq6Oho56rNzKykqa+SjogdA9OSrgW+lmb7gFmlrjNTG3XazcysIk0dCUiaXpp9FzBw5dAdwEJJB0iaDcwB7gPuB+ZImi1pf4qTx3c0X7aZmbXDsEcCklYCpwLTJPUCVwCnSpoLBLAF+ABARGyUtIrihO9e4OKIeCGt5xLgG8AkYHlEbGz7ozEzsxFp5OqgRUM0X1+n/xJgyRDtq4HVI6rOzMxGlT8xbGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWsWFDQNJySTslbSi1/bWkhyU9IOk2SYem9k5Jz0pan26fKy1zsqQHJfVIukaSRuchmZlZoxo5ErgBmD+obQ1wQkT8MvB94PLSfY9FxNx0u6jUvhR4PzAn3Qav08zMxtiwIRARdwO7BrXdGRF70+y9wMx665A0HZgaEfdGRAArgLObK9nMzNqlHecEfh/4eml+tqTvSvo3SW9ObTOA3lKf3tQ2JEmLJXVL6u7v729DiWZmNpSWQkDS/wb2Al9MTduBV0fEicBHgJskTR3peiNiWUR0RURXR0dHKyWamVkdk5tdUNIFwJnAaWmIh4h4DnguTa+V9BhwDNDHvkNGM1ObmZlVqKkjAUnzgY8B74yIZ0rtHZImpenXUpwA3hwR24Hdkk5JVwWdB9zecvVmZtaSYY8EJK0ETgWmSeoFrqC4GugAYE260vPedCXQW4BPSHoeeBG4KCIGTip/kOJKowMpziGUzyOYmVkFhg2BiFg0RPP1NfreCtxa475u4IQRVWdmZqPKnxg2M8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy1lAISFouaaekDaW2wyWtkfRo+ntYapekayT1SHpA0kmlZc5P/R+VdH77H46ZmY1Eo0cCNwDzB7VdBtwVEXOAu9I8wBnAnHRbDCyFIjSAK4DXA/OAKwaCw8zMqtFQCETE3cCuQc1nATem6RuBs0vtK6JwL3CopOnA6cCaiNgVET8E1vDSYDEzszHUyjmBIyJie5r+AXBEmp4BbCv1601ttdpfQtJiSd2Suvv7+1so0czM6mnLieGICCDasa60vmUR0RURXR0dHe1arZmZDdJKCOxIwzykvztTex8wq9RvZmqr1W5mZhVpJQTuAAau8DkfuL3Ufl66SugU4Edp2OgbwNslHZZOCL89tZmZWUUmN9JJ0krgVGCapF6Kq3z+Algl6UJgK3BO6r4aWAD0AM8A7wWIiF2S/hy4P/X7REQMPtlsZmZjqKEQiIhFNe46bYi+AVxcYz3LgeUNV2dmZqPKnxg2M8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy1nQISDpW0vrSbbekSyVdKamv1L6gtMzlknokPSLp9PY8BDMza9bkZheMiEeAuQCSJgF9wG3Ae4FPR8RV5f6SjgMWAscDRwLflHRMRLzQbA1mZtaadg0HnQY8FhFb6/Q5C7g5Ip6LiMeBHmBem7ZvZmZNaFcILARWluYvkfSApOWSDkttM4BtpT69qc3MzCrScghI2h94J/Dl1LQUOIpiqGg7cHUT61wsqVtSd39/f6slmplZDe04EjgDWBcROwAiYkdEvBARLwLX8rMhnz5gVmm5mantJSJiWUR0RURXR0dHG0o0M7OhtCMEFlEaCpI0vXTfu4ANafoOYKGkAyTNBuYA97Vh+2Zm1qSmrw4CkHQQ8DbgA6Xmv5I0Fwhgy8B9EbFR0ipgE7AXuNhXBpmZVaulEIiIHwM/N6jt3Dr9lwBLWtmmmZm1jz8xbGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWMYeAmVnGHAJmZhlzCJiZZcwhYGaWsZZDQNIWSQ9KWi+pO7UdLmmNpEfT38NSuyRdI6lH0gOSTmp1+2Zm1rx2HQn8WkTMjYiuNH8ZcFdEzAHuSvMAZwBz0m0xsLRN2zczsyaM1nDQWcCNafpG4OxS+4oo3AscKmn6KNVgZmbDaEcIBHCnpLWSFqe2IyJie5r+AXBEmp4BbCst25va9iFpsaRuSd39/f1tKNHMzIYyuQ3reFNE9En6eWCNpIfLd0ZESIqRrDAilgHLALq6uka0rJmZNa7lI4GI6Et/dwK3AfOAHQPDPOnvztS9D5hVWnxmajMzswq0FAKSDpI0ZWAaeDuwAbgDOD91Ox+4PU3fAZyXrhI6BfhRadjIzMzGWKvDQUcAt0kaWNdNEfHPku4HVkm6ENgKnJP6rwYWAD3AM8B7W9y+mZm1oKUQiIjNwK8M0f4kcNoQ7QFc3Mo2zcysffyJYTOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGMOATOzjDkEzMwy5hAwM8uYQ8DMLGNNh4CkWZK+JWmTpI2SPpzar5TUJ2l9ui0oLXO5pB5Jj0g6vR0PwMzMmje5hWX3An8UEeskTQHWSlqT7vt0RFxV7izpOGAhcDxwJPBNScdExAst1GBmZi1o+kggIrZHxLo0/TTwEDCjziJnATdHxHMR8TjQA8xrdvtmZta6tpwTkNQJnAh8JzVdIukBScslHZbaZgDbSov1UiM0JC2W1C2pu7+/vx0lmpnZEFoOAUkHA7cCl0bEbmApcBQwF9gOXD3SdUbEsojoioiujo6OVks0M7MaWgoBSftRBMAXI+IrABGxIyJeiIgXgWv52ZBPHzCrtPjM1GZmZhVp5eogAdcDD0XEp0rt00vd3gVsSNN3AAslHSBpNjAHuK/Z7ZuZWetauTrojcC5wIOS1qe2jwOLJM0FAtgCfAAgIjZKWgVsoriy6GJfGWRmVq2mQyAi7gE0xF2r6yyzBFjS7DbNzKy9/IlhM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OMOQTMzDLmEDAzy5hDwMwsYw4BM7OMjXkISJov6RFJPZIuG+vtm5nZz4xpCEiaBHwWOAM4Dlgk6bixrMHMzH5mrI8E5gE9EbE5In4C3AycNcY1mJlZoogYu41J7wbmR8T70vy5wOsj4pJB/RYDi9PsCcCGMSuyvaYBT1RdRAtcf7Vcf7Umcv3HRsSURjpOHu1KmhERy4BlAJK6I6Kr4pKaMpFrB9dfNddfrYlcv6TuRvuO9XBQHzCrND8ztZmZWQXGOgTuB+ZImi1pf2AhcMcY12BmZsmYDgdFxF5JlwDfACYByyNi4zCLLRv9ykbNRK4dXH/VXH+1JnL9Ddc+pieGzcxsfPEnhs3MMuYQMDPL2LgNgYn89RKSlkvaKWlCfr5B0ixJ35K0SdJGSR+uuqaRkPRKSfdJ+l6q/8+qrmmkJE2S9F1JX6u6lpGStEXSg5LWj+RSxfFC0qGSbpH0sKSHJL2h6poaJenYtN8HbrslXVp3mfF4TiB9vcT3gbcBvRRXFS2KiE2VFtYgSW8B9gArIuKEqusZKUnTgekRsU7SFGAtcPYE2v8CDoqIPZL2A+4BPhwR91ZcWsMkfQToAqZGxJlV1zMSkrYAXRExIT9oJelG4P9HxHXpKsZXRcRTVdc1Uul1tI/iA7lba/Ubr0cCE/rrJSLibmBX1XU0KyK2R8S6NP008BAwo9qqGheFPWl2v3Qbf+92apA0E3gHcF3VteRG0iHAW4DrASLiJxMxAJLTgMfqBQCM3xCYAWwrzfcygV6EXk4kdQInAt+ptpKRScMp64GdwJqImEj1/w3wMeDFqgtpUgB3SlqbvgJmIpkN9AP/kIbjrpN0UNVFNWkhsHK4TuM1BGwckHQwcCtwaUTsrrqekYiIFyJiLsWn0udJmhDDcpLOBHZGxNqqa2nBmyLiJIpvC744DY9OFJOBk4ClEXEi8GNgQp2TBEjDWO8Evjxc3/EaAv56iYqlsfRbgS9GxFeqrqdZ6VD+W8D8qmtp0BuBd6Zx9ZuBX5f0hWpLGpmI6Et/dwK3UQzvThS9QG/pyPEWilCYaM4A1kXEjuE6jtcQ8NdLVCidWL0eeCgiPlV1PSMlqUPSoWn6QIoLDB6utqrGRMTlETEzIjop/t3/S0T8XsVlNUzSQeliAtIwytuZQN8CHBE/ALZJOjY1nQZMiAsiBllEA0NBMH6/RbSZr5cYNyStBE4FpknqBa6IiOurrWpE3gicCzyYxtUBPh4RqyusaSSmAzemqyNeAayKiAl3qeUEdQRwW/E+gsnATRHxz9WWNGIfAr6Y3oBuBt5bcT0jksL3bcAHGuo/Hi8RNTOzsTFeh4PMzGwMOATMzDLmEDAzy5hDwMwsYw4BM7OMOQTMzDLmEDAzy9h/A/x5iSWEB6SSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate random integers between 1 and 6 for a sample size of 10000\n",
    "x = np.random.randint(1, 7, 10000)\n",
    "\n",
    "\n",
    "plt.hist(x)\n",
    "plt.title(\"Histogram of simulated dice rolls\")\n",
    "plt.xlim([0,7]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the simulated dice rolls to help gain insight into probability questions, particularly complex ones.  For example, we might be faced with a problem asking to find the probability of rolling two sixes in 10 rolls?  This can be straightforwardly simulated with Python and numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.302\n"
     ]
    }
   ],
   "source": [
    "# Simulate\n",
    "n_sims = 1000\n",
    "count = 0\n",
    "for i in range(n_sims):\n",
    "    # simulate ten rolls\n",
    "    rolls = np.random.randint(1,7,10)\n",
    "    \n",
    "    # Look for two sixes in the 10 rolls\n",
    "    if len(np.where(rolls == 6)[0]) == 2:\n",
    "        count += 1\n",
    "        \n",
    "print(count/n_sims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can simulate trials from a geometric distribution.  Here the geometric distribution provides the probability of success on the nth trial for a Bernoulli process.  In this example, we can determine the probability it takes 5 tries to be successful when the probability of success with each trial is 20%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEICAYAAAC9E5gJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGl9JREFUeJzt3X20XXV95/H3RyKCghAkTTGgQY1abCvQCGhtS8XyqAOdUQaKGi2u2IqOzmgVXU5BlBl0WR/ojHSooEGtSPEBfKiQ4nMVJAiigEqEsJIIJBICRBQFvvPH/l05udzn3HvPvdf3a6277t6//dt7//bZ55zP2b+9zz6pKiRJv90e0e8GSJL6zzCQJBkGkiTDQJKEYSBJwjCQJGEYTKsk1yU5uN/t6Kckf5lkbZItSfabonX8SZIfTdGyT03ysalY9hDrOjjJuulYV79N5WsjyQlJLh1j3WnbvzONYTBJkqxJ8vxBZS9P8s2B8ap6RlV9dZTlLE5SSeZNUVP77T3Aa6pqp6q6eipWUFXfqKqnTcWyxyPJR5K8s9/t6KexPgZjeW0MsewxvVaq6uNVdeh4lv3byDD4LTMDQuaJwHV9boNmiKl+Ps6A5/usYRhMo96jhyQHJFmV5O4ktyd5b6v29fZ/c+tKeXaSRyR5W5JbkmxIcl6SXXqW+7I27Y4k/3PQek5NcmGSjyW5G3h5W/e3k2xOcmuS/5Nk+57lVZJXJ7kxyT1J3pHkyUm+1dp7QW/9Qds4ZFuTPCrJFmA74HtJfjLM/B9o3Uh3J7kqyZ+M8HgemeT61sb1Sd7YyrfqXmmPx98luTbJz5Ock2Rhkn9r8/57kvlDzTt4vw3Rhn9NcluSu5J8PckzWvly4ATgTW0/fq6VPz7Jp5JsTHJzkv/Ws6wd2yfpO5NcDzxruG1v9Q9N8qO27g8m+VqSV/ZM/+skN7TlXZLkiT3TnpPkyjbvlUme0zPtq0ne2fb3liSfS/K4JB9v++XKJIt76j89ycokm1p7jh3lMViT5M1JrgV+nmTeoOfsdknemuQnbf9clWSvIR6CoV4rL0/yH0nel+QO4NQMOkIf63MsyQ7pXjd3pHutXJlk4Uj7ZFarKv8m4Q9YAzx/UNnLgW8OVQf4NvDSNrwTcFAbXgwUMK9nvr8GVgNPanU/DXy0TdsH2AI8F9ierhvm1z3rObWNH0MX/jsCfwQcBMxr67sBeH3P+gq4CHgs8AzgPuCytv5dgOuBZcM8DsO2tWfZTxnhcXwJ8LjWtjcAtwE7DFP3VuBP2vB8YP82fDCwbtDjfjmwEFgEbAC+C+wH7AB8GThlqHmH2G+nAh8btL07A48C3g9c0zPtI8A7e8YfAVwF/H3bV08CbgIOa9PPAL4B7AbsBfxgcFt6lrU7cDfwn9tj9bq2n1/Zph/d9sPvtelvA77Vpu0G3Am8tE07vo0/rk3/apv3yT37+8fA81v984APt7qPAdYCr2jT9gN+Buwz1GPQ83he07ZxxyEe478Dvg88DQjwzIG2DVrOYh7+Wnk5cD/w2taeHXn463DY51jv/gVeBXwOeDTdh5g/Ah7b7/eaqfrzyGByfbZ9gticZDPwwRHq/hp4SpLdq2pLVV0+Qt0TgPdW1U1VtQV4C3BcukPgFwGfq6pvVtWv6N5oBt9w6ttV9dmqerCqflFVV1XV5VV1f1WtAf4f8GeD5nl3Vd1dVdfRvSld2tZ/F/BvdC/68bZ1VFX1saq6o7XtH+jeZIfr//81sE+Sx1bVnVX13REW/Y9VdXtVrad7w72iqq6uql8Cnxlhe0Zr77lVdU9V3Uf3RvLM9By1DfIsYEFVnVZVv6qqm4B/Bo5r048FTq+qTVW1FjhzhFUfCVxXVZ+uqvtb3dt6pv8N8L+r6oY2/X8B+7ajg6OAG6vqo+1x/gTwQ+CFPfN/uKp+0rO/f1JV/96W9a889Hi9AFhTVR9uy7oa+BTw4lEeujOram1V/WKIaa8E3lZVP6rO96rqjlGW1+unVfWPrT0PW/44nmO/pguNp1TVA+11c/c42jGrGAaT65iq2nXgD3j1CHVPBJ4K/LAdfr5ghLqPB27pGb+F7lPNwjZt7cCEqroXGPzCWds7kuSpST7fujfupnuj2H3QPLf3DP9iiPGdJtDWUSV5Y+vauKsF6i5DtG3Af6F7U7yldZE8e4RFT3R7RmrrdknOaN0Zd9N9umWE9j4RePygDwxv5aHHZqt9ydaP42CD93sBvd1bTwQ+0LOeTXSfshfx8H00sK5FPeNjfbyeCBw4aJtOAH53hLbDoOfkIHsBQ3YjjtFIyx7Pc+yjwCXA+Ul+muTdSR65De2a0QyDPqmqG6vqeOB3gHcBFyZ5DA//VA/wU7oX3YAn0B0K307XVbLnwIQkO9J9mtlqdYPGz6L7JLikqh5L94aUiW/NmNs6otZ3+ya6T8jzW6DeNVzbqurKqjqa7jH8LHDBtjUdgJ/TdQsMtGk7YMEwdf+Krjvm+XRvKIsHZhto4qD6a4Gbez8wVNXOVXVkm34r3RvhgCeM0M7B+z29421drxq0rh2r6ls8fB8NrGv9COsbzlrga4PWs1NV/W2bPtxtkUe6XfJaui6q0Yx72eN5jlXVr6vq7VW1D/AcuqOgl42hXbOSYdAnSV6SZEFVPQhsbsUPAhvb/yf1VP8E8N+T7J1kJ7pP8p9sh+wXAi9sJwS3p+uqGO2NfWe6/uYtSZ4O/O0o9cdjpLaOZme64NgIzEvy93TnLR4myfbprh/fpap+Tbc9D05C+38M7JDkqPYp8G103QjDtfc+uiOxR9Nta6/b2Xo/fge4p5083bEdWfx+koETxRcAb0kyP8medP3ew/kC8AdJjmldcCex9afxf2rLGjihvUuSga6bLwJPTfJX7eTtf6U79/T5EdY3nM+3Zb00ySPb37OS/N4wj8FYfAh4R5Il6fxhksEfcGDo18poxvMc+/Mkf9A+ENxN1200Gc+xGckw6J/DgevSXWHzAeC41p9/L3A68B/tsPsg4Fy6Q9avAzcDv6S9UbQ+/dcC59N9WtxCd4L0vhHW/Ua6T7X30PVZf3ISt2vYto7BJcCX6N6Qb2nzjnTI/1JgTeui+Ru67olt0vrIX033hrSe7khhuC9+ndfauZ7uJOvg8z7n0J3T2Jzks1X1AN2ny33pHpuftfUMnGN4e1vezcCldI/jcO38GV2//LvpwmgfYBVtv1fVZ+iOOM9vj88PgCPatDtaO97Q5n0T8IK2zHGpqnuAQ+nOe/yU7rzFu3goQLd6DMa42PfSBeOldG/C59CdCB687qFeK6MZz3Psd+k+bN1Nd5HF1xhhn8x26boaNVe0T+Ob6bqAbu53ezQ9kjyCLrROqKqv9Ls9mn08MpgDkrwwyaPbOYf30F2Wt6a/rdJUS3JYkl2TPIqHzvuMdFWaNCzDYG44mu4Q/afAErouJw/55r5n01118zO6y0KPGeZSTWlUdhNJkjwykCR1XwaasXbfffdavHhxv5shSbPKVVdd9bOqGu77MUOa0WGwePFiVq1a1e9mSNKskmSkb68PyW4iSZJhIEkyDCRJGAaSJMYYBu1bjhcm+WG79euzk+yW7teNbmz/B34pKknOTLI63S9L7d+znGWt/o1Jlk3VRkmSxmesRwYfAL5UVU+n+9WhG4CTgcuqagndr2Cd3OoeQfct2CXAcrrbJZNkN+AU4EDgAOCUgQCRJPXXqGGQ7leb/pTuzoG0X2jaTHcLhBWt2gq6n1WklZ/XfqHocmDXJHsAhwEr26843QmspLtzpySpz8ZyZLA33b2/P5zk6iQfajdEW1hVt7Y6t/HQrzUtYutbwq5rZcOVbyXJ8nQ/FL9q48aN49saSdKEjCUM5gH7A2dV1X5093c/ubdCuynapNzkqKrOrqqlVbV0wYJxfYFOkjRBY/kG8jpgXVVd0cYvpAuD25PsUVW3tm6gDW36erb+6b49W9l64OBB5V+deNO3tvjkL0zWoqbcmjOO6ncTJGkrox4ZVNVtwNokT2tFh9D9qtPFwMAVQcuAi9rwxcDL2lVFBwF3te6kS4BD20/6zaf7daRLJm9TJEkTNdZ7E70W+Hj7jd2bgFfQBckFSU6k+/m4Y1vdLwJHAquBe1tdqmpTkncAV7Z6p1XVpknZCknSNhlTGFTVNcDSISYdMkTdovtx7qGWcy7db+RKkmYQv4EsSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRJjDIMka5J8P8k1SVa1st2SrExyY/s/v5UnyZlJVie5Nsn+PctZ1urfmGTZ1GySJGm8xnNk8OdVtW9VLW3jJwOXVdUS4LI2DnAEsKT9LQfOgi48gFOAA4EDgFMGAkSS1F/b0k10NLCiDa8AjukpP686lwO7JtkDOAxYWVWbqupOYCVw+DasX5I0ScYaBgVcmuSqJMtb2cKqurUN3wYsbMOLgLU9865rZcOVbyXJ8iSrkqzauHHjGJsnSdoW88ZY77lVtT7J7wArk/ywd2JVVZKajAZV1dnA2QBLly6dlGVKkkY2piODqlrf/m8APkPX53976/6h/d/Qqq8H9uqZfc9WNly5JKnPRg2DJI9JsvPAMHAo8APgYmDgiqBlwEVt+GLgZe2qooOAu1p30iXAoUnmtxPHh7YySVKfjaWbaCHwmSQD9f+lqr6U5ErggiQnArcAx7b6XwSOBFYD9wKvAKiqTUneAVzZ6p1WVZsmbUskSRM2ahhU1U3AM4covwM4ZIjyAk4aZlnnAueOv5mSpKnkN5AlSYaBJMkwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgSWIcYZBkuyRXJ/l8G987yRVJVif5ZJLtW/mj2vjqNn1xzzLe0sp/lOSwyd4YSdLEjOfI4HXADT3j7wLeV1VPAe4ETmzlJwJ3tvL3tXok2Qc4DngGcDjwwSTbbVvzJUmTYUxhkGRP4CjgQ208wPOAC1uVFcAxbfjoNk6bfkirfzRwflXdV1U3A6uBAyZjIyRJ22asRwbvB94EPNjGHwdsrqr72/g6YFEbXgSsBWjT72r1f1M+xDy/kWR5klVJVm3cuHEcmyJJmqhRwyDJC4ANVXXVNLSHqjq7qpZW1dIFCxZMxyol6bfevDHU+WPgPyU5EtgBeCzwAWDXJPPap/89gfWt/npgL2BdknnALsAdPeUDeueRJPXRqEcGVfWWqtqzqhbTnQD+clWdAHwFeFGrtgy4qA1f3MZp079cVdXKj2tXG+0NLAG+M2lbIkmasLEcGQznzcD5Sd4JXA2c08rPAT6aZDWwiS5AqKrrklwAXA/cD5xUVQ9sw/olSZNkXGFQVV8FvtqGb2KIq4Gq6pfAi4eZ/3Tg9PE2UpI0tfwGsiTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIYQxgk2SHJd5J8L8l1Sd7eyvdOckWS1Uk+mWT7Vv6oNr66TV/cs6y3tPIfJTlsqjZKkjQ+YzkyuA94XlU9E9gXODzJQcC7gPdV1VOAO4ETW/0TgTtb+ftaPZLsAxwHPAM4HPhgku0mc2MkSRMzahhUZ0sbfWT7K+B5wIWtfAVwTBs+uo3Tph+SJK38/Kq6r6puBlYDB0zKVkiStsmYzhkk2S7JNcAGYCXwE2BzVd3fqqwDFrXhRcBagDb9LuBxveVDzCNJ6qMxhUFVPVBV+wJ70n2af/pUNSjJ8iSrkqzauHHjVK1GktRjXFcTVdVm4CvAs4Fdk8xrk/YE1rfh9cBeAG36LsAdveVDzNO7jrOramlVLV2wYMF4midJmqCxXE20IMmubXhH4C+AG+hC4UWt2jLgojZ8cRunTf9yVVUrP65dbbQ3sAT4zmRtiCRp4uaNXoU9gBXtyp9HABdU1eeTXA+cn+SdwNXAOa3+OcBHk6wGNtFdQURVXZfkAuB64H7gpKp6YHI3R5I0EaOGQVVdC+w3RPlNDHE1UFX9EnjxMMs6HTh9/M2UJE0lv4EsSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkxvazl5pki0/+Qr+bMCZrzjiq302QNE08MpAkGQaSJMNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkxhAGSfZK8pUk1ye5LsnrWvluSVYmubH9n9/Kk+TMJKuTXJtk/55lLWv1b0yybOo2S5I0HmM5MrgfeENV7QMcBJyUZB/gZOCyqloCXNbGAY4AlrS/5cBZ0IUHcApwIHAAcMpAgEiS+mvUMKiqW6vqu234HuAGYBFwNLCiVVsBHNOGjwbOq87lwK5J9gAOA1ZW1aaquhNYCRw+qVsjSZqQcZ0zSLIY2A+4AlhYVbe2SbcBC9vwImBtz2zrWtlw5YPXsTzJqiSrNm7cOJ7mSZImaMxhkGQn4FPA66vq7t5pVVVATUaDqursqlpaVUsXLFgwGYuUJI1iTGGQ5JF0QfDxqvp0K769df/Q/m9o5euBvXpm37OVDVcuSeqzsVxNFOAc4Iaqem/PpIuBgSuClgEX9ZS/rF1VdBBwV+tOugQ4NMn8duL40FYmSeqzsfyewR8DLwW+n+SaVvZW4AzggiQnArcAx7ZpXwSOBFYD9wKvAKiqTUneAVzZ6p1WVZsmZSskSdtk1DCoqm8CGWbyIUPUL+CkYZZ1LnDueBooSZp6fgNZkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiTGEAZJzk2yIckPesp2S7IyyY3t//xWniRnJlmd5Nok+/fMs6zVvzHJsqnZHEnSRIzlyOAjwOGDyk4GLquqJcBlbRzgCGBJ+1sOnAVdeACnAAcCBwCnDASIJKn/Rg2Dqvo6sGlQ8dHAija8Ajimp/y86lwO7JpkD+AwYGVVbaqqO4GVPDxgJEl9Mm+C8y2sqlvb8G3Awja8CFjbU29dKxuu/GGSLKc7quAJT3jCBJunybD45C/0uwljtuaMo/rdBGlW2+YTyFVVQE1CWwaWd3ZVLa2qpQsWLJisxUqSRjDRMLi9df/Q/m9o5euBvXrq7dnKhiuXJM0AEw2Di4GBK4KWARf1lL+sXVV0EHBX6066BDg0yfx24vjQViZJmgFGPWeQ5BPAwcDuSdbRXRV0BnBBkhOBW4BjW/UvAkcCq4F7gVcAVNWmJO8Armz1TquqwSelJUl9MmoYVNXxw0w6ZIi6BZw0zHLOBc4dV+skSdPCbyBLkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJMbwewbSbLD45C/0uwljsuaMo/rdBGlIHhlIkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEn7pTJpWs+XLceAX5H7beGQgSZr+I4MkhwMfALYDPlRVZ0x3GySNbrYcxXgEMzmm9cggyXbA/wWOAPYBjk+yz3S2QZL0cNN9ZHAAsLqqbgJIcj5wNHD9NLdD0hwxW45gZrrpDoNFwNqe8XXAgb0VkiwHlrfR+5L8YJra1g+7Az/rdyOmkNs3u83l7ZvL2wbwtPHOMOOuJqqqs4GzAZKsqqqlfW7SlHH7Zje3b/aay9sG3faNd57pvppoPbBXz/ierUyS1EfTHQZXAkuS7J1ke+A44OJpboMkaZBp7SaqqvuTvAa4hO7S0nOr6roRZjl7elrWN27f7Ob2zV5zedtgAtuXqpqKhkiSZhG/gSxJMgwkSTM4DJIcnuRHSVYnObnf7ZlsSdYk+X6SayZyGdhMk+TcJBt6vxeSZLckK5Pc2P7P72cbJ2qYbTs1yfq2/65JcmQ/27gtkuyV5CtJrk9yXZLXtfK5sv+G2745sQ+T7JDkO0m+17bv7a187yRXtPfQT7aLdoZfzkw8Z9BuW/Fj4C/ovph2JXB8Vc2ZbyonWQMsrao58cWXJH8KbAHOq6rfb2XvBjZV1Rkt0OdX1Zv72c6JGGbbTgW2VNV7+tm2yZBkD2CPqvpukp2Bq4BjgJczN/bfcNt3LHNgHyYJ8Jiq2pLkkcA3gdcB/wP4dFWdn+SfgO9V1VnDLWemHhn85rYVVfUrYOC2FZqhqurrwKZBxUcDK9rwCroX4KwzzLbNGVV1a1V9tw3fA9xAd7eAubL/htu+OaE6W9roI9tfAc8DLmzlo+6/mRoGQ922Ys7svKaAS5Nc1W7BMRctrKpb2/BtwMJ+NmYKvCbJta0baVZ2oQyWZDGwH3AFc3D/Ddo+mCP7MMl2Sa4BNgArgZ8Am6vq/lZl1PfQmRoGvw2eW1X7093B9aTWFTFnVdcfOfP6JCfuLODJwL7ArcA/9Lc52y7JTsCngNdX1d290+bC/hti++bMPqyqB6pqX7q7OhwAPH28y5ipYTDnb1tRVevb/w3AZ+h24Fxze+uvHei33dDn9kyaqrq9vQAfBP6ZWb7/Wl/zp4CPV9WnW/Gc2X9Dbd9c24cAVbUZ+ArwbGDXJANfLB71PXSmhsGcvm1Fkse0E1kkeQxwKDAX7856MbCsDS8DLupjWybVwJtk85fM4v3XTkCeA9xQVe/tmTQn9t9w2zdX9mGSBUl2bcM70l14cwNdKLyoVRt1/83Iq4kA2mVe7+eh21ac3ucmTZokT6I7GoDuliD/Mtu3L8kngIPpbg18O3AK8FngAuAJwC3AsVU1607EDrNtB9N1LxSwBnhVT//6rJLkucA3gO8DD7bit9L1q8+F/Tfc9h3PHNiHSf6Q7gTxdnQf8C+oqtPa+8z5wG7A1cBLquq+YZczU8NAkjR9Zmo3kSRpGhkGkiTDQJJkGEiSMAwkSRgGkiQMA0kS8P8BPUNqkTWh2ywAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate random integers between 1 and 6 for a sample size of 10000\n",
    "y = np.random.geometric(.2, 10000)\n",
    "\n",
    "plt.hist(y)\n",
    "plt.title(\"Histogram of a simulated geometric trials\")\n",
    "plt.xlim([0,30]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "https://docs.scipy.org/doc/numpy/reference/routines.random.html\n",
    "https://matplotlib.org/api/_as_gen/matplotlib.pyplot.html\n",
    "https://stackoverflow.com/questions/6294179/how-to-find-all-occurrences-of-an-element-in-a-list\n",
    "https://en.wikipedia.org/wiki/Geometric_distribution\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}