RiskEngineering notebooks

basic-statistics.ipynb

coins-dice.ipynb

correlations.ipynb

Poisson-processes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Poisson processes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://risk-engineering.org/static/img/logo-RE.png\" width=\"100\" alt=\"\" style=\"float:right;margin:15px;\">\n",
    "\n",
    "This notebook is an element of the [risk-engineering.org courseware](https://risk-engineering.org/). It can be distributed under the terms of the [Creative Commons Attribution-ShareAlike licence](https://creativecommons.org/licenses/by-sa/4.0/).\n",
    "\n",
    "Author: Eric Marsden <eric.marsden@risk-engineering.org>. \n",
    "\n",
    "---\n",
    "\n",
    "This notebook contains an introduction to the use of Python and SciPy to analyze data concerning Poisson processes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import scipy.stats\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"bmh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Earthquakes\n",
    "\n",
    "Suppose we live in an area where there are typically 0.03 earthquakes of intensity 5 or more per year. Assume that earthquake arrival is a [Poisson counting process](https://en.wikipedia.org/wiki/Poisson_point_process) (the interval between earthquakes follows an exponential distribution, and events are independent). \n",
    "\n",
    "Let's simulate the random intervals between the next earthquakes of intensity 5 or greater."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.23658881e+02, 4.66527507e-01, 1.54656670e+01, 7.25553121e+01,\n",
       "       1.20041381e+01, 6.26361707e+01, 3.96577482e-01, 3.04319082e+01,\n",
       "       4.66770724e+01, 2.01411541e+01, 1.70674321e+01, 4.00946655e+01,\n",
       "       5.45694133e+01, 4.23367155e+00, 7.46744783e+01, 1.18744067e+01,\n",
       "       2.78208171e+01, 6.53253854e+00, 7.15911038e+00, 1.18011111e+01,\n",
       "       1.68632189e+01, 3.02152066e+01, 8.36335319e+00, 5.26832457e+01,\n",
       "       2.64991671e+01, 6.75709511e+01, 7.40832911e+00, 1.99123785e+02,\n",
       "       1.30016001e+02, 1.51028636e+01, 4.57147639e+00, 8.98047189e+01,\n",
       "       1.45742938e+01, 1.08856820e+02, 1.51040948e-01, 6.74018364e+00,\n",
       "       8.58840371e-01, 2.10760556e+00, 8.00526734e+01, 2.34880577e+01,\n",
       "       1.00823265e+01, 1.19254111e+01, 4.05922179e+01, 5.82675162e+01,\n",
       "       1.14731748e+01, 2.48160927e+01, 4.77813117e+01, 4.30768342e+01,\n",
       "       1.78342941e+01, 9.39478890e+01])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intervals = scipy.stats.expon(scale=1/0.03).rvs(50)\n",
    "intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do we really have an average rate of 0.03 per year?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026523651792807357"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(intervals) / intervals.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's the probability of seeing more than one earthquake of intensity greater than 5 in the next year? It we call $X$ the random variable \"number of earthquakes next year\", then we are interested in $P(X > 1)$. \n",
    "\n",
    "$P(X > 1)$ is $1 - P(X=0) - P(X=1)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9704455335485082"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# P(X = 0)\n",
    "scipy.stats.poisson(0.03).pmf(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.029113366006455248"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.poisson(0.03).pmf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0004411004450365977"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - scipy.stats.poisson(0.03).pmf(0) - scipy.stats.poisson(0.03).pmf(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's the probability of at least one 5+ earthquake in the next 20 years? It's 1 - the probability of zero earthquakes in the next 20 years, which is 1 minus the probability that the next arrival time is greater than 20. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4547"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Monte Carlo simulation\n",
    "N = 10000\n",
    "count = 0\n",
    "for i in range(N):\n",
    "    if scipy.stats.expon(scale=1/0.03).rvs() > 20:\n",
    "        count += 1\n",
    "1 - count / float(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4579900279842801"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.sum([scipy.stats.expon(scale=1/0.03).pdf(i) for i in range(20)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8111243971624382"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - scipy.stats.poisson(1/(20 * 0.03)).pmf(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cumulative distribution function for the arrival times tells us the probability that the arrival time for the next earthquake does not exceed a given number of years. Obviously, it's an increasing function which asymptotically reaches 1 as the arrival time increases (the more time passes, the more likely it is that an earthquake will occur, but since we're talking about a probability we are necessarily less than or equal to 1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VfWdx/H3l0ASkkASkrBkAQIEwiaIUXBXxL2IU23H3Vor02mttmo7zmi1dtVuto5MK1q12CpVx7HUIpQqVauiBNn3AMGEkAUSQlay3N/8kegTMZgL3OTknvt5PU+e5Nx7uPkcT/Lxl989iznnEBERf+njdQAREQk9lbuIiA+p3EVEfEjlLiLiQyp3EREfUrmLiPhQl+VuZk+aWbmZbTjC82Zmj5hZgZmtM7NpoY8pIiJHI5iR+9PARZ/x/MVATvvHXOA3xx9LRESOR5fl7px7E6j8jFXmAAtcmxVAkpkNC1VAERE5eqGYc88AijosF7c/JiIiHukbgtewTh7r9JoGZjaXtqkb4uPjT8rNzQ3BtxcRiRyrVq3a55xL62q9UJR7MZDVYTkTKOlsRefcfGA+QF5ensvPzw/BtxcRiRxmtjuY9UIxLbMIuKH9qJkZQLVzbm8IXldERI5RlyN3M3sOOAdINbNi4H6gH4Bz7rfAYuASoACoB27qrrAiIhKcLsvdOXd1F8874OshSyQiIsdNZ6iKiPiQyl1ExIdU7iIiPqRyFxHxoVAc5y4iIp1wzlFV38ze6gbKDx6i7GAjpQcbmZk7mBMyk7r1e6vcRUSOgXOOyrom9lY3UnKgoe1zdQOl1Y3srW5kb3UDZQcP0dQS+NS/TU2IUbmLiHihNeAoO9hIcVUDxVX1FFc1sKeqgZLqts97DjRw6LDi7hdlDE2MZVhif07MSmZYYixDBsYytMPntIQYovt2/4y4yl1EIlZ1fTO7K+v4sLKeDyvrKWr/XFzVQMmBBppbP3mZrNSEaDKS+jNu6ABm5g4mPak/6UmxpCf1Z1hif1Lio+nTp7PLbfU8lbuI+NZHUyeF++vYta+e3fvrKNzf9nn3/nqqG5o/sX5KfDSZg+KYnJHIJZOHkZncn8zkODKT+5OR1J/YflEebcnRU7mLSNhrbG5l1746dlbUsaOitu3rfXXsqqjlYGPLx+v1MchMjmNEShyzpwxjxKB4hqfEMXxQHFmD4kiI8U8l+mdLRMT3quubKaiooaC8lu1ltRRU1LKjopbiqgZchxmUjKT+ZKfGc9nUdLJTE8hOjWNkSjyZyXE9Mt/dG6jcRaTXqTvUwvbyWraWHmRraS3by2vYWlpDec2hj9eJ6duHUWkJTM1K5oppmYwZnMCo1ASyU+PpHx0+0yfdReUuIp4JBBxFVfVs3nuQzXtr2FJ6kC2lNezeX//xOrH9+pAzeABn5KQydsgAcgYnkDN4ABnJ/YnqJW9e9kYqdxHpEU0tAbaV1bCxpJpNJQfZWNJW5LWH2ubE+xiMTI1nUnoiV0zLZNzQAeQOHUBWclyvOQIlnKjcRSTkmloCbC2tYd2eA2zYU836PdVsLa35+NDC+OgoJqQP5IppGYwfNpAJ6QMZO2RAWB2N0tup3EXkuAQCjp376lhTdIC1RQdYV3yAzXtraGptO8EnsX8/JmUM5MtnZDMpPZFJGYmMGKTReHdTuYvIUTlQ38TqogOs3l3F6qIDrCk6QE374Ybx0VFMzkzkptNHMjkzkRMyksga1B8zFXlPU7mLyBE55yjcX8/Kwko+2F1F/u4qCsprgbY58nFDBzJ7SjpTs5KYmpXE6LQEvcnZS6jcReRjrQHH5r0HeX9XJe/vqiR/dyX7apuAtumVk0Ykc/nUdKaNSGZKZhLxPjrpx2+0Z0QiWGvAsankICt27mfFzv28X1j58RRLRlJ/zspJI2/kIE4emczotATNk4cRlbtIBHHOsb28lrcL9vHOjrZC/6jMR6XG87kThjE9O4WTsweRkdTf47RyPFTuIj5XfrCRt7bv458FbR8V7Wd5jkiJ49LJwzh1dAozRqUwZGCsx0kllFTuIj5zqKWV/MIq3thWwZvbKthSWgO0XfHwtDGpnDkmldPGpJCZHOdxUulOKncRH9hzoIHlW8r5x9Zy3i7YT0NzK9FRfcgbmczdF+dyZk4q44cO1Jx5BFG5i4Sh1oBjTdEBXttcxmuby9la1jY6z0zuz5UnZXLOuDRmjErR0SwRTHteJEw0Nrfy1vZ9/H1TGa9tKWNfbRNRfYyTRybzX5fkMjN3MKPTEnTCkAAqd5FerbqhmeVbylmyoZQ3tlXQ0NzKgJi+nJM7mFnjB3PO2MEkxvXzOqb0Qip3kV6mqq6JZZvKeHXDXv5ZsI/mVseQgTFceVImF0wcwvTslIi54YQcO5W7SC9QXd/M0k2lvLJuL28X7KM14MhM7s9Np2dz4cShnJiVpDdD5aio3EU8Ut/Uwt83l7NoTQlvbCunubWt0L9yZjafm5zOpIyBmj+XY6ZyF+lBLa0B3t6xn5dX72HpxlLqm1oZOjCWG08dyewp6ZyQmahCl5BQuYv0gE0lB3npg2JeXlPCvtpDDIzty5yp6cyZmsEpIwdpykVCTuUu0k0q65p4efUeXlhVzOa9B+kXZczMHcznp7Udhx7TV3cdku4TVLmb2UXAr4Eo4Ann3IOHPT8c+D2Q1L7O3c65xSHOKtLrtQYcb22v4Pn8IpZtKqO51XFCZiLfnzOR2Sekkxwf7XVEiRBdlruZRQHzgPOBYmClmS1yzm3qsNq9wPPOud+Y2QRgMTCyG/KK9Ep7qxv408oiXsgvZs+BBgbFR3PDqSP5Ql4muUMHeh1PIlAwI/dTgALn3E4AM1sIzAE6lrsDPvoJTgRKQhlSpDcKBBxvbK/gjys+5PUtZTjgjDGp3HPpeGaNH6Jj0cVTwZR7BlDUYbkYmH7YOt8D/mZm3wDigVmdvZCZzQXmAgwfPvxos4r0ClV1TTyfX8Qf3ttNUWUDqQnR/Ps5o7nq5OFkDdKVFqV3CKbcO3sb3x22fDXwtHPuF2Z2KvCMmU1yzgU+8Y+cmw/MB8jLyzv8NUR6tY0l1Tz9diF/XltCU0uAU7IH8Z0Lc7lw4lCN0qXXCabci4GsDsuZfHra5WbgIgDn3LtmFgukAuWhCCnildaAY9mmUp78ZyHvF1bSv18UXzgpk+tPHaG5dOnVgin3lUCOmWUDe4CrgGsOW+dD4DzgaTMbD8QCFaEMKtKTag+18PzKIp56ZxdFlQ1kJvfnnkvG88W8LF2oS8JCl+XunGsxs1uBpbQd5vikc26jmX0fyHfOLQLuBB43s2/RNmXzJeecpl0k7JQdbOSptwv543u7qWlsIW9EMvdcMp7zJwwlSicaSRgJ6jj39mPWFx/22H0dvt4EnB7aaCI9Z0dFLY+9sYP/W72H1oDj4knDuOWsUUzNSvI6msgx0RmqEtHWF1fzP/8oYMnGUqKj+nD1KcO5+YxsRqTEex1N5Lio3CUirdpdySOvFfDGtgoGxPbla+eM5qbTs0lNiPE6mkhIqNwlory/q5Jf/X0b7+zYz6D4aL594TiuP3UEA2P1Jqn4i8pdIkLHUk9NiOHeS8dzzfThxEXrV0D8ST/Z4mtrig7wi79t5a3t+z4u9Wunj6B/tK7IKP6mchdf2lpaw8//tpVlm8pIjuvHf12Sy/UzRqrUJWKo3MVXiqvq+eWybfzf6j0kRPflzvPHctMZ2STE6EddIot+4sUXDtQ38ejrBSx4dzcYzD1zFF89e7Suny4RS+UuYa2xuZUF7xby6OsF1Bxq4cppmdxxwViGJfb3OpqIp1TuEpaccyxeX8qDSzZTVNnAOePSuPviXF3MS6Sdyl3CzoY91Tzwl42sLKwid+gAnrn5FM7MSfM6lkivonKXsLGv9hA/X7qVP+UXMSgumh//y2T+9eQsXdBLpBMqd+n1WloDPLNiN79cto2GplZuPj2b22bl6KxSkc+gcpdebWVhJd99eQNbSms4MyeV+2dPZMzgBK9jifR6KnfplSrrmvjx4s28uKqY9MRYfnvdNC6cOBQzTcGIBEPlLr1KIOB4cVUxP351M7WNLXz17NHcdt4YXQNG5CjpN0Z6jR0VtfznS+t5f1clJ49M5kf/MpmxQwZ4HUskLKncxXNNLQHmv7mDR14vILZvHx78/GS+mJdFHx0FI3LMVO7iqfXF1Xz7xbVsKa3h0snDuP+yCQweEOt1LJGwp3IXTzQ2t/LIa9t57M2dpMRH8/gNeZw/YYjXsUR8Q+UuPW59cTV3PL+G7eW1fOGkTO69dAKJcTpmXSSUVO7SY5paAjy6vIB5ywtIS4jhqZtO5txxg72OJeJLKnfpEQXlNXzzT2vYsOcgnz8xg/tnT9RoXaQbqdylWwUCjgXvFvKTV7cQH9OX3153EhdNGup1LBHfU7lLtymvaeSuF9bx5rYKzh2XxkNXnqAjYUR6iMpdusXyLeXc9cJaag+18IPLJ3Hd9OG6dIBID1K5S0gdamnlwVe38NTbheQOHcDCuTPI0VmmIj1O5S4hU7ivjluf+4ANew7ypdNGcvfFucT2i/I6lkhEUrlLSLyyroS7/3c9fQzmX38SF0zUm6YiXlK5y3E51NLKj/66mQXv7mba8CQeufpEMpPjvI4lEvFU7nLMiqvq+fqzq1lbdIBbzszmOxfl0i+qj9exRASVuxyjN7dVcNvC1bS2On573TQumjTM60gi0kFQwywzu8jMtppZgZndfYR1vmhmm8xso5k9G9qY0ls455i3vIAbn3qfoQNjWfSNM1TsIr1QlyN3M4sC5gHnA8XASjNb5Jzb1GGdHOA/gdOdc1VmpguG+FDtoRbufH4NSzeWMXtKOg9dMVl3SBLppYL5zTwFKHDO7QQws4XAHGBTh3VuAeY556oAnHPloQ4q3tq9v45bFuSzo6KOey8dz81nZOukJJFeLJhyzwCKOiwXA9MPW2csgJm9DUQB33POLQlJQvHcP7fv4+vPfoAZLPjyKZw+JtXrSCLShWDKvbPhmevkdXKAc4BM4C0zm+ScO/CJFzKbC8wFGD58+FGHlZ7lnGPBu7t54C8bGTM4gSduOJnhKTrMUSQcBPOGajGQ1WE5EyjpZJ0/O+eanXO7gK20lf0nOOfmO+fynHN5aWlpx5pZekBza4Dv/nkD9y/ayMzcIbz0tdNV7CJhJJhyXwnkmFm2mUUDVwGLDlvnZeBcADNLpW2aZmcog0rPqW5o5qanVvKHFR/yb2eP4rHrTyIhRm+cioSTLn9jnXMtZnYrsJS2+fQnnXMbzez7QL5zblH7cxeY2SagFfi2c25/dwaX7lFUWc+Xn15J4f46fnblCXwhL6vrfyQivY45d/j0ec/Iy8tz+fn5nnxv6dzaogPc/Pt8mlpaeez6PE4dneJ1JBE5jJmtcs7ldbWe/tYWAJZtKuMbz31A2oAYFs6dzpjBukyvSDhTuQt/fG833315A5MzEnnixpNJGxDjdSQROU4q9wjmnOPhZdt45PUCzh2Xxrxrp+mMUxGf0G9yhGppDXDvyxtYuLKIf83L4kf/Mom+uqKjiG+o3CNQY3Mrty9czdKNZXxj5hjuOH+sLiUg4jMq9whT09jM3AWreHfnfu6fPYGbTs/2OpKIdAOVewSpqmvixqfeZ1PJQX71r1O5/MQMryOJSDdRuUeI8oONXPe79yjcX8/8G05iZu4QryOJSDdSuUeA4qp6rn3iPSpqDvH0TSdz2mhd1VHE71TuPrd7fx3XPP4eNY3N/OEr05k2PNnrSCLSA1TuPrajopZrHl9BU0uAZ2+ZwaSMRK8jiUgPUbn71PayGq554j0CAcdzc2eQO3Sg15FEpAep3H1oe1kNVz++AjNj4dwZ5AzRdWJEIo3K3WcKymu4+vH3Pi720WkJXkcSEQ/ofHMfKSiv5ar572EGz92iYheJZCp3nyjcV8c1j68A2op9zGAVu0gkU7n7wEfHsbcEHM/eMl3FLiIq93BXWt348XHsz9x8CmP15qmIoDdUw9r+2kNc+8QKKuua+MNXpjMxXcexi0gbjdzD1MHGZm586n2Kqxr43Y15TM1K8jqSiPQiKvcw1NDUyleezmfL3hp+e91JTB+lG1mLyCdpWibMNLcG+NofV7FydyWPXHUi5+YO9jqSiPRCGrmHkUDA8R8vrmP51gp+ePkkZk9J9zqSiPRSKvcw8tCSLby0eg93nD+Wa6eP8DqOiPRiKvcw8cRbO3nszZ3ccOoIvjFzjNdxRKSXU7mHgUVrS/jhXzdz6eRh3D97om5mLSJdUrn3cit27ueu59dySvYgfvHFKUT1UbGLSNdU7r3YtrIa5i7IZ3hKHI9fn0dsvyivI4lImFC591LlBxu56amVxPSL4umbTiYxrp/XkUQkjKjce6H6pha+siCfqvomnvrSyWQmx3kdSUTCjMq9lwkEHN/60xo27KnmkatO1H1PReSYqNx7mYeWbGHpxjLuvXQCsyYM8TqOiIQplXsv8vzKIh57cyfXzxjBTaeP9DqOiISxoMrdzC4ys61mVmBmd3/GeleamTOzvNBFjAzv76rknpfXc2ZOKvfPnqBj2UXkuHRZ7mYWBcwDLgYmAFeb2YRO1hsA3Aa8F+qQfldUWc9X/7CKrOQ4Hr16Gn2j9AeViByfYFrkFKDAObfTOdcELATmdLLeD4CfAo0hzOd7dYdauGVBPi2tAZ64MU+HPIpISART7hlAUYfl4vbHPmZmJwJZzrlXQpjN95xz3Pn8WraV1TDv2mmMStO9T0UkNIIp984mf93HT5r1AR4G7uzyhczmmlm+meVXVFQEn9Kn5i0vYMnGUv7rkvGcmZPmdRwR8ZFgyr0YyOqwnAmUdFgeAEwC/mFmhcAMYFFnb6o65+Y75/Kcc3lpaZFdZq9tLuMXy7Zx+dR0bj4j2+s4IuIzwZT7SiDHzLLNLBq4Clj00ZPOuWrnXKpzbqRzbiSwArjMOZffLYl9YGdFLd9cuIaJ6QN58IoTdGSMiIRcl+XunGsBbgWWApuB551zG83s+2Z2WXcH9Jv6pha++odV9Ovbh8d0MTAR6SZB3UPVObcYWHzYY/cdYd1zjj+WPznnuPt/11NQXsuCL08nI6m/15FExKd0QHUPevqdQhatLeHOC8ZxRk6q13FExMdU7j0kv7CSH/11M7PGD+Hfzx7tdRwR8TmVew/YX3uIW59dTXpSf37xxSn00d2URKSbBTXnLscuEHB86/m1VNY38dK/n0Zif52BKiLdTyP3bvabN3bw5rYK7vvcBF2bXUR6jMq9G63YuZ9f/G0rl01J59rpw72OIyIRROXeTSrrmrh94WpGpsTz489P1olKItKjNOfeDZxz3PXCWqrqmnnySyeTEKP/zCLSszRy7wZPvl3I61vKuefS8UxM1zy7iPQ8lXuIrS+u5sFXN3P+hCHccOoIr+OISIRSuYdQfVMLty1cTWpCDD+7UhcEExHvaDI4hH7wyiYK99fx7FdmkBQX7XUcEYlgGrmHyJINpTz3fhFfPXs0p45O8TqOiEQ4lXsIlFY3cvdL65ickci3Zo31Oo6IiMr9eAUCbYc9HmoO8KurphLdV/9JRcR7aqLj9MyK3fyzYB/3XDqe0brBtYj0Eir347CjopafvLqZc8al6fICItKrqNyPUXNrgDv+tIbYflH8VPdBFZFeRodCHqP/Wb6DtcXVzLtmGoMHxnodR0TkEzRyPwYbS6r579e3M2dqOpeeMMzrOCIin6JyP0pNLQHufH4tyfHRPHDZRK/jiIh0StMyR+nR17ezpbSGx2/I01moItJraeR+FDbsqWbeP3bw+RMzOH/CEK/jiIgckco9SE0tAe56YS0p8dHcP1vTMSLSu2laJki/fWMHW0preOKGPBLjdJNrEendNHIPwrayGv779e3MnpLOLE3HiEgYULl3oTXg+M6L60iI6cv3Zk/wOo6ISFA0LdOFp98pZE3RAX591VRSEmK8jiMiEhSN3D9DUWU9P1+6lZm5g7lsSrrXcUREgqZyPwLnHN/98wbM4AeXT9K1Y0QkrKjcj+CVdXv5x9YK7rpgHBlJ/b2OIyJyVFTunaiub+aBv2zkhMxEbjxtpNdxRESOmt5Q7cSDSzZTVd/M7798ClF9NB0jIuEnqJG7mV1kZlvNrMDM7u7k+TvMbJOZrTOz18xsROij9oxVuyt57v0ibj4jm4npiV7HERE5Jl2Wu5lFAfOAi4EJwNVmdvgB36uBPOfcCcCLwE9DHbQntLQGuOf/NpCeGMvt5+V4HUdE5JgFM3I/BShwzu10zjUBC4E5HVdwzi13ztW3L64AMkMbs2c8/U4hW0pruG/2ROJjNGMlIuErmHLPAIo6LBe3P3YkNwOvdvaEmc01s3wzy6+oqAg+ZQ/YW93Aw8u2ce64NC6cqEsMiEh4C6bcO3tH0XW6otl1QB7ws86ed87Nd87lOefy0tLSgk/ZA374ymZaAo4HLtMx7SIS/oKZeygGsjosZwIlh69kZrOAe4CznXOHQhOvZ7y1vYK/rt/LneePZXhKnNdxRESOWzAj95VAjpllm1k0cBWwqOMKZnYi8BhwmXOuPPQxu09TS4DvLdrIiJQ4bjlrlNdxRERCostyd861ALcCS4HNwPPOuY1m9n0zu6x9tZ8BCcALZrbGzBYd4eV6naff2cWOijrunz2B2H5RXscREQmJoA4Jcc4tBhYf9th9Hb6eFeJcPaLsYCO//vt2zssdzMxcvYkqIv4R0Zcf+MnizTS3Ou7TddpFxGcittzzCyt5eU0Jc88axYiUeK/jiIiEVESWeyDgeOAvmxg6MJavnTva6zgiIiEXkeX+4gfFrN9Tzd0X5xIXrTNRRcR/Iq7caxqb+emSrUwbnsScqbq7koj4U8QNW+ct38G+2kP87sY8nYkqIr4VUSP3D/fX8+Q/d3HFtEymZCV5HUdEpNtEVLk/tGQLUX2Mb184zusoIiLdKmLKfdXuSv66fi9zzxrF0MRYr+OIiHSriCh35xw//OtmBg+I4d/O1vVjRMT/IqLcX1m3l9UfHuCuC8bp0EcRiQi+L/dDLa08tGQLuUMHcMVJYXmDKBGRo+b7cn/m3d0UVzVw76UTiOqjQx9FJDL4utyrG5p5dHkBZ41N44ycVK/jiIj0GF+X+2/+sYPqhmb+4yId+igikcW35b63uoGn3t7F5VMzmJie6HUcEZEe5dtyf3jZNpyDO84f63UUEZEe58ty31ZWw4urirn+1BFkDdINr0Uk8viy3H++dCtx0X35+rljvI4iIuIJ35X7mqID/G1TGbecOYpB8dFexxER8YTvyv1nS7eQEh/NzWdmex1FRMQzvir3twv28XbBfr527hgSYnSZARGJXL4pd+ccP126lfTEWK6dPtzrOCIinvJNuS/bVMbaogPcPiuH2H5RXscREfGUL8o9EHD8ctk2slPjuWKaLg4mIuKLcl+ysZQtpTXcfl4OfaN8sUkiIscl7JuwNeB4eNk2xgxOYPaUdK/jiIj0CmFf7q+sK2F7eS3fnJWjS/qKiLQL63JvaQ3w679vJ3foAC6ZNMzrOCIivUZYl/uf15Swc18d35w1lj4atYuIfCxsy72lNcCjywuYMGwgF04c4nUcEZFeJWzL/S/rSti1r47bzsvBTKN2EZGOgip3M7vIzLaaWYGZ3d3J8zFm9qf2598zs5GhDtpRa8Dx368XkDt0ABdM0KhdRORwXZa7mUUB84CLgQnA1WY24bDVbgaqnHNjgIeBh0IdtKNX1pWws6Jt1K65dhGRTwtm5H4KUOCc2+mcawIWAnMOW2cO8Pv2r18EzrNumiv5aNQ+bsgALpo4tDu+hYhI2Aum3DOAog7Lxe2PdbqOc64FqAZSQhHwcIvX76WgvJZvnDdGo3YRkSMIptw7a1B3DOtgZnPNLN/M8isqKoLJ9ynxMVGcP2GIjmsXEfkMwVz0vBjI6rCcCZQcYZ1iM+sLJAKVh7+Qc24+MB8gLy/vU+UfjJm5Q5iZqzdRRUQ+SzAj95VAjpllm1k0cBWw6LB1FgE3tn99JfC6c+6YyltERI5flyN351yLmd0KLAWigCedcxvN7PtAvnNuEfA74BkzK6BtxH5Vd4YWEZHPFtS96Jxzi4HFhz12X4evG4EvhDaaiIgcq7A9Q1VERI5M5S4i4kMqdxERH1K5i4j4kMpdRMSHzKvD0c2sAth9jP88FdgXwjjhQNscGbTNkeF4tnmEcy6tq5U8K/fjYWb5zrk8r3P0JG1zZNA2R4ae2GZNy4iI+JDKXUTEh8K13Od7HcAD2ubIoG2ODN2+zWE55y4iIp8tXEfuIiLyGcKu3Lu6WbcfmFmWmS03s81mttHMbm9/fJCZLTOz7e2fk73OGkpmFmVmq83slfbl7PYbrm9vvwF7tNcZQ8nMkszsRTPb0r6vT42Affyt9p/pDWb2nJnF+m0/m9mTZlZuZhs6PNbpfrU2j7T32TozmxaqHGFV7kHerNsPWoA7nXPjgRnA19u3827gNedcDvBa+7Kf3A5s7rD8EPBw+/ZW0XYjdj/5NbDEOZcLTKFt2327j80sA7gNyHPOTaLtEuJX4b/9/DRw0WGPHWm/XgzktH/MBX4TqhBhVe4Ed7PusOec2+uc+6D96xrafukz+OSNyH8PXO5NwtAzs0zgUuCJ9mUDZtJ2w3Xw3/YOBM6i7V4IOOeanHMH8PE+btcX6N9+x7Y4YC8+28/OuTf59J3ojrRf5wALXJsVQJKZheQeouFW7sHcrNtXzGwkcCLwHjDEObcX2v4HAAz2LlnI/Qr4DhBoX04BDrTfcB38t69HARXAU+1TUU+YWTw+3sfOuT3Az4EPaSv1amAV/t7PHznSfu22Tgu3cg/qRtx+YWYJwP8C33TOHfQ6T3cxs88B5c65VR0f7mRVP+3rvsA04DfOuROBOnw0BdOZ9nnmOUA2kA7E0zYtcTg/7eeudNvPebiVezA36/YFM+tHW7H/0Tn3UvvDZR/9ydb+udyrfCF2OnCZmRXSNtU2k7aRfFL7n+/gv31dDBQ7595rX36RtrL36z4GmAXscs5VOOeagZeA0/D3fv7IkfZrt3VauJV7MDfrDnvt882/AzZbZRaXAAABCUlEQVQ7537Z4amONyK/EfhzT2frDs65/3TOZTrnRtK2T193zl0LLKfthuvgo+0FcM6VAkVmNq79ofOATfh0H7f7EJhhZnHtP+MfbbNv93MHR9qvi4Ab2o+amQFUfzR9c9ycc2H1AVwCbAN2APd4naebtvEM2v40Wwesaf+4hLZ56NeA7e2fB3mdtRu2/RzglfavRwHvAwXAC0CM1/lCvK1Tgfz2/fwykOz3fQw8AGwBNgDPADF+28/Ac7S9p9BM28j85iPtV9qmZea199l62o4kCkkOnaEqIuJD4TYtIyIiQVC5i4j4kMpdRMSHVO4iIj6kchcR8SGVu4iID6ncRUR8SOUuIuJD/w8f33xgRhFYzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "support = numpy.linspace(0, 100, 100)\n",
    "plt.plot(support, scipy.stats.expon(scale=1/0.03).cdf(support));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4511883639059736"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.expon(scale=1/0.03).cdf(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $Y$ denote the amount of time (in years) until the next earthquake, and $N(t)$ the counting function for the number of earthquakes that have arrived by year $t$. Because $Y$ will be greater than $t$ if and only if no events occur within the next $t$ units of time, we have\n",
    "\n",
    "$P(Y>t)=P(N(t)=0) = \\text{scipy.stats.expon}(scale=1/0.03).cdf(Y)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.029113366006455248"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.poisson(0.03).pmf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00043670049009682845"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.poisson(0.03).pmf(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.sum([scipy.stats.poisson(0.03).pmf(i) for i in range(20)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Births in a hospital"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Births in a hospital occur randomly at an average rate of 1.8 births per hour.\n",
    "What is the probability of observing 4 births in a given hour at the hospital?\n",
    "\n",
    "Let X = number of births in a given hour\n",
    "(i) Events occur randomly\n",
    "⇒ X ∼ Poisson(1.8)\n",
    "(ii) Mean rate λ = 1.8\n",
    "\n",
    "We can now use the formula to calculate the probability of observing exactly 4\n",
    "births in a given hour\n",
    "\n",
    "Pr(X = 4) = 0.0723\n",
    "\n",
    "What about the probability of observing more than or equal to 2 births in a given\n",
    "hour at the hospital?\n",
    "\n",
    "We want Pr(X ≥ 2) = Pr(X = 2) + Pr(X = 3) + . . .\n",
    "\n",
    "but (trick) Pr(X ≥ 2) = Pr(X = 2) + Pr(X = 3) + . . . \n",
    "\n",
    "= 1 - Pr(X < 2) = 1 - Pr(X=0) - Pr(X=1) = 0.537"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sum of two Poisson variables"
   ]
  },
  {
   "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

Semmelweis-fever.ipynb

data-analysis-earthquakes.ipynb

data-analysis-speed-light.ipynb

data-analysis-weather.ipynb

dike-case-study.ipynb

distributions.ipynb

elementary-effects.ipynb

exercise-descriptive-statistics.ipynb

exercises-reliability-engineering.ipynb

helmet-data-bootstrapping.ipynb

intervention-analysis-hierarchical-model.ipynb

monte-carlo-hoop-strength.ipynb

monte-carlo-hydrogen-tank-failure-rust.ipynb

monte-carlo-hydrogen-tank-failure.ipynb

monte-carlo-integration.ipynb

monte-carlo-LHS.ipynb

monte-carlo-pi.ipynb

monte-carlo-project-risk.ipynb

monte-carlo-reliability.ipynb

monte-carlo-slope-stability.ipynb

poker.ipynb

regression-CCPP.ipynb

regression-piecewise-linear.ipynb

regression-smoking.ipynb

saint-petersburg.ipynb

sensitivity-analysis.ipynb

simple-descriptive-statistics.ipynb

stock-market.ipynb

uncertainty-propagation.ipynb