Added supporting notebooks subfolder.

Wrote a notebook explaining how to plot PDFs and run monte
carlo simulations and compute their results.

Supporting_Notebooks/Computing_and_plotting_PDFs.ipynb

@@ -0,0 +1,508 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Computing and plotting PDFs of discrete data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "So let's investigate how to compute and plot probability distributions.\n",
+    "\n",
+    "\n",
+    "First, let's make some data according to a normal distribution. We use `numpy.random.normal` for this. The parameters are not well named. `loc` is the mean of the distribution, and `scale` is the standard deviation. We can call this function to create an arbitrary number of data points that are distributed according to that mean and std."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "50000\n",
+      "3.00761397\n",
+      "1.99265941012\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import numpy.random as random\n",
+    "\n",
+    "\n",
+    "mean = 3\n",
+    "std = 2\n",
+    "\n",
+    "data = random.normal(loc=mean, scale=std, size=50000)\n",
+    "print(len(data))\n",
+    "print(data.mean())\n",
+    "print(data.std())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As you can see from the print statements we got 5000 points that have a mean very close to 3, and a standard deviation close to 2.\n",
+    "\n",
+    "We can plot this Gaussian by using `scipy.stats.norm` to create a frozen function that we will then use to compute the pdf (probability distribution function) of the Gaussian."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 64,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGQ1JREFUeJzt3XmQXXWZxvFvdzqdkIVAFraQhYgSAooEDFEWm7AkwTEs\nonEdXAYpHa2aGRwxs0hXSdVYWFM15Vgzo6Il1KjBkS2yBBKghQECBklAhgABAiQEEhL2GTVI5o/3\nNLnd6Zu7nXve857zfKpu9XLv6fvQ3H5z7nt+C4iIiIiIiIiIiIiIiIiIiIiIiIiISIoWAOuAx4GL\nhrj/U8Ba4EHgLuA9DRwrIiI5MAxYD0wHhgNrgMMHPeb9wLjk8wXAqgaOFRGRDHTWuH8OVrA3ADuA\npcCZgx5zD/BK8vm9wMENHCsiIhmoVewnA89WfL0x+V41XwBubPJYERFpk64a9+9s4GedDHweOL6J\nY0VEpI1qFftNwJSKr6dgZ+iDvQf4Idazf6nBY/WPgohIczrS+kFdwBPYRdZuhr7IOhXrzc9t4lhQ\nsU9br3eAJnwM2AL8IzCh4vsjgE8Dm4Fvk+ILuwG9Ds9ZVL3eAQom9dq5EHgUK+hLku9dkNwALgO2\nAQ8kt/tqHDuYin26er0DNOiz2Du+o/fwmAnYQIAfYKO8stSb8fMVWa93gIIJVzvDBc65Xu8ADTgH\na/fNrOOxY4E7gEvbmmh3vRk/X5H1egcomHC1M1zgnOvxDlCnQ4GtwDENHDMBG8p7bjsCVdGT4XMV\nXY93gIIJVzvDBZaWdQP3A19p4thjsH8kpqWaSCSecLUzXGBp2RLgBpq/4PoPwK9aOF6kCMLVznCB\npSXTsQv6h7TwM7qBh8m2nSOSN+FqZ7jA0pKrsSGWrToReBoYmcLPEokoXO0MF1iaNgdbQiOtAr0M\n+KuUfpZINOFqZ7jA0rSb2TU/Iw3vBl4A9k7xZ4pEEa52hgssTTkReBLrt6fpP4G/T/lnikQQrnaG\nCyxNuQlbFTVtRwLPo969lE+42hkusDTsCGx9mxFt+vnXA19s088WyatwtTNcYGnYZdjY+HY5CXiM\n7NfNEfEUrnaGCywN2Q/YzsDVLNPWAawGzmjjc4jkTbjaGS6wNOQb2Jl9u30Bm1UrUhbhame4wFK3\nTmxPg/dl8FyjsZm5WjNHyqKh2llrD1qRVpyKbUa/OoPnegMbhqkLtSI5pTP74rqKdCdR1TILeI7a\n222KFEG42hkusNRlIvAy2c9uvQ/bC1mk6NTGkVxYDNwIvJrx814OnJfxc4pIHXRmX0z34nOGPQG7\nTrCPw3OLZClc7QwXWGo6DJsx69U7/yVwvtNzi2RFbRxx92ng58CbTs//M+DjTs8tIlXozL5YOrCl\nC7IYW1/NXsBLwIGOGUTaTWf24uoorH2Txdj6av4PWxxN2xaKJFTsJW2LgV/g/45taZJFRHLCuyhI\nejqw5RFmewfBNknZBkzxDiLSJmrjiJujsRfgA95BgD8CNwCLvIOI5IGKvaTpLOAa8vNu7TrgTO8Q\nImLyUhikdWuAE7xDVBgDvIYmWEkxhaud4QLLkKYBW8nfblHXA5/wDiHSBurZi4tFWI/8T95BBlEr\nRyQndGZfDCuAs71DDOEAbIJVuzY7F/ESrnaGCyy7GYetbjnGO0gVdwPzvUOIpExtHMncQuBO4HXv\nIFWolSOSAzqzj+9n5Hs7wJnARmzSl0hRhKud4QLLAMOB7cBB3kFqWAcc6x1CJEVq40imTgLWY3u/\n5tl1aDatlJiKvbTqw8Ay7xB1uAntTSviSm2c2NYBx3iHqEM3tl3hRO8gIikJVzvDBZa3TQe2EOcd\n4rVoNq0Uh3r2kpn5wC3AW95B6nQzauWIuNGZfVxXA5/xDtGAGcDz6CRHiiFc7QwXWAAbcvkSsL93\nkAY9iq27LxKd2jiSieOAp4AXvIM0aDlaOkFKSMVemjUf64FHo769iBO1cWL6DdDjHaIJo7ENTfb2\nDiLSonC1M1xgYSI2Zr3bO0iTbsG2UBSJTD17abvTgD5sU++IlqNWjshuFmCzJB8HLhri/pnAPcDv\ngQsH3bcBeBB4ALivys/XmX08PwG+7B2iBUdgF5dFIku1dg7DFrmajg21WwMcPugxk7DVBC9h92L/\nFDC+xnOo2MfSgS169g7vIC3o/2+Y4R1EpAWptnHmYMV+A7ADWMrum0BsBVYn9w9Fa4gXy0zgD8AT\n3kFasBO4DTjFO4hIVmoV+8nAsxVfb0y+V6+dwErsH4PzG4smOdWD9eujuxU41TuESFa6atzfaovl\neGAz1upZgfX+7xzicb0Vn/dRjGJSVCcD13uHSMGtwKXYCU+UtX2k3Hpo43DnudjIhX5LGPoiLcDF\n7N6zr+d+9ezj6MBWuZzqHSQljwFHeYcQaVKqPfvVwDuxC7TdwGKqb1QxuDc/ChibfD4aOB14qJFw\nkjuzsAlJz3gHScmtqG8v8raF2OJR67Eze4ALkhvAAVhf/xVsYaxngDHYSIc1ye13FccOpjP7OL4C\nXOYdIkXnAjd4hxBpUrjaGS5wiV0FfMo7RIomYCcpw72DiDQhXO0MF7ikOoEXaWw0VgS/BU7wDiHS\nBC2XIG1xJLAd2OQdJGUrUd9eSkDFXup1MnC7d4g20EVakYyojRNDUTfrHg28nnwUiSRc7QwXuIQ6\nsRbOgd5B2uTXaBVMiUc9e0ndUdj2g5u9g7RJH/BB7xAi7aRiL/U4mWIvYdFHzF23REJRGyf/lgEf\n8w7RRnthffsx3kFEGhCudoYLXDLDsJnR+3kHabM+bBN1kSjUs5dUHY1t9LHFO0ib9aFWjhSYir3U\n0kMxx9cP1oddmxCRNlEbJ99uAD7iHSIDI7G+/dhaDxTJiXC1M1zgEunCFgqb6B0kI31ovL3EoZ69\npGY28DS2AFoZ9KG+vRSUir3sSVHXw6mmDxV7kbZRGye/lgNne4fIkPr2Ekm42hkucEkMB14FxnsH\nyVgf6ttLDOrZSyqOBZ7AFkArkz7UypECUrGXaoq+Hk41t6NiL9IWauPk0y3AIu8QDkYCr6G+veRf\nuNoZLnAJdGMFbx/vIE5uR317yT/17KVlc4DHgJe9gzjpQ60cKRgVexlKD+UaXz9YH1onRyR1auPk\nz63Ah7xDOOofb7+3dxCRPQhXO8MFLrgRWL9+nHcQZ7cDC71DiOyBevbSkuOAR7AF0MqsD/XtpUBU\n7GWwsq2HU00fKvYiqVIbJ1/6UPsC1LeX/AtXO8MFLrD+jbc1ocioby95pp69NG0u8BB2gVa0dIIU\niIq9VCrrejjV9KFiL5IatXHy407gdO8QOTIC9e0lv8LVznCBC2oUVthGewfJmdso9wQzyS/17KUp\nHwDWAm94B8kZ9e2lEFTspZ/G1w+tD62TI5IKtXHy4S7gFO8QOdS/fERZl3uW/ApXO8MFLqAxWL9+\nlHeQnFoBfNg7hMgg6tlLw44Hfgv8r3eQnOpDrRwJTsVeQOvX16KLtCIpUBvH3yp05ronw4FXgfHe\nQUQqhKud4QIXzFisXz/SO0jO3Qyc5R1CpIJ69tKQE4DfAL/3DpJzauVIaCr2ovVw6nM7anWJtERt\nHF+/AU7yDhHAcGz3roneQUQS4WpnuMAFMg6bMDTCO0gQNwLneIcQSahnL3U7EbgX+IN3kCD6UCtH\nglKxL7d5qF/fCPXtpdAWAOuAx4GLhrh/JnAPNprjwgaPBbVxPK0B3u8dIpAu4GVgP+8gIqRcO4cB\n64Hp2AWqNcDhgx4zCTgWuISBxb6eY1MPLHWbiF1w7PIOEsz1wEe9Q4iQcs9+DlawNwA7gKXAmYMe\nsxVYndzf6LHipwfbmepN5xzRqJUjIdUq9pOBZyu+3ph8rx6tHCvtdwq2C5M0RpOrJKRab+FbabE0\ncmxvxed96KJhFuYB3/cOEdBa4EDgAOB55yxSLj20cKJRq9hvAqZUfD0FO0OvRyPH9tb5MyUdBwMT\ngAe9gwT0J+AO7I9uqW8UKZk+Bp4IX9zIwbXaOKuBd2IXWbuBxcCyKo/taOFYydY8rB3xlneQoNTK\nkUJaCDyKXWxdknzvguQG9nb2WWxkx0vAM9jOR9WOHUyjcbL3E+BL3iECey/2uhbxFK52hgscXAf2\nD/Jh3kEC6wS2AQd5B5FS03IJskeHYv/fH/MOEthbwK9RK0cCUbEvn3nYkEu9o2pNHxpvL9IQFZ1s\nXQl81jtEARwBPOkdQkotXO0MFziwTmALMNU7SAF0AJuBGd5BpLTUs5eqjsRGTT3jHaQAdgIrgVO9\ng4jUQ8W+XPr79ZKOlcBp3iFEolAbJzvLsMltko7J2BDMYd5BpJTC1c5wgYPSWuzt8T/AMd4hpJTU\ns5chHYP16rd4BykY9e0lBBX78tCSxu2hYi9SJ7VxstEHnOEdooDGAa8Be3kHkdIJVzvDBQ5oLFaQ\nRnsHKai70Nm9ZE89e9lND3Af8IZzjqJSK0dyT8W+HOYDt3iHKDAVe5E6qI3Tfo8BR3uHKLBubGby\nRO8gUirhame4wMEcAryA3sW123XAJ7xDSKmoZy8DnA6sQFsQtttyYIF3CJE805l9e10FfMY7RAnM\nAJ5HJ1CSnXC1M1zgQLqwfYEP9A5SEo+iayOSHbVx5G1zsCUSNnsHKYmbUCtHckrFvtjmAzd7hyiR\n5cBC7xAieaU2TvusQuO/s7QXNlN5nHcQKYVwtTNc4CAmYWO/R3gHKZnlwDneIaQU1LMXwHrHtwF/\n8A5SMhqCKVKFzuzbYynwF94hSmgmdlG8wzuIFF642hkucADDsSGXB3kHKaEO4ClglncQKTy1cYQP\nAE8Cz3kHKaGdWCtHewdIrqjYF9OHgBu8Q5TYr4APe4cQyRu1cdL3MDahSnyMxEZCTLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x85bfd30>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy.stats as stats\n",
+    "\n",
+    "def plot_normal(xs, mean, std, **kwargs):\n",
+    "    norm = stats.norm(mean, std)\n",
+    "    plt.plot(xs, norm.pdf(xs), **kwargs)\n",
+    "\n",
+    "xs = np.linspace(-5, 15, num=200)\n",
+    "plot_normal(xs, mean, std, color='k')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "But we really want to plot the PDF of the discrete data, not the idealized function.\n",
+    "\n",
+    "There are a couple of ways of doing that. First, we can take advantage of `matplotlib`'s `hist` method, which computes a histogram of a collection of data. Normally `hist` computes the number of points that fall in a bin, like so:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {
+    "collapsed": false,
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEQ9JREFUeJzt3W+MHGdhx/HvNbYhadQYl9Z2bEvnokTCahtCaRKq0Lht\n6rgIJeENBLVVaHhRKW1J1Ar/CRJc3yDXCIVKbVQJCr1CSQlpFCUVBDsVJ6RWTRDENMQ5bOdP4Ux9\nQUXQwhuCsn3xPNubW8/eeWf/PM/s8/1Io312dnbmsXfmN8888+weSJIkSZIkSZIkSZIkSZIkKWOf\nAJaBpyvztgDHgVPAMWBz5bXDwGlgEdhXmf8rcR2ngb8cY30lSUN6C3A1q4P/KHAglg8CR2J5D3AC\n2AjMAmeAmfjak8A1sfx5YP/YaixJGtosq4N/Edgay9vicwit/YOV5R4DrgO2A89W5t8G/M04KipJ\nWt9PNXjPVkL3D/GxexK4HFiqLLcE7KiZfzbOlyQl0CT4qzpxkiS1xIYG71kmdPGcI3TjvBTnnwV2\nVZbbSWjpn43l6vyzfdbtSUSSmplZf5ELN8v5N3e7ffmHOP/m7iZgN/BcpSJPANfG52vd3M0x+OdS\nV6DGXOoK1JhLXYE+5lJXoMZc6grUmEtdgRpzqStQYy51BfoYKDvXa/HfD9wAvBb4NvABQtA/ALwH\neBF4R1z2ZJx/EvgJcGelMncCfwdcTAj+xwappCRpdNYL/nf1mX9jn/kfilOvrwK/dKGVkiSNz7A3\nd0uwkLoCNRZSV6DGQuoK9LGQugI1FlJXoMZC6grUWEhdgRoLqSswjXLs45ek3A2Unbb4JakwBr8k\nFcbgl6TCGPySVBiDX5IKY/BLUmEMfkkqjMEvNeMv06q1DH5JKozBL0mFMfglqTAGvyQVxuCXpMIY\n/JJUGINfkgpj8EtSYQx+SSqMwS9JhTH4JakwBr8kFcbgl6TCGPzS8PylTrXKhtQVkDJWDfOZZLWQ\nRswWv7QmG/KaPrb4pdVMek09W/zSecx+TTeDX5IKY/BLF8aRO5oaBr+0Yo1gN/M1PQx+SSqMo3pU\nsnGM03fsv7Jni1+FG7oLp2YFdgspbwa/NJiem7yd3tek7Bn80sD65bu5r3YYJvgPA88ATwOfAV4F\nbAGOA6eAY8DmnuVPA4vAviG2K0lKYBZ4nhD2AJ8FbgeOAgfivIPAkVjeA5wANsb3nqH+pGOTSZPU\ngU6366YydWrKdfMuZFlpIgba15q2+P8HeBm4hDAy6BLgO8DNwHxcZh64NZZvAe6P73mREPzXNNy2\nNIw+oWxGqxxNg/97wEeAbxEC//uELp6twHJcZjk+B7gcWKq8fwnY0XDb0pAmGvK2/JWdpuP4Xwfc\nTei2+QHwOeD3epZZb4fv99pcpbwQJ6mlOjicX2OwN06NNA3+NwH/Bvx3fP4Q8GbgHLAtPm4HXoqv\nnwV2Vd6/M86rM9ewTpJUigVWN4o/OMibm3b1LALXARcTmjM3AieBRwk3eYmPD8fyI8BtwCZgN3AF\n8GTDbUuSEjnAynDOecKInS3A49QP57yHcFN3EbipzzrtC9U4jGDUTtNRPY7w0UQMtH/l1vloh6jG\nobN616orD/v6hSzrvq2xGSg7/eauJBXG4Jekwhj8klQYg19Kwxu+Ssbgl6TCGPzS5NjKVxYMfmli\nzHzlweCXpMIY/JJUGINfkgrT9Nc5pVxVO9L9iQSphi1+TSFvokprMfglqTAGvyQVxj5+TQP7dqQB\n2OLXlOib/TmeFHKskwpi8EsTZ+4rLYNfkgpj8EtSYQx+SSqMwa+2W6vD3M50qYbBrylm7kt1DH5J\nKozBL0mFMfil9Kp/ktE/z6ixM/il/Hky0Ej5Wz1qK4NQasgWv1psKrK/06csjY3BLyVl7mvyDH5J\nKozBL0mFMfglqTAGvyQVxuCXpMIY/FKe/NKWxmaY4N8MPAg8C5wErgW2AMeBU8CxuEzXYeA0sAjs\nG2K7UgHMfOVpHrgjljcAlwFHgQNx3kHgSCzvAU4AG4FZ4Az1Jx33dl2oDnS6reI+5fVeH2TZVOvq\n/ls9NrSmiewflwHP18xfBLbG8rb4HEJr/2BluceA62re786ttfSEYa5hbfBr4gbaP5p29ewGvgt8\nEvga8DHgpwmhvxyXWWblJHA5sFR5/xKwo+G2JUlDaBr8G4A3AvfFxx8Bh3qWWa+VYgtGWp/HiUau\n6a9zLsXpK/H5g4TunHOELp5zwHbgpfj6WWBX5f0747w6c5XyQpykQnWAmdSVUH72xmnivgxcGctz\nhBu7R1npyz/E+Td3NxG6iZ6jfm+2daO1dFZPufbLj2tdUl8T2z+uIrT4vw48RLjhuwV4nPrhnPcQ\nRvMsAjf1Wac7t9bSwrA2+DURA+0fuV1DdsivTspHZ+VhhtW7S115vdfbti6PDfU1UHb6zV1JKozB\nL0mFMfil9rCvXyNh8EutYeZrNAx+SSqMwS+1j01/DcXgl6TCNP3JBmmSbOFKI2SLX7nrrHqQNDSD\nX5IKY/BLUmEMfkkqjMEvSYUx+CWpMAa/1E7+bo8aM/ilVupUC54ANBCDX2o1M1+DM/iVI1ux0hgZ\n/JJUGH+rRzmxlS9NgC1+Zcbsl8bN4Jekwhj8ypnNf2kMDH5lzNyXxsGbu9J0qJ4lZ5LVQq1gi1+a\nGl4h6cIY/JJUGINfkgpjH79yYB+FNEG2+JWaf0xdmjCDX5IKY/BLUmHs41cq9u1IidjiV0Jmv5SC\nwS9JhbGrR5o+/nyD1jRsi/8i4Cng0fh8C3AcOAUcAzZXlj0MnAYWgX1DblfSmuxGU3/DBv9dwElW\n9rJDhOC/EviX+BxgD/DO+LgfuG8E25YkNTBM+O4E3gp8nJXLyZuB+VieB26N5VuA+4GXgReBM8A1\nQ2xb7dTBP6QuJTdM8N8LvA94pTJvK7Acy8vxOcDlwFJluSVgxxDbVmuZ+VJqTYP/bcBLhP79fjeP\n1mvZmQCSlEDTUT2/RujWeSvwauBngE8RWvnbgHPAdsLJAeAssKvy/p1xXp25SnkhTpKkFXvjlMwN\nrIzqOQocjOVDwJFY3gOcADYBu4HnqL9S8CpgunWgU+nnryuv9/ogy7ouPKZKMdDnPKpx/N2NHgEe\nAN5DuIn7jjj/ZJx/EvgJcCfukJKURG5f7uiQX500Op3VH3Fdeb3XB1nWdfXw2JpeA2WnY+k1KV7h\nJeN/vVYz+CWpMAa/JBXG4Jekwhj8klQYg1+SCmPwS1JhDH5JKox/gUvj5iDyfHQ/C7/IVThb/JoA\nsz8PnWrB3/EpmMEvFcnML5nBL0mFMfglqTAGvyQVxuCXpMIY/JJUGMfxaxwcL94e1eE9fl6FsMWv\ncXLMYCv4MZXG4NcYGShSjgx+SSqMwS9JhfHmrkbJvh2pBWzxa1Q6qx7URv5wWyEMfkmRmV8Ku3o0\nLNNCahlb/BoBs19qE4NfTdgXLLWYwS9JhTH4JakwBr+GYXeP1EIGv4Zg7kttZPBL6uUZfcoZ/JJU\nGINfkgpj8EtSYQx+SSpM0+DfBXwJeAb4BvDeOH8LcBw4BRwDNlfecxg4DSwC+xpuV5KUyDbgDbF8\nKfBN4PXAUeBAnH8QOBLLe4ATwEZgFjhD/UnH0QR566yeOj2P/crrvT6uZV1X83X5sxwtM9Bn1bTF\nf44Q5AA/BJ4FdgA3A/Nx/jxwayzfAtwPvAy8SAj+axpuW0mZBWXonFfQ9BhFH/8scDXwBLAVWI7z\nl+NzgMuBpcp7lggnCknShA37e/yXAv8E3AX8b89r610q9nttrlJeiJMkacXeOE3cRuCLwN2VeYuE\n/n+A7fE5wKE4dT0GXFuzTi8r85VJP3QO/d+lrav6+StTA302Tbt6ZoC/BU4CH63MfwS4PZZvBx6u\nzL8N2ATsBq4Anmy4bUlSAtcDrxBu8D4Vp/2E4ZyPUz+c8x7CTd1F4KY+67VFka9MWqW5tYZLWFf1\n81emBvpsZsZVi4Y65FenknV3ppmVcvUj6lRe6p03yOvjWtZ1jWZdq3h85mmg7PSbu5LWYUN/2hj8\nklSYYYdzajrZxJOmmC1+9WH2S9PK4JekwtjVowth819dvfuCo3xayBa/utYYp23uq6vT8/j/+407\nSYsY/OrlAawBucu0jcGvHh7E0rQz+CWpMAa/JBXGUT2SRq3aX+ionwzZ4pc0Kp3aorJj8JfNYXga\nNfenFjD4JY2Qud8GBr/Ao1UqisFfpp4uHnNfKonBXyzDXhPhjpYhh3OWxYNQki3+gnRWPUgqlsEv\nSYWxq2f62cRXDvwd/4zY4i+C2a+kHEGWGYN/unmUKQPuhrkx+CWpMPbxTx+bV2qDDvbzJ2OLv938\nO7mSBmaLf3qY9Gqb7j5ry3/CbPFPFbNfbXLe/trBnwqfCFv808EDRW3WM9yz7wWAf9lrRGzxTwVz\nX21Wu//2af27r4+CwS8pU4b8uNjVk6f1Lmk9IiQ1ZvDno+aSdqZ3/szKc4dBS2rGrp6s9GvId9Zb\nQJpmjvQZsUkH/35gETgNHJzwtqeA+75Ktea+33sj2BPFOiYZ/BcBf0UI/z3Au4DXT3D7Te0d8fo6\n1O+oKsZC6gq0xMJaL/YcQxM7hPZOakPjNMngvwY4A7wIvAz8I3DLBLff1N7xrNbum3ItpK5ASyys\n8/q6h07dkNBhvyS2t+H7sjLJ4N8BfLvyfCnOmwb9dqS61n3NIpJGpN9v/9ddHRT7TeFJBv+o/3O7\nI1w6cNGHR7C+3pDuTh9c47W6nWatS9DidjBpstYdINFvXr/jer08aKVJjge8Dpgj9PEDHAZeAf6i\nskxr/yMlKbEsx3dvAJ4DZoFNwAnacXNXkjSE3wG+SbjJezhxXSRJkiSl8CfAs8A3WLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x9bf1588>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(data, bins=200)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "that is not very useful to us - we want the PDF, not bin counts. Fortunately `hist` includes a `density` parameter which will plot the PDF for us."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEMlJREFUeJzt3WuMXOddx/Gv8a6BXJTEYGJwtlqKimoJuS0KrqME4nCJ\nNhGqxau+iIDWJfILwkVE1DWVmkWqEBehogqoorSgAAFzUVKlVUpiS10RRYmDwZcW4mA7NrFdq1iQ\nWqJSVUcML56z3rOzZ/bMmdvznHm+H+loz3X2PztnfvPM85yZBUmSJEmSJEmSJEmSJEmSJI3QAnAK\nOA3sr9j+EHACOAm8BOwobTtfrD8GvDrWKiVJA9sInAHmgVngOLC9a5+7gFuK+QXgldK2c8Dm8ZYo\nSarzHTXbdxLC/jxwDTgI7Ona52XgajF/BLija/uG4UqUJA2rLuy3ARdKyxeLdb18BHiutNwBDgNH\ngYcHKVCSNLyZmu2dBrd1H7AXuLu07m7gMrAFOETo+3+xSYGSpOHVhf0lYK60PEdo3XfbATxB6LN/\nq7T+cvHzCvAMoVuoO+ybvKBIklaMrJt8BjhLGKDdRPUA7TsI/fq7utbfANxczN9IuFLn/orfkWLY\nL8YuoIfF2AVUWIxdQIXF2AVUWIxdQIXF2AX0sBi7gAqLsQuo0Cg761r2bwOPAM8Trsz5HPAasK/Y\n/jjwCeA24DPFumuEFvxW4OnS73kKeKFJcZKk0agLe4AvFVPZ46X5Xyqmbm8A7x2wLknSCNVdjZOr\npdgF9LAUu4AKS7ELqLAUu4AKS7ELqLAUu4AelmIXUGEpdgHTIMU+e0lKXaPstGUvSRkw7CUpA4a9\nJGXAsJekDBj2kpQBw16SMmDYS1IGDHtJyoBhL0kZMOwlKQOGvSRlwLCXpAwY9pKUAcNekjJg2EtS\nBgx7ScqAYS9JGTDsJSkDhr0kZcCwl6QMzMQuQGqR8j943hCtCmkAtuylRjr1u0gJMuwlKQOGvSRl\nwLCXpAwY9pKUAcNekjJg2EtSBgx7ScqAYS9JGTDsJSkDhr0kZcCwl6QM+EVo0lrLX4Cz3ped+aVo\napV+WvYLwCngNLC/YvtDwAngJPASsKPBsVKL+aVomh4bgTPAPDALHAe2d+1zF3BLMb8AvNLgWPAZ\nozR0WDkXy/Nd+3Q6K9s7vfaTJqHRuVfXst9JCOzzwDXgILCna5+XgavF/BHgjgbHSpImoC7stwEX\nSssXi3W9fAR4bsBjJUljUjdA2+Rtwn3AXuDuAY5dLM0vFZMUWz8DtdKk7C6mgdSF/SVgrrQ8R2ih\nd9sBPEHos3+r4bGwOuylRHQw55WQJVY3hB8b5Y3PAGcJg6ybqB5kfQehb37XAMeCA1yKq9M1ldat\nGYB1gFYpGfm59wDwOiHQDxTr9hUTwGeB/waOFdOrNcd288miSesK9s6wYd99vDQJjc63FN6j+l5Z\nk1bui+90nYIbVrZfX79h5bjO6sW1857LmpRG2enXJUir9Wot2WpXqxn20iqVmW7Qq/UMe6mWWa/2\nM+yVs3GkuAO1SpJhr1xUhPC4sn5sNy4NzLCXpAwY9pKUAcNekjJg2GuaVQ2WTrIv3cFaJcOwV2Ym\nnfVSGgx7ScqAYS9JGTDsJSkDhr0kZcCwlybHq3MUjWEvjZ8Br+gMe+XAsFX2DHtlwKyXDHtJyoBh\nL0kZMOw1LbzSRVqHYS9JGTDsJSkDhr0kZcCwV5vZTy/1aSZ2AdIIxPwHJf1KsSZlxJa9pkBVjqaW\nranVo9wY9pKUAcNekjJgn73aaL0+EftLpAq27NVSvTLdrJeqGPaSlAHDXpIyYJ+9FEev/qYNE61C\n2bBlL0XT6TF/fYWfENbI9BP2C8Ap4DSwv2L7u4GXgW8Bj3ZtOw+cBI4Brw5cpRRkFn4Z3VVFtxE4\nA8wDs8BxYHvXPluAO4FPsjbszwGba36HZ7T61YFOZ/XPXvN122PdFnXbK+6rVKnRuVHXst9JCPvz\nwDXgILCna58rwNFiexX7IKVVzG9NXl3YbwMulJYvFuv61QEOE14MHm5WmiRpVOquxhm2CXI3cJnQ\n1XOI0Pf/4pC3KUlqqC7sLwFzpeU5Quu+X5eLn1eAZwjdQlVhv1iaXyomSdKK3cU0FjPAWcIA7Saq\nB2iXLbJ6gPYG4OZi/kbgJeD+iuPswFS/pmSAdt19K+6rVGnk58YDwOuEgdoDxbp9xQSwldCvfxV4\nC3gTuAl4J+HF4Tjw1dKxYy9YU8uwl1Y0OjdSuFKmQxp1KH2dldOlfNpUzddtb7LvKG+rn32XN15f\n5/NDVRplp5+glaQM+N04Sp3dGNII2LJXC5j30rAMeyk9vrpp5Ax7KTlmvUbPsJekDBj2kpQBw16S\nMmDYS1IGDHtJyoBhL6XPy3M0NMNekjJg2EtSBgx7qR38umMNxbBXygy36/xTaDh+66VSZLJJI2bL\nXoky76VRMuwlKQOGvSRlwLCXpAwY9pKUAa/GkdqlPHK9IVoVah1b9lLreKWSmjPsJSkDduMoFTZX\npTGyZa/YSt/5Yt5L42LYKwGGvDRuhr0kZcCwl6QMGPaSlAGvxpHayw9YqW+27KVWc3Bb/THsJSkD\nhr0kZcCwl6QMOECrGOxoliasn5b9AnAKOA3sr9j+buBl4FvAow2PVbbMeyklG4EzwDwwCxwHtnft\nswW4E/gkq8O+n2PBZ32OOtBZ/k6c0nzVukH3TfW2xvV7laFGj3tdy34nIbDPA9eAg8Cern2uAEeL\n7U2PVX4MJimCurDfBlwoLV8s1vVjmGMlSSNUF/bDtMJswUlSIuquxrkEzJWW5wgt9H40OXaxNL9U\nTJL651cnTL/dxTQWM8BZwiDrJnoPskII7PIAbb/H+g4gLwkObk7DAG15nTIx8sf6AeB1wmDrgWLd\nvmIC2Erom78KvAW8Cdy0zrFjL1hJSzwoDXu1RqPHOoW3ex3SqEPj0eOELD/sy/NV6wbdN9XbmsTv\n9fmUiUbZ6dclaAI6XT8lTZphL0kZMOyl6WPfvdYw7KWpY85rLcNeml6mvq4z7DVOho2UCMNekjJg\n2EtSBgx7ScqAYS9JGTDsJSkDhr0kZcCwl6abn6YVYNhrPAyYZHS6F3xcMmXYa0zMFCklhr2UB199\nM2fYS1kw63Nn2EtSBgx7ScqAYS9JGZiJXYCmhp3C7VF+rPzn5JmwZa8RMu/bwccpR4a9JGXAsJek\nDBj2kpQBB2g1DDt/pZawZa8hmfct55ejZcKw1yAMiKnhw5gLw14DMiSkNjHsJSkDDtCqXzblpRaz\nZa8GzHuprQx7ScqAYS9JGTDsJSkDhr0kcEBm6vUT9gvAKeA0sL/HPp8utp8A3ldafx44CRwDXh24\nSknSWG0EzgDzwCxwHNjetc+DwHPF/PuBV0rbzgGba36HLYp26EBn+ZOzpfmqdYPum8NtpXwf1DKN\nHrO6lv1OQtifB64BB4E9Xft8AHiymD8C3ArcXtruf8KRpMjqwn4bcKG0fLFY1+8+HeAwcBR4ePAy\nJUnDqPsEbb9vE3q13u8BvgZsAQ4R+v5frNhvsTS/VExKh2/x87D8OPtuPE27i2kgdWF/CZgrLc8R\nWu7r7XNHsQ5C0ANcAZ4hdAvVhb2kKDoUOW/op2mJ1Q3hx5ocXNeNcxR4F2GAdhPwQeDZrn2eBX6h\nmN8FfAP4OnADcHOx/kbgfuArTYqTFINv5KZRXcv+beAR4HnClTmfA14D9hXbHydcifMgYSD3m8CH\ni21bgadLv+cp4IVRFS5p7Mqpbyu/5VJ4AK+/d1RyOmsXN3TNV60bdN8cbqu198HnaHoaZaefoFUN\n39JL08Cwl6QMGPbq5qcppSlk2EtSBgx7ScqA/4NWy+y6kaaYLXuVmPfStLJlr15MfpX5AauWM+zz\ntk6g+1k3dVvz3TngSdIaduNkzwa8mvKcaSPDXpIyYNhLUgYMe0nKgAO0kkbBQdvE2bKXNIzSdyk5\ncJsyW/b58RmpEfIS3bawZZ8l817KjWEvSRkw7PPh99RrUjzPEmTY58EBNClzDtBON9NdEmDLPgPm\nvSRb9pLGp7ul4TWaEdmyn1426RWT40SJMeynj1fdKAGegqkx7KeSTzRJq9lnPx1Md7XB8nlq330E\nhn37+K8E1VKr/q2hJ+qE2Y3TSp0e85JUzbCXFJsXFUyA3TiSYugV7t39+v5TlBGxZS8pgk7Xz7UL\ntavViGGflrq3s571mmKe3uNk2EtSBuyzj6/Xtcc2c5Q7nwMj1E/LfgE4BZwG9vfY59PF9hPA+xoe\nq6DH9ZSVfZvStOvnhPcqnhHaCJwB5oFZ4DiwvWufB4Hnivn3A680OBbSfLB2j/j2Ol1TxbZOZ/35\nuu1N9h3lbX050bpSu62qv1Pb7kPs2+o5UbE8arvHdLvDaHRf61r2OwmBfR64BhwE9nTt8wHgyWL+\nCHArsLXPY1O1e/Q32eleGOeJOUFLsQtoiaXYBUyBTo/5Neuqnl/DPt92D3FsEurCfhtwobR8sVjX\nzz4/0MexbdYrsHu1PLp2oXqTpKbWeY7VPj+zURf2/f4xJvhhh5nPMNoHqurBf4zVidwjwHu1JAxz\naXJqn2edtbOVz91ez3lYyYTWvkjUXY1zCZgrLc8RWujr7XNHsc9sH8cua/DHe3vA4/pW916xpNdr\n3IaK7VXzddvHte8ob+u3E60rtduq+ju17T609baoWN+zfbpu/1DNulabAc4SBlk3UT9Au4uVAdp+\njpUkJeIB4HXCYOuBYt2+Ylr2x8X2E8CP1hwrSZIkaZr9CvAa8FXg9yLXUvYo8H/A5tiFAH9A+Bud\nAJ4GbolYS4ofmJsDvgz8G+E8+tW45ayyETgGfCF2IYVbgX8gnE//TuiCje0A4bH7CLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x861ca90>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(data, bins=200, normed=True)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "I may not want bars, so I can specify the `histtype` as 'step' to get a line."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGjdJREFUeJzt3XuUFOWZx/FvyUUUuYiCKKDjDUWJimsA19uYqCHqakzi\najSJihqyR41JSFRi1Ik5UZM1J4myG413E5SELHE18QYnO9EYLxC5KIqCgIGRCKJ4RwZ49o+nxu5p\neujume5+q7t+n3PqTHVV18xPnHnmnbfe9y0QEREREREREREREREREREREREREZEyGgcsBBYBl+Y5\nfyYwD5gPPAEckHVuWXx8DvBMRVOKiEindQMWAw1AD2AuMCLnPYcC/eL9ccBTWeeWAgMqG1FERArZ\nqsD50XixXwa0AlOBk3Pe8yTwdrz/NDA053zUtYgiItJVhYr9EGB51usV8bGOnAs8mPXagJnAbOD8\nzgQUEZGu617gvJXwuY4GxgOHZR07DFgJDARm4H3/j5cSUEREuq5QsW8BhmW9Hoa37nMdANyC99m/\nlXV8ZfxxNfAHvFsot9iX8gtFREQyytZN3h14Bb9B25P8N2h3xfv1x+Yc3xboE+/3xkfqHJfnaySx\n2DeFDtCBptAB8mgKHSCPptAB8mgKHSCPptABOtAUOkAeTaED5FFS7SzUst8AXAg8go/MuQ14EZgQ\nn78ZuBLYHvhlfKwVb8EPBqZnfZ0pwKOlhBMRkfIoVOwBHoq3bDdn7Z8Xb7mWAAd1MpeIiJRRodE4\nadUcOkAHmkMHyKM5dIA8mkMHyKM5dIA8mkMH6EBz6AB5NIcOUA+S2GcvIpJ0JdVOtexFRFJAxV5E\nJAVU7EVEUkDFXkQkBVTsRURSQMVeRCQFVOxFRFJAxV5EJAVU7EVEUkDFXkQkBVTsRURSQMVeRCQF\nVOxFRFJAxV5EJAVU7EVEUkDFXkQkBVTsRURSQMVeRCQFVOxFRFJAxV6kKHYk2E3xNiR0GpFSqdiL\nFGc/YH/gJGBA4CwiJeseOoBIDVkA9A8dQqQz1LIXEUkBFXsRkRRQsRcRSQEVexGRFFCxFxFJARV7\nEZEUULEXEUkBFXsRkRRQsRcRSQEVexGRFFCxF2nHdgMbAbbnFt50JJjF22lViybSBVobR6S924BP\nAGuBfbbwvseBlVVJJFIGxbTsxwELgUXApXnOnwnMA+YDTwAHlHCtSALYkWDHgh0UH7g6aByRCijU\nsu8GTAaOAVqAWcD9wItZ71kCHAm8jRf3XwFji7xWJAmmAu+g702pY4Va9qOBxcAyoBX/oTg55z1P\n4oUe4GlgaAnXiiTFjaEDiFRSoWI/BFie9XpFfKwj5wIPdvJakSTpA3Yq2N6hg4iUQ6FuHCvhcx0N\njAcO68S1TVn7zfEmEsp7+P2na4Eb8HtOIqE1xlunFCr2LcCwrNfD8BZ6rgOAW/A++7dKvBbaF3uR\nKrELgG8DO+WcaIHoVLAbAoQS6Ugz7RvCV5VycaFunNnA3kAD0BM4Db/Jmm1XYDrwZbyPvpRrRarM\ndgM7EOwT+CMGHwGGA6s6+Qn7gP0q3k4qW0yRMitU7DcAF+I/EC8Av8VHLEyIN4Arge2BXwJzgGcK\nXCsS0jX4faXH49dvQfQK/v1arEOy9rcGzgB2Bv6lLAlFKqCYSVUPxVu2m7P2z4u3Yq8VCe0a4Ed5\njg8Gttn8sA0E+sQvHsa7JP8MjIiPr8P/ktWMdEkszaAVca8Db8b7a7KOnwP8Ij53B0TfzZyyW4GJ\nOe8XSSQVexEAor8BJ+QcvBP4a7z/e4g25VyT9VetlXSzTKTaVOwlrbYFTsS7YzoQPQs8W8LnHAfW\nH/g+RO92KZ1ImanYSwpYI9CXzLDgD4G27phZZfoij+CLp/0IuAWsBaK3ClwjkiqlTL4S6QSbA/Yc\n2N/ApoCdWcGvtQBsLdi0yn0NEaDE2qnRA5IWN1Xny0T7A+dX52uJFE/FXkQkBdRnL3XK2hbd+yDr\n4Nb4jdlq2CbO8AZEH1Xpa4p0SC17qVdL8ZnbP4lfr8cnTY3Gb9BW0jpgVJzh0Ap/LZGaoRu0UgG2\n3hc6s1viG7SjAmRojkcCiVSCbtCKiEh7KvYiIimgYi8ikgIq9iJVYT8AW+ybSPWp2ItU1q1g04Ed\ngd8AewbOIymlYi/1ri8+vj6E84FLgX3j16sD5RDRpCqpa+/i4+rBx9lXWbQIrEf1v67I5lTspZ5N\nheju0CFEkkDdOCIiKaBiL3XAHgRrBVsYOolIUqnYSz3ohj8LVt2SIh1QsZd6sTF0AJEkU7EXEUkB\n/dkrNcoeAAYBj4dOIlIL1LKXWnUQXui/CRyTOWzXgb0DJGl8+77A10KHEAlN69lLJ9hysD3A+sbb\nvvHaMz8DmxQfi0KnBOuWlbEnmL7fpVxK+l5SN47UslaI3vFdG5R1fF3meGjRRiAri2q9hKFuHBGR\nFFDLXmqM7QOMJP+Dw3vgC5+JSA617KXWnAJcCzTT/sHhG4AIOBZfAE1EEkadmFICu8xH3NQq3aCV\nstEDx0VEpD312YtUnfUFdohfrMnaXwXR+2EySb1Ty16k+r4CzAdejvdfBBYAn/XTNhzs+HjrGSqk\n1Be17KVGWHfgUGB34K3AYcrhbuBofATRk3gLv82XgS8BDcDgnHMinVJMy34csBBYhD9PM9e++Dfr\nOnyZ2WzL8BbMHOCZTqcUgW3wETj74d9Xtezy+GMr8OP4Y65f024ylkjXFGrZdwMm42uPtACzgPvx\nPzvbrAEuAj6X53oDGoE3uxpUBPgAoiNCh+iiA+OPayBqyRy23wdJI6lRqGU/GliMt6RaganAyTnv\nWQ3MJn/rBHzss4gAEM2Pt5bC7xUpn0LFfgiwPOv1ivhYsQyYif8yOL+0aCIiUi6FunG6OgHkMGAl\nMBCYgff9a/1xEZEqK1TsW4BhWa+H4a37Yq2MP64G/oB3C+Ur9k1Z+83xJiIiGY3x1imFiv1sYG98\nCNhrwGn4kLB8cvvmt8Vv8L4L9AaOA37QwbVNhaOKiKRaM+0bwleVcnGhYr8BuBB4BC/ct+EjcSbE\n52/GxwHPwlcb3ARcjA+PGwRMz/o6U4BHSwknkjK/Avs+PuJNpO5oYSgpgvUBq+PVLG0HsFFgr4Nd\nDXYl2Bo/LpKXFkKTemEjwWb5xgmh01RWtAbvKhWpCBV7SbLe+L2f1WQWCxORTlCxl6R7D3gjdIgq\n2groFTqE1B8Ve6kV3UIHqALDF0b7GvBR4CwiZacbtNIBGwP2NNidYBvA1oZOVF26QStbpBu0Um+i\nsyHqDlH/0ElEapWKvYhICqjYi4ikgIq9iEgKqNiLJNs/waaEDiG1T8VeJLl2Bs4hHcNOpcJU7EUS\nK1qPL0Yo0mUq9iLJNwbsdrDhoYNI7VKxl4Sy7wITQ6dIgGeAHwJHADsFziLSJZpBK1lsW7CdweaD\nXQ92WuhEyWCPgx0ROoUkSs3VzpoLLJVk54B9CLYS7MDQaZJDxV42o+USpOZNhWhniOaFDiJSL1Ts\nRURSQMVeRCQFVOxFRFJAxV6kdoyMH0reDUw/u1ISfcOI1IYFwOXAs/is2vPDxpFa0z10ABEpRvT1\nzL7dHC6H1CoVe0kAGw98Pn7xp5BJROqVunEkIIvAegIjgdXAseh7UqQi9IMlIe0GfARcCDwPbAwb\nR6R+qdhLaK9C1BOin4YOIlLPVOwlaT4ROoBIPVKxlyR5CBgKaE0ckTLTaBxJkOgLoRPUkG7+QBMA\nZkJ0T9A0knhq2YvUpsHA2cA2wJiwUaQWqNiL1J5VwOnAS8BTgbNIjVCxF6k50RUQDYdoROgkUjtU\n7EVEUkDFXqrMtgMbAfYXYGnoNCJpoWIv1XYkMBsYBBwN7B42jkg6FFPsxwELgUXApXnO7ws8CawD\nJpZ4raTTX7y/OWqGSA+cF6mCQsW+GzAZL9r7AV8Ccm8KrQEuAq7vxLWSKjYKbxyISJUVKvajgcXA\nMqAVmAqcnPOe1fif5a2duFbSZTrwVXzIoIhUUaEZtEOA5VmvV1D8BI6uXCv16xSIdGNWpMoKtey7\n0p+qvliR6ugF9gzY62DTQoeRZCrUsm8BhmW9Hoa30ItRyrVNWfvN8SYihX0EnBrvX4c/AEbqU2O8\nVUR34BWgAegJzKXjm6xNtB+NU+y1+gsgNWwpmIZaVowdCzYjdAqpmpJqZ6GW/Qb8KUKP4KNrbgNe\nBCbE52/GF2SaBfQFNgEX46Nv3uvgWkkN6wXMwBfrMmCHsHlEJCS17OuWbQv2EdhFYN8BmwjWL3Sq\n+qWWfcqUtWUv0lUbIboxdAiRtNNyCSL1ZX+w28AODh1EkkXFXqR+vABcARwM7Bo4iySMir1I3Yha\nILoNeBUYDfaZ0IkkOVTspUKsF7Bd6BQpNQs4BLgjdBCRbBqNU5fsjngkzurQSdLJdgF7LXQKqaiS\naqda9lJJX4doYOgQIqJiLyKSCir2IiIpoGIvIpICKvYiIimgYi9Sv7qD7ar1iAS0No6Une0CnAsc\nBDwWOEyabQTWAc8DPwZ7Cv95XwnR/KDJJLU0zr6u2CFgK8F+qPVZksB+BHY52FqwF8GmhE4kZaNx\n9hJcC0RXQPRs6CAC+GzarYHJoYNIOCr2IvXt78AaYAqwPnAWCUh99iJ1LZoOTPd9OzNoFAlKLXsR\nkRRQsRcRSQF140gZ2CeBti6Cx0MmkYL6gO0V76+H6B9B00jVqNhLOYwADgcGAb0CZ5GOvQfsBzyM\nj855AxgVNJFUjbpxpFxeAKYBQ1DrPqGi/4VoL984KXQaqS617KWMoomhE4hIfmrZi4ikgFr20knW\nD3/G7LeB4cBbYfOISNJpbZyaZNeBvQNmYNeAqQ+4ptgosFawFrAzQqeRTimpdqplLLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x9a4e780>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(data, bins=200, normed=True, histtype='step')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To be sure it is working, let's also plot the idealized Gaussian in black."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFOWdx/FPz8ExAwKjBARBQFEIAYVVQMAwEg8gKhp1\nSaKJq1nXaA6TuKuSRB2TrDndTVaTeCTZaA6JrtFgggokjhKjIIbDiMgh9yn3NcNcz/7xVEPTztDd\nM9391NP1fb9e/arq7qquL0P3r6ufeuopEBERERERERERERERERERERERERGRLJoILANWALc38/zV\nwGJgCfAKMCzhuTXB4wuB+TlNKSIirVYMrAT6AaXAImBw0jLnAF2C+YnAawnPrQYqchtRRERSKUrx\n/EhssV8D1APTgSlJy7wK7Anm5wEnJT0fa1tEERFpq1TFvjewPuH+huCxlnwGmJlw3wBzgAXADa0J\nKCIibVeS4nmTwWudB1wPjE14bCywGegOzMa2/c/NJKCIiLRdqmK/EeiTcL8Pdu8+2TDgEWyb/a6E\nxzcH0/eAp7HNQsnFPpMvFBEROSJrzeQlwCrsAdp2NH+Ati+2XX900uNlQOdgvhzbU+fCZrahYn9E\nlesAIVLlOkCIVLkOECJVrgOESEa1M9WefQPweeAFbM+cnwNvAzcGzz8E3AV0A34aPFaP3YPvCfw+\nYTu/AWZlEk5ERAqH9uyPqHIdIESqXAcIkSrXAUKkynWAEMmodqbqjSP5Ve06QIhUuw4QItWuA4RI\ntesA0nrasxcRyZz27EVE5Ggq9iIiEaBiLyISASr2IiIRoGIvIhIBKvYiIhGgYi8iEgEq9iIiEaBi\nLyISASr2IiIRoGIvIhIBKvYiIhGgYi8iEgEq9iIiEaBiLyISASr2IiIRoGIvIhIBKvYiIhGgYi8i\nEgEq9iJtYm4C8yCYC1wnETkWFXuRtrkIuBIY5jqIyLGo2Iu03R7XAURSUbEXEYkAFXsRkQhQsRcR\niQAVexGRCFCxFxGJABV7EZEIULEXEYkAFXsRkQhQsRcRiQAVexGRCFCxF0mbKQczGMzAZp68CYwB\nMzvvsUTSoGIvkr6zgCVAdQvPr8pfFJHMpFPsJwLLgBXA7c08fzWwGPsheIWjR/9Lta6IB8zAYAjj\nIcBa12lEcqEYWAn0A0qBRcDgpGXOAboE8xOB1zJYF8BkM7BI9pl7wdSAWQVmJZiNCc89k/D4SjXj\nSB5lVDtT7dmPxBbsNUA9MB2YkrTMqxwZ4nUecFIG64r4YmPqRUTCK1Wx7w2sT7i/IXisJZ8BZrZy\nXREP1AMLOsE374QeJ6VcXCQkSlI8n8nPhPOA64GxrVi3KmG+mpYPgIk4Ug98qyv85DjYXgJ8w94+\ntBMerIFxtY4DSuGrDG45MRp4PuH+NJo/0DoM22RzaivWVZu9hJi5H7bugX+qwb5XDbAaRjUCB+39\nmIFvblebveRZVmtnCbY7WT+gHc0fZO2LLfSjW7EuqNhLKJmhYM6ATS9Ajw2AgZ4N8PgSIGYP2N5y\nPHx6BYe/BL6+A8wiMA8HN3VtllzKeu2cBLyDLejTgsduDG4APwN2AAuD2/wU6yZTsZcQMnugoQ4q\ng733Lltg3RowLwXP14DpaHvj/HQrFAUFf/pmMNuDE6xU7CWXvKud3gWWKDB74K7tgIH2e+HP/w5m\nS1Kx7w3mOTBPwZcW2mU7N8HKdWCaVOwlx7yrnd4FlihYug9Kg731W78N5kIwfwLzXfu8qQn23g2Y\nyUAMes21y1fuhkYVe8k172qnd4Gl4MXgvHrAwCf3gvnn9y9irrCPmy4JD/YEdtv1fqdiL7nmXe30\nLrAUvEsAA10bYWsLxb5Fn7XrDjCw/wEwH8hRRhHvaqd3gaVQmQqonQJdVgEGpj4B5hYwp2XwIiXA\n23b9H+0DMxZM+9zklYjzrnZ6F1gKlRkN0w8BBno1wa9bu1d+hX2N7oegpg7MpCyGFInzrnZ6F1gK\nVd1oGBoU+/trktrjMxHDjgJr4DtLVOwlR7yrnd4FlkJ1w83YPfImOLinDcUe4Br7Wn33Q+3kLAUU\nSeRd7fQusBQac5ztM3/SPOyZsAfB7GtjsS/Fjn1v4I4fgOmUnawih3lXO70LLIXG3A3vNtoxbto1\nwebNdsx6c1wbX/g2wMDkBjBfyEZSkQRZHc9eJCLu32s/OxN3wIknQqw3xPa28UX/F6iD54rhpW5Z\nCCniNe3Zi2MH7oHuDYCBmW9l+cV/Y1934vMplxTJjHe107vAUmj+83HAwKl10PhaysUzM96+9nE7\n0C9pyS7vaqd3gaXQnLkUMHDPdjDZLvZFULHPvj7jsvzaEm3e1U7vAktBqYDiBjtE8drVOSj2wJVv\nAAZOmxFc4GRW9rchEaQDtCIZmAqNxXBuDfRtzM0m/vUdO90yCQ61w17wRyRytGcvLr0EGHhoK5i1\nudmzN32g9za7nembwCzL/jYkgrRnL5KmHsC5UNwIV70HNAGbs7+Z2Ho4J/gS+W3n7L++iB+0Zy+u\nBMMRf+gdMPfkdlNz77bbat8EO5fndlsSEdqzF0nTFXZy0dLcb2rcbhhdC4di8KSGTpC8U7GXiNq2\nGYrPh1gj3PxOfrZ52QE7naFiL3mnYi8R9VQ5NAIVb8CAmvxs84qg2FeXA7qgieSVir1E1NNBsR36\nYv62eWoDDKmDA0XAhPxtV0TFXqKpK7zYzr79b305v5ueHOzdMyW/25WoU7GXCDGDwMyD+1+BemBk\nLVy8K78ZPrY/mJmCPn+SR3qzSZSUAb3htQ/auxcfAB4Bvpqn7Q+AkXXQqwHoCZydp+2KhIL62Uue\nmBFQvxBiOwADf18cjFXzb2A65Hjb7YMrYo2Az+y02+fe3G5TCpz62Yu07KUyMBXAu3DGweDBGojV\n5na7sUPBxVAOwmXxppzLcrtNkSNU7CViZsQvNTjT3dv/woPAHmAw0N9RCIkYFXuJCHMS8BGYHb+I\n+J+CaWn+s7QDmBPcuSj/25coUrGXqDgH1n8F3u4I1GBHu9wKxICDx1wzN+KXKZzoYNsSQSWuA4jk\nz6MbsL1g/oxtp3fZZv5CMJ2A3dWvc5hFIkB79hIhs04MZmY6jWGtB5YCnYFzHGeRCFCxl4jYUwwL\negR3XBf7GJj+0PNv9m6vqbZbpkhhUz97yYNv3Yl9r73tNocZBGYVGAMvNNpMZ5qgr39PMJPBnOs2\no3hC/exFjmb6wML4wGOznUY5ogHGbYAOBhYB87oAY4D/Ax5wG00KUTrFfiKwDFgB3N7M84OAV4Fa\n4Nak59YAS4CFwPxWpxRpm2tg7ZhgPgzFvo+dlNVBZbB39otgCIdcXBZRJLViYCXQD9sfeRH2RJBE\n3YGzgG/x/mK/GqhIsQ0140iOra2Coibs6GeOrwFrOoAZZm8AfBF7acT5YD4WNPEsdplQvJHVZpyR\n2GK/BvtBmc77h2Z9D1gQPN+cWCaBRLLvsQHQFANeA/a5zRKrhdgSewMO/9JYNQganKWSwpeq2PfG\ndhGL2xA8li6DPVNwAXBDZtFEsuUvA4OZOcdczI1l0PUg1HSGR/u4DiOFK1Wxb2sTy1hgODAJ+Byg\nXgbiwJJTg5kwtNcnM3DmJjv7zFC3UaSQpTqDdiOHDyZBML8hg9ePH2x6D3ga2yw0t5nlqhLmq4Ob\nSDb0hR0nQIdDUBvSTgLnboLqU+EfKvZyLJXBLSdKgFXYA7TtaP4AbVwVRx+gLePIwbBy4BXgwmbW\n0wFayaXrAQPDl7sO0rKFvwUMlB6AOh2glXRlvXZOAt7BHqidFjx2Y3ADO9bIeuyQrbuAdUAnYAD2\ny2ER8I+EdXMeWCRBUEhvnuU6SMvMw9AvOMFq7gYVe0mTd7XTu8DijRiwBTAw+2HXYVpmusF5CwED\nd+1QsZc06QxaEeul6UAPOL4eJuxwnaZlsV0wcqmdr+7oNosUKhV7KWAvD7fTcw+F/61+9So7nd8B\nanRuimRd2D8BIm3wclc7raxxmyMdQ/fD6fVQG4PqMtdppPCo2EuhKoLXg2J/QY4vJp4Vh+D8oA12\ndie3UURyQwdoJReGAgZ61kNjA5j/dB0oDZcCBkbsdx1EvOBd7fQusHghGGDs5BB3uXyfrkAjlBhs\n92WRY1FvHBHgPDv54EK3MTKyGzotDQZEG+c6jBQWFXspREXAeDs71adiD/SMD+nwEacxpOCo2Esh\nOgPoBj0PwbVbXIfJzPB5wcyEYy4mkiEVeylEQRPO6F1uY7TG1xbb6wTFRmDb8EWyQsVeClFQ7Mfv\ndhujNc54HXgtOPb2YcdhpICo2EuhKeFwkZzsYbGPNUFpdXCn0mEQKTAq9lJohgPHAavgtEOuw7RO\nRXDNh37XglFTjmSFir0UmngTzl5g4DGXDK1xc6GoAdZWwA9PdJ1GJFt0UpVk00zAwP/uAnMXmJNc\nB2qlasDA0OtcB5HQ8q52ehdYQqsE2AcYWLcVTC/XgdqgCjDQ7RHXQSS0vKud3gWW0BoJGBjY5DpI\nFozHXqpwqesgEloaLkEiq9JOxhfCDsQ86ADUDwZOcB1G/KdiL4WkMpgUwp59LZxVF8yPd5pECoKK\nvRSKEuBcO1sQxR4YHy/25zmNIZIlhfCTW9wL2utjy8F42r8+2V+2AwaKl4LRjpkkU5u9RFKlnRS/\n5DRFVo2aDx2BxsGw4HOu04jfVOylUFTaSYeXnabIqvLJUDPHzj82yG0W8Z2KvRSCEigKDmK+Ptlt\nlKyrtpPXVOzFe2qzlzYqGw0YOKUBzMrCabMHYAz25KpNroNI6HhXO70LLGEz9VnAwDW7CrDYl0K7\neuznpKfrMBIqOkArUbPoFDtde6/bHDlRD0O2BfOVLoOI31TsxXel8O4AO3v160C50Line truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x830f6d8>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.hist(data, bins=200, normed=True, histtype='step')\n",
+    "plt.plot(xs, norm.pdf(xs), color='k', lw=2)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There is another way to get the approximate distribution of a set of data. There is a technique called *kernel density estimate* that uses a kernel to estimate the probability distribution of a set of data. NumPy implements it with the function `gaussian_kde`. Do not be mislead by the name - Gaussian refers to the type of kernel used in the computation. This works for any distribution, not just Gaussians. In this section we have a Gaussian distribution, but soon we will not, and this same function will work."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGlZJREFUeJzt3XvUXXV95/H3zo2QC4QIIRCSBhEUUEQawiWRRuWSoAVc\nyjDKqlqVYVjI1I51GDozJWu1XbaOdk3V1l5kUFqFQSgVEFrDJRohXDLkIkggKJEEAgZzIUAdI3zm\nj72TnDx5nudcnn3Od//2+bzW2us5l73P+XB48n32+e7f/m0wMzMzMzMzMzMzMzMzMzMzMzMzM7MS\nLQTWAuuAKwd5/mJgNbAGuA84oY1tzcysAkYDTwGzgbHAKuDYAeucBhxY3F4IPNDGtmZm1gOjmjw/\nl7xgrwd2AjcA5w9YZzmwvbj9IHBEG9uamVkPNCv2M4ANDfc3Fo8N5RPAHR1ua2ZmXTKmyfNq47Xe\nBXwcmNfBtmZm1kXNiv2zwMyG+zPJ99AHOgH4e/Ke/dY2t/UfBTOzzmRlvdAY4CfkB1nHMfhB1lnk\nvflTO9gWXOzLtjg6QM0sjg5QI4ujA9RMW7Wz2Z79r4FPAf9KPrrmGuBx4NLi+b8F/gg4CPhq8dhO\n8oOzQ21rNgxlwCjIXotOYmbl8p59uRZHB2ifRoEuAN0J2gJ6FfQw6DLQ2OBwi4Pfv04WRweomeRq\nZ3KBK25BdID26AzQatAK0IdA00GTQe8C3QX6Mej4wIALAt+7bhZEB6iZ5GpncoGtDDoM9A+gZ0AX\nFu2bgetkoI+ANud7/mbWILnamVxgGym9ryjgnwNNamH9k0EbQFfnLR8zI8HamVxgGwldBnoOdEqb\n200H3Q+6HrRfd7KZJSW52plcYOuUFoA2gY7scPv9QTeD7m7tG4FZrSVXO5MLbJ3QG4pWzKIRvs5o\n0NdAP3DBtz6XXO1MLrB1Ql8Ffbmk1xpVFPy7KzA00yxKcrUzucDWLh1WjJ+fVuJrjgbdUd4fELPk\nJFc7kwts7dL/BH2pC687BbQWdHH5r21WecnVzuQCWzt0ULFXP6tLr/+boOdBB3fn9c0qK7namVxg\na4cuBX27y+/xF6BvdPc9zConudqZXGBrh5aBzuvye0wCbQTN6e77mFVKcrUzucDWKs0GvQga14P3\nugJ0S/ffx6wykqudyQW2VukPQX/do/favzhh6229eT+zcMnVzuQCW6v0KGhe8/VKe7/Pgr7Vu/cz\nC5Vc7UwusLVCbyr2tHs4cZmmgLZ5ZI71ibZqp2cQtG75beA2yF7v3Vtm24DbgQ/37j3NrFXes68l\n3ZtPZdzz9303+cVQSrsQs1lFJVc7kwtszWgqaHt+0LTn7z0K9HR+spVZrbmNY+HOBe6F7N96/9bZ\n68B1uJVjVjnes68dXQ+6JPD93w76qVs5VnPJ1c7kAttwNBr0C9DMwAxZUexPjMtg1nVu41ioOcAm\nyDbERcgE3AK8Py6DWbW42FvZFgJ3RofAxd6sctzGqRUtB70nOkXRTnoBdFR0ErMuSa52JhfYhqKp\noJdA+0UnyenvQH8QncKsS5KrnckFtqHoQtB3o1PsoUWg+6JTmHVJcrUzucA2FP0d6NPRKfbQfqCt\noOnRScy6ILnamVxgG4p+Cjo+OsXe9K38allmtZNc7UwusA1GR4Geq96JTLoQ9K/RKcy6ILnamVxg\nG4wuBV0XnWJfmlQcND4oOolZyXxSlYU4C7grOsS+speBu4EuXwfXzJrxnn3yNBq0BXR4dJLB6WLQ\nbdEpzEqWXO1MLrANpJNBj0WnGJoOLFo5B0QnMSuR2zjWc2cCS6JDDC3bDvwACLiYilk1uNhbGSra\nr9/LzcAHo0OY9TO3cZKmCaCXQZOjkwxv91QOk6KTmJXEbRzrqXcCKyHbER1keNkW4AFgUXQSswgu\n9jZSZ1L9Fs4uNwEfiA5h1q/cxkmaVoJOj07RGk2LuxC6WemSq53JBbZdNA20DTQ2OknrdA/IFzWx\nOnDP3nrm3cD3IdsZHaQNNwEXRocw60fes0+WrgFdEZ2iPZpefBtxK8dSl1ztTC6wQT67pZ4BvSU6\nSfu0FHRBdAqzEXIbx3riaPLfnyeig3Tg27iVY9Zz3rNPki4HXRudojO7Wznjo5OYjYD37K0nzgES\nvShI9jywivy/wcwKC4G1wDrgykGefwuwHPgl8JkBz60H1gArgYeGeH3v2SdH44rx6gdHJ+mcLgf9\nY3QKsxEotXaOBp4CZgNjyfeGjh2wziHAHOBP2LfYPw1MbfIeLvbJ0Rmgh6NTjIymFxcjdyvHUlVq\nG2cuebFfD+wEbgDOH7DOZmBF8fxgKnZNUivBOcD3okOMTPY8+bfOs6OTmPVCs2I/A9jQcH9j8Vir\nRD5vygrgkvaiWYWdTfLFHsinPfbZtNYXxjR5fqQtlnnAJvJWzxLy3v+yQdZb3HB7abFYJelg4Bjy\n4zSpuwO4CjQKstejw5g1saBYOtKs2D8LzGy4P5N8775Vm4qfm4FbyNtCzYq9VduZ5FMk/Co6yMhl\nT+Vz3HMi8Eh0GrMmlrL3jvDV7WzcrI2zgvzkmdnAOOAi4NYh1h3Ym58A7LqgxUTyr/4/aiecVVJd\nWji73InnuDcD8n8IT5AfqL2qeOzSYgGYTt7X3w5sBZ4BJgFvJB+9swp4tGHbgTwaJxnKQM+Cjo5O\nUh6dDbovOoVZB5KrnckF7l86HvR0XvTrQuOLyxW+ITqJWZt8Bq11TdHCyWr0Bzr7JflxpHdHJzHr\nJhd7a0fCUyQM625c7M26rkZ7iXWm8aAdoCnRScqnt4NSnL3T+ltytTO5wP1JZ4Luj07RHRoFehF0\nRHQSsza4Z29dUdcWDsUJVffiVo7VmIu9tapu4+sHct/erMvcxqk8HVbMENnsjOuE6Zj8MotmyXAb\nx0p3FnAPZL+ODtJF64D9QLOig5h1g4u9teJsatuv3yUT8ENgfnQSs25wsbcmNIp8z35JdJIe+CH5\nTK1m1gXu2Vea3gF6MjpFb2guaHV0CrMWuWdvpeqDFs5uK4E31vPEMet3LvbWTN2HXDbIdgIPA6dF\nJzErm4u9DUMTyS84szQ4SC/dh/v2VkMu9jac3wL+L2Q7ooP0kEfkWC252Ntw+qiFs9tyYA5oXHQQ\ns7rxaJzK0uOgOdEpek+rQKdEpzBrwqNxrAyaBRxMf16I260cqx0XexvKWcBdxYyQ/cbF3qwL3Map\nJN0I+t3oFDE0E/Tzel1r12ooudqZXOD602jQL0AzopPE0c9Ab45OYTYM9+xtxH4T2ATZs9FBAt2P\nT66yGnGxt8GcQ/8NuRxoOS72ZqVyG6dytAx0TnSKWJoLWhOdwmwYydXO5ALXmw4A7QBNiE4SS+NA\nr4AmRycxG4J79jYi7wIegOzV6CCxsl+Rz4I5NzqJWRlc7G2gM4G7okNUxAPAqdEhzMrgYm8DvRu4\nOzpERfggrVmJ3LOvDB0G2pqPszfQ4aAXfXKVVZR79taxdwFLIXstOkg1ZM8BrwJvik5iNlIu9tbo\nPcA90SEqxq0cqwUXe2vkfv2+luODtGalcM++EvRG0Cb3pwfSqaCV0SnMBpFc7UwucD3pk6BvRqeo\nHu1XnFw1MTqJ2QA+QGsdcb9+UNn/A9YAJ0cnMRsJF3ujaN24Xz80H6S15LnYG8DxwCuQrY8OUlE+\nSGvJc7E38F59Mw8Ap/ngtaXMxd4gL/bu1w8p2wD8CjgyOolZyjwaJ5SyYkqAPr4EYSv0bdDF0SnM\nGng0jrXlzcDLfX4JwlYUrRyzNLnY2zzgh9EhEuCDtJY0F3ubB9wXHSIBjwDH+gpeZp1zzz6UngSd\nEJ0iDXoQdEZ0CrOCe/bWKk0DpgGPRSdJxH3k34TMkuNi399OB5Z7/vqWLQPeGR3CrBMu9v3N/fr2\n3Aec7it5WYpaKfYLgbXAOuDKQZ5/C/lIhV8Cn2lzW4s1Hxf7NmQ/B54H3hadxKxso4GngNnAWGAV\ncOyAdQ4B5gB/wt7FvpVtwQdog2h/T93bCf096IroFGaUfIB2LnnBXg/sBG4Azh+wzmZgRfF8u9ta\nnDnAY5C9Eh0kMcvIvxGZJaVZsZ8BbGi4v7F4rBUj2da6z/36zhQHaT0pmqVlTJPnR9JiaWfbxQ23\nlxaLddd84NroEAlaD/waOBp4MjaK9ZkFxdKRZsX+WWBmw/2Z5HvorWhn28UtvqaVQqPIh11+MjpJ\nejKB7iGfKdTF3nppKXvvCF/dzsbN2jgryPdgZgPjgIuAW4dYd+DX2na2td56C7AVsuejgyRqV7E3\nq5VFwBPkB1uvKh67tFgAppP35rcDW4FngEnDbDuQR+P0nP4D6OvRKdKlmaDNxTcksyjJ1c7kAqdP\n/wByC2dEPKeQhfPcONbUfDyt8Ui5lWNJcbHvOzoCmEzeXrPO3QOcFR3CLCVu4/SULgL9c3SK9Gkq\n6CXQ+Ogk1rfcxrFhuYVTimwLsAb4regkZq1wse8/LvbluZN8xJmZtcBtnJ7RgaAdoHHRSepB7wD5\n2IdFSa52Jhc4XVoIujc6RX0oA20CHRWdxPqSe/Y2JLdwSpUJ+BfcyrEEuNj3Fxf78t0BnBsdwiwF\nbuP0hMaBXgYdEJ2kXjSlGIK5f3QS6ztu49igTgLWQfZSdJB6ybYBq/EQTKs4F/v+cQZu4XSLWzlm\nLXAbpyf0PZAvC9kVOhG0LjqF9Z3kamdygdOj/Yvx9QdGJ6knZaBnQcdEJ7G+4p697WM+sAay7dFB\n6ikT8F3gt6OTmA3Fxb4/nAUsiQ5Rc7cB74sOYVZlbuN0nVaC5kWnqDdNKIZgHhSdxPpGcrUzucBp\n0TTQNtDY6CT1p9tBH4pOYX3DPXvby/uAJZDtjA7SB9zKMRuG9+y7St8BXRydoj9oBLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x8531710>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "kde = stats.gaussian_kde(data)\n",
+    "\n",
+    "xs = np.linspace(-5, 15, num=200)\n",
+    "plt.plot(xs, kde(xs))\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Monte Carlo Simulations\n",
+    "\n",
+    "\n",
+    "We (well I) want to do this sort of thing because I want to use monte carlo simulations to compute distributions. It is easy to compute Gaussians when they pass through linear functions, but difficult to impossible to compute them analytically when passed through nonlinear functions. Techniques like particle filtering handle this by taking a large sample of points, passing them through a nonlinear function, and then computing statistics on the transformed points. Let's do that.\n",
+    "\n",
+    "We will start with the linear function $f(x) = 2x + 12$ just to prove to ourselves that the code is working. I will alter the mean and std of the data we are working with to help ensure the numbers that are output are unique It is easy to be fooled, for example, if the formula multipies x by 2, the mean is 2, and the std is 2. If the output of something is 4, is that due to the multication factor, the mean, the std, or a bug? It's hard to tell. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 93,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X24VGW9//H3kgeVB8UHhERwe5RQfCi1EFNzi2Wg5whm\nZY+Umhd1MrW0+Fmmdx3tVJfHnvylqHhOmVdopqQmB8XcR34aKgpiKiYqR0BABQUEeb5/f3zXds8M\ns/fM7Hm415r1eV3XXHtm1qw1X2aG79xzr/v+3iAiIiIiIiIiIiIiIiIiIiIiIiIiIjU0FlgIvAhM\nLrJ9PPA0MA94EhiTs20xsCDe9nhdoxQRkW7rASwCWoBewHzgkILH9M25fnj8+HavAHvWMT4RESnD\nTiW2j8KS92JgCzANa8nnWp9zvR/wZsH2qIr4RESkBkol+yHAkpzbS+P7Ck0AngdmABfk3O+BWcBc\n4LzuhykiItXoWWK7L/M40+PLCcAtwIj4/uOA5cBA4AGs73925WGKiEg1SiX7ZcDQnNtDsdZ9Z2bH\nx9wLWIUleoA3gLuwbqHCZF/uF4qIiOSrWTd5T+Al7ARtb4qfoD0w5wmPih8P0AfoH1/vCzwCnFLk\nOdKe7F3oAKrkQgdQJRc6gCq50AFUwYUOoEoudABVqih3lmrZbwXOB2ZiI3OmYn3zk+LtU4AzgYnY\nCdx3gM/G2wYDd+Y8z63A/ZUEJyIitVEq2YOddJ1RcN+UnOs/iy+FXgY+2M24RESkhkqNxpHS2kIH\nUKW20AFUqS10AFVqCx1AFdpCB1ClttABZE3a++xFREKoKHeqZS8ikgFK9iIiGaBkLyKSAUr2IiIZ\noGQvIpIBSvYiIhmgZC8ikgFK9iIiGaBkLyKSAUr2IiIZoGQvIpIBSvYiIhmgZC8ikgFK9iIiGaBk\nLyKSAUr2IiIZoGQvIpIBSvYiIhmgZC8ikgFK9iIiGaBkLyKSAeUk+7HAQuBFYHKR7eOBp4F5wJPA\nmAr2zRC/F/grwY8KHYmISKEewCKgBegFzAcOKXhM35zrh8ePL3dfAF+zaBPHHwF+SHz9IPAe/NfD\nxiQiTaKi3FmqZT8KS9iLgS3ANKwln2t9zvV+wJsV7Nvsnga+F18/OmQgIpJtpZL9EGBJzu2l8X2F\nJgDPAzOACyrcNyvODR2AiGRXzxLby/2ZMD2+nADcAhxcYRwu53pbfGk2G4AVoYMQkdRqjS/dUirZ\nLwOG5tweirXQOzM7Puae8ePK3deViCPNDgG/f+ggRCT12shvCF9Ry4P3BF7CTrL2pvhJ1gOBKL5+\nVPz4cveF5j5B6+NLG/jp4JfrBK2I1EhFubNUy34rcD4wExtdMxXrm58Ub58CnAlMxE7CvgN8tsS+\nWXQY6sIRkYzLQsu+/aKWvYjUSk2HXoqISBNQshcRyYBSffbSLX4YsFvoKEREkqQJ++z93UX66z34\nNfHfH4SOUERST332CdCrxPbeDYlCRCSmZB9GE/6aEZEkU7Kvvy+EDkBERMm+PraFDkBEJJeSfX2s\nBpbH1zeHDEREBJTsG2E68NHQQYhItinZi4hkgJJ9Y/wD+ArwVOA4RCSjlOwbIloJ0W/pei0AEZG6\nUbKvn6jIfS82PAoREVQbp8b8WdgqXV8Cfosl9+0d26MfxfOp9LqLSOY00WzSvDo4B3bymMvBL40f\nM6Kx8YlIE1FtnAR4BqKXutjeI/5brKtHRKTmlOzDGBw6ABHJFiX7mvF7hY5ARKQzSva1c3qFj9cC\n5CLSMBoVUlvLUT+8iCSQWvYiIhmgZB/OYOAu8AeHDkREml85yX4ssBCbIDS5yPYvAE8DC4BHgCNy\nti2O758HPF5NoE3o+8DBwPOhAxER6QEsAlqwdVXnA4cUPOZYYPf4+lhgTs62V7AZpV1pkklV/mzw\nr4FfDn5BF4+7PJ5Q9fmcCVilXiMRkUI1nVQ1Ckv2i4EtwDRgfMFj/gasia8/BuxXsF0nLPNNByYC\n7+bcp6ULRaSuSiX7IcCSnNtL4/s6cy5wX85tD8wC5gLndSfAFBoMHN755mgBRLcAs4ExwH8Dv2pI\nZCKSWaWGXlbyM+Ek4BzguJz7jsOGIw4EHsD6/mcX2dflXG+LL2nzNeycxafKe3j0JvAQ+IWozr2I\nlNYaX+piNNbybHcpxU/SHoF19xzUxbGuAC4ucn+z9Nl78BPBr7TrZe/3PuvnFxGpSE1zZ0/gJewE\nbW+Kn6AdhiX60QX39wH6x9f7Yq3eU4o8R8qTvY/ADwS/0UomKNmLSENUlDtLdeNsBc4HZmIjc6Zi\nQwUnxdunAJcDewDXxfdtwU7sDgbuzHmeW4H7KwkuJXoDr4cOQkQk6dLest85ZwilWvYi0iiqZx/Y\nttABiIgUUiG02psI7Bo6CBGRXEr2NRfNCh2BiEghdeNUr4sJVCIiyaBkX72zc66n/GSziDQrJfvq\nbQcugCiCaHXoYEREilGyFxHJACX7ZBgMvqtSEyIiVVGyD29T/HdC0ChEpKkp2QcXrQauAT4ROhIR\nkXpK8QgW/59xmYRvVnmcsysrsSAionIJjXRkbQ4T/WdtjiMiUpySfXW21vZw3oPfp7bHFBFRsq+C\nPxQ4ug4H7lGHY4pIxinZd99uoQMQESmXkn11VoQOQESkHEr2IiIZoGSfGFGEfimISJ0o2SfPSPAj\nQgchIs1Fi5ckzyxgObBv6EBEpHmoZV+d7fHfd2t83H41Pp6ISHApLBPgT4wnQK2r8XF9fFlb2+OK\nSBNSuYQG2Cv+W68W+FN1Oq6IZFQ5yX4ssBB4EZhcZPsXgKeBBcAjwBEV7CvFpfDXjoikWQ9gEdAC\n9ALmA4cUPOZYYPf4+lhgTgX7QioTm/9kR5dLTY/b3o3zUG2PKyJNqKbdOKOwhL0Y2AJMA8YXPOZv\nwJr4+mPAfhXsm1YDQwcgIlKJUsl+CLAk5/bS+L7OnAvc18190+T60AGIiFSi1Dj7Sn4mnAScAxzX\njX1dzvW2+JJki7HuKRGRRmmNL91SKtkvA4bm3B6KtdALHQHciPXZv1XhvpCf7EVEZEdt5DeEr6jl\nwXsCL2Gt2N4UP8k6DOubH92NfSGdJ2hfAf8Y+NtqfNw54J/QCVoRKUPNc+c44AUsoV8a3zcpvgDc\nBKwC5sWXx0vsWyityf6AOh37JCV7ESlD6nJn6gJWsheRBNAMWhERyadkXzHfl45JZCIiqaBkX7mv\nAHsA2wLHISJSNiX77vkNRK+GDkJEpFxK9iIiGaBkLyKSAUr2ydQK/uTQQYhI81Cyr4gfB1yNlWyu\nlznAk8CudXwOEckYJfvKnADsgtXqr5PoXWBF/Y4vIlmkZF+ZFM72FRFRshcRyQQl+2TqD5xZv/o7\nIpI1SvbJdAA2U/djgeMQkSahZN89UZ2P/5c6H19EMqbUSlVS3Bv1PXz0dfB6b0SkqaRkhIufCt6D\n/36Dnu9G8Oc15rlEJIUqyp1qPZYvAs4DpoYORESkUuqzr8x2iBr1S2QnbFlHEZGqKdkn1/PAGaGD\nEJHmoGSfXDfYH//xsGGISDNQsk+uTcAi4NjQgYhI+inZJ1a0CfhD6ChEpDko2SdbBBwFvo5VNkUk\nC8pJ9mOBhcCLwOQi2w8G/gZsBC4u2LYYWADMAx7vdpTZtRQYj5VVFhHptlLj7HsA12I1WpYBTwB3\nYyNF2q0CvglMKLK/B1qB1dUGmk3RFPDXhI5CquCIgH3jWytwbAsZjmRXqWQ/CjtJuDi+PQ1raeYm\n+zfiy2mdHKPedWQawJ8AfAL4f6EjkQRyTABagFk4/l6wtSf2C80Dg6h7qQ2R4kp14wwBluTcXhrf\nVy4PzALmYrNP0+q/6GidiRQ6D/gxMDrvXscwYDOwFfsFLBJMqZZ9tbNFjwOWAwOBB7C+/9lFHudy\nrrfFlyQJWVaiD3Am8LuAMUhpa4CP4OgPLMDxYHz/UmAY8HqwyKRZtMaXbimVxJYBQ3NuD8U+vOVa\nHv99A7gL6xYqlewl3xzsS1PJPvnGAWcBj+A4AugNeBxen3CpgTbyG8JXVLJzqW6cucBwrD+yN/ZB\nvruTxxb2zffBVlwC6AucAjxTSXACWJLfHjoIKdta4DDgSuCiwLGIvKdUy34rcD4wExuZMxU7OTsp\n3j4FGIyN0tkNS0oXAiOBfYA7c57nVuD+GsbeSIuxpKtFRaR9hE17Ir+hk0etLXLfKTg2A7NxrKhL\nbCKdKKcvekZ8yTUl5/oK8rt62r0DfLCbcSXNBuBRiFaGDkQS4xrgXfJnOffDPveduRrYAxvRpmQv\nDaV69ulwJPgREL0QOhDJsxb7lTsIG3X2JrZ+8Kgu9nmrAXGJ7EDlEtLhGFQnJ4luAL6BDUeehuNb\nWFfl4KBRiRShZJ8eO4cOQAo4LmfH0Wm3APsDX8W6bEQSQd04JfldsHkCIp3pKFTnWAesw7EBmIid\n72lXqk9fpG6U7Es7Hjgaqy8vWeQ4Fvg08Bw2wix3gMLr2MizjQX7rAJuz7lnCjYE+WPx9sOxIoKL\ncMyrU+Qi71GyL8+DEP010HO3jwAaCX4gRKqt0niHY90yDwP3YQn7LAAcY8s6guOy+O9MYE/gdOBc\n4EaskKBIXSnZd8mfCdxB0GFy0Z3gR2Elog9GhbQay1r1Hye/++VdXF6rvVLXYGWr/7ea0EQqoWTf\ntT7x315BoyB6ArwqboZxCjaPZBbWIq+Vb2EzzN9fw2OKdErJvjxaPCTbZmKzxL9Wh2N/+r3iabAG\nx9Q6PIeIkn2ZNMZdaql9SOYdWJmRb2MngDfh6IvjV8Eik6alcfal/R6iJNTifx/5o0Ck3hx7Arvn\n3HM8trZBNa4CvgS04ViOLen5D6wkyW3AL3EsxTGgyucRyaOWfXpcBFwPfl+IXgsdTLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0xaf517f0>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "mean = 14.01\n",
+      "std  = 2.80\n"
+     ]
+    }
+   ],
+   "source": [
+    "def f(x):\n",
+    "    return 2*x + 12\n",
+    "\n",
+    "mean = 1.\n",
+    "std = 1.4\n",
+    "data = random.normal(loc=mean, scale=std, size=50000)\n",
+    "\n",
+    "d_t = f(data) # transform data through f(x)\n",
+    "\n",
+    "plt.hist(data, bins=200, normed=True, histtype='step')\n",
+    "plt.hist(d_t, bins=200, normed=True, histtype='step')\n",
+    "\n",
+    "plt.ylim(0, .35)\n",
+    "plt.show()\n",
+    "print('mean = {:.2f}'.format(d_t.mean()))\n",
+    "print('std  = {:.2f}'.format(d_t.std()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is what we expected. The input is the Gaussian $\\mathcal{N}(\\mu=1, \\sigma=1.4)$, and the function is $f(x) = 2x+1$. Therefore we expect the mean to be shifted to $f(\\mu) = 2*1+12=14$. We can see from the plot and the print statement that this is what happened. \n",
+    "\n",
+    "Before I go on, can you explain what happened to the standard deviation? You may have thought that the new $\\sigma$ should be passed through $f(x)$ like so $2(1.4) + 12=14.81$. But that is not correct - the standard deviation is only affected by the multiplicative factor, not the shift. If you think about that for a moment you will see it makes sense. We multiply our samples by 2, so they are twice as spread out as before. Standard deviation is a measure of how spread out things are, so it should also double. It doesn't matter if we then shift that distribution 12 places, or 12 million for that matter - the spread is still twice the input data.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Nonlinear Functions\n",
+    "\n",
+    "Now that we believe in our code, lets try it with nonlinear functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 99,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8FNWZ//FPcy8gmwqIIIiCW9wSxYXgFlETRbNgTNQ4\njlnM4iRxyTYx6i9aSV6TzBizmxgnmoz606j5RRwTdxOvUaMgCriigKAgqKi4oCBb/f54qm+f7tvd\nVdVbVXV/369Xv6juqur7AHWfPn3qnOeAiIiIiIiIiIiIiIiIiIiIiIiIiIg00DRgPrAAOLvM/unA\nPGAO8DBwuLNvCfBosG9WU6MUia+eazvsXJFM6QIWAhOA/sBcYLeSY4Y42+8Njs9bDIxoYnwitarn\n2o5yrkiq9AvZPxm7qJcA64FrsdaO621neyjwSsn+XB3xiTRLPdd2lHNFUiUs2Y8DljrPlwWvlToW\neAq4FTjTed0H7gJmA1+sPUyRhqvn2o56rkhqhCV7P+L73Ih9jf0ocJXz+kHAJOBo4KvAIXEDFGmS\nWq9tfVOVTOoO2f8CMN55Ph5rxVRyb/CeI4FXgRXB6yuBGdjX33tLzon6SydSq3IJutZre0RwXJRz\ndW1LszWs8dENLMJuRA2g/I2oHZ0fuE9wPMBgYFiwPQS4HziyzM9oxC+Ep/doy/doRAyVrq96ru0o\n51b72Unzkg6gAi/pACrwkg6ggljXV1jLfgNwOnA7NgLhcqz/8rRg/6XAJ4BPYzeqVgOfCvaNAW5w\nfs7VwB1xghNponqu7UrniqRWWLIHuzF1a8lrlzrbFwaPUs8Ce9cYl0gr1HptVzpXJLXCbtBmRY/e\noy3foxExSF89SQdQQU/SAVTQk3QA7SKt/ZrSHpK8vnRtSzPFur7apWUvIiJVKNmLiHQAJXsRkQ6g\nZC8i0gGU7EVEOoCSvYhIB1CyFxHpAEr2IiIdQMleRKQDKNmLiHQAJXsRkQ6gZC8i0gGU7EVEOoCS\nvYhIB1CyF5E4Nge2SzoIySbV/JZmUj37xvo2cE/SQQiQwesrcwFLpijZN9b12N9rfNKBiBYvEZHm\n2Q94CDgq6UAkHiV7EYlqZPC4Bdg+4VgkJiV7EYlqe+BZ4Dl0kzZzlOxFJKqtgFeA51Gyz5woyX4a\nMB9YAJxdZv90YB4wB3gYODzGuSKSHSOBV1Gyb0tdwEJgAtAfmAvsVnLMEGf7vcHxUc+F9hyxUIG/\nI/h/Bv/kpCPpIBqN0zinA78GBgFrUc9A0ho6GmcylrCXAOuBa7GWvOttZ3so9jUv6rmdZkvgOGCP\npAMRqUG+G2cN8CYwOtlwJI6wZD8OWOo8Xxa8VupY4CngVuDMmOeKSDbku3EAlgNjEoxFYuoO2R/1\na8KNweMQ4Cpg15hxeM52T/AQqcXU4CGNNxJ4INh+Gdg6wVgkprBk/wLFM+XGYy30Su4N3nNEcFzU\nc72QONrFTkkH0AF6KG4sXFDl2GnAz7H7S5cB/1Wy/2SsPEAOeAv4MvBosG8J1pWxEeumnFxX1Nmw\nFYWW/UqU7NtKN7AIu8k6gPI3WXfEfhkA9gmOj3outN9NrCr8leD74P8w6Ug6SKXrK8oAggOALYLt\nacCDzr7FWKOmlp+dVQ9jM2gBfgZ8M8FYJOb1Fday34Ddgb8d++W4HOubPy3YfynwCeDTWOtmNfCp\nkHM73c+TDkCA4gEEUBhA4F6jDzjbM4FtS94jR2cZSWEAhrpxMiYs2YPddL215LVLne0Lg0fUczuU\nvwPWgpR0KDeA4P1Vjv88ViYgzwfuwrpxLgV+1+gAU2g48Hqw/TKwS4KxSExRkr00xveBVdg3Hkle\nnK/AhwGnAgc5rx0ErABGAXdikwfvbVh06ZPDhlavDp6/jP3dJSOU7Fvru9jMw6FJByKRBx+8D2u1\nT8M+rPNWBH+uBGZg3ULlkr3nbPeQ3ZFmg4B1FBor6sZpvalkfKRZu93EKsPfCvzbwP9X8M8F/wrw\nN086qg5R6fqKMoBgO6xff0rJ64OBYcH2EOB+4MgYPzuLtsYSfN5ECvc7JBmZu74yF3B8/ingvwr+\n0eB/Fvy3wD806ag6RLXr62jgaSyhnxO8dhqFAQiXYUMN5wSPWcHrO2AfDnOBx51z4/zsrNkRq3iZ\nNwR4h867SZ0mmbu+MhdwfP4p4F/lPL9Hyb5lVBunMfbCCh663kZdkknSSlUi0nDDKNyczVO/fYYo\n2YtIFEOxWcQuJfsMUbIXkSjUss84JXsRiWIYfVv2qo+TIUr2IhKFO6EqTy37DFGyF5EoyrXslewz\nRMleRKKo1LJXyYSMULIXkSjUss84JfvkjA8/RCQ11GefcUr2yViGTcUXyQoNvcw4JftkfD7pAERi\nGkzfZP8KtlSh8kgG6D9JRKIYDKwpeW0d1o8/vPXhSFxK9iISxWCsymUpdeVkhJK9iEShZJ9xSvYi\nEoWSfcYp2YtIFJWSverjZISSvYhEoZZ9xinZi0iYHOVH44CSfWZESfbTgPnAAuDsMvtPxpYrexRb\nePl9zr4lwevu+p0iki0DgI3A+jL7VB+nTXRhizFPAPpjCyzvVnLMAcAWwfY04EFn32JgRMjPaKd1\nOivoswbtZuCvTS6ejqI1aOs3HHi9wr5DgX+0MBYpaOgatJOxZL8E+1S/FphecswDwBvB9kxg25L9\nWn1eJNsq9deDunEyIyzZjwOWOs+XBa9V8nngFue5D9wFzAa+WEuAIpK4QSjZZ153yP44XxMOA04F\nDnJeOwhYgfXp3Yn1/d9b5lzP2e4JHiK1mBo8pHGqtexXYUXSBmDlEySlwpL9CxSX4h2Pte5LvQ/4\nHdZnv8p5fUXw50pgBtYtFJbsRerRQ3Fj4YJkwmgr1ZL9JgoF0Za3LCKJLawbZzawM3aDdgBwInBT\nyTHbATcA/4r17+cNxj7xAYYARwKP1RduWxkI/o+SDkIkgmrJHqwrZ3SLYpEahbXsNwCnA7djI3Mu\nB54CTgv2Xwqcj92tvyR4bT3Wgh+DfQjkf87VwB2NCjzj1gHnYf9OImkXluyXA2OxIdaSUmHJHuDW\n4OG61Nn+QvAo9Sywd41xtbncJvCfRMlesiEs2S9FK6+lnmbQikiYsGS/jL5DriVllOxFJIySfRtQ\nspdOVk8pkLBz24mSfRtQspdO1QVcjCXt3YGT6FsK5FngA1iS/wHw3zHObSdR+uyV7FNOyT5ZY8Df\nPukgOlQ9pUCinNtOorTsx6PSKKmmZJ+cl7DhasckHUiHqqcUSNxzs65SeeO8t7Bh2lu2JhypRZSh\nl9IUuQfAvznpKDpYPaVA4pzrOds9ZLMUSFjLHgr99qtCjpPaTaWOUiBK9tKp6ikFEvVcaI9SIFGS\nfb7fXrPkm6eHOkqBqBtHOlU9pUCinNtOorbsNbEqxdSybzp/c6zkxDNJRyJF6ikFUuncdhWnG0ek\nonZZzacCfxT4b4E/pcy+S8D/cutj6ihaqap+9wGHhBzzBeD3LYhFCmJdX2rZt8ZayD0YfphIKkXt\ns98OAI+RWF2sl4HH8drmQy/T1GcvImGiJPvFwEQ8vovdwL6L/MxjD80lSQElexEJEyXZP8eRbA98\nHxjovH4AcB8eY5oVnESjZC8iYaqtQWs8juRAupxXlgEbg+1tgSvxNMM2SUr2ydsJ/M2SDkKkiuot\ne49BwC+dV24CJgIfxZYtBPgQcHyT4pMIlOyTtQgbwjc26UBEKshhLftq5RJOxeYcwHpWA5/DYwMe\ntwK/co77Hl5R619aSMk+UbmLsJtZImm1GfAuhRZ6MUve/977fB4z8XjNOeJ7wJvB9q5Ya18SoGSf\nDgeD/2Pwd086EJESYTdnPwjBaJsNvMWdvVVCjccq4DfOK59vcHwSkZJ9OpwAfAvYIelAREqEJfvP\n9W69zl94t+wwy8uc7WPw1G2ZBCX75C0EPpx0ECIVVE72HsOBY3ufv8FvgR3LHLeIQgGvfsApjQ1R\nolCyT95Xkw5ApIpqLfsTKYypf5iruA+rFTS8zLFuKYXPNC48iUrJXkSqGULlkTgfd7avxGq1PEu5\n1j38mcKHxm547NKwCCWSKMleizKLdK4hwNt9XvUYRvFCGjOCP58B3lPm+HeAO51XNCqnxcKSvRZl\nFulsgymX7OFIrJY/wDy83mUaHwf2rPBef3G2lexbLCzZa1Fmkc5WvmUPH3G23ST+OLBHhff6K4Wy\nvAfjMaL+8CSqsGSvRZmb73msHOzcpAMRKaNvsvfoR/EIsr8625Vb9h4vAbOCZ13A0Y0KUsKF1bPX\nosz1+1H13bl3gXngayZtY0yljkWZpY8h9B2NsxcwKth+GXjI2bcIGAMMBVaXeb+bgfcH20cAVzcs\nUqkqLNlrUeb67QOcG/HYz4J/K+Q2hh4plfRQx6LM0ke5PvvDnO2/4xWVUtgIPI3dn3uIvv6OlUEu\nfR9psrBuHC3K3BgPRzjmCuA4NBxW0qVcn72bpO8uc061m7QPOe83AY+J9YUnUYUlFndh5SeB6ygs\nypxfmNldlHkOhT65SudKWbk/Uaj/LZIWxcneoxsbfZfXU+acav3267A1bfPUum+RKGvQ3ho8XJc6\n218IHlHPFZHsKO2znwRsHmwvx+bQlHocOLPKe94NHBVsH4YWKm8JdRmISDWlffZTne27KywmPhe7\nV1VpZaq/O9uHaQWr1lCyF5FqSvvsD3G2eyqc8wKwDiouND4HeCvYHofd95MmU7IXkWoKyd5a4FOc\nffdXOe8hbGJlXx4bsAmYeQfWFaFEomQvItW4ffYTKYyvfwMbYlnJLCole/NPZ/uAmqOTyJTsRaQa\nt8/ebdXPLBlfX+ohYP8q+91kr5Z9CyjZi0g1bp+9m+wfDDlvNnaTttIC4zMpzLLfGLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0xb0126d8>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "mean = -1.60\n",
+      "std  = 2.10\n"
+     ]
+    }
+   ],
+   "source": [
+    "def f2(x):\n",
+    "    return (np.cos((1.5*x + 2.1))) * np.sin(0.3*x) - 1.6*x\n",
+    "\n",
+    "d_t = f2(data)\n",
+    "plt.subplot(121)\n",
+    "plt.hist(d_t, bins=200, normed=True, histtype='step')\n",
+    "\n",
+    "plt.subplot(122)\n",
+    "kde = stats.gaussian_kde(d_t)\n",
+    "xs = np.linspace(-10, 10, 200)\n",
+    "plt.plot(xs, kde(xs), 'k')\n",
+    "plot_normal(xs, d_t.mean(), d_t.std(), color='g', lw=3)\n",
+    "plt.show()\n",
+    "print('mean = {:.2f}'.format(d_t.mean()))\n",
+    "print('std  = {:.2f}'.format(d_t.std()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here I passed the data through the nonlinear function $f(x) = \\cos(1.5x+2.1)\\sin(\\frac{x}{3}) - 1.6x$. That function is quite close to linear, but we can see how much it alters the pdf of the sampled data. \n",
+    "\n",
+    "There is a lot of computation going on behind the scenes to transform 50,000 points and then compute their PDF. The Extended Kalman Filter (EKF) gets around this by linearizing the function at the mean and then passing the Gaussian through the linear equation. We saw above how easy it is to pass a Gaussian through a linear function. So lets try that.\n",
+    "\n",
+    "We can linearize this by taking the derivative of the function at x. We can use sympy to get the derivative. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 75,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-1.5*sin(x/3)*sin(1.5*x + 2.1) + cos(x/3)*cos(1.5*x + 2.1)/3 - 1.6"
+      ]
+     },
+     "execution_count": 75,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import sympy\n",
+    "x = sympy.symbols('x')\n",
+    "f = sympy.cos(1.5*x+2.1) * sympy.sin(x/3) - 1.6*x\n",
+    "dfx = sympy.diff(f, x)\n",
+    "dfx"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now compute the slope of the function by evaluating the derivative at the mean."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 95,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-1.66528051815545"
+      ]
+     },
+     "execution_count": 95,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "m = dfx.subs(x, mean)\n",
+    "m"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The equation of a line is $y=mx+b$, so the new standard deviation should be $~1.67$ times the input std. We can compute the new mean by passing it through the original function because the linearized function is just the slope of f(x) evaluated at the mean. The slope is a tangent that touches the function at $x$, so both will return the same result. So, let's plot this and compare it to the results from the monte carlo simulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 112,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8XPP9x/HX3JtFEmvITwiSNGShiKjYooJWg1pKlVZb\nSlE/1FpELGkttVRLSNHlp61qqZaqPbbYagtJrBHZZJHIIrLJcpfP74/vmTvfmcy9M3Pv3DnznXk/\nH4/DzNnmcydz3/fM95zz/YKIiIiIiIiIiIiIiIiIiIiIiIiIiBTRCGAK8BFwcZblRwKTgYnAm8CB\n3rJZwNvRstfbtUoREWm1WmAa0AfoCEwCBmWs0817vHO0ftJMoHs71iciInmoybF8KC68ZwF1wL24\nI3nfKu/xhsDijOWJNtQnIiJFkCvsewFzvOdzo3mZjgI+AB4HfurNN+BpYAJwauvLFBGRtuiQY7nl\nuZ9/R9N+wN3AgGj+vsB8oAfwFK7t/8XCyxQRkbbIFfbzgG2959viju6b82K0z82BJbigB1gEPIhr\nFsoM+3z/oIiISLqiNZN3AKbjTtB2IvsJ2n7eCw6J1gfoCmwUPe4GvAwcnOU1FPbFNTruAirM6LgL\nqCCj4y6gwhSUnbmO7OuBs4AncVfm/BHXNn96tPxO4Bjgh7gTuCuB46NlPYEHvNe5BxhXSHEiIlI5\ndGRfXKPjLqDCjI67gAoyOu4CKkxB2ZnrahwJz/i4C6gw4+MuoIKMj7sAiZeO7EVEClfUNnsRkXJQ\n7QeFbb7qRmEvIqGo1rvxi/KHTm32IiJVQGEvIlIFFPYiIlVAYS8iUgUU9iIirTcL+AJY4U23AieR\n3g/YxrguY+7HjQ3yJ2BtxnbHlqbk+FT7JVUiklu55sRM0kfnSzqJVNhvBrwB/JXUAfZdwC/yfI3m\nfnbdQSsiUiZ6AM/hhmf9PtAYVyG6zl5EQlfso/5Cr+dvbv3uuC4ixgNnFuF1gleuX89EpHy0lBNW\n5KkQs3Dt7Uu96cfAicByXLv80Czb/QlY7W2zMMfPV8j8shVcwSJScuUa9rna7H+KG8hpcMbykrfZ\nqxlHREJXzs0hY4DOuGFZhwPvectKWrfCXkSkbXKF9o24wH8a2B+Ymsc2RaewFxFpm4eBBu/5U8BD\npDezXE0q8IfTuiaj4FXdDywiBavmnNB19iIikh+FvYhIFVDYi4hUAYW9iEgVUNiLiFSBfMJ+BDAF\n+Ai4OMvyI4HJwETgTdLvJsu1rYiIlIFaYBrQB9cH8yRgUMY63bzHO0fr57stVPclVRXMvgL2r2ja\nMu5qJHjVnBMlufRyKC6wZwF1wL24I3nfKu/xhsDiAraVytUT94d+P9IPCEQkBrnCvhcwx3s+N5qX\n6SjgA+BxXMc/hWwrlWs+sDLuIkQkd3cJ+X5N+Hc07QfcDQwssI7R3uPx0SQiUu5mAVsBW+N6t0ya\nCOyK+3Y7G9fSMRrYGzeAyTTgdlxXx/kaHk2tkivs5wHbes+3xR2hN+fFaJ/do/Xy3XZ0jjokKLY9\n0C/uKkRKwIAZwHeB26J5OwNdSB0s7w2Mw3Vp/H3gM2AIcBGFhf140g+Er2xdydl1AKbj/jp1IvtJ\n1n6kenAbEq2f77ZQ3SdeKpSNA5sONgZsBtiX4q5IgleuOTETGAW87s37FXAp7gi+N/ASbhDy1irZ\n4CWHAB/ivnaMjOadHk3g/jq9i/va8iKwR45tM5XrP6K0mo0DOzh6rLCXYmg2JwysmFOBdc0EDsJd\nYj4QdxXiHGA7XNgPAupxXRu3lkaqknJjG0SXXL6msJciK/ewHwVci7u36Elc6DcC+0T/79+qn7rp\nRyxoflbqz16KqTfwMvAObvxNkWpguAtTXgT6An8h1bS9FBf2W+EGLYmNwl6KbRYkvhJ3EVI9EuUx\nLOFs3InaQ4CTvflfAK8A3waej6GuJuobR4rE+uIuNctmMFjvUlYjEoNTcN3FrPbmGe685knAhcDm\n0fxdgb+XsrhyoDb7imBjwWaB3Z8x/0Gw2WA3x1KWVIpyzYmZpPcHltQBN1ThdtHzPYDHgM9x1+O/\nirsMMx86QSvlxMaCndnMsnMV9tJG1ZwTGpZQRETyo7AXEakCCnsRkSqgsBcRqQIKexGRKqCwFxGp\nArqDVkRCUc2XX7aZwl5EQlAOXSIETc04UirdwQodwUxEikRhL6WwFNgTeC/uQkSqlcJeSiDxZ2DH\nuKsQqWYKexGRKqCwFxGpAgp7EZEqoLAXEakCCnsRkSqgsBcRqQIKexGRKpBP2I8ApgAfARdnWX4C\nMBl4G3gZ2MVbNiuaPxF4vS2FiohI+6kFpgF9gI7AJGBQxjp7A5tEj0fgBtJNmgl0z/Ea6tyoIrQ0\nBi2A1YI1lK4ekYpX1DFoh+LCfhZQB9wLHJmxzivAsujxa8A2GcvVgZGISMxyhX0vYI73fG40rzmn\nAI95zw14GpgAnNqaAkVEpO1ydXFcyNeEA4CTgX29efsC84EewFO4tv8Xs2w72ns8PppERCRleDS1\ni72AJ7znI8l+knYXXHPP9i3s60rggizz1WZfEdRmL1JiRc3ODsB03AnaTmQ/QbsdLuj3ypjfFdgo\netwNd6XOwVleQ2FfEfIKewN7Cezo0tUlUrEKys5czTj1wFnAk7grc/4IfACcHi2/E7gC2Ay4PZpX\nhzux2xN4wHude4BxhRQnFaUR2A/37W6rmGsRkRjoyL4i5DqyL3Q9EcmhqJdeiohIBVDYi4hUAYW9\niEgVUNiLiFQBhb3E4ctgw+MuQqSaKOyl1N4DdsRdtisiJaKwlxJL/BY4Le4qRKqNwl5EpAoo7EVE\nqoDCXtrIasDuA74RdyUiUt7UXULQmjo4Ow5shzy3GQD2YfvWJVLxgsvO4AoWX2u6LlbYixSB+sYR\nEZF0CnsRkSqgsBcRqQIKe4lLV7BhYJvFXYhINVDYSxxWAx8DjwB7xFyLSFVQ2EsMErMhMQx4Pe5K\nRKpFrjFoRUrERgLdgZcg8VDc1YhUGoW9xK0/2AzgVGA+0BFQ2IsUmcJe4vQRcB4wInr+Zoy1iFQ0\ntdlLjBJnAufEXYVINVDYi4hUgXzCfgQwBfeV++Isy08AJgNvAy8DuxSwrYiIlIFaYBrQB3fibBIw\nKGOdvYFNoscjgFcL2BbUEVrgWtMRWtr23wR7xJ2ktTFgNxevNpGKVtSO0IbiAnsWUAfcCxyZsc4r\nwLLo8WvANgVsKwKwP6nPjYi0g1xh3wuY4z2fG81rzinAY63cVqrTC8Aw3J2082KuRaRi5br0spCv\nCQcAJwP7tmLb0d7j8dEkZc+OAX7Stn0kluPO+QB2QFsrEqlgw6OpVXKF/TxgW+/5trgj9Ey7AL/H\ntdkvLXBbSA97Ccd2wArgkCLucxjYLyExsoj7FKkE40k/EL6ymDvvAEzHnWTtRPaTrNvh2ub3asW2\noBO0AbPzwH5TxP3tAXZp2074ilSNomfnIcCHuEBPHm2dHk0AfwCWABOj6fUc22ZS2Aer2GEPbb+6\nR6RqBJedwRUsSaUNe4MOBt0Muhb3NUWCpDFopTIY1Bp8zeBOg3dx/eCvBFYZLDJ4zuACc02JIlLm\ndGQfrPY5sq+lrsHg+wYfGlgeU73B3QY7FLcWkbIWXHYGV7AkFT/sl7LJoJfYp6VgX2WwtpllawxG\nmbtjW6TSBZedwRUsScUNe4PjG+GLjAD/3GCMwXCDjaP1agx6G/zI4Nksof+CQc9i1SVSpoLLzuAK\nlqTihL1BwuCqjMCuM7jWYNM8tt/L4I2M7eca7NzW2kTKWHDZGVzBktT2sI+O0m/1g/oDBpjB4AL3\n08HgUoNGb19LTAOaS+UKLjuDK1iS2hb20RF9WtA3kHhsYz43sPvBRuTey3r7/IbBcm+fywx2bW2N\nImVMl15KMEYBZ3nP753CwCOXs8l3cO3z/QrdYQKeBA7E3ehHtJ8nrRX7EpHi0pF9sFp/ZG9wQkYb\n+73mxkBIrjEW7MxWVwa7Rid3k/ufYqlxF0QqgY7spbwZ7I7rZiPpaeDEBBStm4SE60nzcGBNNGsA\n8Pf0Pygi1UNhLyVlsBnwALBBNOsD4JgErM2yeiewDbLMz0sCXsR1u510CHB5a/cnIm2jZpxgFd6M\nY/C3jGvom7nr1W4GWwv2YZurhF96r9lgbegTXKSMBJedwRUsSYWFvcH3Mtrpv5WxyibAt4HrgYeh\n61swaA1s9Cb0eg763wp8Dddldv5Vuj52nvNed170DUMkZMFlZ3AFS1L+YW+wXcYJ02SbfQI36M1D\nuLGKLY/pc+BOCrhpymBrc52nJV//rny3FSlTwWVncAULgF0D9ko+YR/dOOUfWU832Aj4OjCB/AK+\nuelfwMC8KoajMr5ZfKPAH1qknASXncEVLAA2A+wysH1yrgln+23mL8ChwH1kD+8JwLXAcXDsd+G/\ns+CB12HU03DsLGBWlm3qgJ8DnfOo5V6vltkW9bcjEqDgsjO4ggWisP9SzrVgq+guVjOw1+BvwCLS\nw3o1cDPQP2PrAWDzwd4BuxXsEVyzz1dxR/SZoT+JHN0cG/TIaM75bWE/t0jZCC47gytYoICwvycZ\nrJ/Cok7rB/SfgV7NbL012D3RdGl0dc4CsN7RCrsDL2fsbznu+vqWavpuRnPOV/P7mUXKSnDZGVzB\nAvmEvcGBfqgemB7Kn+Da7PN9vQ3AeoLNznjdGuCnuOv0k/tujOY1V1fC4CGvtskGHfKvRaQsBJed\nwRUskCvsDTqZ66LALDo0JzU9DfxPG173L9EJ4qu8BUOAGaS/ztW4Zp9s9W1jbiCUZOD/pHX1iMQm\nuOwMrmCBPMJ+ZDJIl0WH5G42d9CmkaTspuhvx8NgjRkLewCvkB74v6L5wL/MC/vFpLine truncated
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0xb0c47b8>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "-2.33139272541763"
+      ]
+     },
+     "execution_count": 112,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "plt.hist(d_t, bins=200, normed=True, histtype='step')\n",
+    "plot_normal(xs, f2(mean), abs(float(m)*std), color='k', lw=3, label='EKF')\n",
+    "plot_normal(xs, d_t.mean(), d_t.std(), color='r', lw=3, label='MC')\n",
+    "plt.legend()\n",
+    "plt.show()\n",
+    "m*stdb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see from this that the estimate from the EKF (in red) is not exact, but it is not a bad approximation either. "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.4.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}