joezuntz/gist:4c791f9592d1798577da

## gistfile1.txt
{
 "metadata": {
  "kernelspec": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "display_name": "IPython (Python 2)",
   "language": "python",
   "name": "python2"
  },
  "name": "",
  "signature": "sha256:55c4f9d7a7aae3543fe006737001ce10ac9cc48f7a29a15371c3b9c8fa640653"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Lecture 1 Medium Problem Worked Solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is an ipython notebook, an interactive python session saved as a web page.  You can download this and run it on your own machine - google \"ipython notebook\" for details.\n",
      "\n",
      "These special commands set up the notebook environment, and import numpy, matplotlib, and the factorial function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from math import factorial"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The general approach when you have a probability problem like this one is to write down all the thing you do know as probability distributions, and then use the laws of probability to convert them to what you want. \n",
      "\n",
      "In this case we have some quantities specified as conditional probabilities, like $P(n)$ and the $P(e_i)$ for a single photon, and we want $P(e)$, the overall distribution for the sum of all the energies.\n",
      "\n",
      "First define the constants of the problem (the word \"lambda\" is a python keyword so we can't use that - we call it \"lam\" instead)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mu = 1.0     #The mean of the Gaussian for each photon's energy, in keV\n",
      "sigma = 0.1  #The standard deviation of the Gaussian, in keV\n",
      "lam = 1.0    #The mean arrival rate in the Poisson distribution, in s^-1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first probability we have is the easiest one, the probability of getting a given number of photons. This is a Poisson distribution, which is discrete since you can't have half a photon.\n",
      "\n",
      "A Poisson distribution is:\n",
      "\n",
      "$\n",
      "P(n) = \\frac{\\lambda^n \\exp(-\\lambda)} {n!}\n",
      "$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def P_n(n):\n",
      "    return lam**n * exp(-lam) / factorial(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are told that the probability of the energy of each photon is Gaussian distributed.  We need the distribution of the total energy of the n photons that we are assuming are arriving in a second.  \n",
      "\n",
      "There's a handy formula for the distribution of the sum of Gaussians, which is just that if $X \\sim N(\\mu_1, \\sigma_1^2)$ and $Y \\sim N(\\mu_2, \\sigma_2^2)$, then $X+Y \\sim N(\\mu_1+\\mu_2, \\sigma_1^2+\\sigma_2^2)$.  \n",
      "\n",
      "We can generalize this to give us the distribution of the sum of n identical distributions as $N(n\\mu, n \\sigma^2)$. So:\n",
      "\n",
      "$\n",
      "P(e|n) = (2 \\pi n \\sigma^2)^{-\\frac{1}{2}}\\exp{[-(e-n \\mu)^2 /  2 n \\sigma^2]}\n",
      "$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def P_e_given_n(e,n):\n",
      "    return exp(-0.5*(e-n*mu)**2/(n*sigma**2)) / sqrt(2*pi*n*sigma**2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can use the marginalization rule, which says we can sum over all possible values of n to get the total energy.\n",
      "\n",
      "We need to know a reasonable upper limit for the number of photons, otherwise the sum is infinite (we could also try doing this analytically and doing the infinite sum; in general this sometimes works but not always).  We can look at the Poisson distribution to find a reasonable value, which in this case we take as $10 \\lambda$\n",
      "\n",
      "$\n",
      "P(e) = \\sum_{n=1}^{n=\\infty} P(e|n) P(n)\n",
      "$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def P_e(e):\n",
      "    N = arange(1,int(10*lam))\n",
      "    return sum(P_e_given_n(e,n)*P_n(n) for n in N)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can now plot the distribution as a function of $e$.  We can say waht we expect it to look like - peaks for each possible photon count with a Gaussian around each one for the corresponding energy."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "e = arange(0.01, 6.0, 0.01)\n",
      "plot(e,P_e(e))\n",
      "xlabel(\"Energy\")\n",
      "ylabel(\"Probability\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<matplotlib.text.Text at 0x1098ce610>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEPCAYAAACk43iMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGW95/HPlywCSSBAECEEQRZFh10IikrELXIVFBcG\ntwE3hisuV3FcZq6Jw3UcvYrrFbkoiF4FHESFK5sILctLgWgAhQTDaoLIDrIJBH7zx3MaKp3u6qqu\n85xzuur7fr361VXVp079Gjr1rWc5z6OIwMzMBtM6dRdgZmb1cQiYmQ0wh4CZ2QBzCJiZDTCHgJnZ\nAHMImJkNsKwhIOkESbdL+kObYxZIWirpj5KGctZjZmZrUs7rBCS9FHgQ+H5E7DTKz2cDlwKviYhV\nkuZExF3ZCjIzszVkbQlExMXAvW0OeRvwk4hYVRzvADAzq1DdYwLbAxtLulDSEknvrLkeM7OBMrXm\n158G7A68Algf+I2k30bEinrLMjMbDHWHwErgroh4BHhE0kXALsAaISDJCxyZmU1ARKjdz+sOgZ8D\n35Q0BXgGMB84ZrQDx/tFJjNJiyNicd115NLPv18//27g32+y6+QDdNYQkHQysC8wR9JKYBGpC4iI\nOC4ilks6B7gaeBI4PiKuzVmTmZk9LWsIRMQhHRzzJeBLOeswM7PR1T07yJKhugvIbKjuAjIaqruA\nzIbqLiCzoboLqFvWi8XKIin6eUzAzCyHTt473RIwMxtgDgEzswHmEDAzG2AOATOzAeYQMDMbYA4B\nM7MB5hAwMxtgDgEzswHmEDAzG2AOATOzAeYQMDMbYA6Bmkm8RGKDuusws8HkEKiRxObAecDimksx\nswHlEKjX64GlwMF1F2Jmg8khUK+dgZ8AMyXm1F2MmQ0eh0C9dgGuJG2vuXPNtZjZAHII1Ou5wDLg\n2uK2mVmlsoaApBMk3S7pD+Mct6ek1ZIOyllPk0isC2wA3A6sBObVW5GZDaLcLYETgYXtDpA0BfgC\ncA4wSFtIbgn8JYIngVU4BMysBllDICIuBu4d57APAqcBd+aspYHmkd78wS0BM6tJrWMCkuYCBwLH\nFg81f9f78swjvfmDQ8DMalL3wPBXgU9GRJC6ggapO+hZwG3F7b8Am9dYi5kNqKk1v/4ewCmSAOYA\nr5X0eEScMfJASYtb7g5FxFAlFeazCXB3cfshQBIzInioxprMbBKTtABY0M1zag2BiHjO8G1JJwJn\njhYAxbGLq6qrIpsANwBEEBJ3AJuCQ8DMJqb4cDw0fF/SovGekzUEJJ0M7AvMkbQSWARMA4iI43K+\n9iTQ2hKANDD+TODmWqoxs4GUNQQi4pAujj0sZy0NNFoIbFpTLWY2oOoeGB5kI0NguDvIzKwyDoH6\nuCVgZrVzCNRAQsDGrBkC9wAb1VORmQ0qh0A9ZgGPRfBoy2P3AbNrqsfMBpRDoB4ju4LAIWBmNXAI\n1GO0ELgfh4CZVcwhUA+3BMysERwC9XAImFkjOATq4RAws0ZwCNRjrBDYsIZazGyAOQTqMZv0pt/q\nYWC6xPQa6jGzAeUQqMdM4IHWByII3Bows4o5BOoxE3hwlMc9TdTMKuUQqMcsRrQECh4cNrNKOQTq\nMVZLwCFgZpVyCNTDIWBmjeAQqIe7g8ysERwC9XBLwMwawSFQj3azgzxF1MwqkzUEJJ0g6XZJfxjj\n52+XdJWkqyVdKmnnnPU0QbGhzAzGbgl4Yxkzq0zulsCJwMI2P78ReFlE7AwcDfx75nqaYH3g0Qie\nGOVnfyONF5iZVSJrCETExcC9bX7+m4i4v7h7GbBlznoaYqyuIEiDxQ4BM6tMk8YE3gOcVXcRFRhr\nZhCkcHAImFllptZdAICklwPvBvZpc8zilrtDETGUuaxc3BIwsywkLQAWdPOc2kOgGAw+HlgYEe26\njhZXVlReDgEzy6L4cDw0fF/SovGeU2t3kKStgNOBd0TE9XXWUqF23UEPkELCzKwSWVsCkk4G9gXm\nSFoJLAKmAUTEccBnSFMij5UE8HhE7JWzpgZwS8DMGkMRUXcN45IUEaG66yiDxGHAvhEcOsrPppM2\nl5lW7C9gZjZhnbx3Nml20KAYszsogseAJ4B1K63IzAaWQ6B67bqDwF1CZlYhh0D1OgkBDw6bWSUc\nAtWbhVsCZtYQDoHqrbXJ/Ai+atjMKuMQqJ7HBMysMRwC1XMImFljOASq1+6KYfDAsJlVyCFQPbcE\nzKwxHALVGy8EPDBsZpVxCFSvk+4gh4CZVcIhUD13B5lZYzgEKiQxBViPtEjcWDwwbGaVcQhUawbw\nUARPtjnGLQEzq4xDoFrjdQWBB4bNrEIOgWp1GgIzKqjFzMwhULHxZgZBCgGPCZhZJRwC1eq0JeAQ\nMLNKOASq5RAws0bJGgKSTpB0u6Q/tDnm65JWSLpK0m4562kAdweZWaPkbgmcCCwc64eS9ge2i4jt\ngfcDx2aup26dtAQeAaZLTK2gHjMbcFlDICIuBu5tc8gBwEnFsZcBsyVtlrOmmo0bAhEE8BCeIWRm\nFah7TGAusLLl/ipgy5pqqUIn3UHgLiEzq0gTuhw04n6MepC0uOXuUEQM5Sooo5m0bxkNcwiYWdck\nLQAWdPOcukPgVmBey/0ti8fWEhGLqygos5ms2fIZi0PAzLpWfDgeGr4vadF4z6m7O+gM4F0AkvYG\n7ouI2+stKatuuoM8JmBm2WVtCUg6GdgXmCNpJbAImAYQEcdFxFmS9pd0PWkw9LCc9TRAJ7ODwC0B\nM6vIuCEg6QDgPyOi3cqXo4qIQzo45shuzzuJOQTMrFE66Q46GLhe0hclPS93QX3Os4PMrFHGDYGI\neDuwG3Aj8D1Jv5H0fkle7rh7bgmYWaN0NDAcEfcDpwGnAlsAbwSWSvpQxtr6Uach8BAOATOrwLgh\nIOlAST8lTTuaBuwZEa8FdgY+mre8vuPuIDNrlE5mBx0EfCUiLmp9MCIelvTePGX1rW66g/p5+Qwz\na4hOuoNuHxkAkr4AEBHnZ6mqD0lMI4Xu3zs43C0BM6tEJyHwqlEe27/sQgbATODBYoG48TgEzKwS\nY3YHSToC+Edg2xH7AcwCLs1dWB/qtCsIHAJmVpF2YwI/As4G/i/wCZ5e6O2BiLg7d2F9yCFgZo3T\nLgQiIm6W9AFGrOwpaeOIuCdvaX2n05lB4BAws4q0C4GTgX8Afsfoyztvk6Wi/uWWgJk1zpghEBH/\nUHzfurJq+lu3IeBVRM0su3YDw7u3e2JE/L78cvqau4PMrHHadQcdwxi7fBVeXnIt/a6blsBDwEwJ\ndTil1MxsQtp1By2osI5B0HEIRLBa4jFgXeCRrFWZ2UBr1x20X0RcIOlNjNIiiIjTs1bWf7rpDoKn\nu4QcAmaWTbvuoH2BC4DXM3q3kEOgOzOBbrbOHF5J9M485ZiZte8OWlR8P7SyavrbTOCGLo734LCZ\nZdfJUtJzJH1D0lJJv5f0NUmbVFFcn5lod5CZWTadLCB3CnAHaUnpN5O6J07t5OSSFkpaLmmFpE+M\n8vM5ks6RdKWkP0o6tIvaJ5tuZgeBQ8DMKtBJCDwrIo6OiJsi4saI+Bc6WOte0hTgm8BC4PnAIZJ2\nHHHYkcDSiNgVWAB8WVInexxMRg4BM2ucTkLgPEmHSFqn+DoYOK+D5+0FXB8RN0fE46QWxYEjjrkN\n2KC4vQFwd0Ss7rT4SWYm7g4ys4ZpN0X0QZ6eFfQR4AfF7XVIM1c+Ns655wIrW+6vAuaPOOZ44AJJ\nfyH1mb+1s7InpVm4JWBmDdNudlCvb0CdXOn6aeDKiFggaVvgl5J2iYi1PjFLWtxydygihnqsr2ru\nDjKzrCQtIHWtd6yj/ndJGwHbk65gBWDklpOjuBWY13J/Hqk10OrFwOeK890g6SbgucCSkSeLiMWd\n1NpgDgEzy6r4cDw0fF/SovGe08kU0fcBF5HGAT4LnAss7qCeJcD2kraWNB04GDhjxDHLgVcWr7MZ\nKQBu7ODck4qE8BRRM2ugTgaGP0wa5L05Il4O7AbcP96TigHeI0mhcS1wakQsk3S4pMOLw/4P8EJJ\nVwHnA/+jTzermQ5EBI918RwvJ21m2XXSHfT3iHhEEpLWjYjlkp7byckj4mzSFpWtjx3Xcvsu0rIU\n/a7briBwS8DMKtBJCKwsxgR+Rhq4vRe4OWtV/afbriBwCJhZBcYNgYh4Y3FzsaQh0nz+c3IW1Yfc\nEjCzRup0dtAewEtI0z4viYhu+rbNIWBmDdXJ7KDPAN8DNgbmACdK+ufMdfWbiXQHDS8lbWaWTSct\ngXcAO0fE3wEkfR64Cjg6Z2F9xi0BM2ukTqaI3gqs13J/Xda+6MvacwiYWSO1WzvoG8XN+4FrJA0v\nGvcq4PLchfUZzw4ys0Zq1x30O9JA8BLS9NDhtYCG6GxdIHvaRFoCfwemSUyL4PEMNZmZtV1A7nvD\ntyU9A9ihuLu8WBraOtd1CEQQ0lNXDd+XpSozG3jjDgwXq9KdBNxSPLSVpP8WEb/OWVifmcWay2p3\narhLyCFgZll0MjvoGODVEXEdgKQdSBvE7J6zsD4zke4g8LiAmWXWyeygqcMBABARf6LDi8zsKQ4B\nM2ukTt7MfyfpO8B/AALezijr/VtbE5kdBA4BM8uskxD476QloT9U3L8Y+Fa2ivpTLy0BLydtZtm0\nDQFJU4GrIuJ5wJerKakvuTvIzBqp7ZhAsTHMdZKeXVE9/crdQWbWSJ10B21MumL4ctKiZgAREQfk\nK6vvuCVgZo3USQj8r+K7Wh7zFcPdmWgIeCVRM8uq3dpB65EGhbcDrgZO8JXC3Ss2me+lJTCr3IrM\nzJ7WbkzgJGAPUgDsD3yp25NLWihpuaQVkj4xxjELJC2V9Mdi57J+sz7wWASrJ/BcdweZWVbtuoN2\njIidACR9F7iimxNLmgJ8E3glaTnqKySdERHLWo6ZDfwb8JqIWCVpTre/wCQwC/jbBJ/rEDCzrNq1\nBJ765FrMEurWXsD1EXFz0Y10CnDgiGPeBvwkIlYVr3PXBF6n6SY6MwgcAmaWWbsQ2FnSA8NfwE4t\n9zv5ZDuXNRdNW1U81mp7YGNJF0paIumd3ZU/KWyAQ8DMGqrdUtJTejx3JzOIppEWonsFqe/8N5J+\nGxErRh4oaXHL3aGIGOqxvqq4O8jMKlGs+rygm+fkXAjuVmBey/15rL0t5Urgroh4BHhE0kXALsBa\nIRARizPVmZu7g8ysEsWH46Hh+5IWjfecTlYRnaglwPaStpY0HTgYOGPEMT8HXiJpiqT1gfnAtRlr\nqoO7g8yssbK1BCJitaQjgXOBKcB3I2KZpMOLnx8XEcslnUOahvokcHxE9FsIuDvIzBpLEc2/+FdS\nRITGP7J5JD4ObBbBURN47gbAqgg2KL8yM+t3nbx35uwOsqSX7qCHgBnFVcdmZqVzCOQ34e6gCJ4A\nHiXNnDIzK51DIL9eZgdRPNfrB5lZFg6B/HrpDoLUivCYgJll4RDIr5fZQeAQMLOMHAL59doddD+w\nYUm1VE5i3bprMLOxOQTyG8juIIl1JI4DHpD4icT0umsys7U5BPLrtSUwKUMAOJS0H8WmpAvePl5r\nNWY2KodAfr2OCUy67iCJqcBngA9EcB9wBPBPkmc5mTWNQyCj4iKvQWwJvBq4LYLLACK4Efg1af0o\nM2sQh0BezwCI4NEeznE/ky8E3kLaRKjVqcCba6jFzNpwCOTVa1cQxfMnTXdQMQB8AHD6iB+dBbxY\nYuPqqzKzsTgE8up1ZhBMvu6gfYEVEWvsKkcEDwLnA6+rpSozG5VDIK8Ngft6PMdk6w7aj7R8+GjO\nK35uZg3hEMhrNr2HwKTqDiJtbTc0xs+G6HLrOzPLyyGQV1khMClaAhIzgZ2A345xyHXAuhJbV1aU\nmbXlEMhrQ1J3Ti8mU3fQPsDvInhktB9GEKTWwL5VFmVmY3MI5DVQLQFSCFwyzjGXFMeZWQM4BPIq\nIwQeADaYJLuL7QXpArE2rgBeWEEtZtaBrCEgaaGk5ZJWSPpEm+P2lLRa0kE566lBz91BETzOJNhd\nrAipPUlv8u1cDTzPq4uaNUO2EJA0BfgmsBB4PnCIpB3HOO4LwDkwKT7tdqOMlgBMji6h5wAPR3Bb\nu4OK8YLrgF0qqcrM2srZEtgLuD4ibo6Ix0nLCBw4ynEfBE4D7sxYS13KDIGmTxPdC7i8w2PdJWTW\nEDlDYC6scdXoquKxp0iaSwqGY4uHImM9dSjjYjGAe0mB0mQvZPyuoGFLcAiYNcLUjOfu5A39q8An\nIyIkiTbdQZIWt9wdioih3sqrxGx6nyIKKQSavubOrsC/dnjsFcCHMtZiNpAkLaDLCzJzhsCtwLyW\n+/NIrYFWewCnpPd/5gCvlfR4RJwx8mQRsThTnTmV1R10L7BRCefJohgU3gW4qsOnXANsIzGzWFPI\nzEpQfDgeGr4vadF4z8nZHbQE2F7S1pKmk9aSX+PNPSKeExHbRMQ2pHGBI0YLgEmsrO6ge2h2S2Au\n8ATw104OjuAxUhDsmrMoMxtfthCIiNXAkaTFxK4FTo2IZZIOl3R4rtdtCol1SCHQ61LS0PCWAEUr\noLgiuFNLcQiY1S5ndxARcTZw9ojHjhvj2MNy1lKDmaQpk6tLONc9wFYlnCeXbrqChi0lXVdgZjXy\nFcP5lNUVBM0fGJ5ICFyJWwJmtXMI5FPWoDCklkCTu4N2pfsQuBrYUWJahnrMrEMOgXzKmh4KDW4J\nSMwgzfxa3s3zIngYuAVY6ypyM6uOQyCfTYC7SjpXkweG/wuwvFjjqFtLgd1KrsfMuuAQyGcOcHdJ\n52ryFNGJjAcM87iAWc0cAvnMYTBaAhMZDxjmloBZzRwC+WxCSS2BYuXNkFivjPOVrJeWwFXArpNk\nrwSzvuQQyKfMlgA0cHC4uCBuJyYYAhHcATwE3nPYrC4OgXxKawkUmtgltA1wXwT39HAOXzlsViOH\nQD5ltwSaODjcy3jAsCvxuIBZbRwC+ZQdAneTWhdNsgvpTbwXbgmY1cghkE/Z3UF3AM8s8Xxl2JXe\nQ8AtAbMaOQQykJhCWjvo3hJP29QQWNrjOW4CNpCYU0I9ZtYlh0AeGwH3R/BEiedsVAhIbEIKupt7\nOU8ET1JMFS2hLDPrkkMgj7LHAwBup0EhQDEoXLyJ98rjAmY1cQjkkSME7gA2K/mcvShjPGCYxwXM\nauIQyKPsQWFoWHcQ5YaAWwJmNXEI5PFM4M6Sz9nPIXAtaeP59Us6n5l1KHsISFooabmkFZI+McrP\n3y7pKklXS7pU0s65a6rAFsCtJZ/zbmBDKe+WoJ2QWBfYjvTm3bNi4/nrSMtSm1mFsoaApCnAN4GF\nwPOBQySN3ETkRuBlEbEzcDTw7zlrqsgWwG1lnrAYgL0b2LTM807QC4AVEfy9xHN6XMCsBrlbAnsB\n10fEzRHxOHAKcGDrARHxm4gY3oHrMmDLzDVVYQvgLxnO25QuoTK7goZ5bwGzGuQOgbnAypb7q4rH\nxvIe4KysFVVjc/KEQFOmie5O+SHgwWGzGuTuX45OD5T0cuDdwD75yqlM6d1BhaZME50P/Kjkc14F\n7CQxpeSL7MysjdwhcCtpE/Jh80itgTUUg8HHAwsjYtSlFiQtbrk7FBFD5ZVZnmLJiE1Jn9rLdhsp\nYGpTbGzzfOD3ZZ43gvslbge2p8tN680skbQAWNDNc3KHwBJge0lbk7pHDgYOaT1A0lbA6cA7IuL6\nsU4UEYuzVVmuZwJ3T3Dj9fH8mfQmWafdgGXFbmdlGx4cdgiYTUDx4Xho+L6kReM9J+uYQESsBo4E\nziVNJzw1IpZJOlzS4cVhnyGttXOspKWSLs9ZUwVyDQpDCoGtMp27U/NJA/g5eFzArGLZ55xHxNnA\n2SMeO67l9nuB9+auo0K5xgOgOSFw9rhHTcyVwIczndvMRuErhsuXa2YQpBB4dqZzd2pv8rUEfgvM\nb8IFcWaDwiFQvq0YZfC7JHcD60rMynT+tiQ2Iy0f/acc54/gLlLQuUvIrCIOgfJtC9yQ48QRBOlN\nct54x2byYuCykpaPHstFwMsynt/MWjgEyrcdMOYspxLUOS7wcuDCzK9xEbBv5tfomcQciQMl3iWx\nr8Qz6q7JbCIcAuXL1hIo3EJ9IbCAlulnmVwEvFRq5t+mxKYS3wFWAEcArwK+CNwi8dHiOhGzSaOR\n/9AmK4mNgSmUv6FMqz8DW2c8/6iKPYCfDfwu5+tE8BfS2McLcr7OREi8GLgauA/YNoKFEbwzgvnA\nK4DXAUPeL9kmE4dAubYFbij67nNZDjwv4/nHsgC4NILVFbzWRaSup8aQeAPwc+CwCI6K4J7Wn0dw\nDfBK4BLgUoltaijTrGsOgXLlHg+AdNFdHZ+S9wfOqei1zgFeW9FrjUviJaQlzhdGjP3fIIInI/gU\n8C3gfInNq6rRbKIcAuXKPR4AqS96q2Jjl0oU/fP7A7+o6CXPA14iMaOi1xuTxLbAacA7IzrrCovg\na8CJwLkSG+asz6xXDoFyPY9Mc+iHFbtw3QTskPN1RtgduDcie8ABaTE54ApgvypebywSG5GC77MR\nnNvl0z9H6hr6kQeLrckcAuXahTRwmNs1VNsl9AbgPyt8PYrXe0PFr/kUiWmkFsDZERzb7fOLcaEP\nA+sBny+5PLPSOARKUswT3470Bp1bZSEgIdLKr6dU8Xotfgy8sY7598Xv/C3gYeCoiZ6nWEn2LcCb\nJN5ZUnlmpXIIlGcn4PqS990dy9VUtx/vXsBqSt4/YDwRrCL9nvtX+bqFj5F+77f1usFNBHcDBwDH\nSMwvozizMjkEyvMi4DcVvdYlwD4V9TW/C/hh5mmvY/lh8fqVKaaCfgR4XQQPlHHOYvrou4HTpdpX\ngTVbg0OgPJWFQAR3kFYqzbrQmsRsUlfQd3K+ThsnAy+Tqrk4TmJv0g53b4hYY2/snkVwJvBl4Bee\nMWRN4hAoQTGFcj/yL6nQ6tfkX2PnPcA5xVW8lYvgQdJUyyNzv5bEc4GfAYdGsCTTy3yF9P/tNInp\nmV7DrCsOgXLsAvwtgpsqfM0hMl5VKzGT1Df+5Vyv0aGvA4dJPCvXC0jMI22U8+mIfNdCtMwYehA4\n1UFgTeAQKMebSZ8iq3QOaaG1XOvUfBT4dacXSOUSwZ+B7wGfzXH+YnmHi4BvRnBCjtdoVQw0HwwI\n+H9efdTq5hDoUbEL1ttJ/deVKS6oOrN47VJJvAD4IPDpss89QZ8DXi+V2/1VzNa5GPhSBMeUee52\nigv+3kqadXWOF5yzOmUNAUkLJS2XtELSJ8Y45uvFz6+SVNW0xzIdBNwWwdIaXvsE4PAyt2OU2IAU\naJ+suHtrTMVibe8BfiCxRa/nk1hH4h9JIXpEBP/W6zm71RIElwGXS+xTdQ0AEhtK7CdxlMS3Jc6T\nWCaxSuI+iQckbpNYIXGJxPclFkscJLFVcU2FTWKKyDPzT9IU4DrSyoq3kpYBOCQilrUcsz9wZETs\nL2k+8LWI2HuUc0VENO6PTWJ94CrgyAksK9ByHi2IiKEJvL6A84EzI/jqRF+/5XyzSStlXgN8oKxp\noRP9/dY+D58C3gG8priOYCLn2I00QLsuaRB4eW819f67FdNSvwX8FDg6gr/2cr42rzONNKNsb2A+\nsCcwF7iStET4dcCNwErgfuBB2O7FcP1SYBbwLOA5pIsidymePwVYAlxO2iP68uLaiEmhrL/Npurk\nvTPnht57AddHxM1FMacABwLLWo45ADgJICIukzRb0mYRcXvGukpRfPr+LumPfsIBUFjABGYWRRAS\nR5CWLl4awa8nWoDEfsBxwFnAP5V8XcACSpg5FcHnJZ4ArpA4Cji5k60ui/9XrwHeT3rj+hzw7V4v\nBCssoMffLYKfSVwEfAa4VuJU0r+Lyye6lWdxDck2wM6kN/0XkS4wvIk0lflC0mY417ZbHly6Ya8I\nzgJuI62LdVHLawjYgvTfdC/S1dV7StxOauH8tvh+VdHyaaIFVDurr3FyhsBcWGOu9SpY64rJ0Y7Z\nEmhsCBSrd74CWETa/OSgOuuJ4E8SB5OmHR4PfLeThd6Kzep3IO3n+xbSp7yPRfDTrAX3KIIvSlwK\nHAP8i8RPSW80t5A+vU4DZpA2wNmOtC/yPqQluE8EDong4Tpqb6fo8vqIxBeBw0i1bipxMbCU9An9\nL8BDxZeA9Yuv2aR/S3NJ+0/vSFrM8E7gj6RP6f+bFCr3l1hzkFr5t1JMjCjCZ0fSv/X5pODdXuIW\nUoisIK20+1fgDtK/9btJS3Q82smHjyJ8ppJac+sV30d+PQOYPs4XcMCLJD4CPFF8PVl8fwx4vPg+\n0dura7rIsis5Q6DTX35kU2XU50n8ojh2+Hh1cX8izxntHBsBm5Oazl8FTsm86XpHIrhAYg/gU8Al\nRTfVStI/rMdIf9TrtXxtVHy/kTQw+nnSQmlVbBjTswguLS7s2p20vePBpDfA2aTf9xHSDmw3kFpr\n74rIuttbaYprMj4HfE5iLimkX0BaPmNzUsDNIL1ZPUz6Xf9GeiNeBfwS+AbpE34pVzx3Wf8TpOD5\nI+m/PRLrkbqRdgC2J7VOXgU8s/iaQ/p7nC7xCOl3Ctb+9ydSyK9L+v3/3vL1SMvtR4vvj7X5ejxV\nPHMDUotpCmmMdErxNY0UFMPfu709DZgqrREO7QLj8Zb/jOO9d3X6PjdEB3KOCewNLI6IhcX9TwFP\nRsQXWo75NjAUEacU95cD+47sDpLU+DQ1M2uiOscElgDbS9qa1Iw9mLQEQaszSFeDnlKExn2jjQc0\ncVDYzKwfZAuBiFgt6UjgXFLz6rsRsUzS4cXPj4uIsyTtL+l6Uj/nYbnqMTOztWXrDjIzs+Zr9BXD\nnVxsNplJOkHS7ZL+UHctZZM0T9KFkq6R9EdJH6q7pjJJWlfSZZKulHStpL7cPUzSFElLJZ1Zdy1l\nk3SzpKuL3+/yuuspUzHd/jRJy4q/z7Wuv3rq2Ka2BDq52Gyyk/RS0mJi34+Inequp0ySngU8KyKu\nlDSTNKPqDX32/2/9iHhY0lTSHg9HRcQldddVJkkfBfYAZkXEAXXXUyZJNwF7RMQ9dddSNkknAb+O\niBOKv88dnoBvAAAEEElEQVQZETHq9OAmtwSeutgsIh4nbW94YM01lSoiLgburbuOHCLirxFxZXH7\nQdJFgj0v+dAkETF8vcF00rhXX72ZSNqSNDX1O6w9lbtf9N3vJWlD4KURcQKk8dmxAgCaHQKjXUg2\nt6ZarAfFDLHdSBd19Q1J60i6knTB04URcW3dNZXsK8DHof5rYTIJ4HxJSyS9r+5iSrQNcKekEyX9\nXtLxktYf6+Amh0Az+6msK0VX0GnAh4sWQd+IiCcjYlfSVe4vk7Sg5pJKI+l1wB0RsZQ+/LRc2Cci\ndgNeC3yg6J7tB1NJF1J+KyJ2J828/ORYBzc5BG4lXQI/bB5MbNEwq4ekacBPgP+IiKr3W6hM0dT+\nBfDCumsp0YuBA4p+85OB/SR9v+aaShURtxXf7yQt3rdXvRWVZhWwKiKuKO6fRgqFUTU5BJ662EzS\ndNLFZmfUXJN1SJJIywZcGxE9r3DaNJLmSJpd3F6PtAxCHcuJZxERn46IeRGxDfBfgQsi4l1111UW\nSetLmlXcngG8GuiLWXoR8VdgpaQdiodeSVoZeFQ5rxjuyVgXm9VcVqkknUzaJ3gTSSuBz0TEiTWX\nVZZ9SMs+Xy1p+M3xUxFxTo01lWlz4CRJ65A+TP0gIn5Vc0059Vv37GbAT9NnFaYCP4yI8+otqVQf\nBH5YfIC+gTYX4jZ2iqiZmeXX5O4gMzPLzCFgZjbAHAJmZgPMIWBmNsAcAmZmA8whYGY2wBp7nYBZ\nbpKeAK5ueejkiPhiXfWY1cHXCdjAkvRARMwq+ZxTI2J1mec0y8ndQWYjFJuNLJb0u2LTkecWj88o\nNgK6rFid8YDi8UMlnSHpV8AvJa0n6cfFhjqnS/qtpD0kHSbpKy2v8z5Jx9T0a5oBDgEbbOsVu0oN\nf72leDyAOyNiD+BY4Kji8f8J/Coi5gP7Af/askTvbsCbIuLlwAeAuyPiBcA/kzZlCeDHwOuLDZMA\nDiWtr2RWG48J2CB7pFhKeDSnF99/DxxU3H416U18OBSeAWxFeoP/ZUTcVzy+D/BVgIi4RtLVxe2H\nJF1QnGM5MC0ixlzYy6wKDgGz0T1afH+CNf+dHBQRK1oPlDSftGb7Gg+Pcd7vkFoUy4ATSqjTrCfu\nDjLr3LnAh4bvSBpuRYx8w78UeGtxzPOBp/aPjojLSZvQvI20Tr9ZrdwSsEG2Xssy1wBnR8SnRxwT\nPL2M8tHAV4vunXWAG4EDRhwD8C3SMtPXAMtJa7m37vH6Y2CXdvu+mlXFU0TNSlbsMTAtIh6VtC3w\nS2CH4amjks4EjomIC+us0wzcEjDLYQZwQbG9poAjik2SZgOXAVc6AKwp3BIwMxtgHhg2MxtgDgEz\nswHmEDAzG2AOATOzAeYQMDMbYA4BM7MB9v8B+JufWYBp2mAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1098c5b50>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've been a little careless about the n=0 case, which is quite likely since $\\lambda$ is small.  That looks like a $\\delta$ function spike at Energy=0 with total probability $1/e$"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}
	{
	"metadata": {
	"kernelspec": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 2
	},
	"display_name": "IPython (Python 2)",
	"language": "python",
	"name": "python2"
	},
	"name": "",
	"signature": "sha256:55c4f9d7a7aae3543fe006737001ce10ac9cc48f7a29a15371c3b9c8fa640653"
	},
	"nbformat": 3,
	"nbformat_minor": 0,
	"worksheets": [
	{
	"cells": [
	{
	"cell_type": "heading",
	"level": 1,
	"metadata": {},
	"source": [
	"Lecture 1 Medium Problem Worked Solution"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"This is an ipython notebook, an interactive python session saved as a web page. You can download this and run it on your own machine - google \"ipython notebook\" for details.\n",
	"\n",
	"These special commands set up the notebook environment, and import numpy, matplotlib, and the factorial function."
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"%pylab inline\n",
	"from math import factorial"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": [
	"Populating the interactive namespace from numpy and matplotlib\n"
	]
	}
	],
	"prompt_number": 2
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"The general approach when you have a probability problem like this one is to write down all the thing you do know as probability distributions, and then use the laws of probability to convert them to what you want. \n",
	"\n",
	"In this case we have some quantities specified as conditional probabilities, like $P(n)$ and the $P(e_i)$ for a single photon, and we want $P(e)$, the overall distribution for the sum of all the energies.\n",
	"\n",
	"First define the constants of the problem (the word \"lambda\" is a python keyword so we can't use that - we call it \"lam\" instead)."
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"mu = 1.0 #The mean of the Gaussian for each photon's energy, in keV\n",
	"sigma = 0.1 #The standard deviation of the Gaussian, in keV\n",
	"lam = 1.0 #The mean arrival rate in the Poisson distribution, in s^-1"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 3
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"The first probability we have is the easiest one, the probability of getting a given number of photons. This is a Poisson distribution, which is discrete since you can't have half a photon.\n",
	"\n",
	"A Poisson distribution is:\n",
	"\n",
	"$\n",
	"P(n) = \\frac{\\lambda^n \\exp(-\\lambda)} {n!}\n",
	"$"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"def P_n(n):\n",
	" return lam*n exp(-lam) / factorial(n)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 4
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We are told that the probability of the energy of each photon is Gaussian distributed. We need the distribution of the total energy of the n photons that we are assuming are arriving in a second. \n",
	"\n",
	"There's a handy formula for the distribution of the sum of Gaussians, which is just that if $X \\sim N(\\mu_1, \\sigma_1^2)$ and $Y \\sim N(\\mu_2, \\sigma_2^2)$, then $X+Y \\sim N(\\mu_1+\\mu_2, \\sigma_1^2+\\sigma_2^2)$. \n",
	"\n",
	"We can generalize this to give us the distribution of the sum of n identical distributions as $N(n\\mu, n \\sigma^2)$. So:\n",
	"\n",
	"$\n",
	"P(e\|n) = (2 \\pi n \\sigma^2)^{-\\frac{1}{2}}\\exp{[-(e-n \\mu)^2 / 2 n \\sigma^2]}\n",
	"$\n"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"def P_e_given_n(e,n):\n",
	" return exp(-0.5(e-nmu)*2/(nsigma*2)) / sqrt(2pinsigma**2)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 5
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"Now we can use the marginalization rule, which says we can sum over all possible values of n to get the total energy.\n",
	"\n",
	"We need to know a reasonable upper limit for the number of photons, otherwise the sum is infinite (we could also try doing this analytically and doing the infinite sum; in general this sometimes works but not always). We can look at the Poisson distribution to find a reasonable value, which in this case we take as $10 \\lambda$\n",
	"\n",
	"$\n",
	"P(e) = \\sum_{n=1}^{n=\\infty} P(e\|n) P(n)\n",
	"$"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"def P_e(e):\n",
	" N = arange(1,int(10*lam))\n",
	" return sum(P_e_given_n(e,n)*P_n(n) for n in N)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 6
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We can now plot the distribution as a function of $e$. We can say waht we expect it to look like - peaks for each possible photon count with a Gaussian around each one for the corresponding energy."
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"e = arange(0.01, 6.0, 0.01)\n",
	"plot(e,P_e(e))\n",
	"xlabel(\"Energy\")\n",
	"ylabel(\"Probability\")\n"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"metadata": {},
	"output_type": "pyout",
	"prompt_number": 7,
	"text": [
	"<matplotlib.text.Text at 0x1098ce610>"
	]
	},
	{
	"metadata": {},
	"output_type": "display_data",
	"png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEPCAYAAACk43iMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGW95/HPlywCSSBAECEEQRZFh10IikrELXIVFBcG\ntwE3hisuV3FcZq6Jw3UcvYrrFbkoiF4FHESFK5sILctLgWgAhQTDaoLIDrIJBH7zx3MaKp3u6qqu\n85xzuur7fr361VXVp079Gjr1rWc5z6OIwMzMBtM6dRdgZmb1cQiYmQ0wh4CZ2QBzCJiZDTCHgJnZ\nAHMImJkNsKwhIOkESbdL+kObYxZIWirpj5KGctZjZmZrUs7rBCS9FHgQ+H5E7DTKz2cDlwKviYhV\nkuZExF3ZCjIzszVkbQlExMXAvW0OeRvwk4hYVRzvADAzq1DdYwLbAxtLulDSEknvrLkeM7OBMrXm\n158G7A68Algf+I2k30bEinrLMjMbDHWHwErgroh4BHhE0kXALsAaISDJCxyZmU1ARKjdz+sOgZ8D\n35Q0BXgGMB84ZrQDx/tFJjNJiyNicd115NLPv18//27g32+y6+QDdNYQkHQysC8wR9JKYBGpC4iI\nOC4ilks6B7gaeBI4PiKuzVmTmZk9LWsIRMQhHRzzJeBLOeswM7PR1T07yJKhugvIbKjuAjIaqruA\nzIbqLiCzoboLqFvWi8XKIin6eUzAzCyHTt473RIwMxtgDgEzswHmEDAzG2AOATOzAeYQMDMbYA4B\nM7MB5hAwMxtgDgEzswHmEDAzG2AOATOzAeYQMDMbYA6Bmkm8RGKDuusws8HkEKiRxObAecDimksx\nswHlEKjX64GlwMF1F2Jmg8khUK+dgZ8AMyXm1F2MmQ0eh0C9dgGuJG2vuXPNtZjZAHII1Ou5wDLg\n2uK2mVmlsoaApBMk3S7pD+Mct6ek1ZIOyllPk0isC2wA3A6sBObVW5GZDaLcLYETgYXtDpA0BfgC\ncA4wSFtIbgn8JYIngVU4BMysBllDICIuBu4d57APAqcBd+aspYHmkd78wS0BM6tJrWMCkuYCBwLH\nFg81f9f78swjvfmDQ8DMalL3wPBXgU9GRJC6ggapO+hZwG3F7b8Am9dYi5kNqKk1v/4ewCmSAOYA\nr5X0eEScMfJASYtb7g5FxFAlFeazCXB3cfshQBIzInioxprMbBKTtABY0M1zag2BiHjO8G1JJwJn\njhYAxbGLq6qrIpsANwBEEBJ3AJuCQ8DMJqb4cDw0fF/SovGekzUEJJ0M7AvMkbQSWARMA4iI43K+\n9iTQ2hKANDD+TODmWqoxs4GUNQQi4pAujj0sZy0NNFoIbFpTLWY2oOoeGB5kI0NguDvIzKwyDoH6\nuCVgZrVzCNRAQsDGrBkC9wAb1VORmQ0qh0A9ZgGPRfBoy2P3AbNrqsfMBpRDoB4ju4LAIWBmNXAI\n1GO0ELgfh4CZVcwhUA+3BMysERwC9XAImFkjOATq4RAws0ZwCNRjrBDYsIZazGyAOQTqMZv0pt/q\nYWC6xPQa6jGzAeUQqMdM4IHWByII3Bows4o5BOoxE3hwlMc9TdTMKuUQqMcsRrQECh4cNrNKOQTq\nMVZLwCFgZpVyCNTDIWBmjeAQqIe7g8ysERwC9XBLwMwawSFQj3azgzxF1MwqkzUEJJ0g6XZJfxjj\n52+XdJWkqyVdKmnnnPU0QbGhzAzGbgl4Yxkzq0zulsCJwMI2P78ReFlE7AwcDfx75nqaYH3g0Qie\nGOVnfyONF5iZVSJrCETExcC9bX7+m4i4v7h7GbBlznoaYqyuIEiDxQ4BM6tMk8YE3gOcVXcRFRhr\nZhCkcHAImFllptZdAICklwPvBvZpc8zilrtDETGUuaxc3BIwsywkLQAWdPOc2kOgGAw+HlgYEe26\njhZXVlReDgEzy6L4cDw0fF/SovGeU2t3kKStgNOBd0TE9XXWUqF23UEPkELCzKwSWVsCkk4G9gXm\nSFoJLAKmAUTEccBnSFMij5UE8HhE7JWzpgZwS8DMGkMRUXcN45IUEaG66yiDxGHAvhEcOsrPppM2\nl5lW7C9gZjZhnbx3Nml20KAYszsogseAJ4B1K63IzAaWQ6B67bqDwF1CZlYhh0D1OgkBDw6bWSUc\nAtWbhVsCZtYQDoHqrbXJ/Ai+atjMKuMQqJ7HBMysMRwC1XMImFljOASq1+6KYfDAsJlVyCFQPbcE\nzKwxHALVGy8EPDBsZpVxCFSvk+4gh4CZVcIhUD13B5lZYzgEKiQxBViPtEjcWDwwbGaVcQhUawbw\nUARPtjnGLQEzq4xDoFrjdQWBB4bNrEIOgWp1GgIzKqjFzMwhULHxZgZBCgGPCZhZJRwC1eq0JeAQ\nMLNKOASq5RAws0bJGgKSTpB0u6Q/tDnm65JWSLpK0m4562kAdweZWaPkbgmcCCwc64eS9ge2i4jt\ngfcDx2aup26dtAQeAaZLTK2gHjMbcFlDICIuBu5tc8gBwEnFsZcBsyVtlrOmmo0bAhEE8BCeIWRm\nFah7TGAusLLl/ipgy5pqqUIn3UHgLiEzq0gTuhw04n6MepC0uOXuUEQM5Sooo5m0bxkNcwiYWdck\nLQAWdPOcukPgVmBey/0ti8fWEhGLqygos5ms2fIZi0PAzLpWfDgeGr4vadF4z6m7O+gM4F0AkvYG\n7ouI2+stKatuuoM8JmBm2WVtCUg6GdgXmCNpJbAImAYQEcdFxFmS9pd0PWkw9LCc9TRAJ7ODwC0B\nM6vIuCEg6QDgPyOi3cqXo4qIQzo45shuzzuJOQTMrFE66Q46GLhe0hclPS93QX3Os4PMrFHGDYGI\neDuwG3Aj8D1Jv5H0fkle7rh7bgmYWaN0NDAcEfcDpwGnAlsAbwSWSvpQxtr6Uach8BAOATOrwLgh\nIOlAST8lTTuaBuwZEa8FdgY+mre8vuPuIDNrlE5mBx0EfCUiLmp9MCIelvTePGX1rW66g/p5+Qwz\na4hOuoNuHxkAkr4AEBHnZ6mqD0lMI4Xu3zs43C0BM6tEJyHwqlEe27/sQgbATODBYoG48TgEzKwS\nY3YHSToC+Edg2xH7AcwCLs1dWB/qtCsIHAJmVpF2YwI/As4G/i/wCZ5e6O2BiLg7d2F9yCFgZo3T\nLgQiIm6W9AFGrOwpaeOIuCdvaX2n05lB4BAws4q0C4GTgX8Afsfoyztvk6Wi/uWWgJk1zpghEBH/\nUHzfurJq+lu3IeBVRM0su3YDw7u3e2JE/L78cvqau4PMrHHadQcdwxi7fBVeXnIt/a6blsBDwEwJ\ndTil1MxsQtp1By2osI5B0HEIRLBa4jFgXeCRrFWZ2UBr1x20X0RcIOlNjNIiiIjTs1bWf7rpDoKn\nu4QcAmaWTbvuoH2BC4DXM3q3kEOgOzOBbrbOHF5J9M485ZiZte8OWlR8P7SyavrbTOCGLo734LCZ\nZdfJUtJzJH1D0lJJv5f0NUmbVFFcn5lod5CZWTadLCB3CnAHaUnpN5O6J07t5OSSFkpaLmmFpE+M\n8vM5ks6RdKWkP0o6tIvaJ5tuZgeBQ8DMKtBJCDwrIo6OiJsi4saI+Bc6WOte0hTgm8BC4PnAIZJ2\nHHHYkcDSiNgVWAB8WVInexxMRg4BM2ucTkLgPEmHSFqn+DoYOK+D5+0FXB8RN0fE46QWxYEjjrkN\n2KC4vQFwd0Ss7rT4SWYm7g4ys4ZpN0X0QZ6eFfQR4AfF7XVIM1c+Ns655wIrW+6vAuaPOOZ44AJJ\nfyH1mb+1s7InpVm4JWBmDdNudlCvb0CdXOn6aeDKiFggaVvgl5J2iYi1PjFLWtxydygihnqsr2ru\nDjKzrCQtIHWtd6yj/ndJGwHbk65gBWDklpOjuBWY13J/Hqk10OrFwOeK890g6SbgucCSkSeLiMWd\n1NpgDgEzy6r4cDw0fF/SovGe08kU0fcBF5HGAT4LnAss7qCeJcD2kraWNB04GDhjxDHLgVcWr7MZ\nKQBu7ODck4qE8BRRM2ugTgaGP0wa5L05Il4O7AbcP96TigHeI0mhcS1wakQsk3S4pMOLw/4P8EJJ\nVwHnA/+jTzermQ5EBI918RwvJ21m2XXSHfT3iHhEEpLWjYjlkp7byckj4mzSFpWtjx3Xcvsu0rIU\n/a7briBwS8DMKtBJCKwsxgR+Rhq4vRe4OWtV/afbriBwCJhZBcYNgYh4Y3FzsaQh0nz+c3IW1Yfc\nEjCzRup0dtAewEtI0z4viYhu+rbNIWBmDdXJ7KDPAN8DNgbmACdK+ufMdfWbiXQHDS8lbWaWTSct\ngXcAO0fE3wEkfR64Cjg6Z2F9xi0BM2ukTqaI3gqs13J/Xda+6MvacwiYWSO1WzvoG8XN+4FrJA0v\nGvcq4PLchfUZzw4ys0Zq1x30O9JA8BLS9NDhtYCG6GxdIHvaRFoCfwemSUyL4PEMNZmZtV1A7nvD\ntyU9A9ihuLu8WBraOtd1CEQQ0lNXDd+XpSozG3jjDgwXq9KdBNxSPLSVpP8WEb/OWVifmcWay2p3\narhLyCFgZll0MjvoGODVEXEdgKQdSBvE7J6zsD4zke4g8LiAmWXWyeygqcMBABARf6LDi8zsKQ4B\nM2ukTt7MfyfpO8B/AALezijr/VtbE5kdBA4BM8uskxD476QloT9U3L8Y+Fa2ivpTLy0BLydtZtm0\nDQFJU4GrIuJ5wJerKakvuTvIzBqp7ZhAsTHMdZKeXVE9/crdQWbWSJ10B21MumL4ctKiZgAREQfk\nK6vvuCVgZo3USQj8r+K7Wh7zFcPdmWgIeCVRM8uq3dpB65EGhbcDrgZO8JXC3Ss2me+lJTCr3IrM\nzJ7WbkzgJGAPUgDsD3yp25NLWihpuaQVkj4xxjELJC2V9Mdi57J+sz7wWASrJ/BcdweZWVbtuoN2\njIidACR9F7iimxNLmgJ8E3glaTnqKySdERHLWo6ZDfwb8JqIWCVpTre/wCQwC/jbBJ/rEDCzrNq1\nBJ765FrMEurWXsD1EXFz0Y10CnDgiGPeBvwkIlYVr3PXBF6n6SY6MwgcAmaWWbsQ2FnSA8NfwE4t\n9zv5ZDuXNRdNW1U81mp7YGNJF0paIumd3ZU/KWyAQ8DMGqrdUtJTejx3JzOIppEWonsFqe/8N5J+\nGxErRh4oaXHL3aGIGOqxvqq4O8jMKlGs+rygm+fkXAjuVmBey/15rL0t5Urgroh4BHhE0kXALsBa\nIRARizPVmZu7g8ysEsWH46Hh+5IWjfecTlYRnaglwPaStpY0HTgYOGPEMT8HXiJpiqT1gfnAtRlr\nqoO7g8yssbK1BCJitaQjgXOBKcB3I2KZpMOLnx8XEcslnUOahvokcHxE9FsIuDvIzBpLEc2/+FdS\nRITGP7J5JD4ObBbBURN47gbAqgg2KL8yM+t3nbx35uwOsqSX7qCHgBnFVcdmZqVzCOQ34e6gCJ4A\nHiXNnDIzK51DIL9eZgdRPNfrB5lZFg6B/HrpDoLUivCYgJll4RDIr5fZQeAQMLOMHAL59doddD+w\nYUm1VE5i3bprMLOxOQTyG8juIIl1JI4DHpD4icT0umsys7U5BPLrtSUwKUMAOJS0H8WmpAvePl5r\nNWY2KodAfr2OCUy67iCJqcBngA9EcB9wBPBPkmc5mTWNQyCj4iKvQWwJvBq4LYLLACK4Efg1af0o\nM2sQh0BezwCI4NEeznE/ky8E3kLaRKjVqcCba6jFzNpwCOTVa1cQxfMnTXdQMQB8AHD6iB+dBbxY\nYuPqqzKzsTgE8up1ZhBMvu6gfYEVEWvsKkcEDwLnA6+rpSozG5VDIK8Ngft6PMdk6w7aj7R8+GjO\nK35uZg3hEMhrNr2HwKTqDiJtbTc0xs+G6HLrOzPLyyGQV1khMClaAhIzgZ2A345xyHXAuhJbV1aU\nmbXlEMhrQ1J3Ti8mU3fQPsDvInhktB9GEKTWwL5VFmVmY3MI5DVQLQFSCFwyzjGXFMeZWQM4BPIq\nIwQeADaYJLuL7QXpArE2rgBeWEEtZtaBrCEgaaGk5ZJWSPpEm+P2lLRa0kE566lBz91BETzOJNhd\nrAipPUlv8u1cDTzPq4uaNUO2EJA0BfgmsBB4PnCIpB3HOO4LwDkwKT7tdqOMlgBMji6h5wAPR3Bb\nu4OK8YLrgF0qqcrM2srZEtgLuD4ibo6Ix0nLCBw4ynEfBE4D7sxYS13KDIGmTxPdC7i8w2PdJWTW\nEDlDYC6scdXoquKxp0iaSwqGY4uHImM9dSjjYjGAe0mB0mQvZPyuoGFLcAiYNcLUjOfu5A39q8An\nIyIkiTbdQZIWt9wdioih3sqrxGx6nyIKKQSavubOrsC/dnjsFcCHMtZiNpAkLaDLCzJzhsCtwLyW\n+/NIrYFWewCnpPd/5gCvlfR4RJwx8mQRsThTnTmV1R10L7BRCefJohgU3gW4qsOnXANsIzGzWFPI\nzEpQfDgeGr4vadF4z8nZHbQE2F7S1pKmk9aSX+PNPSKeExHbRMQ2pHGBI0YLgEmsrO6ge2h2S2Au\n8ATw104OjuAxUhDsmrMoMxtfthCIiNXAkaTFxK4FTo2IZZIOl3R4rtdtCol1SCHQ61LS0PCWAEUr\noLgiuFNLcQiY1S5ndxARcTZw9ojHjhvj2MNy1lKDmaQpk6tLONc9wFYlnCeXbrqChi0lXVdgZjXy\nFcP5lNUVBM0fGJ5ICFyJWwJmtXMI5FPWoDCklkCTu4N2pfsQuBrYUWJahnrMrEMOgXzKmh4KDW4J\nSMwgzfxa3s3zIngYuAVY6ypyM6uOQyCfTYC7SjpXkweG/wuwvFjjqFtLgd1KrsfMuuAQyGcOcHdJ\n52ryFNGJjAcM87iAWc0cAvnMYTBaAhMZDxjmloBZzRwC+WxCSS2BYuXNkFivjPOVrJeWwFXArpNk\nrwSzvuQQyKfMlgA0cHC4uCBuJyYYAhHcATwE3nPYrC4OgXxKawkUmtgltA1wXwT39HAOXzlsViOH\nQD5ltwSaODjcy3jAsCvxuIBZbRwC+ZQdAneTWhdNsgvpTbwXbgmY1cghkE/Z3UF3AM8s8Xxl2JXe\nQ8AtAbMaOQQykJhCWjvo3hJP29QQWNrjOW4CNpCYU0I9ZtYlh0AeGwH3R/BEiedsVAhIbEIKupt7\nOU8ET1JMFS2hLDPrkkMgj7LHAwBup0EhQDEoXLyJ98rjAmY1cQjkkSME7gA2K/mcvShjPGCYxwXM\nauIQyKPsQWFoWHcQ5YaAWwJmNXEI5PFM4M6Sz9nPIXAtaeP59Us6n5l1KHsISFooabmkFZI+McrP\n3y7pKklXS7pU0s65a6rAFsCtJZ/zbmBDKe+WoJ2QWBfYjvTm3bNi4/nrSMtSm1mFsoaApCnAN4GF\nwPOBQySN3ETkRuBlEbEzcDTw7zlrqsgWwG1lnrAYgL0b2LTM807QC4AVEfy9xHN6XMCsBrlbAnsB\n10fEzRHxOHAKcGDrARHxm4gY3oHrMmDLzDVVYQvgLxnO25QuoTK7goZ5bwGzGuQOgbnAypb7q4rH\nxvIe4KysFVVjc/KEQFOmie5O+SHgwWGzGuTuX45OD5T0cuDdwD75yqlM6d1BhaZME50P/Kjkc14F\n7CQxpeSL7MysjdwhcCtpE/Jh80itgTUUg8HHAwsjYtSlFiQtbrk7FBFD5ZVZnmLJiE1Jn9rLdhsp\nYGpTbGzzfOD3ZZ43gvslbge2p8tN680skbQAWNDNc3KHwBJge0lbk7pHDgYOaT1A0lbA6cA7IuL6\nsU4UEYuzVVmuZwJ3T3Dj9fH8mfQmWafdgGXFbmdlGx4cdgiYTUDx4Xho+L6kReM9J+uYQESsBo4E\nziVNJzw1IpZJOlzS4cVhnyGttXOspKWSLs9ZUwVyDQpDCoGtMp27U/NJA/g5eFzArGLZ55xHxNnA\n2SMeO67l9nuB9+auo0K5xgOgOSFw9rhHTcyVwIczndvMRuErhsuXa2YQpBB4dqZzd2pv8rUEfgvM\nb8IFcWaDwiFQvq0YZfC7JHcD60rMynT+tiQ2Iy0f/acc54/gLlLQuUvIrCIOgfJtC9yQ48QRBOlN\nct54x2byYuCykpaPHstFwMsynt/MWjgEyrcdMOYspxLUOS7wcuDCzK9xEbBv5tfomcQciQMl3iWx\nr8Qz6q7JbCIcAuXL1hIo3EJ9IbCAlulnmVwEvFRq5t+mxKYS3wFWAEcArwK+CNwi8dHiOhGzSaOR\n/9AmK4mNgSmUv6FMqz8DW2c8/6iKPYCfDfwu5+tE8BfS2McLcr7OREi8GLgauA/YNoKFEbwzgvnA\nK4DXAUPeL9kmE4dAubYFbij67nNZDjwv4/nHsgC4NILVFbzWRaSup8aQeAPwc+CwCI6K4J7Wn0dw\nDfBK4BLgUoltaijTrGsOgXLlHg+AdNFdHZ+S9wfOqei1zgFeW9FrjUviJaQlzhdGjP3fIIInI/gU\n8C3gfInNq6rRbKIcAuXKPR4AqS96q2Jjl0oU/fP7A7+o6CXPA14iMaOi1xuTxLbAacA7IzrrCovg\na8CJwLkSG+asz6xXDoFyPY9Mc+iHFbtw3QTskPN1RtgduDcie8ABaTE54ApgvypebywSG5GC77MR\nnNvl0z9H6hr6kQeLrckcAuXahTRwmNs1VNsl9AbgPyt8PYrXe0PFr/kUiWmkFsDZERzb7fOLcaEP\nA+sBny+5PLPSOARKUswT3470Bp1bZSEgIdLKr6dU8Xotfgy8sY7598Xv/C3gYeCoiZ6nWEn2LcCb\nJN5ZUnlmpXIIlGcn4PqS990dy9VUtx/vXsBqSt4/YDwRrCL9nvtX+bqFj5F+77f1usFNBHcDBwDH\nSMwvozizMjkEyvMi4DcVvdYlwD4V9TW/C/hh5mmvY/lh8fqVKaaCfgR4XQQPlHHOYvrou4HTpdpX\ngTVbg0OgPJWFQAR3kFYqzbrQmsRsUlfQd3K+ThsnAy+Tqrk4TmJv0g53b4hYY2/snkVwJvBl4Bee\nMWRN4hAoQTGFcj/yL6nQ6tfkX2PnPcA5xVW8lYvgQdJUyyNzv5bEc4GfAYdGsCTTy3yF9P/tNInp\nmV7DrCsOgXLsAvwtgpsqfM0hMl5VKzGT1Df+5Vyv0aGvA4dJPCvXC0jMI22U8+mIfNdCtMwYehA4\n1UFgTeAQKMebSZ8iq3QOaaG1XOvUfBT4dacXSOUSwZ+B7wGfzXH+YnmHi4BvRnBCjtdoVQw0HwwI\n+H9efdTq5hDoUbEL1ttJ/deVKS6oOrN47VJJvAD4IPDpss89QZ8DXi+V2/1VzNa5GPhSBMeUee52\nigv+3kqadXWOF5yzOmUNAUkLJS2XtELSJ8Y45uvFz6+SVNW0xzIdBNwWwdIaXvsE4PAyt2OU2IAU\naJ+suHtrTMVibe8BfiCxRa/nk1hH4h9JIXpEBP/W6zm71RIElwGXS+xTdQ0AEhtK7CdxlMS3Jc6T\nWCaxSuI+iQckbpNYIXGJxPclFkscJLFVcU2FTWKKyDPzT9IU4DrSyoq3kpYBOCQilrUcsz9wZETs\nL2k+8LWI2HuUc0VENO6PTWJ94CrgyAksK9ByHi2IiKEJvL6A84EzI/jqRF+/5XyzSStlXgN8oKxp\noRP9/dY+D58C3gG8priOYCLn2I00QLsuaRB4eW819f67FdNSvwX8FDg6gr/2cr42rzONNKNsb2A+\nsCcwF7iStET4dcCNwErgfuBB2O7FcP1SYBbwLOA5pIsidymePwVYAlxO2iP68uLaiEmhrL/Npurk\nvTPnht57AddHxM1FMacABwLLWo45ADgJICIukzRb0mYRcXvGukpRfPr+LumPfsIBUFjABGYWRRAS\nR5CWLl4awa8nWoDEfsBxwFnAP5V8XcACSpg5FcHnJZ4ArpA4Cji5k60ui/9XrwHeT3rj+hzw7V4v\nBCssoMffLYKfSVwEfAa4VuJU0r+Lyye6lWdxDck2wM6kN/0XkS4wvIk0lflC0mY417ZbHly6Ya8I\nzgJuI62LdVHLawjYgvTfdC/S1dV7StxOauH8tvh+VdHyaaIFVDurr3FyhsBcWGOu9SpY64rJ0Y7Z\nEmhsCBSrd74CWETa/OSgOuuJ4E8SB5OmHR4PfLeThd6Kzep3IO3n+xbSp7yPRfDTrAX3KIIvSlwK\nHAP8i8RPSW80t5A+vU4DZpA2wNmOtC/yPqQluE8EDong4Tpqb6fo8vqIxBeBw0i1bipxMbCU9An9\nL8BDxZeA9Yuv2aR/S3NJ+0/vSFrM8E7gj6RP6f+bFCr3l1hzkFr5t1JMjCjCZ0fSv/X5pODdXuIW\nUoisIK20+1fgDtK/9btJS3Q82smHjyJ8ppJac+sV30d+PQOYPs4XcMCLJD4CPFF8PVl8fwx4vPg+\n0dura7rIsis5Q6DTX35kU2XU50n8ojh2+Hh1cX8izxntHBsBm5Oazl8FTsm86XpHIrhAYg/gU8Al\nRTfVStI/rMdIf9TrtXxtVHy/kTQw+nnSQmlVbBjTswguLS7s2p20vePBpDfA2aTf9xHSDmw3kFpr\n74rIuttbaYprMj4HfE5iLimkX0BaPmNzUsDNIL1ZPUz6Xf9GeiNeBfwS+AbpE34pVzx3Wf8TpOD5\nI+m/PRLrkbqRdgC2J7VOXgU8s/iaQ/p7nC7xCOl3Ctb+9ydSyK9L+v3/3vL1SMvtR4vvj7X5ejxV\nPHMDUotpCmmMdErxNY0UFMPfu709DZgqrREO7QLj8Zb/jOO9d3X6PjdEB3KOCewNLI6IhcX9TwFP\nRsQXWo75NjAUEacU95cD+47sDpLU+DQ1M2uiOscElgDbS9qa1Iw9mLQEQaszSFeDnlKExn2jjQc0\ncVDYzKwfZAuBiFgt6UjgXFLz6rsRsUzS4cXPj4uIsyTtL+l6Uj/nYbnqMTOztWXrDjIzs+Zr9BXD\nnVxsNplJOkHS7ZL+UHctZZM0T9KFkq6R9EdJH6q7pjJJWlfSZZKulHStpL7cPUzSFElLJZ1Zdy1l\nk3SzpKuL3+/yuuspUzHd/jRJy4q/z7Wuv3rq2Ka2BDq52Gyyk/RS0mJi34+Inequp0ySngU8KyKu\nlDSTNKPqDX32/2/9iHhY0lTSHg9HRcQldddVJkkfBfYAZkXEAXXXUyZJNwF7RMQ9dddSNkknAb+O\niBOKv88dnoBvAAAEEElEQVQZETHq9OAmtwSeutgsIh4nbW94YM01lSoiLgburbuOHCLirxFxZXH7\nQdJFgj0v+dAkETF8vcF00rhXX72ZSNqSNDX1O6w9lbtf9N3vJWlD4KURcQKk8dmxAgCaHQKjXUg2\nt6ZarAfFDLHdSBd19Q1J60i6knTB04URcW3dNZXsK8DHof5rYTIJ4HxJSyS9r+5iSrQNcKekEyX9\nXtLxktYf6+Amh0Az+6msK0VX0GnAh4sWQd+IiCcjYlfSVe4vk7Sg5pJKI+l1wB0RsZQ+/LRc2Cci\ndgNeC3yg6J7tB1NJF1J+KyJ2J828/ORYBzc5BG4lXQI/bB5MbNEwq4ekacBPgP+IiKr3W6hM0dT+\nBfDCumsp0YuBA4p+85OB/SR9v+aaShURtxXf7yQt3rdXvRWVZhWwKiKuKO6fRgqFUTU5BJ662EzS\ndNLFZmfUXJN1SJJIywZcGxE9r3DaNJLmSJpd3F6PtAxCHcuJZxERn46IeRGxDfBfgQsi4l1111UW\nSetLmlXcngG8GuiLWXoR8VdgpaQdiodeSVoZeFQ5rxjuyVgXm9VcVqkknUzaJ3gTSSuBz0TEiTWX\nVZZ9SMs+Xy1p+M3xUxFxTo01lWlz4CRJ65A+TP0gIn5Vc0059Vv37GbAT9NnFaYCP4yI8+otqVQf\nBH5YfIC+gTYX4jZ2iqiZmeXX5O4gMzPLzCFgZjbAHAJmZgPMIWBmNsAcAmZmA8whYGY2wBp7nYBZ\nbpKeAK5ueejkiPhiXfWY1cHXCdjAkvRARMwq+ZxTI2J1mec0y8ndQWYjFJuNLJb0u2LTkecWj88o\nNgK6rFid8YDi8UMlnSHpV8AvJa0n6cfFhjqnS/qtpD0kHSbpKy2v8z5Jx9T0a5oBDgEbbOsVu0oN\nf72leDyAOyNiD+BY4Kji8f8J/Coi5gP7Af/askTvbsCbIuLlwAeAuyPiBcA/kzZlCeDHwOuLDZMA\nDiWtr2RWG48J2CB7pFhKeDSnF99/DxxU3H416U18OBSeAWxFeoP/ZUTcVzy+D/BVgIi4RtLVxe2H\nJF1QnGM5MC0ixlzYy6wKDgGz0T1afH+CNf+dHBQRK1oPlDSftGb7Gg+Pcd7vkFoUy4ATSqjTrCfu\nDjLr3LnAh4bvSBpuRYx8w78UeGtxzPOBp/aPjojLSZvQvI20Tr9ZrdwSsEG2Xssy1wBnR8SnRxwT\nPL2M8tHAV4vunXWAG4EDRhwD8C3SMtPXAMtJa7m37vH6Y2CXdvu+mlXFU0TNSlbsMTAtIh6VtC3w\nS2CH4amjks4EjomIC+us0wzcEjDLYQZwQbG9poAjik2SZgOXAVc6AKwp3BIwMxtgHhg2MxtgDgEz\nswHmEDAzG2AOATOzAeYQMDMbYA4BM7MB9v8B+JufWYBp2mAAAAAASUVORK5CYII=\n",
	"text": [
	"<matplotlib.figure.Figure at 0x1098c5b50>"
	]
	}
	],
	"prompt_number": 7
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We've been a little careless about the n=0 case, which is quite likely since $\\lambda$ is small. That looks like a $\\delta$ function spike at Energy=0 with total probability $1/e$"
	]
	}
	],
	"metadata": {}
	}
	]
	}