From d280a9f34422fb7b99121f1db334b91e0d6b8d2f Mon Sep 17 00:00:00 2001 From: Roger Labbe Date: Tue, 16 Jun 2015 21:56:21 -0700 Subject: [PATCH] Added supporting notebooks subfolder. Wrote a notebook explaining how to plot PDFs and run monte carlo simulations and compute their results. --- .../Computing_and_plotting_PDFs.ipynb | 508 ++++++++++++++++++ 1 file changed, 508 insertions(+) create mode 100644 Supporting_Notebooks/Computing_and_plotting_PDFs.ipynb diff --git a/Supporting_Notebooks/Computing_and_plotting_PDFs.ipynb b/Supporting_Notebooks/Computing_and_plotting_PDFs.ipynb new file mode 100644 index 0000000..bb875be --- /dev/null +++ b/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\nAE8Cz3kHKaGdWCtHewdIrqjYF9OHgBu8Q5TYr4APe4cQyRu1cdL3MDahSnyMxEZCTfAOIoUWrnaG\nC5xzWuUyH65BaxJJe6lnX3IfAm5Eq1x6W4ZaOSID6Mw+XTcB53qHEPYDXkaT2qR9wtXOcIFzbDTw\nKpqunxd3Y/sJiLSD2jgldjq2sfgr3kEEUCtHZACd2afnCuDL3iHkbbOAp9FsWmmPcLUzXOCc6ga2\nA5O9g8jbOoAngKO8g0ghqY1TUj3YTkmbnHPILjtRK0dyQsW+OM4GrvYOIbtZBpzlHUIkD9TGad0w\nbOvBQ72DyG66gC3YZDeRNKmNU0Jzga3Aeu8gsps3sdm0mvsgrlTsi+Ec1MLJs18CH/UOIeJNbZzW\ndGDLGWvER351Ye+8pnkHkUJRG6dkjsL+pz/oHUSqehO4FrVyxJGKfXwfwXrCeoeUb2rlSOmpSDWv\nA3gcONY7iNQ0HHgRmOodRApDbZwSORb7H36/dxCpaQdwHfZOTCRzKvaxfRL4GXp3FIVaOZJrC4B1\nWLvgoiqP+W5y/1rg6Irvb8AuHD6ArcY4FBWq5vRPpDrMO4jUrRvYhlo5ko5Ua+cwbKLOdKznuAY4\nfNBjzsB2RgI4DlhVcd9TwPgaz6Fi35xTgNXeIaRh3weWeIeQQki1Zz8HK/YbsJ7jUuDMQY9ZBFye\nfH4vsA+wf8X9Wt61PfpbOBLLFdjetPq7kEzVKvaTgWcrvt7I7kvo7ukxO4GV2Bno+c3HlEFGYguf\nXekdRBp2N7ZV4THeQaRcumrcX+/bhGpnKScAzwGTgBVY7//OIR7XW/F5X3KT6hZi10e0nHE8O7Gz\n+z9HbThpTE9ya4u5wPKKr5ew+0Xa/wA+XvH1Oga2cfpdDFw4xPfVs2/cf6F3SpHNwFbC7PYOIqGl\nWju7sJ12pmMvzFoXaOey6wLtKGBs8vlo4C6G3nxZxb4x44GXqX3hW/Ltv7HrXSLNSr12LsR2QFrP\nrlEEFyS3ft9L7l8LzE6+NwP7x2EN8Duqj0BQsW/MV9GF2SL4IjbuXqRZ4WpnuMCOOrB5C/O8g0jL\n9gVeST6KNCNc7QwX2FH/UFjNfC6GXwBf8g4hYYWrneECO/oBmpBTJPOx2eUacy/NCFc7wwV2MgbY\nDhzoHURS04m9U5vrHURC0qqXBbUYuANbD0eK4S1s6PKXvYOIZEFn9vW5B/gz7xCSugnYUNoJ3kEk\nnHC1M1xgB+/GlqGoNeNZYroc+Jp3CAknXO0MF9jBD4FveoeQtpmLRllJ48LVznCBMzYJeCn5KMXU\ngY3Kme8dREIJVzvDBc7YN7Ehl1Js52PbForUK1ztDBc4QyOB54FZ3kGk7UZji6O9yzuIhBGudoYL\nnKHPATd5h5DM9GLXZ0TqEa52hguckQ7gIeA07yCSmYnYxLmDvINICOFqZ7jAGTkNK/aaSl8u/wJc\n6h1CQghXO8MFzsgKrI0j5TIV2IZWw5TawtXOcIEzcBLwJDDcO4i4uBz4O+8Qknvhame4wBm4Hfis\ndwhxcwQ2Cmsv7yCSa+FqZ7jAbTYPeBwtjVB21zD0ns0i/cLVznCB26gD25v0095BxN0sbNz9Pt5B\nJLfC1c5wgdvodGAdMMw7iOTCZcA/eYeQ3ApXO8MFbpMOYBXwce8gkhsHYyNzJnsHkVwKVzvDBW6T\njwJr0Vm9DPRtNKtWhhaudoYL3AYjgaeAk72DSO7sC2wFDvcOIrkTrnaGC9wG3wCu9Q4huXUhcAOa\nTS0Dhaud4QKnbCrwInCodxDJrW7gEeBM7yCSK+FqZ7jAKbsGuNg7hOTePGADMMo5h+RHuNoZLnCK\nFgGPYT17kVp+CnzHO4TkRrjaGS5wSvbFNhHvcc4hcUzCllE4zjuI5EK42hkucEouB/7VO4SEsxh4\nGK2bIwFrZ7jAKTgXWI9tRSfSiA7gSuC73kHEXbjaGS5wi6Zha568zzuIhLUv8AxwhncQcRWudoYL\n3IJu4G7g695BJLyTsP79VO8g4iZc7QwXuAXfA64DOr2DSCF8HbgXGOEdRFyEq53hAjfp89gwy3He\nQaQwOoCrgJ+g2bVlFK52hgvchFOBF4CZ3kGkcEYD9wNLvINI5sLVznCBG/Re7ILsB72DSGEdBDyN\ntrIsm3C1M1zgBswENmNDLUXaqf+1do53EMlMuNoZLnCdZmEzZM/zDiKlMRsbofMp7yCSiXC1M1zg\nOszGzrK0l6xkbRbwLPAV7yDSduFqZ7jANZyBbTaht9PiZTrwOPBNNEqnyMLVznCBq+jERkRsBt7v\nnEXkAGA18Atgb+cs0h7hame4wEPYH7gZuBPbJFokD0YC/4bN7zjKOYukL1ztDBd4kAXAJuASoMs5\ni8hQPoG1Fv8avUaLJFztDBc4MQ24GngCmzQlkmfvBG4F1qA2Y1GEq53RAu+NXfjalnzULlMSRQfw\nSeA54MfAIb5xpEXRameYwBOBb2Gbg/8U/aFIXOOwtuM2bF2dd7mmkWZFqZ1vy3PgDuBE4EfAduD7\nwDtcE4mkZ1/s3elW4CZspne3ayJpROq1cwGwDhu3e1GVx3w3uX8tcHSDx+at2Hdie3x+C+vHPwx8\nDRvKJlJEo7AJgLdh6zj9OzZfRC3KfEu1dg7Dts+bDgzHLu4cPugxZwA3Jp8fB6xq4NjUAzehEzgS\n+CJwBfZifxi4FNtNKtqklB7vAAXT4x0gYzOAvwXuAF4BliVff4DW183vafF4Gaih2llrGNYcrGBv\nSL5eCpwJPFLxmEXY5tlgGynsg50FH1LHsVkajo2Bn4YtGnVkcnsv9jb2buAu7G3tBp+IqegB+pwz\nFEkP5fp9Pgl8J7lNwEaaHY9tvHMY9k79IeBB7G/56eT2Rh0/u4dy/S5zpVaxn4yts9FvI3b2Xusx\nk7FlV2sd24wubA3vMcnHys/HYxdSJyQfJ7KrwO+HzW59GngU+B1wDfai3ZpCLpGi2YZtbn5l8vUY\n7ATp3cltPva3NQ14Dfvb2oj9Pb2Y3LYCLwGvYzVhZvL5G8nHHdn8p0itYl/v24RWWx33Y2feXcnH\nPX3egb1Q+l8sb1TctrPrRfYI9mLdiG3OvAl4s8WcImX2OtamXTXo+53YydR07ERvIjAJ2x93NnYS\nNhob6/9Bdp2gjcH+nndU3P446Ov+25tYPRrq9tYe7qt2f6vy8jPqVqvYbwKmVHw9BSuee3rMwclj\nhtdxbL/ZNZMONDa5ydAu9g5QMPp9pufAIb43Au2j664LG5EyHRuSVesC7Vx2/atfz7EiIpITC7Ee\n93p27XN5QXLr973k/rUMPEsf6lgRERERESmqXqyf/0ByW+CaJqZ6JrBJ/TZgI7UeAO7zjRLSj4EX\nsGGa/cYDK7All2/BhmlLfYb6ffYSsG5eDPyNd4jA6p3AJvV7CitO0pwTsdn0lcXpUuDryecXAd/O\nOlRgQ/0+G6qbnWknakG0map5Ujn5bQe7JrBJa/SabN6d2Pj6SpUTMC8Hzso0UWxD/T6hgddonor9\nV7ELvD9Cb+8aVW1imzRvJ7AS29rvfOcsRbE/1oog+bi/Y5aiqLtuZlnsV2BvQQbfFmELLx2CLV2w\nGfjnDHMVgff6QkV0PPa2eSHwl9jbaElPWpObyqyhupnlFmWn1fm4y4BftTNIAdUz+U0aszn5uBVb\nVmMO9lZamvcCtm7W89jkqi2+ccKr/P3VrJt5aeNUzqo7m4EXIaS21dhU9OnYBLbF2GqF0pxR7Jqh\nPRo4Hb0m07AMOC/5/DzgWscsRRCybl6BDXNbi70A1MtrnCawpecQbETTGmzBPP0+G/dzbPvDP2LX\nkz6HjW5aiYZeNmPw7/PzqG6KiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiFT3/9pdDcltiJd4AAAA\nAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "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+SbjJezhxXSRJkiSl8CfAs8A3WN3/n4M/I9yX\n2JK6IsCHCf9PXwceAi5LWJfcvpi3C/gS8AxhP3pv2uqschHwFPBo6opUbAYeJOxPJwn341I7TPj8\nngY+A7wqQR0+ASzHOnRtAY4Dp4BjhP+71HXKKQsa+Q3Cf+rG+PznEtal1y7gMeAF8gj+32ZlVNaR\nOKVwEaHrbpbwueVw72Yb8IZYvpTQvZi6Tl1/CvwD8EjqilTMA3fE8gbSB8cs8DwrYf9Z4PYE9XgL\ncDWrQ/YocCCWDzL5466uTrlkQWMPAL+ZuhJ9fA74ZfIJ/qq3A59OtO03E06IXYfilJOHgd9KXQlg\nJ/A4oYGTS4v/MkLI5mQL4WT9GsKJ6FHgxkR1mWV1yC4CW2N5W3w+abOsrlPVBWVBbj/SdgXw68C/\nE76296aktVlxC+ELZ/+RuiJ93AF8PtG2c/9i3iyhhfRE4noA3Au8j9BdmIvdwHeBTwJfAz4GXJK0\nRvA94CPAt4DvAN8nnDBzsJXQ1UJ83LrGsilcUBak+Fnm44QzZa/3E+rzGkIf468SrgB+IYN6HQb2\nVeZNarxsvzrdw0qL8f3Ajwn9oCnk/N2LSwl913cBP0xcl7cBLxH69/emrcoqG4A3An8MfAX4KOGK\n7QMJ6/Q64G7CSfsHhKvt3yV0keUkty9xpc6Cxr4A3FB5fgb42UR16fpFwpn9hTi9TPi9oZ9PWKeu\ndwP/Crw6YR2uY3VXz2HyuMG7EfgiIUBy8CHCldELwH8BPwL+PmmNgm2EOnVdD/xzorp0vRP4eOX5\n7wN/nagus5zf1dNtjG0nn66ed5M+Cxr7Q+DPY/lKwqVebnLp499PGPXw2sT1yPGLeTOEUL03cT36\nuYF8+vgBvkw43iB8uz71aLqrCKOxLiZ8lvPAHyWqyyzn39ztNmwOkeZG6iyr65RLFjS2EfgU4R/1\nVfK6JO56njyC/zTwn4Sug6eA+xLWJbcv5l1P6Ec/wcr/z/413zFZN5DXqJ6rCN08OQ0HPMDKcM55\nVkb6TdL9hHsMPyZcrf0B4dh/nHTDOXvrdAd5ZYEkSZIkSZIkSZIkSZIkSZIkSZIkle3/AIfkyzOW\nVo4xAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "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//TuiCje0A4bH7CvDXwHdGqOHP\ngK8XNSzbDBwC/gN4gfC3S6GulPKgb/cR/pizxfKWiLWUzQH/CJwjjbD/GVauoPrdYoqh3w/MTdpW\n4L3F/E2ELsQU6gL4DeAp4NnYhRSeBPYW8zPED4p54A1WAv5vgV+MUMePE74FoByqvw98tJjfT5zn\nXVVdqeRBI38H/GTsIir8PbCDdMK+7OeAv4r0u+8ivAgu+1gxpebzwE/FLoJwhdphQqMmhZb9LYRg\nTclmwovzbYQXny8APx2plnlWh+op4PZifmuxHMM8q+sqq82DVL718l3ATxCu5FkC7oxaTbCHcKno\nydiF9LCXlaugJq2fD9vFNk9oCR2JXAfAp4DfJHQHpuAHgSvAnwP/CjwB3BC1Ivgf4A+BN4GvAd8g\nvECm4HZCFwrFz9vX2TeW2jyY5LdeHiK8Knb7eFHHbYR+wx8jtPTfGbmmA8D9pXWT+uBYr5p+i5VW\n4ceBbxP6NWNI/Rrjmwj90b8G/G/kWn4W+C9Cf/3uuKVcN0O4au4R4J+BPyK8M/tExJp+CPh1wov0\nVcK76ocIXV8pSfFDVbHzoJEvAfeWls8A3xOpFoAfIbyCnyuma4Tv+Pm+iDUt+xDwEvBdEWvYxepu\nnAOkM0g7CzxPCI4U/A7hXdA54DLwTeAvolYUGhPnSsv3AF+MVMuyDwKfLS3/PPAnkWqZZ203znID\n7PtJqxvnQ8TPg0b2sfIRwx8mvJVLSSp99guEqxW+N3IdqX5gbgMhSD8Vu5Ae7iWNPnuAfyI81wAW\niX8F3HsIV1B9N+FxfBL45Ui1zLN2gHa5MfMx4g2EzrO6rlTyoJFZ4C8Jd+RfSOft7rI3SCPsTwP/\nSegSOAb8acRaUvzA3D2EfvHjrPyNFqJWtNq9pHM1znsIXTgpXbb3UVYuvXySlavzJulvCGMG3ya8\nI/sw4bl/mLiXXnbXtZe08kCSJEmSJEmSJEmSJEmSJEmSJEmS0vL/VmDxdj504cEAAAAASUVORK5C\nYII=\n", + "text/plain": [ + "" + ] + }, + "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+4ptgosFawFrAzQqeRTimpdqplLyWyKN6JgGsg\nui5kGum0BUADcD2wbdgoUg3qs5dSDQU2AZeEDiJdEa2HqAV4P3QSqQ617KUzVkA0LHQIESmeWvZS\nBDsY7H7fGBM6jYiUTi17KcZAYBjQCuwUOIuIdIKKvRRrFfBR6BAi0jnqxpFSfSt0ABEpnVr2Uorv\nAQNQC19EOkGTqhLPPgP2SOgUUgl2K9gbYAtDJ5GS6YHjIlK0ScARwO6hg0hlFVPsxwELgUXApR28\n54b4/Dzar4+9DJgPzAGe6XRKEamQaDWwJHQKCa8bsBifVt0DmIs/qCLb8cCD8f4Y4Kmsc0vxPt4t\nUTdO4qkbp77Z1mC6D1N7ytqNMxov9svwMdZTgZNz3nMScFe8/zTQn/ZjsSNERCSoQsV+CLA86/WK\n+Fix7zFgJjAbOL/zMSUcGwnsHTqFVFxPsAVgHXXVSo0rNPSy2D8TOmq9Hw68hs/AnIH3/ed7ZF1T\n1n5zvEky3IcPtWwOnEMqpxUYCZwFXAd2LbA3RK+EjSU5GuOtIsbiDyduM4nNb9LeBJye9Xoh+afU\nXwXke2yd+uwTzRaD7RU6hVSDRWBbgS0B2zN0GimorH32s/E/4RuAnsBpwP0577kf+Gq8PxZYC7yO\nr5HdJz7eGzgOeK6UcCJSTZFBtAkvIjuAfTreBoZOJtXxWeAl/EbtpPjYhHhrMzk+Pw84OD62Bz56\nZy7wfNa1udSyTyQ72yfaWKta9mljr4BdAPZmPOHq30InkrxqrnbWXOB0sG+B3QW2L1jP0Gmkmj4u\n9g+DPaBin1h6LKGUzZsQaRq9SB3QcgmSxa4HWwn2j9BJJLhrQweQ8lKxl2z9gBuBwaGDSFCHAvsA\nZ4QOIuWjbhzJ9XboABJatCqzr1tq9UItewHsQLAPgfGhk4hIZajYC/gM6EXAdsAdgbOISAWo2Eub\nTRB9CGzAVzi9InAeSY4pYAY2FSx3IUSpEUlYkdJIRo6UsePJjLj4OXAxRAf5lPmPH2TxNkRrgsST\nhLDB+NpWo/EVbwcDa/BlzZ+DqDlcttQrqXYmociq2AdhZwBfAXYG+gLveLEX6YiNAQ4AvoAvh7IU\norPCZkq1mqudut0fhJ0Bdk98c/ZwsFGFrxFpY2f5DGsJSDNopRTRvNAJRKTydINWRCQF1LIXkS6y\nHYFx8Ysn9dCTZFLLXkQ66wywtcAXgV8AVwNHhI0kHVGxTxW7BGyFb/QPnUZq2r34kMwl+DMvlgCP\nBU0kW6Riny598R/SPlTwWZaSBtF6iNYCtwMvAhqZk3Dqs0+ft4FbgG3wx06KdEE0ObNvdwaLITVB\n4+wrznYA+zPYq2DfD51G6pXdCfYnsB+GTpISNVc7ay5wbbHe8aMF3wT7FFhD6ERSr+wEsCvB3g+d\nJCVqrnbWXODaYP3BDgV7DOwDsJdCJ5I0sN4q9lVTc7Wz5gLXBjsG7G2wJ8GOCp1G0kLFvoq0XIJ8\nbBZEx4QOIWlk++CLprX5J0SPh0ojKvZ1yp4ERqDRNhJGT/xZxg3APHxl1bcBFfuUUzdO2dg4sHvB\n3gX7V7DtQyeStLGeYFfH22fiYyf4KB0pM3XjpNhwYADwNeB5iN4JnEdSJ1oPXBk6hWxOxb7m2SeB\nvYFTgB3xpwfdGzaTyGZ2ATsdeAyi10KHSSMV+5piI4FT4xe/Bt4BLgRGAhuBqcDcMNlEOrQSeAn4\nGTAR7FGI3gicSQJQn33R7DSwBfFM2AvANoKtAjs7dDKRwuxesDVgS0MnqRMl1U4thFZ7ngdejfdf\ngGgQRHcGzCNSpOhLwCHtj9lssOe0xELlqRundp1a+C0iidQA1tYqNeA6YAjYrsDWwDqIloP1An4c\nv28mRA9UP2r9ULGvTTcCuwBrQgcRKdEyoAewGzAoPjYSOIHMX6xPA2Px8fpfB2YAHwAq9l2gYp8Y\nNgboDayF6Nk857cDevl+NK2KwUTKKDJgA/BKvAH2ITAt3hYAl2Zd8BHwBP4sBukCFfvkuB1v8bxG\n/geLLMPXoP+f6kUSqYZoLh+PIrOxQaPUMRX7oD7ut/xl/PEXwKlg1+KPehuC92kOjM8P1JA1SQfb\nBe+qlDJRsQ/vIuBEMn+mDgSOwvsn1wNrgWfiTSQN9gNa4v08XZoAdi6wF7AGouurE6u2FVPsxwE/\nB7oBt5K5O57tBrwl+gFwNjCnhGvrnPUAts06sAmid7NeLwE+xBctWwq8HB9/GKInqpNRJDFexGeD\ngw8tXgk2Cfgi2Iis9x0GxAv+WV/gvPj4bhC1Vi9u/egGLMZXr+uB96uNyHnP8cCD8f4Y4KkSroVk\nTqpqLN+nsi+CbYjXlrd4mw32erwfhclVNo2hA+TRGDpAHo2hA+TRGDpABxrbv7ThYJ/Lsx0Gtgjs\n52BXxT9n3wL7RrzF3Z+2VWYrV6ZEKOukqtF4wV4GtOLT8U/Oec9JZJ4s/zTQHxhc5LVJ1bjl0zYA\nbP9465V1vBfYOfFzOF8C+zuwA3AfRP3wpV4n478QLwf2LW+uIBpDB8ijMXSAPBpDB8ijMXSADjS2\nfxm9DNF9m2+sAnYFzsCXUJ4M7I5371wJHAm2B/AbfDmRjWB/zPxstv0y6Eym2lOoG2cIsDzr9Qq8\n9V7oPUPwmyuFrg3MRuPfLL2BZnytmd4wtB+sGIQPAftX/Bvpfvwb6ct4Ae+B96n/O9jhwKfxoZEj\ngLvxPvaLge/ga3oD0T/xPnoR6brl+M8nwAqIXs+csqH4+lHbAJvwlWAPj0/OxteU+l7cOI6AP+Jd\nSB8CfbK+xnvA1nFN2C0+9gZENfc0rkLFvtg/E0rpiugim4j/1u6L/w9qMxLvMloFbIf3k7edPzH+\n+D5e2LO9S/v/ucTdf9/EJy39DvgP/BtlLnAP8Cf8BtI04L74osuBN/Glhf8aZ52J/1JYXfp/p4hs\nWbQO+HsH5z4P1g2fkQuwHqJbMudtKN54643XjotzPsE/8IZgbPxG4POZY7Yq5/0D/WswA2/sjoqP\nr8IbfqPJTCJ7Fh9ifRz568MgYB0fz6vhNojOo8LGAg9nvZ5E+wkPADcBp2e9XgjsVOS14L9QtGnT\npk1b6VvZdMdnuTXgU5cL3aAdS+YGbTHXiohIQnwWX4t6Md46B5gQb20mx+fnAQcXuFZEREREROrZ\nRfjd8OdJ1uSrifjd/AGhgwD/if8bzQOmA/0CZhmH359ZRP57MSEMA/4PX0zreeAbYeO00w2fbJiU\nlRv7A7/Hv59ewLtgQ5uE/797Dh8IsfWW314RtwOvxxnaDMBvvL4MPIr/2yUhV5LqQdGOxv8xe8Sv\nSxn/WknD8JvMS0lGsT+WzNyI6+IthGInzFXbYOCgeH87vAsxCbkAvg1MwYfwJsFdwPh4vzvhC0UD\nPpu8rcD/FjgrQI4j8JE02UX1J8Al8f6lhPm5y5crKfWgJL8DPhU6RB7TgANITrHPdgo+WSSEQ2k/\n0uqyeEua+/D5D6ENBWbijZoktOz74YU1SQbgv5y3x3/5PAAcEyhLA+2LatsIQ/AGxcJqB4o10D5X\ntoL1ICmPJdwbOBIfydPMZo8uC+JkfCLY/NBBOjCezCioautoIl2SNOAtoacD5wB/0PZ38e7AJNgd\nH9t9Bz7m+xbar98UwpvAT/Hx7a/hCwDODJooYye8C4X4405beG8oBetBNVe9nIH/Vsx1eZxje7zf\n8JN4S3+PwJkm4ZMe2lRr4lhHmb5HplV4OT6B454qZcpV1vG9FbAd3h99MT4DMqQT8Yk1c0jOlPvu\n+Ki5C4FZ+GKFl+GTFUPZE5/I2IDPWJ8GnIl3fSVJ2ce3l0HoelCSh/BlfdssxpckCGUk/ht8aby1\n4mv8DNrCNdVyNv7knl4F3ldJxU6YC6EH8AheOJLgGvyvoKXASnwW991BE3ljYmnW68NpPxs9hNPw\nlXHbfAX4r0BZGti8G6etAbYzyerGOZvw9aAkE4AfxPvD8T/lkiQpffbj8NEKOwbOkdQJcxFeSH8W\nOkgHjiIZffYAj+E/awBNhB8BdyA+gmob/P/jXcAFgbI0sPkN2rbGzGWEuxHaQPtcSakHJemBL1r0\nHL7WRWPQNJtbQjKK/SL8ocxz4u2/A2ZJ4oS5w/F+8blk/o3GBU3U3lEkZzTOgXgXTpKG7V1CZujl\nXWRG51XTvfg9g/X4X2Tn4D/7Mwk79DI313iSVQ9ERERERERERERERERERERERERERERERERERJLl\n/wHgMWtmnHhV2wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "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/bHDlRD0O2BfOVLoOI31TsxXel8O4AO3v160C50zQ5MTp+aUX1\ntxevqRlH2mIUYKDHTjAng/kjmGdch8quP/wR+zl5x3USCRXvaqd3gSVUbgcMXLDEdZDcqXkUOgTt\n9q/e7jqNhIba7CVSKu3k3PVOU+RUBwNjG+z8X6e4zSK+UrEXn5UC4+zs5RucJsmtrTBuv5196QNu\no4i0nppxpLWC9vqKbWD+y3WYHDsbMNBzn+sgEhpqxpHIqLSTQSudpsiPhVBaA1s6AX1chxH/qNiL\nx8on2unnz3GbIy8aoPeyYL7SZRCR1lIzjrRGKXRsAAxU3wwm5jpQ7l38a+zn5Reuk0goZL12TgSW\nASuw3dySDQJeBWqBWzNcF1TspXWC9vryAu6Fk+wn07Cfl3ddJ5FQyGrtLAZWAv2wPR8WAYOTlukO\nnAV8i6OLfTrrZj2wRMWwHwEGTv+T6yT5c/CTUF6H/cyc7DqNOJfVA7QjsQV7DVAPTAeS+/m+BywI\nns90XZFWOv4GOx3wqtsc+dTRwDCNkyOtkqrY9wYSfyZvCB5LR1vWFTmWUpjfwc4+F7H267FbgxmN\nkyMZSVXs29LEouYZyZURcCAG7VYDERvn/coT7LTX1UAEDkpLtpSkeH4jR/fp7YPdQ09HJutWJcxX\nc/jqPCLNCvZqu8yzrYiRUQdndYbjGmFTCfZ42GrHmSR/Kslh810JsAr7pmpHywdZwRbsxAO06a6r\nXwCSqdmAgbO/6DqIGx3io2Be7zqJOJX12jkJO7TqSmBa8NiNwQ3s1XPWA3uAXcA6oNMx1s15YClU\npg/s+xu0N4CBh892nciN7ndgPze/cp1EnPKudnoXWFwxA2HmduwYMRvBdHWdyI3xY7Gfm/Wo3T7K\nvKud3gUWV8xAuGUn9j3zA9dp3Nl4PFQ0Yf8Op7pOI85oIDQpZNVlwcwcpzGc6mWOHKfrdpHLJOIP\nFXvxyJPHwZL22JP05rpO49ABaP9bO1tb6TSJSAbUjCNpGncz9v1S7ThIGAwBDHTcBbu7uQ4jTqgZ\nRwqRKYF144M7f3YaJRyWwnGHoKYrPDLBdRiRdGjPXtJgJsJpQZdLojB+fTqeAAyM+47rIOKE9uyl\nED1zAiwH2Ae87jhMWAS/cFZE9HwD8Y327CUNl9+Hfa/8wXWSEOkPGCjeS+qhT6TwaM9eCtHiM4MZ\ntdcfsRr6HILGzthrSoi0SMVefBCDjcOD+Qj3r2/O+D3BzIVOY4ikQc04kkrQzfD4WjQ8QJLHlmM/\nQ391nUTyTs04UnDOt5PRO9DOQZKP7oVYEzAa6OI6jYSXir144Pipdjphh9scYVSxOeiJWgz/8mUw\nQ10nEmmJ9tTkWDpCcXCR7XlfcB0mnG5cAhi4sQHMo67TSN54Vzu9Cyx5NREw0G9byiUj69lZgIGT\n61TsI0Vt9lJQJtnJmLVuY4TZ+U9BWR2sLYUFnVIvL+KG9uzlWILeJo//n+sgIRcMnfAf81wHkbzx\nrnZ6F1jy5hTAQPuDUPtj12FC7lrAwIiNroNI3qgZRwrGRDsZuCy47qy07DnAwJs9gXLXYSR8VOwl\nxE79Nzv9RKPbHF7YBj1WQ30RfPhqMP1dBxJJpj02aU4HaN8AGFi/DcwDrgOF3/jHAQPX7gOzxXUa\nyTnvaqd3gSUvLgIMnLIHzLNgPus6kAfOwo6CuREaVewLn3e107vAkhc/Bgyc87jrIB4pArYABv6+\n3XUYyTkdoBXvxYBL7ezl850m8UsTMNPOPtveaRKRZmjPXpJcEjTh9GiA2ktdp/HMFYCB0XWug0jO\neVc7vQssuXZnI2Dg4rfBDHSdxjPHAYcgZmDTZjBlrgNJznhXO70LLLliYvZ2ZvzC4pNcJ/JTbCZg\n4KE6MOpzX7i8q53eBZZcMX1hbbzQ7wc6OA7kqxsAAxMbVOwLmne107vAkiumL/xgJ/Y98ZTrNB7r\nATRBOwNf7uE6jOSMd7XTu8CSbeYiMDPA/ArG12DfE592ncpzc7GXcvyU6yCSM97VTu8CS7aZG+wZ\nn9vehRIDNAAnuE7lua8ABkqfcB1Eckb97MVLB+CJclvneRHQSUFt84yd1E8E1OdeVOwlTJ6N781r\nb7Tt3oWhTUBnYILrMOKeir2ExKY3YZYBGoGnXacpDJfFRwu90mkMCQUVewmJB7uBKQb+gppwsmRq\ng50WXYWaciJPxV5C4tkBwcyTTmMUlCFDYPAhaOrM4QvBiLRsIrAMWAHc3sIy/xM8vxgYnvD4GmAJ\nsBBoaUAr9caJvDe/BMVNqBdODnxtE/YzNt11Esm6rNbOYmAl0A8oBRYBg5OWmczhkfYYBbyW8Nxq\noCLFNlTsI+/6X2LfB7McBylAy5di/7YHgU6Ow0h2ZbXr5UhssV8D1GP3DqYkLXMp8GgwPw/oij2D\nLy6WSSCJojmjgxn1wsm6gfXQeRHQkcPDRksUpSr2vYH1Cfc3BI+lu4wB5gALsON1iCT50ihYdzqU\nNqL2+hy586Rg5hNOY4hTJSmeT/dnQkt77+OATUB3YDa27X9uM8tVJcxXBzeJhO3ft9MzlsKCPW6z\nFKQr4J9/Dnd8GJgETRXATtehpFUqg1tOjAaeT7g/jfcfpH0Q+HjC/WUc3YwTdzdwazOPq80+umLQ\nPT7wmXqL5IwpgiGrsX/nm12nkazJau0sAVZhD9C2I/UB2tEcOUBbhj17D6AceAW4sJltqNhH1zmA\ngU57sZ0BJGfumQEY6LICzAjXaSQrsl47JwHvYA/UTgseuzG4xT0QPL8YiL+RBmC/HBYB/0hYN5mK\nfWRd9CZ23PVq10kK3+77oEtwBbCXdF3fwuBd7fQusGRFByivAwycP8Z1mMJnvgf/uhsw8IkNrtNI\nVnhXO70LLFnxSex46ytdB4kG8z14fR222awe2xVT/KYhjiXszCoY+Us7XznbaZRIOas3DDsE+0uA\nj7lOI9GjPfvIeeMgYKBzI6y62nWaaDDHgTkRLn8Q+5l70XUiaTPvaqd3gaWtrq8HDNywG8wVrtNE\nyxuXQIfgQC1DXKeRNvGudnoXWFrL9ICdNVBmAAOL16nY55uZAFM3YD93D7tOI22iNnsJs5812DG5\nxuyHYfWu00TTbfHeOJ9Co4xKHmnPPjK29YRTGgADj64G06A9+3wzE8AYmBz8uir5uutE0mre1U7v\nAkumzKVgFsPDC7H/32vh3ZPBDACjYXfzynQEMxS+cR9goKIG9t0Dpq/rZJIx72qnd4ElU+ZaaFwP\nI5qw/9+fd51IiEHFRsDAz3eAudh1IMmYd7XTu8CSKXMt/OavgIHSndhxk8S9GwAD/fdC7SWuw0jG\nvKud3gWWTJlrYehm7P91S2MkSf61x16LwsC/fNt1GMmYd7XTu8CSqR99A/v/vAfo4jiMHO0m7GiY\na1DvPN94Vzu9CyyZ6rcE+/98r+sk8j7toXsNYGDcLWCusT12xAPe1U7vAktGPgIY6FiHvWKZhM49\nwZfxwFpo2A/mOdeJJC3e1U7vAku66h+AYcEwxje94TqNtGR6N6jYBxh4aKuKvTe8q53eBZZ0ffcV\nwMCJjbDr167TyLFc9UPAQNcDsHOW6zSSFu9qp3eBJS0doPt+wMDX/wzm064DyTEVAfOxA9TpGgN+\n8K52ehdY0nHdPMBAt3Xo+rK+GA0YaNcEf1FXzPDzrnZ6F1hSGXs2lAZny47SRTK8UvF77HWB97tO\nIil5Vzu9CywtMSVwcAyMCQY7O+lp14kkY70gFjS/cZXrMHJM3tVO7wJLS0w3+EmwRx/bCnRznUha\n40N3Yj+X24GejsNIy7yrnd4FlpZMGgGdg6FzudJ1Gmmt2oEwPrh05MkLoOEqMCr64eNd7fQusDRn\n+SgYGjTf8CQQcxxIWs0MhLW7oUvwxf3jPWAudJ1K3se72uldYElmxsBVwYBaxavR+DeeM8eD+Qbc\nERysLW2Eb37OdSp5H+9qp3eBJdl192Pb6euAEa7TSNbEgIcAAx3eQ+33YeNd7fQusMSZUfCLJ6Gk\nEfv/eJPrRJJ17WDYLsBA2RvYYZElHLyrnd4Fji7TDsw9wa0vzHkeOgWFfugM1E5foN68F06sxzbp\nzABKXScSwMPa6V3g6DJlYOrB7Ie31kKP+GUGf4fGQi9wM34JXeM9rX4HlLjNI3hYO70LHF2mDMxB\nePU9OD7e82YO+mkfEY/+EDocAgxM2Q01f3edKOK8q53eBY4uUwazaqE82KPvPg/o5DqV5Iu5Ff5m\noFOwh1/ZCI8Nt101TZntnml+BOY610kjwrva6V3giCqC8jugJPign/kGLOznOpTkkykC8wG47Xoo\n3Y694ImBtw2Y8WC+CmY7mCdcJ40I72qnd4ELnxljL01nKoIHjofYH7H/Vwb4Pmqjj7je/aDXJmy3\nzEY477+h7mtgVqrY5413tdO7wIXPrAFTAzXnA9cD72EvbGHgyT1us0mIdAJ+xeGdgNNWwNK1KvZ5\n413t9C5wYTKngzFgGqBxDbywET60m8Mf5PEGlr/tOKSE0iNvQfu9gLHNfFeuQNcbzgfvaqd3gQuT\nOR3q34U/NMHIoE81Bipq4Uu/hPqPglE/emmGeQm2GLjOQCze1LcfeABm/zeYb4M52XHIQpT12jkR\nWAasAG5vYZn/CZ5fDAzPcF0VeydMOZgusKsL9B8FfR6BD8S7Uxoor4XPvg1LLnedVMLO/BOYj4Dp\nBY/cChfUcuT4joHx9fCDtbBts+3RI1mS1dpZDKwE+mHPmlsEDE5aZjIwM5gfBbyWwbpZD+y5yjxs\nowg4He5bBJ9pgN4JH0oMsBz4Cu67VFY63n6YVLoOkD7TG8xl8NMvwpBZUBT0y8dAsYGhO+Gjz8Hd\nVUCXoIdPUQa/GitzldxDWa2d5wDPJ9y/I7glehCYmnB/GXbApHTWBRX7RFXZe6lnB8HXPgqDPw58\nATr+BE5+K+EsyIRblwNw3t9hypWEZ8iDKtcBQqTKdYA26Ab/tRIubIKi+BnXCbf+TXCFgbsa4MHN\n8P17gb60fN3iqvzE9kJGtTPVKc+9gfUJ9zdg995TLdMb6JXGunLEKGAA9pdSKdAumCbdLhgOxWXw\n7h7YEQPTGcp6wslnwJ5i2AIc6ASXJBXtGmBtMN+jyW5u92Pw8v2wZyG8qC9dyYVd8OWL4cvl8P0d\ncNsw6DcVulwFb5fC6hisBp4qxu4kTrO3mIHYbijeDcW7oGQ3FO2EQwOhay9o3Afbt2Pf2PFbY9Kt\noZnHGoH4l05cS/PHeq6t6+RdqmKfbriw7A16xNwC/BDYA8yFvhNgfRnwqWOvN7uZx3Ziv0sTFdfB\nCXUwtBP0qYGmv8K+OfD738KW3kBH4E2I7WjzP0XkmGLLEu6sgTUzwFwDv+sMd58OyweBGQKDPgyx\nQbCjC2wrAtMNmrpBff+jX2/rGflMXyhSFfuNQJ+E+314f1VJXuakYJnSNNaNi+Be5eHvxy7Axdl/\n/cZ2sLUdbAVb2C8Ibt/16Lv5btcBQqTA/hbNvQeXNfOY5EsJsAp7kLUdqQ/QjubIAdp01hURkZCY\nBLyD7VkzLXjsxuAW90Dw/GKOvlJRc+uKiIiIiEghuQp4C3t0PPm6pdOwJ2ItA6J2Vfsq7LGNhcFt\notM0bqRzMl5UrAGWYN8L891GybtfYA86vZnwWAW2l8JyYBbQ1UEuF5r7W1ThSa0YBJwGvMjRxf6D\n2Pb9Umx7/0qiNcLi3diTmqIq3ZPxomI1tsBF0bnYM/ITC9z3gNuC+duB7+Q7lCPN/S0yqhUui+gy\n7LdzsinA40A9dq9mJTAyf7FCwZvuMjkwEvt/vgb7HpiOfU9EWVTfD3OBXUmPXQo8Gsw/ClyW10Tu\nNPe3gAzeG2HcY+7F0V004ydpRckXsAe7f050fqbGtXSSXlTFL/24ALjBcZYw6EHQnziY9nCYJQzS\nrhW5LvazsT87km+XZPg6hdYPv6W/y6XAT4H+wJnAZuA+RxldKbT/67Yai/35Pgn4HPbnvFjxYRei\nKqNakesrxF/QinWaO0lrY3bihEa6f5efAc/mMkgIpXMiX5RsDqbvAU9jm7nmuovj3FbssApbgBOB\nbW7jOJX4b09ZK8LSjJPY7jQD+Dj2RKz+wECi1QvhxIT5yzn6gEwULMD+n/fDvgemYt8TUVQGdA7m\ny7E906L2fkg2A7g2mL8WeMZhFte8qRWXY9tma7Df0s8lPPdV7EG6ZcBF+Y/m1GPYrnaLsW/kKLZJ\n6mQ8qz8E0JuHAAAAQklEQVS2N9Ii4B9E72/xOLAJqMPWiuuwPZPmEL2ul8l/i+tRrRARERERERER\nERERERERERERERERERERERERadn/A38U8Gwenr9vAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "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+y7St8BXRydoj9oBugXoDHRSawv\nJFc7kwucDk10H7nX9AjIZ9NaL7iNY7udBTwM2dboIH3ErRyzIXjPvmt0Leg/RafoL5oDWhudwvpC\ncrUzucBp0DjQi6BZ0Un6i0YVFzQ5OjqJ1Z7bOAbA2cDjkD0THaS/ZK8Dt+NWjlWMi319fRi4PjpE\nn7oNT51gtg+3cUqnSaDtoEOik/Sn3aOgpkQnsVpzG8c4D7gPss3RQfpT9gqwDDgnOonZLi729fQR\n4JvRIfqcWzlmA7iNUyrNAG3JJ+ayOJrps2mty9zG6XO/A9wE2avRQfpbtgF4BjgtOokZuNjXjDLg\nY8DXY3NY4XbcyjHbzW2c0uhU0JNF0bdwmgv6cXQKqy23cfrY7wDXFZfJs3grgKmgo6KDmFWBC1Mp\ndk+PMDs6iTXSNaBPR6ewWkqudiYXuJp0Aej70SlsIF0Auis6hdVScrUzucDVpH8CfSI6hQ2kSaAd\noAOjk1jtJFc7kwtcPZpaTI/gglJJuhP076JTWO34AG0fugi4E7Lt0UFsUL6giRnesy+B7ge9NzqF\nDUWzioPno6OTWK0kVzuTC1wtOhr0AmhsdBIbjlaD5kensFpxG6fPXAxcD9nO6CA2rFuB86NDmEXy\nnn3HlIEez8+ctWrTiaCf+uxmK1FytTO5wNWht4KecQFJgTLQT0DviE5iteE2Th+5kHyGS//BrLxM\nwE3AB6OTmEVxoeqYHgN5Ct1kaC5orb+JWUmSq53JBa4GHQfaCPK3s2QoK9pux0cnsVoovY2zEFgL\nrAOuHGKdLxXPrwYae5LrgTXASuChdoJZU7taOK9HB7FWZQJuxq0cq6DRwFPAbGAssAo4dsA65wJ3\nFLdPAR5oeO5pYGqT9/CefUf0KGhedAprl+aD1kSnsFoodc9+LnmxXw/sBG5g37HC5wHfKG4/CEwB\nDm143v3J0ulY8s95eXQSa9v9wCGgY6KDWH9pVuxnABsa7m8sHmt1HQF3kV/E4ZLOY9oAFwI3u4WT\noux14J+AD0Qnsf4ypsnzrX5NGGrvfT7wHHAIsIS8979skPUWN9xeWiw2tAuBy6JDWMduBr4AfC46\niCVlQbF0xanAvzTcv4p9D9L+DfDvG+6vZe82zi5XA58Z5HH37Nui40DPehROyjS6+H94XHQSS1qp\nPfsVwNHkB2jHkU+le+uAdW4FPlLcPhXYBrwATAAmF49PBM4GftROOBvURcCNbuGkLHsN+Cb5NYPN\nKmMR8AT5gdqriscuLZZdvlI8vxo4qXjsjeSjd1YBjzZsO5D37FumrDgpx3PhJE9vBW3wtMc2AsnV\nzuQCx9HbQet9BmZdaCXoPdEpLFmeG6fGLgL+j+fCqY2vAx+PDmHWKy5cLVFWTJF7UvN1LQ2aCtoK\nOjg6iSXJe/Y1NQf4NfnUE1YL2RbgO8DHgoOY9YT37FuiL4D+ODqFlU2ngZ70UFrrQHK1M7nAvadR\nxWyJb41OYmVTBloFOic6iSUnudqZXODe07x87nqrJ30U9L3oFJac5GpncoF7T9eAPhudwrpF40DP\ngU6ITmJJSa52Jhe4t3RAMWJjsCkorDb0h6Bro1NYUpKrnckF7i1dAro5OoV1m6aCtoAOi05iyUiu\ndiYXuLf0EGhRdArrBX0F9KfRKSwZydXO5AL3jk4BPe35U/qF3gTaDJoYncSSkFztTC5w7+gG0Kej\nU1gv6RbQ5dEpLAnJ1c7kAveGZoF+kR+gtf6h04tpMZpdWMgsudqZXODe0BdAX4xOYRH0A9CHo1NY\n5SVXO5ML3H2aVozMmBmdxCLoXNAaT2VtTSRXO5ML3H36POivolNYFGWg1XnRNxtScrUzucDdtXuv\n/ojoJBZJH8rbOWZDSq52Jhe4u/Rl0JeiU1g0jSkO1M6LTmKVlVztTC5w9+jNxThrX8zCAF0GujU6\nhVVWcrUzucDdo+94wjPbQ/uDnvfU1jaE5GpncoG7Q+cUX9vHRyexKtFVoOuiU1glJVc7kwtcPo0H\nrfPoC9uXphQn1/1GdBKrnORqZ3KBy6c/yk+TNxuM/hz019EprHKSq53JBS6X3lYclJ0VncSqSoeA\nXgQdFZ3EKiW52plc4PJoLOgR0Mejk1jV6X+AvhWdwioludqZXODyaDHouz4t3prTJPJLF86JTmKV\nkVztTC5wOXQS6AXQ4dFJLBX6GOhhX9/ACsnVzuQCj5zGg34Eujg6iaVEGej7oE9FJ7FKSK52Jhd4\n5PQV0I1u31j7dGxxQP/I6CQWLrnamVzgkdEHQD8BHRidxFKlPwD90Bc46XvJ1c7kAndOR4J+Djo5\nOomlTKNAS/ID/NbHkqudyQXujMaBHsTXlLVS6DDQRtB7o5NYmORqZ3KBO6Mv5zMYuk9vZdHpxTfF\nY6KTWIjkamdygdunT4LWuk9v5dMnQE+BpkUnsZ5LrnYmF7g9mu+9L+su/XHRIpwcncR6KrnamVzg\n1mlWcdbjwugkVmfKQH8L+gFoYnQa65nkamdygVujKaBV+TA5s27TKNA1oGWgg6LTWE8kVzuTC9yc\nJoHuB/0vH5C13tEo0BdBP/ZJV30hudqZXODhaUpxwsvX8n98Zr2mTxXzLp0XncS6KrnamVzgoWkm\naHWxR+9Cb4F0KuhnoC+Axkansa5IrnYmF3hwOqM4GPtZt26sGvQG0O2gBzwarJaSq53JBd6bxhZD\n314AnROdxmxvGlW0dTaDrvA3zlpJrnYmF3gPnQFaA7ojP33drKp0DGg56G7v5ddGcrUzucCgucXX\n45+BPui2jaVBo0GfIb+e7edBB0QnshFJrnYmEliHgv4j6KGiyF8OGh+dyqx9mg66tjjGdBlo/+hE\n1pFEauceFQ6sw4t+572gbaDrQe/Dl4WzWtDJ5JPzPQ/676Cp0YmsLaXXzoXAWmAdcOUQ63ypeH41\n8I42t61QsVcGOhH034qToraCrsvHK3sv3upKxxV7+ttA/wh6jw/kJqHU2jkaeAqYDYwFVgHHDljn\nXOCO4vYpwANtbFt64PbpEND7yecW2UA+g+Bfgs4G7RebrSMLogPUzILoAL2jg0G/Rz7Nx8+KfweL\nQBNKeoMFJb2O5dqqnc0uazaXvGCvL+7fAJwPPN6wznnAN4rbDwJTgOnAkS1s20PKgGnAW4HjgbcD\n88izLgfuAv4CeBKyCn3baNsCYGlwhjpZQN98ntmLwF/mi04A3gv8V+BG0CPk39yfAZ4HXgA2A9uA\n7cBLkL3W5A0W0DefZfU0K/YzgA0N9zeS7703W2cGcHgL27ZIGTABmDTMMnGQxyYDhwJHFJn+DXis\nWFaQt58ebeGX1KzPZGuANcDnyK/DMId8B2kmcBL5v6uDgQPJd/Amg15mT/Hf1rC8BLwMH5wLN/0+\n8AqwA9gCbC1+bsm387/FbmlW7Fvdwx3h0EM9Rt7qGdPwc9ftscB44JfAyw3LKwPuNy4bG26/ADyb\nL9mOkeU060fZduDuYhmCRpPvXE1hzx+AXbcPBCbC6DHAb+S3OQA4CJja8HMyaAfwKvA68NqA5fVi\nURtLu+sP9xoV017ZbVbsnyX/S77LTPJCOtw6RxTrjG1h20J2XPOoTCgWX5GnuaujA9SMP8+mWi08\nN76zyQq7/khYj40BfkJ+kHUczQ/QnsqeA7StbGtmZhWxCHiC/GDrVcVjlxbLLl8pnl9N3s8bblsz\nMzMzM6urxeT9/JXF4mu2tq+VE9isdevJR6OsBB6KjZKk/00+OOJHDY9NBZYATwLfw735dgz2eS4m\nwbp5NfCfo0MkrNUT2Kx1T5MXJ+vMO8nPpm8sTp8H/ktx+0rgz3odKmGDfZ5t1c0qnRLtmSM713jy\n2072nMBmI+Pfyc4tIx9D36jxBMxvABf0NFHaBvs8oY3f0SoV+yvID/Beg7/etWuoE9uscyI/q3oF\ncElwlro4lLwVQfHz0MAsddFy3exlsV9C/hVk4HIe8FXy6RVOBDYBX+xhrjqo4AkfyZtH/rV5EXA5\n+ddoK8+uE5asc23VzWYnVZXprBbX+xpwWzeD1FArJ79ZezYVPzcDt5C3ypbFxamFF8jnonoeOAz4\neWyc5DV+fk3rZlXaOI2X9Hs/ex+EsOZWAEez5wS2i4BbIwMlbgL5qf+Qn9p/Nv6dLMOtwEeL2x8F\n/jkwSx0kWTevIx/mtpr8F8C9vPb5BLbyHEk+omkV8Cj+PDtxPfAc8Cvy40m/Sz666S489LITAz/P\nj+O6aWZmZmZmZmZmZmZmZmZmZmZmZmZmZmZmZmY2tP8PpVUtjJBiTpMAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "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+IXgsdTEb8GBt1cyXw\nCNbP/mOqWefB8XDB7b8Cx8bXI+zL4FmaYoU3SRIl+/R4HZuWvwwlgvpy3IsV8dsCXIrj+vj+BXV+\nXg8sxamktdSeunFEdrQPtn5DsRPzKoEgqaSWfWpEj4MfBjwaOpKmZl0pEbCtyNaXsf72jUW2iSSa\nkr1IviVYd9kKrHRxB8dbwB8DxCRSNSX7rn2I3CJXkhVDsdb7LsDbgWKIcERxP75I1ZTsu3YB8PPQ\nQUgAjjcDR7AKOAzHc0r4UgtK9l17Hfhp6CCkzmw8/X3xrfC/5Bx74HgO+DtWdO3BwBFJE9BoHMk2\nR1/gI8AIul5OsLEcI4G/Aj/EMQfHF0OHJOmmZC9ZdxBWd/5RqpksVT8HAsOxGdQi3aZkX5Q/A/yb\n2HhraX4v4jgtdBCdGARsDh2EpJ/67IsbA+yFLR69OnAskl1nYWsPfw8Ax8HAadgs29sCxiUpVE7L\nfiywEHgRmFxk+8FYMaeNwMUV7ptU7T/nF0C0JWgkUh+OQTiOAw4t2DKo2MODcLyJYxm2ju0I4F+B\nH8V/RSpSqmXfA7gWGxGwDKvpfTfwfM5jVmHLqk3oxr5J5YFXsXilOZ2GVRF9g47Vv4bFfwOuTNap\nU7HVsmZivzpFKlIq2Y8CFgGL49vTgPHkJ+z2/yyFfZ7l7JtU24FfQDQ7dCBSV/lj6a0VnWR/AaZj\nrXuRipTqxhmCTR9vtzS+rxzV7Cud2w/8fqGDSDXHhdg6sCKZUaplX81QtEr2dTnX2+KL7Kh9VMYS\nVOa4Gv2A3wFzgOtJxyLu/UIHIMG1xpduKZXsl2F1QtoNxVro5ahkX1fmMTMuWgm+F6q6WAvrgKeA\nG4DlgWMp5WHsy13ditnWRn5D+IpKdi6V7OdiEzpagNewoWCf6+SxhS3NSvYVaTzHQuCS0GGU5LgH\nuCe+fmLYYCStSiX7rcD52AiAHsBU7ATrpHj7FGAwNtJmN+zE5oXASGzkQLF9pXo9wHtgIEShC3al\ni+MqrMCdFvWWTClnUtWM+JIrd+HrFeR315TaV2pHM6ArtztW3E7JXjJFyUKywfE+HPsDfYG3cawN\nHVIVjsLxp9BBSLoo2UtW3IktKzg2dCBVmgecDRweOhBJFyV7yZKVoQOomv0iWYCtZNUjvmgYrpSk\nZJ8+HviP0EFIcMOwARRb0Rh8KYOqXu7ARyRhtaJORduAS8BPDB1JKjg+gzVq9o7vSU6hs+57C/g1\nNt/iwsCxSEok4eefJxlxAH5n4Dashs+3IPpF4IC64DcAH4bo2dCRJJpjE1bQbhu2nvCrwFoc64PG\nVSuOddhouDVaqzZzEpQ7y5OQD6j/hI1d9x78D0JHU5r34O8IHUXiOTbh2Dl0GHXj4iTvGINjl9Dh\nSENVlDvVZ98hJyFE/xYujLJ9GRgAvnfoQCQgx+7YFPoZwDFhg5EkU7Lf0U9DB1Cm9cDJwCGhA5HA\nHCdhRd0+h9PnQYrTCdp890D0f0IHUZ7oT+AXhI4ikWwoYoJPstfFNOxkbRsqSyJFqGUvzWgksAV4\nPXQgDeO4DpgPDIy7dkTyKNlLs1oB9CJbn/E3gKuAc0MHIsmTpf8IpdwC7Bk6CKmZrVhLdw6JGfFV\nZ45vAjeFDkOSSX32HdYCF4cOQmpmLY4TQgcRyCgc3wWew3Fv6GAkGdSy77ARWB06iArtib6gZEcn\nYqsYfTZ0IJIcSvbp9nvgS+A/EjqQxLCVnB4OHUZg/bBfqjvhMjcqSTqhbpx0+wnwaWyVsGxzfBlb\nGW0L8CzwybABBfMQdr5iP6xlPwibjyEZp2QPgB8OHBQ6ispFa8AvCh1FQgwG/gnoDczBkc3lGvPX\nq50KfC9oPJIY6sYxY+K/64JG0T19ge+HDiIh5gHXAbeHDiRBRuH4Y+ggJDy17DvcANGK0EF0wy3A\n5aGDSIgncUwOHUSCPIGNub85roi5EjgAx7thw5IQ1LJPv0eANaGDkASyVa0eBe7Ffu0MCBuQhKSW\nvaSf42SgFVuuT3I5ltI+BFMt+kxTyx7fApyEvvjS7F+wEtX/EzoQkaQqJ8GNBX6BVRG8ieIlgH8F\njAM2AF/BTpQBLMbG+27DhsSNqira+jgHOAu4OXQgVRgJfj5EHwwdSMM4hmKf3wHYMMO7cdwXNiiR\n5CqV7HsA1wIfA5ZhJ3zuJr+E6qnYsMXh2OIJ1wGj420e+3md5JmpHnAQ/TB0IFX6QOgAGuwh4ABg\nFfb5ui1sOCLJVqobZxSwCGuhb8FqZo8veMzpwG/j649hLa3cRZ0TvEai3wloCR1FlV7ASvoCfknG\nVq56GngHuETDC0W6VqplPwRYknN7KTsufVbsMUOwYV4emIV140wBbqwm2NryEfbFNBG4KHAwVYi2\ngl8Z39iPRH+51txncGhSWWW+hmM18JfMTjzLqFLJvtzSsJ0lmOOB14CBwAPAQmB2kce5nOtt8aXe\nhmOtYiD6ZQOer57eAc4A7godiCTaNdjn/lPYL3Yl+3RpjS/dUirZLwOG5tweirXcu3rMfvF9YIke\nbFGFu7BuoVLJXioWbQamg98cOhJJMBfPtHZ8APgZjl2Aq3H8IWhcUq428hvCV1Syc6k++7lYS6AF\nqzlyFnaCNtfdWFcI2InZt7EunD5A//j+vsApwDOVBNcAK8k/vyCSFf8E7I/96pYMKNWy34pVEpyJ\njcyZio3EmRRvnwLch43IWQSsB86Otw0G7sx5nluB+2sVeI2shaiZ1intDWwEvz9Er4YORhJtd6xh\nJhlRzjj7GfEl15SC2+cX2e9lIDvjvpOlOSfLOXpi54H6YzXbpXt+jjXGTg0diDSOZo02pwnYRLj0\nc4wDrsTq/xyN1e7/CzbOPo1VSsNz3BH/fT9wFY7+OK4KG5TUm5J9cxpDsyR72AMYhnU7vAwMx9FM\nXW8hXYaN5Po+jsHxguXSpLKc7L9Nc56cehubANdMHsDx+dBBNB3H2rhFvwI4KnQ4Ul/N2bdbkv8Q\ndpL5vwIHUmtzgedCB1Ezji+BuhfqyrEeq18lTS6jyZ5+2KzftNfDKRB9GJs40yz2wer1XxY6kAwY\ni+N2HB8KHYjUR1aTPcDLEGnoWVLZyJuewEocL4cOp8n9FatoeyTprv4qXchqsv88sEvoIOrok+AP\nDh1ElX4C/BirqyT15HgVx2XAmUALjt+HDklqL6vJ/qNY/3Yz+hs2ue0q8EeGDqZKk3F8N3QQGfK/\nwDewNSykyWQw2fuRwAia7+RsLHoNW0DmQGxKfLo4vojjIrJXnz88xxraJ1A6TsTRGtfPkSaQxaGX\nHwA2YsPNmlR0B/jPhY6im87Fho8u4b2qpNJge2AFt7YDn8LxFvB3lUROtywmewfMgKiwemez2RVo\nBT8HomWlHhyMnYj9bHzrQqzExik4HgoXVKa9DXw4vv5nOupbnQFMDxKR1ETGkr0/GXg/8J3QkTTA\nYGxd4BewpSWTqje20ll7l+Lx2PKXEoJjK/BUfP0crNFwcciQpDaSsKqRp+5x+AhLKuuBV4EjIHqn\nvs8Zmh+Hrcv6PYiSm+wdfbBFND6CJfzncGwMG5TkcUwH3gf0Am7C8ZvAEYmpKHdmpWU/CFgeXx/f\n/IkeIJoB/hZgAvg1wO0QbQodVacc80OHIF0aijWYzsGxBsetoQOSymRgNI7vT8dSXucDfw8XS8O9\nDZwM/A5bQCY5HAPiejfqm0+H/sA9wD+AkwLHIt2QgWTPQdC+7Fr0fyEqd13dZvBD4OnQQXTicGy2\nZk/g04Fjka79CTvvMwWbbSsplIE+e38MMMeuR0n49wbg1wLDIVoZOhIAHEcDPwD2wnFC6HCkAo6v\nAv+GjdS5GtiG45WwQWVWRbmzyVv2fjzvJfpMF9PqD6wAH77uiWMM8O9YjfpmKtqWFQ9gv5S/AMzC\nZmxLCiShpVunlr2/GGt5AJwMUYZ/fvpzsPWDAY6B6PGGh+AYBHwNq0v0GvATHDMbHofUjmMwdg7s\nZuz80K9xWj2sgSrKnc2c7Nv75k8E5kC0ufbPkRZ+GHAplmxpWHeW4yxgAPBxrMgWWP/vD3E805AY\npH4cA7GW/YHABmA4Nox2J2BrPGZf6kfJHnxPOlZrGgnR87U9fhr5ccB98Y3LIKrPoiCOAcC/YlPt\nv4NNvY+AnwHX4lhSl+eVsByvYWPx203CcUOocDIi68neDwFewSaAABwC0cLaHT+t/D7AKVgLf2RN\nW/cdtef3As4CfgT0wb5w98WxqmbPJcnk2Atb/a0/Ngv6ReCPWB2qN/VLri5qnuzHYotX98AWOPhp\nkcf8CpuavwGruDivgn1rkOz9+Vg/8IPA/2DFzkZD9Fh1x21Gfh+gfVSOA24AVndrwpUtG3gsdrL1\ntJwtvwbav2Bv1ozYjHGcC/wcS/xg5cTvAtbh+HWwuJpPTZN9D6y2yseAZVjNks8Bud0ip2KTlU4F\njgF+CYwuc9+KA96RHwC8BWylY0bwGRA1qmhTK1YhMCX8Ttj7EC9Q0UY85+yjEM0u6xCO/bFzIV+l\noxTxvdgJ8Q24hlarbCVVr/8OWklv/K10FrtjBFaH6oPYSLj2/58vA7cAv4zXvw2plfS+9lDjcgmj\ngEXA4vj2NGA8+Qn7dKyQFcBj2Am5wcABZezbDX5nrIvmIGwt2fYRHT2xL5XDGrzcYCup+sBE24Fb\nwU8FFkHboXGyfxj8A8DOXLT/bxjw6hi7j7uxL+/dsJmT38g52FPYez8AeBD33i+6RmolVa//DlpJ\nb/ytdJ7sXwBewPEEMB/YhBVV+yS2iPwP4sqmz2DndbZjv8x3BVbjuLfOsUO6X/uKlUr2QyDvhNpS\nrPVe6jFDgH3L2LcEPxRbPnAb9kthFBRdMu1Y7Itlq9aV7YT1qR6DLeK9HqKv89YBg7h302cY94dn\nuOeGiRx9w8c59HbosfWjAGzc7Wx6vbuZHlv6An9m/cDhrBt8J4OfuQRYotEWUpJjBVZmod2fcVyN\nnT/6d6z7dzv263xSzn7bsVE9ZwCfwXLA3tj/8xnAO8Cz2DmBfQGPY1EnMfTAGojbcWR2VF6pZF9u\naYFq+9y30THB617gn4HNWOGlzjjsp+FS4LGmKINgQxUvwurY7IS9ru2XwtvF7uvq9t47PN8er0yh\nD5vYm4l8cexidtq6js39Z/GnWw/ijUP/hV1X9wJ6sWVXWDVifLzXSKx19jxEI+vxMkiTs5O1z+C4\nCUvC72AjeQ7DFhWai/1qbwPOxj7Hj2KfvX2B32C/BnYrOC7AbCx39IL3Bg4cjTUWweWUD7mXXfhn\nrsORjJnlgY0G/jvn9qXA5ILHXE/H4hNgJ+YGlbkv2BeKLrrooosulV9qpifwEtCCtbLnA4cUPOZU\nOsZvj6ajPEE5+4qISEK0r3a0CGudg/WtTcp5zLXx9qeBo0rsKyIiIiIizcphJ1rnxZexQaMp31js\nHMWLFD8fkXSLgQXYa9744miVuRmbDJY7E3NPrArjP4D7sSGgSVUsfkd6PvdDsYVmnsWKn10Q35+W\n96Cz+B3Jfw92wYa1zweew0YxQXpe+zxXAN8OHUSFemDdUy3Ymf80npN4BfvApMEJwJHkJ8ufAd+N\nr08GftLooCpQLP40fe4HYxOkwEbKvIB93tPyHnQWf1regz7x357YedHjqfC1T1I9+yTU6alE7oSz\nLXRMGkubtLzus7Gx2LlyJ/T9FpjQ0IgqUyx+SM/rvwLeWyf4HWxy5BDS8x50Fj+k4z3YEP/tjTU0\n36LC1z5Jyf6b2AneqaTj50hnk8nSxGMLUMwFzgscS3cMoqPOz8r4dtqk7XMP9mv2SKxrIY3vQQsW\nf/vIwTS8BzthX1Yr6eiOSuxr/wDEkynyL6djszrbJ/9cScdCG0l2JnBjzu0vQuqKPLWXpB2IfZCS\nvkRgC/ndIIUt5dWNC6VbWsiPP42f+37Ak3S0ItP2HvTDGjft8aftPdgd+5I6ifS99jtogVSUQy13\n0lhaXAFcHDqIElrI/2wsxPpiwb64kl7KuoXOP9tdbUuKXlgtqoty7kvTe1As/lwtJP89AFu/+RIq\nfO2T0o2Tu+jBGaTjBZ+LrczTgvWjnYUVDUuLPnSUoO2L1SpJw+ue627gy/H1LwONqnRaK2n63EdY\nq/c5rGx5u7S8B53Fn4b3YG86upd2xVZ+m0d6Xvs8v8OGAD6NBZyYvqcS0jxp7ACs62Y+NhQt6fH/\nAVuzYDMSjhfZAAAAVUlEQVR2ruRsbCTRLNIx9Kww/nNI1+f+eKxg2Xzyhymm5T0oFv840vEeHI5V\nmJ2Pxfqd+P60vPYiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIo33/wGc2Z6GsRGFgQAAAABJRU5E\nrkJggg==\n", + "text/plain": [ + "" + ] + }, + "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+wfDDlvNnaTttIC4zMpzLLfG48hNUcokSjZ\np895SQcg4nCTvdvd8kDIeauAFVSqdOvxOvBE8KyL6t8CpAGU7NPlhxSPdpDmCltvYVcsqa0Fvlmy\nbwm2hoM7kbAdWZ+9xyCsYF9elL9zWL+9+4GhrpwmU7JPl38kHUAHibLewqvAGcBFZc73sTHnk6ie\n0LIsR6HPfh8KkzCfxuO1COeHJXu3K2i/miKUyJTspVNFWW9hJdb3vL7Ce7T7ZKCBWNmTDcTrr88L\nu0k729lWsm8yJfv06QK/f9JBdIB611vwgbuwhPXFBsaVJrXenM2bi31bGlRh/5NYFxnAeDxGx45Q\nIotSG0daZxNwMPB14MKEY2l3cdZbKOcg7AbkKGxt1fnAvWWO85ztHrK1VoM7xr6WZL8Gq5OzP+W6\nKD024DGHwo3ffSksfiR9TaWOtRqU7FMldzf4PwX2Bf8AyIWNeJDaxVlvoZwVwZ8rscW2JxOe7LPG\n+us9RlFYbnQtlsCj+gc26KDS/ajZFJL9fijZV9NDHWs1qBsnnU7AWovSPHHWWyjtmx8MDAu2h2CL\nbz/W+BATl+/GcUfhPBaUO4jqXopLIpdSv32LqGXfVP7Hqb3IkyaZNJe73kIXcDmFtRrAyniPwW4y\nbo51sZ2FjdzZGluwB+x36GrgjlYF3kL5ZD/JeW1OzPe4D7gS+3cq9yHhrmalZN9ESvbNNRm72JeG\nHSiJCFur4UWKu3ryVlPc2m1X5Vr2c2O+x6vY9b83xa34vGewf8+hwDZ4jMVjeQ2xSgh14zTfTMit\nTDoIkRrkb9DW07KHQr99Xx4bgUecV9S6bxIl+/SZQ+/iDv4vEo1EOt1gNmMt8J7g+SZs1nBc9wBH\nVNmvfvsWULJPndwfsRt+XwO2SDgY6WxDmMhmFG5QP41XtERhVHdgN2krjbd3k71q5DSJkn0q5TZi\n9cL/Bfz/AH9w0hFJRxrCdr3rzUL8/vq8Vdg31kqt+4ed7UkVjpE6Kdmn1yJsxMe5wGYJxyKdaQhj\nGOE8r6W/Pu864JQK+xZSmKk7Go8xdfwcqUDJPrVy90LuU1irSCQJgxnZuzIV1JfsrwGOgjIlEWwR\nlHnOK50w0qnllOxFpLwuhjK0qJVdazcOwOvAH4DvVdjvvreSfRNESfaq+S3SicYzln69c3GW4fFK\nne/4A+AYrJx0KSX7JgubVJWv+f1BrJbIQ9iU8qecY/I1v48tc36+5neU2tdS3nDgLPBvU60caalx\njHWe1dOFk/ca8GHgb1htoR5nn5vsdZO2CcJa9qr5nQ7nYy0ikdbZmq2cZ41I9mA1hE4Crgf2dF5/\nHNgYbO+Mx9AG/TwJhCV71fwW6VRbFY3Eqae/vtTfsHkkt5AvR+GxBusuBmsgvreBP08I78ZRze/0\n+Br4f4Gc7n1UN5U6an5LwCPHu72VPaFxLfu8a4CxWG2iQ7BRZ3OBPYL9exO+qLnEEJbsVfO7Zv5Y\nrO/xnga82UrgOawUr1TXQx01v6XXdgykK9h+Hbv+Gu0nWGXR+4AzWMfjzhWufvsGC+vGUc3v2u0C\n7EVjvv6+CbzbgPcRicpNtnPx6v6WX44P/DvwQ+AnXOsMy3ybw9H9voYKa9mr5nd9noZcaQndemhe\nhLRKvZUuo/Kx3HA1xzAaKysNA9mBLi5hI1+m/u5kIVo9e9X8Tgcf+3/QoibSfD6TnHZ1M5N9wcW8\nhMcyYFu6ybE1h7CCfwMuacnPb3NqKWbHtKQDkA7it6xlX6rQ7XkUV2Izbrds4c9vW0r2IlLMYyT9\nehcYfxd4uoU/vZDsJ7AVNkyzUgE1iUHJPpP8+8C/K+kopG25rfrH8CpOmGyG0rIJ/43m6DSEkn02\nrCl5fhBwBPg/SSIYaXvuvbZWduFAabLfnXuwSpkTWhxH21GyT79vYOORXcuwmkVdfQ/P8/cA/3jw\nRzYvNGlTxcMuW2sxNtQYYCtOYBtsNODRLY6j7SjZp17uZ5D7n+DJIPA9bIjrBmAC+KMqnHg8Vn9k\nl2ZHKG0nqZuz5WrbTwJuw+bpSB2U7LNjLVYiNl8z5DngYLRAszSSx2DyC4zb6PZaFhivV2m//T+B\nKWiSVV2U7DMjt5HiX7wZ2GS2KeDPAv9k67LxJ4E/PJkYpQ28l3xeWM0qvN7lAluptNzxc1iiLzef\nRyKKMqlK0uU4KPoFPD/4cwy2psBlWE0iJXypRaEL5w2WJxSD23W0N/YdYybWun8+kYjagFr22fS/\n2EIQ52DrC1wFfAv4brB/GwqLlE8NbtaebDdsRaoqjMRZxeKEYngSuycFsAMeWwAPo+JodVGyz6Z7\nIPcW5OZC7iasKuYYYPsyx/4Qq/z4O+B68B8G/4AWxirZUkioK1s6marA413gCeeVvbEiinuWP0Gi\nULLPlg3AW/RdFeznFJaF/D1wYMn+ccCgYHsfKKpTLmI8uoH39T5/NtEqtaVdOY+hBU3qomSfKbn/\nhdzmkPtDyetLsa+5pV7DarvrxlZ507AFdRYAZ5fZvyu2gMZa4Jsxz82iXch3/73NOpYl2j/uJvtJ\n2Pj7UVh1XamBkn17+xXkDsNWCZNiXdjEtGlYSe6TgN1KjnkVOAO4qIZzs6jQhfMSa7HVo5JSOiJn\nI1Zeffdkwsk+JfvO8BXg/UkHkTKTgYXAEqxb7FrsZrdrJbaAT2m3WZRzs6iQ7G0cTlqS/e54DMS+\nRe2UUDyZp2TfPF9OOoCC3LvAOxRGOIjdx1jqPF8WvNbsc9PMTfYDsOUIk+HxJrAoeNaNrU27ECX7\nmmmcffNMBz7bwp+3ERuD/Gr53bnHgf7g3wH8Fa1nW8/qR3HO9ZztHorXx00PjxzFyb4/NhggSXOA\nHYPtSViy/1By4SRuavCoiZJ9c93Yuh+VW0Hv0Et/J+AKyhex+ih2w7HTvUDxjevxWAu90ed6sSNL\nxvbkJ+L5vM7rbMJqMCVpLvDJYHsScA2p+sbccj0UNxYuiHOyunHaUm4h5D4LuXIfNutaHU1KzQZ2\nxkrnDgBOBG6qcGxpTZY452bFPr1b63iSJLtwCkpH5Kgbpw5K9h3L3wD+WUlHkaANwOlY+dwngeuw\n0R6nBQ+wiWpLga8D/wfrJhta5dwsKyT7t1hEsjdn89xkvxfH8wr24apSIDVQN07n6kJVBG8NHq5L\nne0XqTxHody5WVZI9i+yDCu5kSyPFXi8hC1eMoQ92Ik/sRDrx5+dbHDZo5Z9U/jTSfcH6fVJByAp\nYjdn9+19/hSvAq8kFk+x0vH26sqpUZoTUpadCtyCjZBJmZwPnAj+iqQjkdTYBtg62F7NU2zA5hik\nwRzgqGB7b5TsaxalZa8p5bW5HHKtXKhZpFb7ONtz2cRWpKdlr5u0DRKW7DWlXKT9ucn+EUhxsu+n\nZF+rsGSvKeUi7a802Y8iPd04i4DVwfbWHMI7wMQE48mssGSvKeWx+UcC+ycdhUgM6W3Zly5Afghj\ngJEUSnZLRGE3aDWlPL6dsIUXsjA0rAv8bsi1U82cqdQxpbzjeIyiMLx0LTZfID3J3swBDgKgm0nY\nfIftsfuBElFYsteU8tosgNwLSQcRYhNwIXAm5Ve4yqoe6phS3oHcpf4exWMD6erGgeKG0/7As8AO\nKNnHEtaNoynlbSv3DewXRjpbaRdOP2AEFQvqJeIhZ3t/bCET9dvHFNayd6eFdwGXU5hSDjbbcAz2\nn7E51lo8Cxt9s7rCuSKSHqXJfjj2u5umYcNPYzENBcYwlldYroZKXFEmVWlKuUj72tfZnoNNsErX\nhDuPjXg8AnwAgL0ZwHK17ONSuQSRTmU3Z/Mt5HXYot7pS/am0JWzA1uhLsjYlOxFOpc7RHgOHu9i\n3bIvJhRPNYVkvyUTsT77Ti/kF4uSvUjnctclnhn8mf6WfTd7Y/cHRyQWTQYp2Yt0riwl+8XAa8H2\ncEaxHI3IiUXJXqQTWVnjyc4rs4I/tyGN3TgePu54+x1ZjfrtY1GyF+lMO1FY8elVrAYNWJ99Glv2\n4HblTKQfatnHomQv0pmKW/Veb3mTccDyBOKJopDsxzACtexjUbJvKH86cELSUcS0HfgXJx2EtJzb\nX5/vwumHzZl5vvXhRPJg79YwxtOtZB+Hkn1j/Q+wH/BcwnFEtRSrjTMg6UCk5crdnB2FzVR9p/Xh\nRGDr0Vp3Uz/6M07rY8ShZN9YrwCTIPdfSQcSTW4T8G7SUUiLeQzElvjLy7fs09yqz/tn79Y4RmOl\nWCQCJXuRzrMPhW9zi/B6i55tR/qT/f29WxNYD2ybXCjZomQvANPBPyzpIKRlDnG273O2s5DsCy37\nbekmpxE5USnZN4y/FdA/6ShqcBNWjVSzETuHm+zvdba3o3h1uTR6AngDgMH0Z8eiUUVShZJ943jA\nMKwsdIbkXsTGWe8Gvvo/251HP+Bg5xU32e9IYbx9OtkyhQ/0Pt+eA5MLJluU7BvrAsgtTjqIGswH\nvo/W9ewEewJbBtsvAQucfTsBC1seUXyFrpxR7J5gHJmiZC9A7jzSOtxOGq24C6cwmaofNkkp3S17\nU7hJO4pxCcaRKUr2Ip2lUn/9OKzQ2NutDacms4CNAIxgMJ7uN0WhZC+unZMOoMWmYV1YC4CzKxzz\ny2D/PIoX514CPIqt7jSr72kpZMXPKiX7rHThgMdq4GHAKtq/xocTjScjlOwlbyFwDfi7JB1Ii3QB\nF2MJf3fgJOgzI/MYLAnuDHwJuMTZ5wNTsQ+ArIwImQiMDbbfxD6s8nYDnml5RLX7m7M9PbEoMkTJ\nXvL+BVuT+MSkA2mRydgH3BJsce1r6Zs0PgZcEWzPxG5sjnb2Z22lJLdVfz9e0BVi9sSWJcyKv/du\nDeKABOPIDCV7CeSexBJepxhH8ZjyZcFrUY/xgbuwGutfbFKMjfZBZ/sfJfv2BB5vYSz1up98qY9B\njMXr/cYiFSjZN4Q/DE1Kyho//BCgcuv9YKwL52jgqxS3mtPH+uvdZH+Xs50D3kuWWvYea3CHYMLh\nSYWSFd1JB9AmpmJf+W9OOA6J7gWs8FfeeKzlXu2YbYPXoFDzfSUwA+sWcm945nnOdk/wSMIe2MIk\nYKNu5jj7tsVayStbHVSd/g7ky3wcDvzfBGNphanBoyZRWvadNWIhNr8b+A5wN+SuTjoaiWw2duN1\nAlYU7ESsdITrJuDTwfYU4HVsItJgbLY0wBDgSCq3ij3n0VN/2DX7kLP9t5L++v1xFwbJjsJNWp8j\ngm8v7ayH4usplrCWfX7EwgexFs1DFGqp5LkjFt6PjViYEuzLj1h4jfaVw/6+hyYdiMSyATgduB27\nzi/HruvTgv2XArdg1/dCbPz554J9Y4Abgu1u4GrgjpZEXTs32d9Vsm8y2WyMzcbq7w8lx3ZYuYds\nDB9NQFiyd0csQGHEgpvsK41YeCl4rd0/bQE2Qe6+8MMkZW4NHq5LS56fXua8ZymuB59uHptR3Bi5\ns+SIyUBG1mBweKzH4x7oHWd/FEr2FYV143TiiIVOty/4o5IOQhrqMKzrCWABHm79pgHY6mpZbNmD\nffvK+0hiUWRAWMu+ESMWlmPLnd2J9f2n+SZWp3sE+DI2+SZrN+vyplLHTaw25SbBv5Ts2x9rDa9q\nXTgNdTPw62D7cDyGBjNspURYsk9ixEJG+AOxv/d3k46kcXIzwD8n6Sjq1ENxY+GCZMJICbtpWS3Z\nTwXublk8jebxHOeygAHsjH1LOQL434SjSqWwbpxWjVjIol2xEUifxkbjtJOjkw5AGua92KIkYIt+\n3F+y/0jc2ahZ1I8ZzjN15VQQluzdEQtPAtdRGLGQH7VwC3bDaiF2c+srwetjsFb8XOzG7V9J/4iF\nWs0JPyQzbgPOSjoIaZiPOdu34bHeeb4lth5tdlv2AN1FLfmP4WkR8nKiTKrqjBELkvdzyv9/Sjad\n4GyXdm8chTXIsr6WwUzWsIZBDAK2xmYz9yQbUvqoXELtTgs/RCRBHrth3TgAa7Fv164Tgf/X0pia\nwWMjL/GI88oJFY/tYEr2tTsOm1/wTTS2V9LpeGf7Fjzecp5vid3MvIF2sLzoQ+sT6srpS8m+Pt+B\n3E8h93zSgYgUsVE4bgv3+pIjPo711b/espiaaQ5XsLp3qPjWaEZ7H0r2Iu1pElb8DGANfbtwTgKu\naWlEzbSSVSwsmhvymcRiSSkl+9j8vcD/D4oXsWg3w8G/MukgpC6nOtt/xitaW3YnbBRO6QdAtj3N\nbc6zT+IsE5MRAAANJklEQVSxeWKxpJCSfXx7AOcmHUQTrQI+CwwCfwD4nVDbqL1YLZyTnVd+X3LE\nWcB/k/1ROMXmcy2v9X6oDUY3aoso2cfibw6cg80QPhVbx7PN5HxgHfBJrMb5mOrHSwp9HLsBCzYE\n+h5n33Dsg+DiVgfVdD73MZv+ziufTyyWFFKyj6cbm018DOT+ALn2ahmVd37SAUhsZzjb/4PHppJ9\nf6VQyqSdvMVjPIrPhuD5FDz2TzSiFFGyj2cgsB5yjyYdSJPdB3wCuBA4MOFYJA6P90PvAtzrgcuc\nveOBM2nnD/C3+AvPFQ2F/lpisaSMkn1k/rex1tAWSUfSfLmlkLuBdhqt0Tm+7mxfg8cK5/mPsQqR\nS1oaUWv9iTsZ6Tw/AY9tE4smRZTsI/GnYLMNL6JQF1wkXTx2xu615P3M2T4Jq1ufvUVK4nmKF3iR\nNcwLnncD30oyoLRQsg/l7wA8gFX/nAG5dxMOqNV2Af/MpIOQSC6A3pmjf8PrTXgfwNaJPp52G4FT\n3vX0FC2o9G9q3SvZR3ER9rV3NOT+mXAsrbYA6/Pt+F+U1PPYA/gX55X8Ogv7Y/VvTqK9qrNWcz0z\n2Q+f2cHzgbT3cOlIlOyL+DnwDwP/A+D/J/jXYq2ib0JuTdLRtV7uHWzJycPA/3HS0UgFVhrhIgor\nxt2MxwPYtXszNgSxdJHxdvYM8Az3F5VU/xIeeyYVUBoo2RfLYQs53AOcjY1XHln1jPbXA/yJ4ht/\nki7HAdOCbR9r1X8U+DPW2i9dnaoT/Iq7OAK/t9RxF3Bx8MHYkZTsK/saNuzwAmzBlg6VmwX8FOgH\n/ieSjkZKeAwHftH7fCO/xePTwG+whUs6qUXvmgEM5Wauh95x94cCn0supGQp2QPgTwD/FgoLF28C\nfMg9DLnvQ66Dkz1grcUb6eBflFSyVurvgHEAbOJVLuJgbOLf3tjAgk61ETib2XyDDUWzhX8ZjFrq\nOB2e7P1Pgr8GWIytu/pvwCOQ64LcL5ONLU1yG7H6KoeCf1LS0Uivr2CT38yf6c8afo1167yaVFAp\ncjMwi4vYBp/5wWtDgOvxGJJgXIno4GTv98e6aW6BotKo7bZ4eKP8A/vl+R74ZwRF0jr4+kmYx0ew\n4ZRmHq/xBAdgS4b6lU7rQF9gLTtwHffi966/uzfwx05b4KSDfln9A8G/GPxJ4H8dWIbddPwnsBtW\n2GwG8HKCQaZY7k2s/3cAlmTeBc4Af5tEw+pEHkfh8yfyv79vsII17AQ8mWhc6bQG+Bjzmcy9POi8\n/lHgKryiwmltLQ13pn2aGof/baxS5SPA4c6OXwA/g9xzzfvZ7cy/EZgePPlP4GLIvZBgQJU0+fpq\n6c/egun8F+/jS3QF77uRpXQxGY8XG/hz2tEw4CqO4UAmM8p5/VbgZDxWJRRXPWJdX22Y7P0tgR2w\nvrk3sdrs+WJIV2I1rh8DDgj6oqUmRckebIjfBmBWyj5As57sRwPTGcIJHMkH2KuoJfo88EE8FtT5\nMzpFjn6cwoe5hH2Lyp4sBk7F6x2mmRUNT/bTgJ9j41Qvo3xtjV9iNzjfwZLrnBjnlgTsTwG2gNzt\nlUPyDwBOwcbEH4H1u78B/IjyVRovDGKZEcS6IeiWkJr578FWO3offe9z/A1b83Qp8LKNakpMtV+I\nFl/boYZh/6b7BY996c9oDmUhU9iJ7qKVl54EpuEVlQWQKIYwmGOZwc4cWbLnWuAHeJnpDmtosu8C\nngY+iC3Y8RA27dodingMcHrw5/ux7pEpEc8tE7B/P5aw7wnOuQl4DzAXG2XwKWAecGzhnB5gav7J\nqiCWocHz/sCTEVqbU4M3qkeHvoc/iuJ7HTOAfaFnO+f/5VzI/ah5MVQOjvLXeQLXNmOxxskIrL+9\nH1aoazdgMrAd8CgDmMNk3mQyOzCMD5FjeMn7XgOchsfqiP8GpaZS/79rM0yllXGdz3HAlfQrGZnz\nDvPo4lIGcmNQNbS1cUUXK9l3h+yfDCykUBL1Wuyru3tRfwy4Itieia2QMwaYGOHcEv4e2NfWt7EJ\nEIdSvmLdROAG4MPAE3DxGzB1AfCVOrpmppK5JJua93gFu5Zy2AW4E7AbXHE+TJ0UHHMB+D8Mtpdg\nLeTvYuOhBwP5NQI+CWxmzwceAeuixhBXK6/tYcD3sG8GdwEr6MJnSwayBQMZx2vsyo1swzD6sQc2\nn2FQmfdZii0peCNeXSNuppLO5DWVVsb1fW7AYyb2Da1QLXQwe2GT0n7Dt1jF7bzNgdxMP2bjs5Ah\nLGMYLwJv1/n/0FJhyX4cFH1NXIa1cMKOGYe1YsLOdfijsfoe64D9bSKTvxc2+uM32Nfaz2J315+B\n3FzwtwA2wp+/BTkv5O8iTZPzsaSd97Q9Dj4RmITdIN8Sm6x2KPbNbQJWjsK1CNgRS6Lj4NgH4brv\nYF1FG7Fvd8uxb3pgLeoc1i13G/G05tr+FktYwzg24x2G8Cb9OBT7MIuzGPYS4CfA5Xh0YI2mJvJ4\nATgej8nAecBHcEcpDmU4IxjONpwGnFZ07nrgbDawifXAWnzWAmvJsZZ+rKcf6+liHf1YRz820jtZ\ns8+ftcYeS1iyjxpII26AXYt9sn+4MGM1ly/RWmFpsdwbDfi50jSzF2EVM9+AXJkuB38bYFfs29wv\nsV+yL0DucvC/C93HYMvojQ1OmI51z91N8cgqgkv1OMjNiBhca67toWwfdChuTrwEvwi4E7ga+GfJ\n0oLSaB6zgOl4jMG65I7FuuwGVDynP9CfbiyPlvsmlilTALfFdA5WIMz1W6wfPW8+9ssb5VywXzo9\n9Gjmoxxd23q0w6NhurEWxgTsE24udjPJdQw2Ggbsl+DBGOeKJEXXtkiJo7E+2IVYCwbo0391cbB/\nHjZ0rNq5Immha1tERESkEY4HnsBGWexTsu8cbEm8+dBn4kMlk4FZ2JC+h6h4UzfUGdgQusepfXHm\nb2J32kfUeP6PgxjmYUNMt4h43jTs32wB5fuQw4zHbn4+gf3961l7tgv7v6h14YwtseX0nsImEE2p\n4T3Owf4uj2Fj0wfWGEscjb6um8HDRhDNCR7Tqh7dfPVet820BBsWPAfLL0n5PfASdi3njcBu4j8D\n3IH9zqTSrsAuWHJxfyl2x/pA+2N9oguJVrCtBzgq2D46eN+4DsP+8fJT0kdVObaS8djNu8XUnuw/\nROHv/J/BI0wX9m81AYu/ln7kMVhFQLBRL0/X8B5538BGktxU4/lXYMXpwPrIo37g5U0AnqWQ4K8D\nPlNjLHE0+rpuhguw/580aMR120z1/B430iHYMGY32V8IfDvYPpuQPJFk1cv52CdSqenAH7FRrEuw\nC2FyhPdbQSEhbInNbIzry1jJhXwp1JVVjq3kpxT+A2p1J/QOtZtJtAW/3UlC6ylM9InjReyXDWA1\n1qoeW/nwirbFbm5eRm1DF7fALu7fB883YOUw4ngT+3cYjH1YDKa2ayKuRl/XzZKGuljQmOu22dLw\nb3Uv9CnW5k76u4KiqgJ9pbHE8VjsK2ZefiJLmO9gE0+ex7pBarlptjO2SPOD2DeF/WKePx2L99Gw\nA2M4lcKIkGoqTQCq1QSsJTGzhnN/Bvw71Dw2fCL2QfsHrFrp76CocFUUr1G4HpYDr5PsEn21XtfN\ncgbWTXg5yX79b/R122g+dt3MBr6YcCylRmNdOwR/jq52cNikqnrdiXUNlDqXeH25+fGkld7vPKx/\n+UysLsvxWKvwQzFiOg/79xiO9Q/vjxXz2iHG+edQ3BdbrUUQ5d/mPGxG8TVV3ievkWNuh2L95WdB\n7PorH8Hq5MzBKYwTUzfWBXI6dv/l59iH+fkx3mNHrNrpBOxbwZ+Ak7GupXo1+rpuhmrX6SXA94Pn\nP8A+FD/fxFiqaea/QSMchPUajML+Tedjrey0afi4+2Yo7dv8DsVVFG+japmFXm4Vyxzxv/aD1bY+\n1Hm+EBgZ8dw9sU/XxcEj/3V96xriACsNcT82tT6KqBN9wvQHbqdQFjquH2IttcXYL8nbWGnpOMYE\n5+cdDPw15nuciHUj5Z1CYY3hVmjUdd1sEyjuB261Rl23rXABNvgiKRMo/r+aT+EDfZvgeardDezr\nPM/fyBqAfZ1fRLQ+s0coJOojsBZhXKdhBavAbrI9X8N75NVzY2caNqJjqxjnNGKiTw5LzD+LeV4l\nh1L7aJx/YP8HYKNH4o6M2gsbUTQI+3tdAXy1xlhq0ajruhnc1cW+TrRvjs2S5glqg7EidmDrY9xP\nsqOoJtD3Bm3+g/E7RBvIkYiPYy3ANdiNwVudfedirer5FEbYhNkP61+eCzyA9TfH1R+4CvsHfZja\nuyHARoLUmuwXAM9RGBr3m4jn1TvR52Csn30ujRmWdyi1j8bZC/vAjjv81PVtCkMvr4CWLEHX6Ou6\nGa7E7ivNA24kpK+3BdI6QW0i9rswF2s4JBnbH7F7T+uw6+tzWH65iwwMvRQRERERERERERERERER\nERERERERERERERERabn/D5wtRkgPiur7AAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "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/vFpmvvJSzBZWdw\nBQu0FPYGPRthZTJIz04F77U0E7yteP2aLGEP0AV4lPTAv7SZOrsYzPQC/5bi1CZSEsFlZ3AFC7QU\n9g1wRzJA3warLXrQQwthD+4SzAdJD/wfZ90LHOOFfb3BjsWrUaRdFT07RwBTgI+Ai7MsH4j76rwG\nuCBj2SzgbWAi8Hoz+1fYByl72BsMqvdGixrRLkEPOcIeXOA/QyrsG4AjstSbsPQbvp4obp0i7aao\n2VkLTAP64NpZJwGDMtbpAXwFdzIsM+xnAt1zvIbCPkjZw34qvJUMzqfBEq5ppMhBD6mwtwPAzgE7\nMctKGwNvkQr81cCe6+0JdjHXQVoy8A8sfr0iRVfU7Nyb9COdS6IpmyvJHvab53gNhX2Q1g/738BP\n/UstD4AXaLf+45vCfizYZJrvHXNL3AFLMvDnRvPS9wZ/9Gp/zdrlD5RIURWUnbn6s+8FzPGez6XZ\nG2CaLeZp3C3wpxZSmISlM/TeD25KPn8Alj4Hh1HEAUla8HwLyz7FNUUujZ73wnXVkHld/c9J9ak/\nFDiqmAWKxC3XjSRtPereF5iPa+p5Ctf2/2KW9UZ7j8dHk4Sj0/Ewbvfo87QGmOH6v1lRwhr6gc0D\nToDE+Ixl04DvAY/hjtj3x3Wj3PRNNOH6yRkLnB/Nusbg4QTUt3vlIvkZTjuOxbAX6c04I8l+khay\nN+Pks1zNOEGx3mCngS1ONuPUwq/e9ppvPoB7SlCH34xzLlgvsOfBWurJ8jLSr9A5Pm2PsIXBcq85\n50ft+AOItFVRs7MDMB13grYT2U/QJo0mPcy74rqxBeiG68Pk4CzbKeyDYt/EdVnwO7AewBHHeUG/\nFtaa+ybX3nX4YR8NTG5P5gj7Gly/+cmwXw70TdsrXO6F/WxLDZ8oUm6Knp2HAB/ivgqPjOadHk0A\nPXHt+stw7aKzgQ2BL+H+OEwC3vW2bfeCpT3ZN3G9TwJs1wE++8AL+0a4pkR1NBf2h7hlzdqE9BO2\nL+E1ZxpsaLDQC/wz2+9nEGmT4LIzuIKrW1PY1wDjf+gFfYNrAsl1qW2x6sgW9k9EpTyaY+OhuLb4\nZOCnDUJucK4X9nMsj37yRWIQXHYGV3B1awr7czuATfPC3tx5mVLVkSXsAeywPMIeYBSpsK/Hu/4+\n6kZhgfdznVHk4kWKoaiXXopk8dqGwC9PBPqlZn6GG4AkFNfhmnDA3QtwD675kYS7+eoGb92ROrqX\n0CnspUBLa+HEwbWwQUbvYjcm3HmbMmM1YDOjyb/LtgH4Ae4kLbi/W/75hjuAhdHjbYGT2rtSkUqn\nZpygDLsbsO+nN98ssdSVVyViNdHLv9VyM07Ten8Bqwd7N2NHPyTVnNMIDGvaEi70fsaPzV2RJlIu\ngsvO4AquYjtBoq4G7P30sL8896bFZgmwkdG0uzc/W9g3gnUE25n1u1VIkN4l8lRcN8kYdPOvzJnE\nLqPAjnWTSOyCy87gCq5StUQDg3w7PeiXGWwad3EpdhjYElw/bFuR1jumDQCbF12e6d8LsA2uCSoZ\n+E3t9QYXJX/WpWyypCNrp9Fyb5sipRJcdgZXcJU6iygMJ6WHfYmuq8+XbQH2NdzA5F/KCPttwB7H\n3f2befPVj0mFfQPu8szkdfeLkj/vNYx8RmEvZSK47Ayu4Cq0La6fG/tmetCvKs3dsq2R7JUzW7/3\n/p229iDY61B/Pq7Tvijwu34Ii74CtrHBJcmfeSFbLOvIWoW9lIPgsjO4gqtMAniYKATfpGatF/Y3\n5dg2RjYD7FSwO5sJ+9U03Qtm/wC7Gdd1whc0Bf6Na8Au7MGnf1xBtzXJn/sk/k9hL+UguOwMruAq\ncwxR+H0t/ah+jcHWcRfXPJsBNgbsWbCfZCwbDHYg2EG4DtQujNb9FfzoBZrCvks9vHsv2LtXMHpV\n8mefQn/bnEUDwIaDRTdj2Xa4DuJOA9u4xD+sVKfgsjO4gqtIN1xfRwbY+zDPC/uxMdeWg82Ijurz\nuNHLzgVb5H6s1f+FrZfQFPiHrwa7ezOWzFxFl6ZvNZdw7RPRyd6p0T4OA5sDtgLsMrA33STSboLL\nzuAKriLXEIXegbA0GXR11DYa9I67uJbZ22DLwa7PY90twXaPptvg7n+SOllrwGFg0//K915LvgdT\n2X5VLXVjwD4BOxLsZLBHwaZH3xLudauKtJvgPl/BFVwldsCN3GS4w+R3kkH3b46YHXNt7chuc2F9\nvH8p5kxYMaM3M3+3ms7rku/DGM66HuwhsIVgt2aE/RiFvbSz4D5fwRVcBRLA40Rhd4AX9I3QsBtv\nPhtzfe0oGfbzpwFec865S8HG3M8xE7ymrDcMErjeNusV9lJiBX2+1DeOZHMkbtxWALvPhR4A0+n3\n4kSGfBFPWaXUsxH4Wer5bZvCs90v4+qJuJEXAb6CG5DnKGBj3MlsaLdB1kXCpqOf8tIVmEV0RPt1\nuM87krW9efk5UoOXVCC7DexP0YnXGlzPmNHR/aC50HCLwRjvPXnJHd03bf8hWJ07Mawje2lXwX2+\ngiu4wv2CVFv14lXwQDLYnuGAT8F+C3ZEzDW2I7sK7GOw8dGML5M20MlVjxtsY7DOC/zhWfaToGmk\nRvvfUlUvVSW47Ayu4ArWD9dEYYAdBaMMGpOhdgDPvAKWbRzhSncTTWHfZSmwscGdXtg/vf4mlgDr\nDHY72FklrleqQ3DZGVzBFewRUkf1rzXAXV6gPQ42rkrDfmPgE1Lvza8NvmRQ770/e2ff1G5T2Es7\nCS47gyu4Qh1OKswaT4dvGtR5YTasisMe4Luk3p96YGeDP3vvTzPnMZJX91h9NGV2wCbSWsFlZ3AF\nV6AuwAxSYXanwVgvyJ53q1V12CeAZ0m9Ry+shYF+M5fBbutvZrdFl2Z+FL1/I9ZfR6RVdOmlFOwi\nXCdgAJ/dDrcApyQXfsDAW8DOB/rEUFu5MFw3z/XR8/06u0sv7/fWGZVlu7dwfeU/ihsJS6Rq6cg+\nXn1xA2wnj1hPN7jRO1p9vRNrBoItBfs12IB4y43dDaTeqwWvwT7+pakGOzW/qT2hI3spoqJn5whg\nCvARcHGW5QNxIxitAS4ocFtQ2MftIVLhNWEG9DBY6YXXkWADwabEXGe52AiYS+o9u9ngIe/9+mvz\nmyrspaiKmp21wDTc1/eOwCRgUMY6PXBfZ68mPezz2bboBUtBDiUVWgbsafALL7jeMahR2K/nO6Te\ns4bb4QTvPWsw169QFgp7KaqittkPxQX2LKAOuBd3K71vETAhWl7othKfDYAx3vM/GnwAnO3Nuzah\nduZs7geeiR7XnAFnGIxLPgcuiacskeblCvtewBzv+dxoXj7asq20vwtxN1EBfA6MBP6X1ODh04B/\ngJ0C/Kb05ZW15Mna5AHOvj+D173lP7Sy7wJaqk2HHMvb0sRSyLajvcfjo0naT2/gUu/5KINVwPne\nvOsSbuDtIcB84NYS1heCKcCvic5F3QSnXQsvd4J9cb9XFwFnZtmuFiz5e9cACTVjSr6Gk7VrjuLY\nC3jCez6S5k+0Xkl6m32+2+rDXnoPkGpzfguoNfip1+4826CTW9XGgmULLYENcd9eDbCT4UHvPcwy\nbKM95t1cZWD9suxTJF9Fzc4OwHTcSdZONH+SFdzRuR/2+W6rsC+tb5B+UnZvg84Gc7yg8m7vV9jn\n8G2i9zIBDcvhbe99bGFAdpuusJc2Knp2HgJ8iGvDHRnNOz2aAHrijm6WAUtxY5Zu2MK27V6wNKsz\nMJVU0N+Fe3CmF1CfmrujNqKwzyGBOzlrgH0P3vfey1XmrlbLQmEvbRZcdgZXcMBGkgr6z4H/MdjA\n0gcSPy99E4V9HgYA64je24WuGSz5fl6TfROFvbRZcNkZXMGB6g18QSrsz8Y98Nvq56cf1YPCPm/X\nEr23x8My7z1dbrDZ+qsr7KXNgsvO4AoO1L9JBf0koINBlyjgk8H00/U3U9jnqRuuCdMSYHPhM+99\nvXz91RX20mbBZWdwBQfI777YiPpeNzjXC6R55m60yqCwL8DRRO/x99N7w1xirpsFj8Je2iy47Ayu\n4MB0BWaSCvrf4x50NViQ/Qocn8K+AAnc5cZWCzbbXX6ZfH8zLjtW2EubBZedwRUcmKtJBf1iYHPc\nkwu8IJqT/ageFPYF2wFYC9gpqffXDBabG/EqorCXNgsuO4MrOCBpV4kQ9VFv0M1goRdEZzS/C4V9\nK1wNWEewWa5jtOT7fEVqFYW9tFlw2RlcwYFI4AbCTgb9f4n6QjK4yAug2eauv2+Gwr4VuuI6ALQf\npB/dL7fom5XCXooguOwMruBAHE8q6BuAXXFPNjZY5AXQaetvaoPcOUY7D+yfCvtWORKwGrD30gP/\nBrfYpoOtBrsx1iolZMFlZ3AFB6A78CmpsG/qtdLS+6ufZU194PjsPLCpYCvBXgc7tlSFV5AEbihC\nOzo97L8w2BqsK9hFYLfEXagEK7jsDK7gANxFKujnEZ0YNOhp7hb+ZPD8IPvmdh6YujVuu364Edzs\njfTAH+sW2zkKe2mDgrJTA45XnoOBk7znZwDLo8dX4NqTAd4G/la6sqrSdOB6gMvS559msH0M9YjE\nSkf2xbMh6dfU35dcYNDfoN47umxheDwd2RdRF6J/k/HpR/f3R0f288H+G3ONEqbgsjO4gsvYb0gF\n/WfAlskFBv/wguZZc23KWdi2YFcr7IvqMMCGpoe9PcMBR4Dt656KFCy4z01wBZepPXHjxSbD/sTk\nAoOhlh40ezS/G5sJNhcsS38u0gZ3A/a39H+HV67kykT0dBjY1jn3IpISXHYGV3AZ6gS8SyronyQ6\ncjc37t0rXsD8o+Vd2Uywvu1cbzXqDszvDbbGC/wGEseCvRQ15zTTZYVIVsFlZ3AFl6ErSQX9Stzo\nYOBmnOAF/VpLDTLeDIV9OzoCsOvTj+6nu5va7DaFvRQouOwMruAysztQRyrsz0kuMNjQ0gcmuS73\n7hT27eyvm4AtTg/8nynspRWCy87gCi4jXYD3SQX9y0BtcqHBVV6gzLf1utnNRmHfzjYHFpydHvYr\nd+Kdu8BuAPty3AVKMILLzuAKLiP+1Tcr8ZpoDPpYehe7P8pvlwr7EjiqA9g7XuBPYMgkXM8Kn4Ad\nDvY/cRcpZS+47Ayu4DJxIKmgN1IDwIOb8S8v6CdY3jfQKexL5G9fzbgU8wWGHQv2MNgisBbugxAB\nAszO4AouA5sRDYEXTY/iXTdvcHjGpZb75r9rhX2JbA7Mu8f7d2qE9ww6gj2hsJc8BJedwRUcswTw\nAKmgXwJslVwYnZSd7QX9/+W3W9sIbEeweQr7kjmoFzQuT//DfIHCXvIUXHYGV3DMziS9+eZof6HB\nTV5wLLKm/tNzsW/ierl8H2ybYhctzbrhQi/s62H1TrzzvMJe8lD07BwBTAE+Yr1xNJuMiZZPBnbz\n5s/Cdbg1EXi9mW0V9vkbTDTkXTSN9RcaDLH0kZF+mN9u7XKwl8EeKXbBklOnzvCWf7L2VYYsqaVO\nYS+5FDU7a4FpuJt0OgKTgEEZ6xwKPBY93hN41Vs2E3fnYEsU9vnZEPiQVNBPwhs31qCDwRte0D9j\nzfZ/k8nG4frDGdYOdUtuA/aC1fVe4D/IEb+Ouygpe0XNzr2BJ7znl0ST7w7gOO/5FFIdcM0kdzOC\nwj63BHAv6ZdZ9vdXMBjlBf0acwNf58nGgR1cxHqlcCfc6IX9KmrWGqivHGlJQdmZ63K8XsAc7/nc\naF6+6xhuHNQJwKmFFCZpLiT9D+oZwNTkE3NNZ6O95b9IuGY1Ccc918Dvk/9oXWns9B5bPJn/tzOR\nlnXIsTzfvxzNfSCHAZ8APYCncEf9L2ZZb7T3eHw0ifN10rs5uB3XgyIA5ppy7ib1b/kqTeOcSkg+\nh7PPhf0ehYEAO7H4y9Nc9xc3x1yalIfh0dQu9iK9GWck65+kvQM3uHWS34zjuxK4IMt8NeM0ry/u\n0spk881LZIwZa3Cj13yzKnfzjSXAto+mncGeAftMzThlo/ed3p3Pq6DhXdg57qKkLBU1Ozvghlbr\ngwuZXCdo9yJ1grYrqb5YuuH6bckWKAr77DbBXcmUDPp5QE9/BYPhBo1e2J+Re7dWG62+HGwM2Byw\ng8DyvERT2ltfOPx9r/1+KomlH2UdGF6qXNGz8xDcVSDTcEf24G7N92/Pvy1aPhkYEs37Eu6PwyRc\nX+sjyU4nqxpuAAAM5UlEQVRhv76OwDhSQb8W94e0icFW5jo3S4bCE/m171otWAPYb6Own9IO9Usb\nHUHfP/j93j/sDpZEfMFlZ3AFt7ME7q5X/8aptOvlo8ssx3tBvzD/KzcU9oFI3AIv+d1eXBMNXi4S\nCS47gyu4nV1OetBfkbmCwS+9EGg0OCi/Xdv+YDd7YT9JYV++doSOj7Np001yn4F9H06Ouy4pG8Fl\nZ3AFt6OfkB70d5HRNGPrd3J2Wf67t/PAXoj+fxDY+WB5dn0scdiKdz7+BOqS/97vQeNg+FrcdUlZ\nCC47gyu4nZxIetA/hWu7b2Kwi8FyL+gfs7y7LoYo5H9TxJql3dn069j/ZL/9/j+wrpuu0JEAszO4\ngtvBcUADqaB/DdjYX8Fga4M5XtB/bHl3cta0F4V9cGw6WL8J8LOM9vuVwE5xVyexCi47gyu4yL4N\n1JMK+om4/uqbmOu2+C3vl3255X1kZ4NcNtg1YFcp7ENj08HWgd00De7xA/8MWAHsGneFEpvgsjO4\ngovoZNKP6N/D3W3cxKCjwSPeL3mdubtq82SHu5OwthTsA4V9aKwz2IVgY07jji4LSTR1dlcP9h0X\n+ENy7kYqUXDZGVzBRXI+6W30U/AGIYGmSyzvyzghe0r+L2HXRVfcPAx2rDsZa3sW8WeQkrBzostk\nbSOWNc6k+4Lk52EN2Ndc4H817iql5ILLzuAKbqMa4FrSg/4tIG2A6Sjo/5YR9NcU9lI2Dmwk2G65\n15XyZedEzTkGduvWzH1vLt1XJD8Xy8CGwTrgu3FXKiUVXHYGV3AbdAP+SXrQv4jrGqGJQa3B3RlB\nf4vl3z99DVgXXL836vMmeLY92NHRNBjsrO2ZOmOtu5nODDfE2IHu83Qx6imzWgSXncEV3Erb4k6+\n+kH/KK4PoSYGGxjcnxH0YwsI+mvB1kabrnbf8qWyWH+wqQY71bmhJ5v+sQ9zn6s/k/G5kooUXHYG\nV3ArHAQsID3obyGji2mDzS3jFnmDO3NfS2/dwHqBjY3a6EeBbdXyNhIuF/YABv3rYW7y87IO7Afu\n8zUZ2D7eOqWdBZedwRVcgI649vlGUiFfB5yWuaJBX4MpWZpucgX9l8Fuig7sDOwitdFXOusffXt7\nA2wfgz4NMN3/7FwPVgvLgGPirlbaTXDZGVzBeeoHvEL60fynwP6ZKxocarAkI+jPy+9lbCbYe2AX\nFq90KW+2AdgeYK+AjQAw6NUI7/qfoUfANnKfu7vJuHdDKkJw2RlcwTl0AC4CviA96Mexfn/0HQyu\nywj5NQbHRmtsC3YA2HCwTcGG4DowGwt2A9ghUdj3LeUPKOXCngB7IPpmd9NZjOlj8B//8zQVbChN\n4yEcEnfFUlTBZWdwBbdgD9Y/CVuHC/+05hiDHbK0z8812Mdb63ywebj7Z36L662yDuxB3GWVH4F9\nrrCvVnYM2AXRtARse4Mag2v9z1Ud2OWuWcdwV4P1jrtyKYrgsjO4grPoA/yV9JA3XPDv7q9o7o7Y\nS8wbei6aHjfYIn23dj7uZpoGsDvA/tdb1h/simjSV/SqZ9Oib3sHARgcb7DM/4y9Arar+1yuxg0T\nqit2whZcdgZXsGdL4EZgDekhvxp3NJ/Za+VXDSZlhHy9wUh3ItauAXsMbDTYsCjgk2HfAJbHsINS\nnexisGfBbmmaA30MXvA/b/Vgt7mjg+Q5pPOALrGVLW0RXHYGVzBuIPCxuFDPPJp/AHdytom5rokf\nzQh5e5PdGocwYQ3uSppJYE+B/QfXvcFMsJcV8JI/OwdsMdjkpjnuBr1LDNb5n71F0dfGru4zOx+4\nANg0ttKlNYLLzlAKrgG+AfyL9F4qk9OrwDB/A4Oh5vq28QcFtwYSqx9jxF86sK4Bd6frbl57/Biw\nCWAL1RYvhbEeuDtsLfqWOBLsLLBOs9lm5wYST2UecHyK62Wtm/sMr8QdxAyK+QeR/ISSnU3KveD+\nuKECZ7B+wBswATia6A7Xhzmsy0Vc9/gU+n+S+YvVQMLu5Tsrt2NW8hv1P91LWIfoF7UH2D64Kywe\nANuy9D+uhM0S0WfnYVxvp4a716o+QcNjbzH4x+voMDfzs/k52K1gO6UfvJxJwWMmSAmVe3aup9wK\nTgCDgUuBSWQPeAOeBb4OdcNrqH93GC+seovB4+qp+TzzF8nAHuHQxr15+VWwPtHUPZafTqqUHYo7\nHzStM6tX3MgFz33GpsuyfVZfBjsbbCv3OV+Hu2z4TGC7eH8GyVBu2ZlT3AUngB2Ak4A/4dovmwn4\n2mWwxe9h1LdG8NhXDL71PgPHLab7umy/NA0k6v7ED20wb9VHv2wds5cg0t7s0OgAfm3UVPhqR9ba\nhdwwfRGbL8z++XUDFl8Ctjuudz3gXeBW3LdZHfXHq6DszKdzrRHAzUAt8Afg+izrjMHdsPEFLjQn\nFrCttVyHbUlqcIbxkFidR83NWHMuTNgDHvwCbloAWx8EKwfD8hauRuhMghFsz/7sxhbsyxvsx4vs\nymRqmnmv59Jr3TbMu34BW/5uKxZ8Ec1eAYm61tcu0ha2GTAgetIV118TQJcEjWeO4IlF5/Pr2v15\nvmdH6rPuYQmui9bXomkCsAJmAW9600RgMfEfxFWDHNlZmFpgGu468o64Zo3MkzeHAo9Fj/fEtfXl\nu22y4BbYEdHJyjVgd4L9mNQdpHtH6xzq5tcdDGP7wrf3gf7fgpvmw5mLYdAbsPFU2MAfFWq9adPo\nCOY4NrRf0G/Fn9i3/g0Gr2kgsSrbkU/GtMjg9wb7WUGDgBfd8BhfuxINj7uA9mc7JD/G/fjoD6dx\nhz3FQdZAojHX534O7u6+W8B+grvdewAs2xDewHXTcAXwPVw3ISegyzyLqahH9nvjbr4YET2/JPr/\ndd46dwDPAfdFz6fgfkH65rFtsuAaoDOwAe7D0MU9Hr0f9BsOtX3gnYnQcwRYX1gKfLwA1mwEi+pg\n9aY1zG3szAI6s65mA29HG+GuJ9vEmzbFdRSyJW5oqK2hsSeJRBeskL+Sjbg/YI8DjwBvJNwQg3Eb\nHU1SHKOp+PfTEqQOUJqOFlfSbfNufPGNKQw4uTcf79qFNQXdvLcYmAt8AnwG3IM7GlwBX6ygdvky\nOq1aRe3nxsp5a+GztbBkFSxczK61i9inxxp6d4Jpb8Ima6DrMrjqbWAV7pLnOsrj9y1OBR3Zd8ix\nvBcwx3s+F/fvlWudXsDWeWwLwPPQWIv7tNV6Uw2j/ed7ZSzv2YHUX4jObTuarsnjj+QCXLi/ifs2\n+0oClrfhNUXKRMLIEpwbwkLgbvjwbnOhMhDY618c/fPtmdZrR95PdKS+2bDZIpoGR8+n0fRXsys0\ndHWZDWTcZe56Z55MHe5uxTrctc4N0f8zpwYSJO86BGs018BqgBm1Cdx8q6GhwUjQSCJhLiOthoaG\nBNYY/fZb9B9r2ihjPnkeTec7+ER+q2U3LPcqaXKFfb7FtKndqIwGz1yNu8RymjdNB95JuLAXqUoJ\nlwUfuOlfDwOdd2dCwwT22ATXPDsIGGTQtwH61MBWNe64rNU6knELerMsKpHoG0qjt6wx27pVKVfY\nz8ONsJS0Le4IvaV1tonW6ZjHtkBZjaHWBdgpmkJ2ZdwFVBi9n2ncb+ybFP67+/Oi1yLF0gF3ZNsH\n6ETuE7R7kTpBm8+2IiJSJg4BPsQ1aYyM5p0eTUm3Rcsnk7pMsrltRURERESkkhwLvIc7yT4kY9lI\n4CPcZZwHl7iuSjAad35kYjSNaHFtyWYE7vP3EXBxzLVUglnA27jP4+vxlhKk/8N1Sf2ON6878BQw\nFdelRdn2WjoQ18nYc6SH/Y649v2OuPb+acR7k1KIrgTOj7uIgOV7Q6DkbyYunKR19gN2Iz3sb8CN\nmwHugCTzHqY0cYboFNxfpExHAn/HXV47C/dLN7R0ZVWMMrrIKThDcZ+7WbjP4b24z6W0jT6Trfci\n7m5S3xHAn6PHfwaOamkH5XjEvDXpl2gmb9KSwpyNO2H+R8r4612Zau5GQWk9A57Gdalzasy1VIot\ncU07RP9vsUv0XNfZt9VTQM8s8y8FHi5gP9V7J0TzmntvRwG3A7+Inl8F3AScUqK6KoE+b8W3L65H\n2R64z+4U3NGqFEfOu3vbO+y/3optst2kNa845VSUfN/bP1DYH1bJ72ZCKcz86P+LgAdxTWUK+7b5\nFHfAtwDXzdfCllYul2Ycvy3vP8DxuBux+uL6mtfZ+8Js5T3+FukndSS3CbjPXR/c5/A43OdSWqcr\nrk9CgG64K+z0mWy7/wAnRo9PBP4dYy0t+hauXXQ17i/T496yS3EnyKbgxn2VwvwFd5nbZNwHQMMb\nFk43BBZPX9wVTZNwg5/o/Szc33EdiK7D5eaPcFc3PU0Al16KiIiIiIiIiIiIiIiIiIiIiIiIiIiI\niIiIiJTc/wOhqrQ4DxElWQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "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 +}