ingle/poisson_process_with_dead_time.ipynb

## poisson_process_with_dead_time.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from random import expovariate,uniform\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 1: using exponential increments with random dead times\n",
    "The model is\n",
    "$$\n",
    "T_n - T_{n-1} = X_n + U_n \n",
    "$$\n",
    "where $T_n$ is the $n$th detection event, $X_n$ are iid Exp$(\\lambda)$ and $U_n$ are iid uniform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def N(t_obs=0.001, lam=1.0e9, taud=100e-7, sigd=10e-7):\n",
    "    \"\"\"Simulates a Geiger counter using random increments\n",
    "    Dead time is simulated using a uniform random variable\n",
    "    \n",
    "    Inputs:\n",
    "    t_obs: observation time i.e. total time to simulate for\n",
    "    lam:   rate of incident Poisson process\n",
    "    taud:  mean dead time\n",
    "    sigd:  standard deviation of dead time distribution\n",
    "    Output:\n",
    "    Nt:    total detections at the end of observation period\"\"\"\n",
    "    Nt = 0\n",
    "    t_cur = 0.0\n",
    "    while t_cur<=t_obs:\n",
    "        t_cur += expovariate(lam) + uniform(taud-np.sqrt(3)*sigd,taud+np.sqrt(3)*sigd)\n",
    "        Nt += 1\n",
    "    return Nt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 2 involves explicitly modeling the Poisson process and then feeding it to the counter\n",
    "\n",
    "The Poisson process event times are given by\n",
    "$$\n",
    "V_m = V_{m-1} + X_m\n",
    "$$\n",
    "where $X_m$ are iid Exp$(\\lambda)$.\n",
    "The detection times are then modeled as:\n",
    "$$\n",
    "T_n = V_{\\min\\{i:V_i \\geq T_{n-1}+U_n\\}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Nalt(t_obs = 0.001, lam=1.0e9, taud = 100e-7, sigd=10e-7):\n",
    "    \"\"\"Simulates a Geiger counter by first explicitly modeling\n",
    "    the incident Poisson process with Exp interarrival times. Dead\n",
    "    time is a uniform random variable, unif. random is sampled every\n",
    "    particle.\n",
    "    \n",
    "    Inputs:\n",
    "    t_obs: observation time\n",
    "    lam:   rate of incident Poisson process\n",
    "    taud:  mean dead time\n",
    "    sigd:  standard deviation of dead time\n",
    "    Output:\n",
    "    particle_count: total detections at the end of t_obs\"\"\"\n",
    "    last_particle_time = 0\n",
    "    last_detection_time = -np.inf\n",
    "    particle_count = 0\n",
    "\n",
    "    while last_particle_time<=t_obs:\n",
    "        last_particle_time = last_particle_time + expovariate(lam)\n",
    "\n",
    "        fuzzy_dead_time = uniform(taud-np.sqrt(3)*sigd, taud+np.sqrt(3)*sigd)\n",
    "\n",
    "        if last_particle_time >= last_detection_time+fuzzy_dead_time:\n",
    "            last_detection_time = last_particle_time\n",
    "            particle_count += 1\n",
    "    \n",
    "    return particle_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Nalt_alt(t_obs = 0.001, lam=1.0e9, taud = 100e-7, sigd=10e-7):\n",
    "    \"\"\"Simulates a Geiger counter by first explicitly modeling\n",
    "    the incident Poisson process with Exp interarrival times. Dead\n",
    "    time is a uniform random variable \n",
    "    ****unif random is sampled once per\n",
    "    detection.***\n",
    "    \n",
    "    Inputs:\n",
    "    t_obs: observation time\n",
    "    lam:   rate of incident Poisson process\n",
    "    taud:  mean dead time\n",
    "    sigd:  standard deviation of dead time\n",
    "    Output:\n",
    "    particle_count: total detections at the end of t_obs\"\"\"\n",
    "    last_particle_time = 0\n",
    "    last_detection_time = -np.inf\n",
    "    particle_count = 0\n",
    "\n",
    "    fuzzy_dead_time = uniform(taud-np.sqrt(3)*sigd, taud+np.sqrt(3)*sigd)\n",
    "    \n",
    "    while last_particle_time<=t_obs:\n",
    "        last_particle_time = last_particle_time + expovariate(lam)\n",
    "\n",
    "        if last_particle_time >= last_detection_time+fuzzy_dead_time:\n",
    "            fuzzy_dead_time = uniform(taud-np.sqrt(3)*sigd, taud+np.sqrt(3)*sigd)\n",
    "            last_detection_time = last_particle_time\n",
    "            particle_count += 1\n",
    "    \n",
    "    return particle_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run 100 simulations of both methods for an observation duration of 0.001 with $\\lambda=1\\times 10^9$, and uniformly distributed dead times with mean $10^{-5}$ and standard deviation of $\\sqrt{3}\\times 10^{-6}.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method 1: Exponential increments\n",
      "mean count= 100.54 stdev= 0.984073168011\n",
      "Method 2: Explicit\n",
      "mean count= 120.0 stdev= 0.0\n",
      "Method 2': Explicit, unif only sampled once per detection\n",
      "mean count= 100.42 stdev= 1.02156742313\n"
     ]
    }
   ],
   "source": [
    "if __name__=='__main__':\n",
    "    Nsim = 100\n",
    "\n",
    "    # Method 1\n",
    "    Nt_nounc = np.zeros(Nsim)\n",
    "    for i in range(Nsim):\n",
    "        Nt_nounc[i] = N()\n",
    "    \n",
    "    # Method 2\n",
    "    Nalt_nounc = np.zeros(Nsim)\n",
    "    for i in range(Nsim):\n",
    "        Nalt_nounc[i] = Nalt()\n",
    "        \n",
    "    # Method 2'\n",
    "    Nalt_alt_nounc = np.zeros(Nsim)\n",
    "    for i in range(Nsim):\n",
    "        Nalt_alt_nounc[i] = Nalt_alt()\n",
    "        \n",
    "    print 'Method 1: Exponential increments'\n",
    "    print 'mean count=',np.mean(Nt_nounc), 'stdev=',np.std(Nt_nounc)\n",
    "    print 'Method 2: Explicit'\n",
    "    print 'mean count=',np.mean(Nalt_nounc), 'stdev=',np.std(Nalt_nounc)\n",
    "    print 'Method 2\\': Explicit, unif only sampled once per detection'\n",
    "    print 'mean count=',np.mean(Nalt_alt_nounc), 'stdev=',np.std(Nalt_alt_nounc)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFvBJREFUeJzt3Xuc1XW97/HXRzHIsqPI4DExGTuiAnI7I7hTS6T98FKp\nKIo7H+mUiT2U2pRpXjLJth0z8oIeKxUbdmHYQcXB487bQc1HW2CQUREx8dooImJeENhAfs8fs6QZ\nmAFmrZlZzJfX8/FYj/mt7/pdPus7ize/+a7f+q5IKSFJytcO5S5AktSxDHpJypxBL0mZM+glKXMG\nvSRlzqCXpMwZ9JKUOYNekjJn0EtS5rqVuwCAXr16pb59+5a7DEnqUubPn/9WSqliS+ttE0Hft29f\n6urqyl2GJHUpEfHK1qzn0I0kZc6gl6TMGfSSlLltYoy+JevWraOhoYE1a9aUuxR1sh49etCnTx92\n2mmncpciZWGbDfqGhgZ22WUX+vbtS0SUuxx1kpQSK1asoKGhgcrKynKXI2Vhi0M3EXFrRLwZEQub\ntPWMiAci4vnCz90K7RERkyNiSUQ8FRHDii1szZo17L777ob8diYi2H333f1LTmpHWzNGXwMcvVHb\nhcBDKaX9gIcK9wGOAfYr3MYBvyylOEN+++TvXWpfWwz6lNKjwNsbNR8PTC0sTwVOaNL+76nR48Cu\nEbFnexUrSWq7Ysfo90gpLS0svwHsUVjeC/hrk/UaCm1LKdE1D/yl1F00891/7teu+9tYfX09r7/+\nOsceeywAEydO5JOf/CTf//73i9pfa9s/+uijTJgwgaeeeorp06czZsyYFrePCL73ve/xi1/8AoBJ\nkyaxcuVKJk6cWFQ9krqOki+vTI3fLt7mbxiPiHERURcRdcuXLy+1jG1OfX099957b4cf5zOf+Qw1\nNTV89atf3ex63bt358477+Stt97q8JokNTH7fzXeyqjYoF/20ZBM4eebhfbXgL2brNen0LaJlNJN\nKaWqlFJVRcUWp2rodC+//DIHHHAA1dXV9OvXj9NOO40HH3yQQw89lP3224+5c+cC8MEHH/CNb3yD\n4cOHM3ToUO6++27Wrl3Lj370I26//XaGDBnC7bffDsCiRYs44ogj2HfffZk8efKGY1199dUMHDiQ\ngQMHcu21125ov+KKK+jXrx+HHXYYzz33XIt19u3bl0GDBrHDDpv/VXbr1o1x48ZxzTXXlNo1krqY\nYoO+FjijsHwGcHeT9tMLV98cArzbZIiny1myZAnnnXceixcvZvHixdx222089thjTJo0iZ/+9KdA\nYxgfeeSRzJ07l9mzZ3P++eezbt06Lr/8csaOHUt9fT1jx44FYPHixdx3333MnTuXH//4x6xbt475\n8+fzm9/8hjlz5vD4449z8803s2DBAubPn8/06dM3/GUwb968kp/Pueeey7Rp03j33XdL3pekrmOL\nY/QR8XvgCKBXRDQAlwFXAn+IiDOBV4BTCqvfCxwLLAFWAV/vgJo7TWVlJQcddBAAAwYMYNSoUUQE\nBx10EC+//DIA999/P7W1tUyaNAlovCz01VdfbXF/X/rSl+jevTvdu3end+/eLFu2jMcee4zRo0fz\niU98AoATTzyRP/3pT3z44YeMHj2anXfeGYDjjjuu5OfzqU99itNPP53Jkyfz8Y9/vOT9Seoathj0\nKaV/aeWhUS2sm4BzSy1qW9G9e/cNyzvssMOG+zvssAPr168HGj/gc8cdd7D//vs323bOnDmb3d+O\nO+64YR+dacKECQwbNoyvf71L/x8sqQ2c66ZERx11FNdffz2N/8fBggULANhll114//33t7j94Ycf\nzsyZM1m1ahUffPABd911F4cffjif//znmTlzJqtXr+b9999n1qxZ7VJvz549OeWUU5gyZUq77E/S\ntm+bnQJhYx19OWSxLr30UiZMmMCgQYP48MMPqays5J577mHkyJFceeWVDBkyhIsuuqjV7YcNG0Z1\ndTXDhw8H4Jvf/CZDhw4FYOzYsQwePJjevXtz8MEHt7j9vHnzGD16NH/729+YNWsWl112Gc8888xm\naz7vvPO44YYbinzGkrqa+OhMtJyqqqrSxl888uyzz3LggQeWqSKVm79/ZeOjSytHtn7CV6yImJ9S\nqtrSeg7dSFLmDHpJypxBL0mZM+glKXMGvSRlzqCXpMx1mevo2332tw641Kmpzpqm+Oqrr+aWW26h\nW7duVFRUcOutt7LPPvtssn0x0xTX1NRQV1fHDTfcwMyZM+nXrx/9+/cvqn5J5eMZfQfprGmKhw4d\nSl1dHU899RRjxozhggsuaHG9UqcpnjlzJosWLSqlVEllYtC3oqtMUzxy5MgNE58dcsghNDQ0tLje\n5qYpnjVrFiNGjGDo0KF88YtfZNmyZc0e//Of/0xtbS3nn38+Q4YM4YUXXmhDT0oqN4N+M7raNMVT\npkzhmGOOafXx1qYpPuyww3j88cdZsGABp556KldddVWzxz/3uc9x3HHH8fOf/5z6+no++9nPtrUr\nJZVR1xmjL4OuNE3x7373O+rq6njkkUdaXae1aYobGhoYO3YsS5cuZe3atVRWVm5dB0nqEjyj34y2\nTFNcX19PfX09r776aqtztHTUNMUPPvggV1xxBbW1tc2O0ZIJEyYwZcoUPvjggw1t3/72txk/fjxP\nP/00v/71r1mzZk271CVp22DQl6jc0xQvWLCAs88+m9raWnr37r3F47U0TfG7777LXnvtBcDUqVNb\n3G5rn4+kbU/XGbrp4Mshi1XuaYrPP/98Vq5cycknnww0fll4bW3tZmveeJriiRMncvLJJ7Pbbrtx\n5JFH8tJLL22yzamnnspZZ53F5MmTmTFjhuP0UhfiNMXaJvn7VzacpliS1NEMeknKnEEvSZkz6CUp\ncwa9JGXOoJekzHWZ6+hvrL+xXfd3zpBz2nV/krSt8oy+g2w8TfHEiRM3zIdTjNa2v/rqq+nfvz+D\nBg1i1KhRvPLKK0Dj7JtHHHHEJuu//PLLRATXX3/9hrbx48dTU1Oz1cevqanh9ddfL+p5VFdXM2PG\nDACuvfZaVq1aVdR+JG09g76DbGvz0TfVu3dvrrvuOtauXVvUMUsJ+qYMeqlzGPSt6Orz0e+44470\n7NmzxW0qKioYNWpUi/Pa3HzzzRx88MEMHjyYk046aZMgnjFjBnV1dZx22mkMGTKE1atXt3iMyy+/\nnIMPPpiBAwcybtw4Nv4E9uTJk3n99dcZOXIkI0eObHEfktqHQb8ZXXk++r333ps777yz1XV/8IMf\nMGnSJP7+9783az/xxBOZN28eTz75JAceeGCzyc8AxowZQ1VVFdOmTaO+vr7ZdMdNjR8/nnnz5rFw\n4UJWr17NPffc0+zx73znO3z6059m9uzZzJ49e4vPTVLxusybseWQ23z0Te27776MGDGC2267rVn7\nwoUL+eEPf8g777zDypUrOeqoo7ZqfxubPXs2V111FatWreLtt99mwIABfOUrXylqX5JKY9BvRlvm\no99///2bbTtnzpzN7q8j5qN/5JFHtjgffVMXX3wxY8aM4Qtf+MKGturqambOnMngwYOpqanh4Ycf\nbnM9a9as4ZxzzqGuro69996biRMnOse9VEZdJui31cshP5qP/vrrryciWLBgAUOHDm3TfPTV1dVc\neOGFpJS46667+O1vf0tKierqai666CLWr1/PrFmzOPvsszfZ/qP56P/4xz9u1Xz0TR1wwAH079+f\nWbNmbZgG+f3332fPPfdk3bp1TJs2bcM89U1t6bl9FOq9evVi5cqVzJgxgzFjxrS6n169erWpbklt\nU9IYfUR8NyKeiYiFEfH7iOgREZURMScilkTE7RHxsfYqdlt06aWXsm7dOgYNGsSAAQO49NJLgcY3\nSRctWtTszdiWNJ2PfsSIERvmox82bNiG+eiPOeaYrZqPfsiQIVsc4tnYJZdc0uwLxX/yk58wYsQI\nDj30UA444IAWt6muruZb3/pWq2/G7rrrrpx11lkMHDiQo446qtXax40bx9FHH+2bsVIHK3o++ojY\nC3gM6J9SWh0RfwDuBY4F7kwpTY+IXwFPppR+ubl9OR+9NubvX9nIYD76bsDHI6IbsDOwFDgSmFF4\nfCpwQonHkCSVoOgx+pTSaxExCXgVWA3cD8wH3kkpffQuYwOw6SCvsjF69OhNvnrwZz/7WdFX60hq\nf0UHfUTsBhwPVALvAP8HOLoN248DxkHj95y2JKVERBRbojrBXXfd1e773Ba+3lLKSSlDN18EXkop\nLU8prQPuBA4Fdi0M5QD0AV5raeOU0k0ppaqUUlVFRcUmj/fo0YMVK1b4j347k1JixYoV9OjRo9yl\nSNko5fLKV4FDImJnGoduRgF1wGxgDDAdOAO4u5id9+nTh4aGBpYvX15CieqKevToQZ8+fcpdhpSN\nUsbo50TEDOAJYD2wALgJ+L/A9Ij4t0LblNb30rqddtqJysrKYsuTJBWU9IGplNJlwGUbNb8IDC9l\nv5Kk9uOkZpKUOYNekjJn0EtS5gx6ScqcQS9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYM\neknKnEEvSZkz6CUpcwa9JGXOoJekzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCX\npMwZ9JKUOYNekjJn0EtS5gx6ScqcQS9JmTPoJSlzBr0kZc6gl6TMGfSSlLmSgj4ido2IGRGxOCKe\njYh/ioieEfFARDxf+LlbexUrSWq7Us/orwP+mFI6ABgMPAtcCDyUUtoPeKhwX5JUJkUHfUT8N+Dz\nwBSAlNLalNI7wPHA1MJqU4ETSi1SklS8Us7oK4HlwG8iYkFE3BIRnwD2SCktLazzBrBHqUVKkopX\nStB3A4YBv0wpDQU+YKNhmpRSAlJLG0fEuIioi4i65cuXl1CGJGlzSgn6BqAhpTSncH8GjcG/LCL2\nBCj8fLOljVNKN6WUqlJKVRUVFSWUIUnanKKDPqX0BvDXiNi/0DQKWATUAmcU2s4A7i6pQklSSbqV\nuP23gWkR8THgReDrNP7n8YeIOBN4BTilxGNIkkpQUtCnlOqBqhYeGlXKfiVJ7cdPxkpS5gx6Scqc\nQS9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYMeknKnEEvSZkz6CUpcwa9JGXOoJekzBn0\nkpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUOYNekjJn0EtS5gx6ScqcQS9J\nmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIy163UHUTEjkAd8FpK6csRUQlMB3YH5gNfSymtLfU424Ib\n628sartzhpzTzpVI0tZrjzP6fwWebXL/Z8A1KaX/AfwNOLMdjiFJKlJJQR8RfYAvAbcU7gdwJDCj\nsMpU4IRSjiFJKk2pZ/TXAhcAHxbu7w68k1JaX7jfAOxV4jEkSSUoOugj4svAmyml+UVuPy4i6iKi\nbvny5cWWIUnaglLO6A8FjouIl2l88/VI4Dpg14j46E3ePsBrLW2cUroppVSVUqqqqKgooQxJ0uYU\nHfQppYtSSn1SSn2BU4H/l1I6DZgNjCmsdgZwd8lVSpKK1hHX0f8A+F5ELKFxzH5KBxxDkrSVSr6O\nHiCl9DDwcGH5RWB4e+xXklQ6PxkrSZkz6CUpcwa9JGWuXcbotxf/+cKKorY7Z0g7FyJJbeAZvSRl\nzqCXpMwZ9JKUOcfoO8E1D/ylzdt895/7dUAlkrZHntFLUuYMeknKnEEvSZkz6CUpcwa9JGXOoJek\nzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUOYNekjJn0EtS5gx6Scqc\nQS9JmTPoJSlzfjl4J3jivduL2OrSdq9D0vbJM3pJypxBL0mZM+glKXMGvSRlruigj4i9I2J2RCyK\niGci4l8L7T0j4oGIeL7wc7f2K1eS1FalnNGvB85LKfUHDgHOjYj+wIXAQyml/YCHCvclSWVSdNCn\nlJamlJ4oLL8PPAvsBRwPTC2sNhU4odQiJUnFa5cx+ojoCwwF5gB7pJSWFh56A9ijlW3GRURdRNQt\nX768PcqQJLWg5KCPiE8CdwATUkrvNX0spZSA1NJ2KaWbUkpVKaWqioqKUsuQJLWipKCPiJ1oDPlp\nKaU7C83LImLPwuN7Am+WVqIkqRSlXHUTwBTg2ZTS1U0eqgXOKCyfAdxdfHmSpFKVMtfNocDXgKcj\nor7QdjFwJfCHiDgTeAU4pbQSJUmlKDroU0qPAdHKw6OK3a8kqX35yVhJypzTFG+jvnbHT4ra7rcn\nOb2xpOY8o5ekzG23Z/Q31t9Y7hIkqVN4Ri9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYM\neknKnEEvSZnbbj8Z+58vrCh3CZLUKTyjl6TMGfSSlDmDXpIyZ9BLUuYMeknKnEEvSZkz6CUpcwa9\nJGXOoJekzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUuQ75KsGIOBq4\nDtgRuCWldGVHHAfgmgf+0lG7lqQstPsZfUTsCPxv4BigP/AvEdG/vY8jSdo6HXFGPxxYklJ6ESAi\npgPHA4s64Fg88d7tHbFbScpGR4zR7wX8tcn9hkKbJKkMOmSMfmtExDhgXOHuyoh4rly1tEEv4K1y\nF7E5v+NHnXm4bb4/Opn98Q/2RXO94OKO6I99tmaljgj614C9m9zvU2hrJqV0E3BTBxy/w0REXUqp\nqtx1bCvsj+bsj3+wL5ord390xNDNPGC/iKiMiI8BpwK1HXAcSdJWaPcz+pTS+ogYD9xH4+WVt6aU\nnmnv40iStk6HjNGnlO4F7u2IfZdZlxpq6gT2R3P2xz/YF82VtT8ipVTO40uSOphTIEhS5gz6JiLi\n1oh4MyIWNmnrGREPRMTzhZ+7FdojIiZHxJKIeCoihpWv8o7Rxv44IiLejYj6wq1Tr/PsaK30xckR\n8UxEfBgRVRutf1HhtfFcRBzV+RV3rLb0R0T0jYjVTV4bvypP1R2nlf74eUQsLuTDXRGxa5PHOvX1\nYdA3VwMcvVHbhcBDKaX9gIcK96Fxiof9CrdxwC87qcbOVMPW9wfAn1JKQwq3yzupxs5Sw6Z9sRA4\nEXi0aWNhyo9TgQGFbW4sTA2Skxq2sj8KXmjy2vhWRxdXBjVs2h8PAANTSoOAvwAXQXleHwZ9Eyml\nR4G3N2o+HphaWJ4KnNCk/d9To8eBXSNiz86ptHO0sT+y1lJfpJSeTSm19EG/44HpKaX/Sim9BCyh\ncWqQbLSxP7LXSn/cn1JaX7j7OI2fKYIyvD4M+i3bI6W0tLD8BrBHYXl7neqhtf4A+KeIeDIi/iMi\nBpShtm3F9vra2JzKiFgQEY9ExOHlLqYMvgH8R2G5018fZZsCoStKKaWI8DKlgo364wlgn5TSyog4\nFphJ47CWtBT4TEppRUT8T2BmRAxIKb1X7sI6Q0RcAqwHppWrBs/ot2zZR0MyhZ9vFtq3aqqHDLXY\nHyml91JKKwvL9wI7RUSv8pVZVtvra6NFhSGKFYXl+cALQL/yVtU5IqIa+DJwWvrHteyd/vow6Les\nFjijsHwGcHeT9tMLV98cArzbZEgjZy32R0T894iIwvJwGl9bK8pSYfnVAqdGRPeIqKTxL5u5Za6p\nbCKi4qM3GyNiXxr748XyVtXxCl/AdAFwXEppVZOHOv/1kVLyVrgBv6fxz8x1NI6bnQnsTuPVJc8D\nDwI9C+sGjV+w8gLwNFBV7vrL3B/jgWeAJ2l84+lz5a6/E/pidGH5v4BlwH1N1r+k8Np4Djim3PWX\nsz+AkwqvjXoah/i+Uu76O6k/ltA4Fl9fuP2qXK8PPxkrSZlz6EaSMmfQS1LmDHpJypxBL0mZM+gl\nKXMGvSRlzqCXpMwZ9JKUuf8PbDUJf+4w/dEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7eed227550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_=plt.hist(Nt_nounc,bins=5,alpha=0.5)\n",
    "_=plt.hist(Nalt_nounc, bins=5, alpha=0.5)\n",
    "_=plt.hist(Nalt_alt_nounc, bins=5, alpha=0.5)\n",
    "_=plt.legend(('method 1 N', 'method 2 Nalt', 'method 2\\' Nalt_alt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note how the two methods produce wildly different counts!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 2,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"from random import expovariate,uniform\n",
	"import numpy as np\n",
	"import matplotlib.pyplot as plt\n",
	"%matplotlib inline"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"## Method 1: using exponential increments with random dead times\n",
	"The model is\n",
	"$$\n",
	"T_n - T_{n-1} = X_n + U_n \n",
	"$$\n",
	"where $T_n$ is the $n$th detection event, $X_n$ are iid Exp$(\\lambda)$ and $U_n$ are iid uniform."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 3,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"def N(t_obs=0.001, lam=1.0e9, taud=100e-7, sigd=10e-7):\n",
	" \"\"\"Simulates a Geiger counter using random increments\n",
	" Dead time is simulated using a uniform random variable\n",
	" \n",
	" Inputs:\n",
	" t_obs: observation time i.e. total time to simulate for\n",
	" lam: rate of incident Poisson process\n",
	" taud: mean dead time\n",
	" sigd: standard deviation of dead time distribution\n",
	" Output:\n",
	" Nt: total detections at the end of observation period\"\"\"\n",
	" Nt = 0\n",
	" t_cur = 0.0\n",
	" while t_cur<=t_obs:\n",
	" t_cur += expovariate(lam) + uniform(taud-np.sqrt(3)sigd,taud+np.sqrt(3)sigd)\n",
	" Nt += 1\n",
	" return Nt"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"## Method 2 involves explicitly modeling the Poisson process and then feeding it to the counter\n",
	"\n",
	"The Poisson process event times are given by\n",
	"$$\n",
	"V_m = V_{m-1} + X_m\n",
	"$$\n",
	"where $X_m$ are iid Exp$(\\lambda)$.\n",
	"The detection times are then modeled as:\n",
	"$$\n",
	"T_n = V_{\\min\\{i:V_i \\geq T_{n-1}+U_n\\}}\n",
	"$$"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"def Nalt(t_obs = 0.001, lam=1.0e9, taud = 100e-7, sigd=10e-7):\n",
	" \"\"\"Simulates a Geiger counter by first explicitly modeling\n",
	" the incident Poisson process with Exp interarrival times. Dead\n",
	" time is a uniform random variable, unif. random is sampled every\n",
	" particle.\n",
	" \n",
	" Inputs:\n",
	" t_obs: observation time\n",
	" lam: rate of incident Poisson process\n",
	" taud: mean dead time\n",
	" sigd: standard deviation of dead time\n",
	" Output:\n",
	" particle_count: total detections at the end of t_obs\"\"\"\n",
	" last_particle_time = 0\n",
	" last_detection_time = -np.inf\n",
	" particle_count = 0\n",
	"\n",
	" while last_particle_time<=t_obs:\n",
	" last_particle_time = last_particle_time + expovariate(lam)\n",
	"\n",
	" fuzzy_dead_time = uniform(taud-np.sqrt(3)sigd, taud+np.sqrt(3)sigd)\n",
	"\n",
	" if last_particle_time >= last_detection_time+fuzzy_dead_time:\n",
	" last_detection_time = last_particle_time\n",
	" particle_count += 1\n",
	" \n",
	" return particle_count"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"def Nalt_alt(t_obs = 0.001, lam=1.0e9, taud = 100e-7, sigd=10e-7):\n",
	" \"\"\"Simulates a Geiger counter by first explicitly modeling\n",
	" the incident Poisson process with Exp interarrival times. Dead\n",
	" time is a uniform random variable \n",
	" ****unif random is sampled once per\n",
	" detection.***\n",
	" \n",
	" Inputs:\n",
	" t_obs: observation time\n",
	" lam: rate of incident Poisson process\n",
	" taud: mean dead time\n",
	" sigd: standard deviation of dead time\n",
	" Output:\n",
	" particle_count: total detections at the end of t_obs\"\"\"\n",
	" last_particle_time = 0\n",
	" last_detection_time = -np.inf\n",
	" particle_count = 0\n",
	"\n",
	" fuzzy_dead_time = uniform(taud-np.sqrt(3)sigd, taud+np.sqrt(3)sigd)\n",
	" \n",
	" while last_particle_time<=t_obs:\n",
	" last_particle_time = last_particle_time + expovariate(lam)\n",
	"\n",
	" if last_particle_time >= last_detection_time+fuzzy_dead_time:\n",
	" fuzzy_dead_time = uniform(taud-np.sqrt(3)sigd, taud+np.sqrt(3)sigd)\n",
	" last_detection_time = last_particle_time\n",
	" particle_count += 1\n",
	" \n",
	" return particle_count"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We run 100 simulations of both methods for an observation duration of 0.001 with $\\lambda=1\\times 10^9$, and uniformly distributed dead times with mean $10^{-5}$ and standard deviation of $\\sqrt{3}\\times 10^{-6}.$"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 5,
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"Method 1: Exponential increments\n",
	"mean count= 100.54 stdev= 0.984073168011\n",
	"Method 2: Explicit\n",
	"mean count= 120.0 stdev= 0.0\n",
	"Method 2': Explicit, unif only sampled once per detection\n",
	"mean count= 100.42 stdev= 1.02156742313\n"
	]
	}
	],
	"source": [
	"if __name__=='__main__':\n",
	" Nsim = 100\n",
	"\n",
	" # Method 1\n",
	" Nt_nounc = np.zeros(Nsim)\n",
	" for i in range(Nsim):\n",
	" Nt_nounc[i] = N()\n",
	" \n",
	" # Method 2\n",
	" Nalt_nounc = np.zeros(Nsim)\n",
	" for i in range(Nsim):\n",
	" Nalt_nounc[i] = Nalt()\n",
	" \n",
	" # Method 2'\n",
	" Nalt_alt_nounc = np.zeros(Nsim)\n",
	" for i in range(Nsim):\n",
	" Nalt_alt_nounc[i] = Nalt_alt()\n",
	" \n",
	" print 'Method 1: Exponential increments'\n",
	" print 'mean count=',np.mean(Nt_nounc), 'stdev=',np.std(Nt_nounc)\n",
	" print 'Method 2: Explicit'\n",
	" print 'mean count=',np.mean(Nalt_nounc), 'stdev=',np.std(Nalt_nounc)\n",
	" print 'Method 2\\': Explicit, unif only sampled once per detection'\n",
	" print 'mean count=',np.mean(Nalt_alt_nounc), 'stdev=',np.std(Nalt_alt_nounc)\n",
	" "
	]
	},
	{
	"cell_type": "code",
	"execution_count": 9,
	"metadata": {
	"scrolled": true
	},
	"outputs": [
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFvBJREFUeJzt3Xuc1XW97/HXRzHIsqPI4DExGTuiAnI7I7hTS6T98FKp\nKIo7H+mUiT2U2pRpXjLJth0z8oIeKxUbdmHYQcXB487bQc1HW2CQUREx8dooImJeENhAfs8fs6QZ\nmAFmrZlZzJfX8/FYj/mt7/pdPus7ize/+a7f+q5IKSFJytcO5S5AktSxDHpJypxBL0mZM+glKXMG\nvSRlzqCXpMwZ9JKUOYNekjJn0EtS5rqVuwCAXr16pb59+5a7DEnqUubPn/9WSqliS+ttE0Hft29f\n6urqyl2GJHUpEfHK1qzn0I0kZc6gl6TMGfSSlLltYoy+JevWraOhoYE1a9aUuxR1sh49etCnTx92\n2mmncpciZWGbDfqGhgZ22WUX+vbtS0SUuxx1kpQSK1asoKGhgcrKynKXI2Vhi0M3EXFrRLwZEQub\ntPWMiAci4vnCz90K7RERkyNiSUQ8FRHDii1szZo17L777ob8diYi2H333f1LTmpHWzNGXwMcvVHb\nhcBDKaX9gIcK9wGOAfYr3MYBvyylOEN+++TvXWpfWwz6lNKjwNsbNR8PTC0sTwVOaNL+76nR48Cu\nEbFnexUrSWq7Ysfo90gpLS0svwHsUVjeC/hrk/UaCm1LKdE1D/yl1F00891/7teu+9tYfX09r7/+\nOsceeywAEydO5JOf/CTf//73i9pfa9s/+uijTJgwgaeeeorp06czZsyYFrePCL73ve/xi1/8AoBJ\nkyaxcuVKJk6cWFQ9krqOki+vTI3fLt7mbxiPiHERURcRdcuXLy+1jG1OfX099957b4cf5zOf+Qw1\nNTV89atf3ex63bt358477+Stt97q8JokNTH7fzXeyqjYoF/20ZBM4eebhfbXgL2brNen0LaJlNJN\nKaWqlFJVRcUWp2rodC+//DIHHHAA1dXV9OvXj9NOO40HH3yQQw89lP3224+5c+cC8MEHH/CNb3yD\n4cOHM3ToUO6++27Wrl3Lj370I26//XaGDBnC7bffDsCiRYs44ogj2HfffZk8efKGY1199dUMHDiQ\ngQMHcu21125ov+KKK+jXrx+HHXYYzz33XIt19u3bl0GDBrHDDpv/VXbr1o1x48ZxzTXXlNo1krqY\nYoO+FjijsHwGcHeT9tMLV98cArzbZIiny1myZAnnnXceixcvZvHixdx222089thjTJo0iZ/+9KdA\nYxgfeeSRzJ07l9mzZ3P++eezbt06Lr/8csaOHUt9fT1jx44FYPHixdx3333MnTuXH//4x6xbt475\n8+fzm9/8hjlz5vD4449z8803s2DBAubPn8/06dM3/GUwb968kp/Pueeey7Rp03j33XdL3pekrmOL\nY/QR8XvgCKBXRDQAlwFXAn+IiDOBV4BTCqvfCxwLLAFWAV/vgJo7TWVlJQcddBAAAwYMYNSoUUQE\nBx10EC+//DIA999/P7W1tUyaNAlovCz01VdfbXF/X/rSl+jevTvdu3end+/eLFu2jMcee4zRo0fz\niU98AoATTzyRP/3pT3z44YeMHj2anXfeGYDjjjuu5OfzqU99itNPP53Jkyfz8Y9/vOT9Seoathj0\nKaV/aeWhUS2sm4BzSy1qW9G9e/cNyzvssMOG+zvssAPr168HGj/gc8cdd7D//vs323bOnDmb3d+O\nO+64YR+dacKECQwbNoyvf71L/x8sqQ2c66ZERx11FNdffz2N/8fBggULANhll114//33t7j94Ycf\nzsyZM1m1ahUffPABd911F4cffjif//znmTlzJqtXr+b9999n1qxZ7VJvz549OeWUU5gyZUq77E/S\ntm+bnQJhYx19OWSxLr30UiZMmMCgQYP48MMPqays5J577mHkyJFceeWVDBkyhIsuuqjV7YcNG0Z1\ndTXDhw8H4Jvf/CZDhw4FYOzYsQwePJjevXtz8MEHt7j9vHnzGD16NH/729+YNWsWl112Gc8888xm\naz7vvPO44YYbinzGkrqa+OhMtJyqqqrSxl888uyzz3LggQeWqSKVm79/ZeOjSytHtn7CV6yImJ9S\nqtrSeg7dSFLmDHpJypxBL0mZM+glKXMGvSRlzqCXpMx1mevo2332tw641Kmpzpqm+Oqrr+aWW26h\nW7duVFRUcOutt7LPPvtssn0x0xTX1NRQV1fHDTfcwMyZM+nXrx/9+/cvqn5J5eMZfQfprGmKhw4d\nSl1dHU899RRjxozhggsuaHG9UqcpnjlzJosWLSqlVEllYtC3oqtMUzxy5MgNE58dcsghNDQ0tLje\n5qYpnjVrFiNGjGDo0KF88YtfZNmyZc0e//Of/0xtbS3nn38+Q4YM4YUXXmhDT0oqN4N+M7raNMVT\npkzhmGOOafXx1qYpPuyww3j88cdZsGABp556KldddVWzxz/3uc9x3HHH8fOf/5z6+no++9nPtrUr\nJZVR1xmjL4OuNE3x7373O+rq6njkkUdaXae1aYobGhoYO3YsS5cuZe3atVRWVm5dB0nqEjyj34y2\nTFNcX19PfX09r776aqtztHTUNMUPPvggV1xxBbW1tc2O0ZIJEyYwZcoUPvjggw1t3/72txk/fjxP\nP/00v/71r1mzZk271CVp22DQl6jc0xQvWLCAs88+m9raWnr37r3F47U0TfG7777LXnvtBcDUqVNb\n3G5rn4+kbU/XGbrp4Mshi1XuaYrPP/98Vq5cycknnww0fll4bW3tZmveeJriiRMncvLJJ7Pbbrtx\n5JFH8tJLL22yzamnnspZZ53F5MmTmTFjhuP0UhfiNMXaJvn7VzacpliS1NEMeknKnEEvSZkz6CUp\ncwa9JGXOoJekzHWZ6+hvrL+xXfd3zpBz2nV/krSt8oy+g2w8TfHEiRM3zIdTjNa2v/rqq+nfvz+D\nBg1i1KhRvPLKK0Dj7JtHHHHEJuu//PLLRATXX3/9hrbx48dTU1Oz1cevqanh9ddfL+p5VFdXM2PG\nDACuvfZaVq1aVdR+JG09g76DbGvz0TfVu3dvrrvuOtauXVvUMUsJ+qYMeqlzGPSt6Orz0e+44470\n7NmzxW0qKioYNWpUi/Pa3HzzzRx88MEMHjyYk046aZMgnjFjBnV1dZx22mkMGTKE1atXt3iMyy+/\nnIMPPpiBAwcybtw4Nv4E9uTJk3n99dcZOXIkI0eObHEfktqHQb8ZXXk++r333ps777yz1XV/8IMf\nMGnSJP7+9783az/xxBOZN28eTz75JAceeGCzyc8AxowZQ1VVFdOmTaO+vr7ZdMdNjR8/nnnz5rFw\n4UJWr17NPffc0+zx73znO3z6059m9uzZzJ49e4vPTVLxusybseWQ23z0Te27776MGDGC2267rVn7\nwoUL+eEPf8g777zDypUrOeqoo7ZqfxubPXs2V111FatWreLtt99mwIABfOUrXylqX5JKY9BvRlvm\no99///2bbTtnzpzN7q8j5qN/5JFHtjgffVMXX3wxY8aM4Qtf+MKGturqambOnMngwYOpqanh4Ycf\nbnM9a9as4ZxzzqGuro69996biRMnOse9VEZdJui31cshP5qP/vrrryciWLBgAUOHDm3TfPTV1dVc\neOGFpJS46667+O1vf0tKierqai666CLWr1/PrFmzOPvsszfZ/qP56P/4xz9u1Xz0TR1wwAH079+f\nWbNmbZgG+f3332fPPfdk3bp1TJs2bcM89U1t6bl9FOq9evVi5cqVzJgxgzFjxrS6n169erWpbklt\nU9IYfUR8NyKeiYiFEfH7iOgREZURMScilkTE7RHxsfYqdlt06aWXsm7dOgYNGsSAAQO49NJLgcY3\nSRctWtTszdiWNJ2PfsSIERvmox82bNiG+eiPOeaYrZqPfsiQIVsc4tnYJZdc0uwLxX/yk58wYsQI\nDj30UA444IAWt6muruZb3/pWq2/G7rrrrpx11lkMHDiQo446qtXax40bx9FHH+2bsVIHK3o++ojY\nC3gM6J9SWh0RfwDuBY4F7kwpTY+IXwFPppR+ubl9OR+9NubvX9nIYD76bsDHI6IbsDOwFDgSmFF4\nfCpwQonHkCSVoOgx+pTSaxExCXgVWA3cD8wH3kkpffQuYwOw6SCvsjF69OhNvnrwZz/7WdFX60hq\nf0UHfUTsBhwPVALvAP8HOLoN248DxkHj95y2JKVERBRbojrBXXfd1e773Ba+3lLKSSlDN18EXkop\nLU8prQPuBA4Fdi0M5QD0AV5raeOU0k0ppaqUUlVFRcUmj/fo0YMVK1b4j347k1JixYoV9OjRo9yl\nSNko5fLKV4FDImJnGoduRgF1wGxgDDAdOAO4u5id9+nTh4aGBpYvX15CieqKevToQZ8+fcpdhpSN\nUsbo50TEDOAJYD2wALgJ+L/A9Ij4t0LblNb30rqddtqJysrKYsuTJBWU9IGplNJlwGUbNb8IDC9l\nv5Kk9uOkZpKUOYNekjJn0EtS5gx6ScqcQS9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYM\neknKnEEvSZkz6CUpcwa9JGXOoJekzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCX\npMwZ9JKUOYNekjJn0EtS5gx6ScqcQS9JmTPoJSlzBr0kZc6gl6TMGfSSlLmSgj4ido2IGRGxOCKe\njYh/ioieEfFARDxf+LlbexUrSWq7Us/orwP+mFI6ABgMPAtcCDyUUtoPeKhwX5JUJkUHfUT8N+Dz\nwBSAlNLalNI7wPHA1MJqU4ETSi1SklS8Us7oK4HlwG8iYkFE3BIRnwD2SCktLazzBrBHqUVKkopX\nStB3A4YBv0wpDQU+YKNhmpRSAlJLG0fEuIioi4i65cuXl1CGJGlzSgn6BqAhpTSncH8GjcG/LCL2\nBCj8fLOljVNKN6WUqlJKVRUVFSWUIUnanKKDPqX0BvDXiNi/0DQKWATUAmcU2s4A7i6pQklSSbqV\nuP23gWkR8THgReDrNP7n8YeIOBN4BTilxGNIkkpQUtCnlOqBqhYeGlXKfiVJ7cdPxkpS5gx6Scqc\nQS9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYMeknKnEEvSZkz6CUpcwa9JGXOoJekzBn0\nkpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUOYNekjJn0EtS5gx6ScqcQS9J\nmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIy163UHUTEjkAd8FpK6csRUQlMB3YH5gNfSymtLfU424Ib\n628sartzhpzTzpVI0tZrjzP6fwWebXL/Z8A1KaX/AfwNOLMdjiFJKlJJQR8RfYAvAbcU7gdwJDCj\nsMpU4IRSjiFJKk2pZ/TXAhcAHxbu7w68k1JaX7jfAOxV4jEkSSUoOugj4svAmyml+UVuPy4i6iKi\nbvny5cWWIUnaglLO6A8FjouIl2l88/VI4Dpg14j46E3ePsBrLW2cUroppVSVUqqqqKgooQxJ0uYU\nHfQppYtSSn1SSn2BU4H/l1I6DZgNjCmsdgZwd8lVSpKK1hHX0f8A+F5ELKFxzH5KBxxDkrSVSr6O\nHiCl9DDwcGH5RWB4e+xXklQ6PxkrSZkz6CUpcwa9JGWuXcbotxf/+cKKorY7Z0g7FyJJbeAZvSRl\nzqCXpMwZ9JKUOcfoO8E1D/ylzdt895/7dUAlkrZHntFLUuYMeknKnEEvSZkz6CUpcwa9JGXOoJek\nzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUOYNekjJn0EtS5gx6Scqc\nQS9JmTPoJSlzfjl4J3jivduL2OrSdq9D0vbJM3pJypxBL0mZM+glKXMGvSRlruigj4i9I2J2RCyK\niGci4l8L7T0j4oGIeL7wc7f2K1eS1FalnNGvB85LKfUHDgHOjYj+wIXAQyml/YCHCvclSWVSdNCn\nlJamlJ4oLL8PPAvsBRwPTC2sNhU4odQiJUnFa5cx+ojoCwwF5gB7pJSWFh56A9ijlW3GRURdRNQt\nX768PcqQJLWg5KCPiE8CdwATUkrvNX0spZSA1NJ2KaWbUkpVKaWqioqKUsuQJLWipKCPiJ1oDPlp\nKaU7C83LImLPwuN7Am+WVqIkqRSlXHUTwBTg2ZTS1U0eqgXOKCyfAdxdfHmSpFKVMtfNocDXgKcj\nor7QdjFwJfCHiDgTeAU4pbQSJUmlKDroU0qPAdHKw6OK3a8kqX35yVhJypzTFG+jvnbHT4ra7rcn\nOb2xpOY8o5ekzG23Z/Q31t9Y7hIkqVN4Ri9JmTPoJSlzBr0kZc6gl6TMGfSSlDmDXpIyZ9BLUuYM\neknKnEEvSZnbbj8Z+58vrCh3CZLUKTyjl6TMGfSSlDmDXpIyZ9BLUuYMeknKnEEvSZkz6CUpcwa9\nJGXOoJekzBn0kpQ5g16SMmfQS1LmDHpJypxBL0mZM+glKXMGvSRlzqCXpMwZ9JKUuQ75KsGIOBq4\nDtgRuCWldGVHHAfgmgf+0lG7lqQstPsZfUTsCPxv4BigP/AvEdG/vY8jSdo6HXFGPxxYklJ6ESAi\npgPHA4s64Fg88d7tHbFbScpGR4zR7wX8tcn9hkKbJKkMOmSMfmtExDhgXOHuyoh4rly1tEEv4K1y\nF7E5v+NHnXm4bb4/Opn98Q/2RXO94OKO6I99tmaljgj614C9m9zvU2hrJqV0E3BTBxy/w0REXUqp\nqtx1bCvsj+bsj3+wL5ord390xNDNPGC/iKiMiI8BpwK1HXAcSdJWaPcz+pTS+ogYD9xH4+WVt6aU\nnmnv40iStk6HjNGnlO4F7u2IfZdZlxpq6gT2R3P2xz/YF82VtT8ipVTO40uSOphTIEhS5gz6JiLi\n1oh4MyIWNmnrGREPRMTzhZ+7FdojIiZHxJKIeCoihpWv8o7Rxv44IiLejYj6wq1Tr/PsaK30xckR\n8UxEfBgRVRutf1HhtfFcRBzV+RV3rLb0R0T0jYjVTV4bvypP1R2nlf74eUQsLuTDXRGxa5PHOvX1\nYdA3VwMcvVHbhcBDKaX9gIcK96Fxiof9CrdxwC87qcbOVMPW9wfAn1JKQwq3yzupxs5Sw6Z9sRA4\nEXi0aWNhyo9TgQGFbW4sTA2Skxq2sj8KXmjy2vhWRxdXBjVs2h8PAANTSoOAvwAXQXleHwZ9Eyml\nR4G3N2o+HphaWJ4KnNCk/d9To8eBXSNiz86ptHO0sT+y1lJfpJSeTSm19EG/44HpKaX/Sim9BCyh\ncWqQbLSxP7LXSn/cn1JaX7j7OI2fKYIyvD4M+i3bI6W0tLD8BrBHYXl7neqhtf4A+KeIeDIi/iMi\nBpShtm3F9vra2JzKiFgQEY9ExOHlLqYMvgH8R2G5018fZZsCoStKKaWI8DKlgo364wlgn5TSyog4\nFphJ47CWtBT4TEppRUT8T2BmRAxIKb1X7sI6Q0RcAqwHppWrBs/ot2zZR0MyhZ9vFtq3aqqHDLXY\nHyml91JKKwvL9wI7RUSv8pVZVtvra6NFhSGKFYXl+cALQL/yVtU5IqIa+DJwWvrHteyd/vow6Les\nFjijsHwGcHeT9tMLV98cArzbZEgjZy32R0T894iIwvJwGl9bK8pSYfnVAqdGRPeIqKTxL5u5Za6p\nbCKi4qM3GyNiXxr748XyVtXxCl/AdAFwXEppVZOHOv/1kVLyVrgBv6fxz8x1NI6bnQnsTuPVJc8D\nDwI9C+sGjV+w8gLwNFBV7vrL3B/jgWeAJ2l84+lz5a6/E/pidGH5v4BlwH1N1r+k8Np4Djim3PWX\nsz+AkwqvjXoah/i+Uu76O6k/ltA4Fl9fuP2qXK8PPxkrSZlz6EaSMmfQS1LmDHpJypxBL0mZM+gl\nKXMGvSRlzqCXpMwZ9JKUuf8PbDUJf+4w/dEAAAAASUVORK5CYII=\n",
	"text/plain": [
	"<matplotlib.figure.Figure at 0x7f7eed227550>"
	]
	},
	"metadata": {},
	"output_type": "display_data"
	}
	],
	"source": [
	"_=plt.hist(Nt_nounc,bins=5,alpha=0.5)\n",
	"_=plt.hist(Nalt_nounc, bins=5, alpha=0.5)\n",
	"_=plt.hist(Nalt_alt_nounc, bins=5, alpha=0.5)\n",
	"_=plt.legend(('method 1 N', 'method 2 Nalt', 'method 2\\' Nalt_alt'))"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"## Note how the two methods produce wildly different counts!"
	]
	}
	],
	"metadata": {
	"kernelspec": {
	"display_name": "Python 2",
	"language": "python",
	"name": "python2"
	},
	"language_info": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 2
	},
	"file_extension": ".py",
	"mimetype": "text/x-python",
	"name": "python",
	"nbconvert_exporter": "python",
	"pygments_lexer": "ipython2",
	"version": "2.7.12"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}