{
 "metadata": {
  "name": "",
  "signature": "sha256:bb3ff2a11114777b7230e90bf632d48347a4fd67d53da451a5dd4b031a99528a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Practice Problem - Monte-Carlo Error Propagation"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Part 1"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You have likely encountered the concept of propagation of uncertainty before (see [the usual rules here](http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas)). The idea is that given measurements with uncertainties, we can find the uncertainty on the final result of an equation.\n",
      "\n",
      "For example, let us consider the following equation:\n",
      "\n",
      "$$F = \\frac{G~M_1~M_2}{r^2}$$\n",
      "\n",
      "which gives the gravitational force between two masses $M_1$ and $M_2$ separated by a distance $r$.\n",
      "\n",
      "Let us now imagine that we have two masses:\n",
      "\n",
      "$$M_1=40\\times10^4\\pm0.05\\times10^4\\rm{kg}$$\n",
      "\n",
      "and\n",
      "\n",
      "$$M_2=30\\times10^4\\pm0.1\\times10^4\\rm{kg}$$\n",
      "\n",
      "separated by a distance:\n",
      "\n",
      "$$r=3.2\\pm0.01~\\rm{m}$$\n",
      "\n",
      "where the uncertaintes are the standard deviations of Gaussian distributions which could be e.g. measurement errors.\n",
      "\n",
      "We also know:\n",
      "\n",
      "$$G = 6.67384\\times10^{-11}~\\rm{m}^3~\\rm{kg}^{-1}~\\rm{s}^{-2}$$\n",
      "\n",
      "(exact value, no uncertainty)\n",
      "\n",
      "Use the [standard error propagation rules](http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas) to determine the resulting force and uncertainty in your script (you can just derive the equation by hand and implement it in a single line in your code).\n",
      "\n",
      "Now, we can try using a **Monte-Carlo** technique instead. The idea behind Monte-Carlo techniques is to generate many possible solutions using random numbers and using these to look at the overall results. In the above case, you can propagate uncertainties with a Monte-Carlo method by doing the following:\n",
      "\n",
      "* randomly sample values of $M_1$, $M_2$, and $r$, 1000000 times, using the means and standard deviations given above\n",
      "\n",
      "* compute the gravitational force for each set of values\n",
      "\n",
      "You should do this with Numpy arrays, and **without any loops**. You should then get an array of 1000000 different values for the forces.\n",
      "\n",
      "Make a plot of the normalized histogram of these values of the force, and then overplot a Gaussian function with the mean and standard deviation derived with the standard error propagation rules. Make sure that you pick the range of x values in the plot wisely, so that the two distributions can be seen. Make sure there are also a sensible number of bins in the histogram so that you can compare the shape of the histogram and the Gaussian function. The two distributions should agree pretty well."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from matplotlib import pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "def gaussian(x, amp, mu, sigma):\n",
      "    return amp/np.sqrt(2*np.pi*sigma**2)*np.exp(\n",
      "        -(x-mu)**2/2/sigma**2)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "G_CONST = 6.67e-11\n",
      "\n",
      "def compute_force(m1, m2, distance):\n",
      "    return G_CONST*m1*m2/distance**2\n",
      "\n",
      "def make_plot(m1, sig_m1, m2, sig_m2, r, sig_r, gauss_scale):\n",
      "    force_est = compute_force(m1, m2, r)\n",
      "    force_error_est = np.sqrt(\n",
      "        (sig_m1/m1)**2\n",
      "        +(sig_m2/m2)**2\n",
      "        +(sig_r/r)**2)\n",
      "    print(\"Est. error: {}\".format(force_err_est))\n",
      "\n",
      "    force_dist = compute_force(\n",
      "        np.random.normal(m1, sig_m1, 1000000),\n",
      "        np.random.normal(m2, sig_m2, 1000000),\n",
      "        np.random.normal(r, sig_r, 1000000))\n",
      "\n",
      "    plt.hist(force_dist, bins=50) and None\n",
      "    eval_at = np.linspace(min(force_dist), \n",
      "                          max(force_dist), 200)\n",
      "    plt.plot(eval_at, \n",
      "         gaussian(eval_at, gauss_scale, force_est, force_error_est),\n",
      "         color=\"red\",\n",
      "         linewidth=3)\n",
      "\n",
      "make_plot(4e5, 0.05e4,\n",
      "  3e5, 0.1e4,\n",
      "  3.2, 0.01,\n",
      "  1000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Est. error: 0.0021650635094610966\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucFNWZ//HPF4arcglGQUCR6KCgoohK4i2jRkWTVTdx\nFU1cskuiG2Kie/sF9ibu/pLI7i9rzCaaaDSiq0QSjZpEETRMTDRcVFBgRAYNynAZDSqoCHJ5fn/U\nmZ7qYZgZmO6uqu7n/XrNa06dPlX91DDM0+ecqlMyM5xzzrk96ZJ0AM4559LNE4Vzzrk2eaJwzjnX\nJk8Uzjnn2uSJwjnnXJs8UTjnnGtTu4lC0rWSlkpaJunaUDdA0lxJKyXNkdQ/1n6qpHpJKySdG6sf\nG45TL+nmWH0PSfeH+vmShhX6JJ1zzu27NhOFpGOALwEnAccBn5F0ODAFmGtmI4AnwzaSRgGXAaOA\n8cAtkhQOdyswycyqgWpJ40P9JGBjqL8JmF7A83POOddJ7fUojgIWmNlWM9sJ/Bb4HHAhMCO0mQFc\nHMoXATPNbLuZrQZWAeMkHQz0MbOFod3dsX3ix3oAOLtzp+Scc66Q2ksUy4DTw1BTb+ACYCgw0Mwa\nQ5tGYGAoDwYaYvs3AENaqV8b6gnf1wCY2Q5gk6QB+3Y6zjnnCq2qrRfNbIWk6cAc4H1gCbCzRRuT\n5OuAOOdcmWozUQCY2Z3AnQCSvknUM2iUNMjMNoRhpTdC87XAIbHdh4b2a0O5ZX3TPocC6yRVAf3M\n7K2WcXgycs65fWNmar/VnnXkqqeDwvdDgc8C9wGPABNDk4nAQ6H8CDBBUndJw4FqYKGZbQA2SxoX\nJrevBB6O7dN0rEuIJsdbZWaZ/br++usTj6FS489y7B5/8l9Zj78Q2u1RAD+XdACwHZhsZpsk3QjM\nkjQJWA1cGv6Q10maBdQBO0L7pkgnA3cBvYBHzWx2qL8DuEdSPbARmFCQM3POOVcQHRl6OqOVureA\nT+2h/beAb7VS/xxwbCv12wiJxjnnXPr4ndklUlNTk3QInZLl+LMcO3j8Sct6/IWgQo1hFZsky0qs\nzjmXFpKwYk9mO+ecq2yeKJxzzrXJE4Vzzrk2eaJwzjnXJk8Uzjnn2uSJwjnnXJs6cme2c+Xp1Vfh\nl7+ELVtg9Gg4/3zo4p+dnGvJE4WrPDt3wr/8C3znO7B9e3P92LFwzz0wcmRysTmXQn7DnassZvA3\nfwO33db664MGwdNPw8c+Vtq4nCsSv+HOub313e/mJ4kTT4QrroCePaPtDRuiIagPPkgmPudSyBOF\nqxyvvQb//M/N25//PCxYAPfeC48/3pwsVq6Eb34zmRidSyEfenKV4+KL4eHwGJTRo2HRIujevfn1\n226Dq6+Oyt26wdKlcOSRpY/TuQLyoSfnOuq555qTBMCtt+YnCYAvfQlOOSUqb98O39pttXznKpIn\nClcZpk9vLk+Y0JwQ4rp0gf/+7+bt++6Lhqucq3AdeRTqVEnLJS2VdJ+kHpIGSJoraaWkOZL6t2hf\nL2mFpHNj9WPDMeol3Ryr7yHp/lA/X9Kwwp+mq2irVsEDDzRvf+Mbe247bhyceWZU3rEjP3E4V6Ha\nTBSSDgO+DJxgZscCXYkeVToFmGtmI4iecT0ltB8FXAaMAsYDt4RnZAPcCkwys2qgWtL4UD8J2Bjq\nbwJiH/2cK4DbboNdu6LyeefB8ce33X7KlObyjBl+BZSreO31KDYTPSu7t6QqoDewDrgQmBHazAAu\nDuWLgJlmtt3MVgOrgHGSDgb6mNnC0O7u2D7xYz0AnN2pM3IubseO6Ca6Jtdc0/4+55wDRxwRlTdt\ngoceKk5szmVEm4kiPBv7O8DrRAniHTObCww0s8bQrBEYGMqDgYbYIRqAIa3Urw31hO9rwvvtADZJ\nGrCvJ+Rcnscfj+6NgOhmuvHj224PIMEXv9i8feedRQnNuaxocwkPSYcD1wGHAZuAn0n6QryNmZmk\nkly3Om3atFy5pqbGn2Xr2jdjRnP5yiuhqoOr1kycCP/6r9Gd3E8+CQ0NMHRocWJ0roBqa2upra0t\n6DHbvI9C0mXAOWb2pbB9JfBx4CzgTDPbEIaV5pnZUZKmAJjZjaH9bOB64LXQZmSovxw4w8y+EtpM\nM7P5YXhrvZkd2Eosfh+F2ztbtsCBB0bfAZYtg6OPzmvSPIWWz8yiIagnnogqvvc9+NrXihmtc0VR\nivsoVgAfl9QrTEp/CqgDfglMDG0mAk2DuI8AEyR1lzQcqAYWmtkGYLOkceE4VwIPx/ZpOtYlRJPj\nznXe44/nkkQdoGOOQVLeV8RafAWXXNJcjl815VyFaW+O4gWiiedngRdD9W3AjcA5klYS9S5uDO3r\ngFlE/y8fAybHugGTgR8D9cAqM5sd6u8ADpBUTzTMFbvkxLlOePDB5iKwe0Jop4d60UXRfAXA734H\nb7xRjCidSz1fwsOVpw8/hIMOiq5aAk4AFreaGMTuCUPkftfOOCNKEgA/+hFcdVWRAnauOHwJD+f2\n5Omnc0mCYcNYvK/H+exnm8uPPtrZqJzLJE8UrjzNnt1cvuCCfT/O+ec3l598MuqpOFdhPFG48hRP\nFPE/9h2Um/A+6ihebap8772op+JchfFE4crPunXwYrj2olu35rWb9krzhPdjTG6ufuyxAgToXLZ4\nonDl5/HHm8unnw7779+pwz1GrEcSP7ZzFcIThSs/8btSzz13j806fDhq2N60iMGLL8Kbb3b6mM5l\niScKV36eeqq5XIBlXt5nfxZycnPFb3/b6WM6lyWeKFx5ef11WL06KvfuDSecUJDD/oazYhu/Kcgx\nncsKTxSuvMR6E3O3bEHdu+9xPae9MY/YhPi8eZ0+nnNZ4onClZdYoniKf6dDS3V0wB/4BNuaNlas\ngPXrO31M57LCE4UrL3mJ4oyCHXYrvZgfr/D7KVwF8UThykdjI7z8MgDb6J4/AV0AeanhmWcKemzn\n0swThSsfsd7EAsaxlV4FPXxeavjDHwp6bOfSzBOFKx9FGnZqkjf09NxzsHVrwd/DuTTyROHKR5ET\nxUaAESOije3b4fnnC/4ezqWRJwpXHt5+G5YuBWAH0VVKRfGJ2HF9nsJViHYThaQjJS2OfW2S9HVJ\nAyTNlbRS0hxJ/WP7TJVUL2mFpHNj9WMlLQ2v3Ryr7yHp/lA/X9Kwwp+qK2sLFkB42NAS4D36FOVt\nrpoxI1d+8B//sSD3aDiXdu0mCjN72czGmNkYYCywBfgF0SNL55rZCKLnXE8BkDQKuAwYBYwHblHz\n/6ZbgUlmVg1USxof6icBG0P9TcD0Qp2gqxALF+aK89to1lnPsDRX/gSDivhOzqXH3g49fYroeddr\ngAuBpo9XM4CLQ/kiYKaZbTez1cAqYJykg4E+Ztb0P/ru2D7xYz0AnL23J+IqXCxRLGyjWWfVMYpN\n9AXgYDZwWBHfy7m02NtEMQGYGcoDzawxlBuBgaE8GGiI7dMADGmlfm2oJ3xfA2BmO4BNkgbsZWyu\nUpmVLFEYXZjPx3PbRZoJcS5VqjraUFJ34M+Ab7R8zcxMUufXSWjHtGnTcuWamhpqCrAyqCsDr7/e\nvPR3376s3Ly5qG/3DKdwHnMAOKWo7+Tc3qutraU2vtR+AXQ4UQDnA8+ZWdNi/I2SBpnZhjCs9Eao\nXwscEttvKFFPYm0ot6xv2udQYJ2kKqCfmb3VMoB4onAuJ9ab4KSTsCefLOrbxa+o8kTh0qblh+gb\nbrih08fcm6Gny2kedgJ4BJgYyhOBh2L1EyR1lzQcqAYWmtkGYLOkcWFy+0rg4VaOdQnR5Lhz7ZLE\nf116aW77W0VOEhDd9d3kWPAb71zZ61CikLQf0UT2g7HqG4FzJK0EzgrbmFkdMAuoAx4DJptZ07DU\nZODHQD3RpPjsUH8HcICkeuA6whVUznXEybGb6xbyi6K/32b68TLRjXfdoPn53M6VKTX/DU83SZaV\nWF3pVElsojf7sQWAwaxlPUPYfWlxtVK3p/r2297LFVzR1MH+wQ9g8uR9PAPniksSZtapG378zmyX\naSMhlyQaGMJ6BpfkfZ9jbGzjuZK8p3NJ8UThMi2+kHihlxVvy7Oc2LzhicKVOU8ULtPiqWERJ5Xs\nfRczpnlj2TL44IOSvbdzpeaJwmVaUj2Kd+mbm9Bm506f0HZlzROFy65t2zgmtpk3b1ACPk/hKoUn\nCpddy5dHl6cCr/AxNtG/zeaFljdP8eyzJX1v50rJE4XLrsWLm4vxOYMS8R6FqxSeKFx2xZ4w9zwn\nlPzt85LT8uU+oe3KlicKl10J9yjepS8rmjZ8QtuVMU8ULpt27oQXXshtJpEoAPIGnHyewpUpTxQu\nm+rrYUt0R/Z6BtGY0NPm4onizmuuQZI/HtWVHU8ULpsSHnZqfu9mx3ECra8R5Vy2eaJw2ZSSRPFC\nrHw0y+nKjsRica5YPFG4bIoliiSueGryNvB6eE5XT7ZxJC8nFotzxeKJwmWPWWp6FAAvcFyufDxL\nEozEueLwROGyZ80a2LgRgHeAPzI80XDiieK4vMEo58pDR59w11/SzyW9JKkuPM50gKS5klZKmiOp\nf6z9VEn1klZIOjdWP1bS0vDazbH6HpLuD/XzJQ0r7Gm6shLrTUSf35O9ymgJx+fKnihcOepoj+Jm\n4FEzGwmMBlYQPa50rpmNIHrG9RQASaOAy4BRwHjgFjVfL3grMMnMqoFqSeND/SRgY6i/CZje6TNz\n5Stv2Cl53qNw5a7dRCGpH3C6md0JYGY7zGwTcCEwIzSbAVwcyhcBM81su5mtBlYB4yQdDPQxs4Wh\n3d2xfeLHegA4u1Nn5cpO0/0JkvjFDTfk6tOQKF7hcN5jPwAG0cjAhONxrtA60qMYDrwp6SeSnpd0\nu6T9gIFm1hjaNELu/8dgoCG2fwMwpJX6taGe8H0NRIkI2CRpwL6ckCtnBhij+ViuJg2f340uLOXY\n3PZxbbR1LouqOtjmBOAaM1sk6buEYaYmZmaSin6n0bRp03Llmpoaampqiv2WLmX2510O51UAtlPF\nipTct7CE4/kE8wFPFC5ZtbW11NbWFvSYHUkUDUCDmS0K2z8HpgIbJA0ysw1hWOmN8PpaCBeWR4aG\nY6wN5Zb1TfscCqyTVAX0M7O3WgYSTxSuMh3L0lx5BUfxIcsSjKZZ/iWyziWn5YfoG2JDtfuq3aEn\nM9sArJEUnvvIp4DlwC+BiaFuIvBQKD8CTJDUXdJwoBpYGI6zOVwxJeBK4OHYPk3HuoRocty53Yym\neYXWFxmdYCT58ie0nSsvHelRAHwNuFdSd+AV4K+ArsAsSZOA1cClAGZWJ2kWUAfsACabWdOw1GTg\nLqAX0VVUs0P9HcA9kuqBjcCETp6XK1PxRBH9cb4vuWBilnIsuxBdMI4E2LoVevZMOiznCkLNf8PT\nTZJlJVZXeFEn1Pg9p3IqzwAwnsd4nPPZfSE+dbCusG1fZgQjqI+qn30Wxpb2Gd7OtUYSZtapm438\nzmyXGWJX3hxFmoaeIH/4iSW+lIcrH54oXGYM4zX68i4Ab/JR1nNwwhHly0sUL6Thwl3nCsMThcuM\n3Sey0/WAIE8Urlx5onCZEV8eI23DTpC/5hMvvBCtcutcGfBE4TIjrZfGNmlgKG/xkWhj0yZ47bVk\nA3KuQDxRuMzY/dLYtJEPP7my5InCZUJv4AhWAbCTLtQxKtmA9sAThStHnihcJhwDdMndr3Ak20jn\nzWzxeYoHr78+b9Xb5tX2ncsWTxQuE+IzEukcdorE505GczhNK962fgOfc9ngicJlQjxRpHEiu0kd\no3Lr2R7BK+zHe4nG41wheKJwmZCVRLGNnrwc247fSe5cVnmicOlnljfYlOZEAfkPU4pfqeVcVnmi\ncOm3Zg39Q/EtPkJD3mNN0ieeGjxRuHLgicKl34vpXrqjpXhqOC4VD2t1rnM8Ubj0eyHdS3e0tPvQ\nk1/x5LKtQ4lC0mpJL0paLGlhqBsgaa6klZLmSOofaz9VUr2kFZLOjdWPlbQ0vHZzrL6HpPtD/XxJ\nwwp5ki7jXkz7Hdn51gEbGQBAX95lGL6Uh8u2jvYoDKgxszFmdnKomwLMNbMRRI8unQIgaRRwGTAK\nGA/couY7jW4FJplZNVAtaXyonwRsDPU3AdM7eV6unOw29JR+8Th9+Mll3d4MPbUcGL4QmBHKM4CL\nQ/kiYKaZbTez1cAqYJykg4E+ZrYwtLs7tk/8WA8AZ+9FXK7MxO9k7iWxc8UKAHYhlnN0wtF1TLzn\n4xPaLuv2pkfxhKRnJX051A00s8ZQbgQGhvJgoCG2bwMwpJX6taGe8H0NgJntADZJGrA3J+LKTXQ3\n89EsomuoqaeaD+idZFAdln+HticKl21VHWx3qpmtl3QgMFfSiviLZmaSfMbOFVz6V4xtnQ89uXLS\noURhZuvD9zcl/QI4GWiUNMjMNoRhpTdC87XAIbHdhxL1JNaGcsv6pn0OBdZJqgL6mdlbLeOYNm1a\nrlxTU0NNTU1HwncZlvZnUOzJco5mJ13oyi6OYBW9eZ8tSQflKkJtbS21tbUFPaasnadwSeoNdDWz\ndyXtB8wBbgA+RTQBPV3SFKC/mU0Jk9n3ESWTIcATwBGh17EA+DqwEPg18D0zmy1pMnCsmX1F0gTg\nYjOb0CIOay9WVx6iax+if+snOYuzmAfAhTzML7mwZWt2v/y0o3XFahvV1TGSkUSd75NZwCLG4b/D\nrtQkYWaduvmoIz2KgcAvwoVLVcC9ZjZH0rPALEmTgNXApQBmVidpFlAH7AAmx/7CTwbuAnoBj5rZ\n7FB/B3CPpHpgI5CXJFylsrxhmywNPUEUb1OiOI4XWJRwPM7tq3Z7FGnhPYrK0dSjGEIDDWEUcxPQ\nn13sfvFdensUU/kW3+KfAfgfruHrfN97FK7kCtGj8DuzXWrlz09A2pfuaMmvfHLlwhOFS638Yafs\niQ+V+ZVPLss8UbjUynqiaGAob4d1b/uzKe9SQOeyxBOFS614osjmwI1aDD85l02eKFwq9eQDRrAS\niJbuWJZwPPsqf/jJuWzyROFS6WiW05VdQLR0R1ZvVvMehSsHnihcKmX5/om4lokivuBh05dzaeeJ\nwqVSVpfuaGkZx7ArXNY7AujJFpoWPPQHGrms8EThUqlcehQf0Jt6qgHoSjSk5lzWeKJwqVQuiQL8\nfgqXfZ4oXOocAnyEdwB4m/6syfgdCH6Htss6TxQudeIzEtEf2WxP+HqicFnnicKlTnygKevDTtDa\n0JNPYrts8UThUqfcEsXrHMom+gIwgLcZwtqEI3Ju73iicKkTTw1ZvjS2WculPHz4yWWLJwqXLlu2\nhItJYSddWMYxiYZTKH7lk8uyDiUKSV0lLZb0y7A9QNJcSSslzZHUP9Z2qqR6SSsknRurHytpaXjt\n5lh9D0n3h/r5koYV8gRdxixblvulXMkIttIr0XAKxXsULss62qO4lujRpk2zcFOAuWY2AngybBOe\nl30ZMAoYD9yi5jUKbgUmmVk1UC1pfKifRPTs7WrgJmB6507JZdoL8RVjy2HYKeKJwmVZu4lC0lDg\nAuDHNF+neCEwI5RnABeH8kXATDPbbmargVXAOEkHA33MbGFod3dsn/ixHgDO3uezcZnTct2j7191\nVe61cpjIbhIt5RE5kpfpwdZE43Fub3SkR3ET8I+Q+z0HGGhmjaHcCAwM5cFAQ6xdAzCklfq1oZ7w\nfQ2Ame0ANkkasBfn4DKvee2j0Zyeqy2nRPE++/NKKFexk1HUJRqPc3ujzUQh6TPAG2a2mD3c9WTR\n0+L9wnBXAJY3LFNOiQLyH77kE9ouS6raef0U4EJJFwA9gb6S7gEaJQ0ysw1hWOmN0H4t5K23MJSo\nJ7E2lFvWN+1zKLBOUhXQz8zeai2YadOm5co1NTXU1NS0e4IuOw5jNf3ZBMBbfIS1uU5neXgB+Fwo\nH8+SJENxZay2tpba2tqCHlNRh6ADDaVPAv9gZn8m6T+JJqCnS5oC9DezKWEy+z7gZKIhpSeAI8zM\nJC0Avg4sBH4NfM/MZkuaDBxrZl+RNAG42MwmtPL+1tFYXXZE1zpE/65/zoM8GP6UPslZfIon4y1p\nvePaWn1H64rVtvX9P434VSg/xel8kqcA4b/XrpgkYWadWgenvR5FS02/0TcCsyRNAlYDlwKYWZ2k\nWURXSO0AJsf+uk8G7gJ6AY+a2exQfwdwj6R6YCOwW5JwleEEns+Vn+eEBCMpjsWx8hgWI3b5mK3L\nhA73KJLmPYryFO9R/IpP82keBeAK7mUmV8RbkvUeBYhGDuQg3gSgmpWsYoT3KFxRFaJH4Xdmu9QY\nE/vMvZgxCUZSPPGe0pi8PoZz6eWJwqXCQDYwmPUAvE90V3Y5iifA+FCbc2nmicKlQvzT9QvALrom\nF0wRxROF9yhcVniicKkQ/6NZzp+zPVG4LPJE4VIhPgxTzn8+X+FwNtMHgIN4s8zuFHHlyhOFS4VK\n6VEYXVjC8bntMey+3lXTl3Np4YnCJa4f73A4rwLwId1YnnA8xZY//ATxta6av5xLD08ULnHx5SyW\nczTbE4ylFOKXyJbfbYWuHHmicInLH3Yq/z+du/conEs3TxQucfkT2eX/p/MlRrKVHgAMAwawMdmA\nnGuHJwqXuLE8lytXQqLYQTeWcmxuO37+zqWRJwqXqP2Bo1gBwA66VkSiAFjESbnySSxKMBLn2ueJ\nwiVqLNAlXOWznKP5gN7JBlQinihclniicIk6KVZelLdV3jxRuCzxROESVamJ4iVG8h77ATCEdRzM\nuoQjcm7PPFG4RFVqothF17xLgb1X4dKszUQhqaekBZKWSKqT9O1QP0DSXEkrJc2R1D+2z1RJ9ZJW\nSDo3Vj9W0tLw2s2x+h6S7g/18yUNK8aJumS1tkTFgRLDw+tb6ZF3JVAl8OEnlxVtJgoz2wqcaWbH\nA6OBMyWdBkwB5prZCODJsE14ZvZlwChgPHCLmhetuRWYZGbVQLWk8aF+EtHzt6uBm4DphTxBlyb5\ny1ScGJ5mB9FlsTvollRgifBE4bKi3aEnM9sSit2BrsDbwIXAjFA/A7g4lC8CZprZdjNbDawCxkk6\nGOhjZgtDu7tj+8SP9QBw9j6fjcuUk1mYK1fSsFOT+DmfyLP4Gk8urdpNFJK6SFoCNALzzGw5MNDM\nGkOTRmBgKA8GGmK7NwBDWqlfG+oJ39cAmNkOYJOkAft2Oi5L4p+iKzFRvMrHcvdkH8BbfCwsjOhc\n2lS118DMdgHHS+oHPC7pzBavm6SSfBSaNm1arlxTU0NNTU0p3tYVhVV8ogDxLHBe2DqJRbzK4UkG\n5MpAbW0ttbW1BT1mu4miiZltkvRronukGiUNMrMNYVjpjdBsLXBIbLehRD2JtaHcsr5pn0OBdZKq\ngH5m9lZrMcQThcu2Q1jDwPBrs4m+ZfuM7PYsIj9R3M+EJMNxZaDlh+gbbrih08ds76qnjzZd0SSp\nF3AO0QPIHgEmhmYTgYdC+RFggqTukoYD1cBCM9sAbJY0LkxuXwk8HNun6ViXEE2OuzIX7008x1is\nQq/Ujk9h+4S2S6v2ehQHAzMkdSFKKveY2ZOSFgOzJE0CVgOXAphZnaRZQB2wA5hsZk3DUpOBu4Be\nwKNmNjvU3wHcI6ke2Aj+kaoSjGNBrlyZw06RhbHyiTxLFdsr7uovl35q/juebpIsK7G63UUdyeZ/\nv99xGqfxNAAX8wsezl0EB5Dfds91e9O2s/sX771e5TCGsxqAE1nEc5wICP99d4UgCTPr1LN1K7O/\n7xLVnW3hctDIM5ySYDTJi5//qSF5OpcmnihcyZ3A8/RkGwD1wJsclGxACYsnilN4Jldu7W5255LQ\n4auenCuU+B/DZ9poVyn2lChaH9JyrvS8R+FKLv7H0AdaYCnH8i77A3AIDQyN7j91LjU8UbgSs7xx\neO9RwE6qWMC43PYp/lNxKeOJwpXUcP7IIKLVXzbRl7qE40kLn9B2aeaJwpVU/NPyH/iEL4MX7Hme\nwrnkeaJwJRX/tPw0pyYYSbrM5+O58vEsqZAnh7us8EThSir/iqfKvn8ibhP9WcbRAFSxs4LvVXdp\n5InClUwfNnMsSwHYSRcWcnLCEaVL/jyFc+nhicIV3J5uFDuVp+kSZiVeZDTv0SfJMFPn95yWK38y\nwTica8kThSsSa/EFZ/Gb3Ku11CQRVKrFfyanAd34MLFYnIvzROFK5kzm5cq/4awEI0mnNRzKK3wM\ngN7kr7DrXJI8UbiS6A+MYTEQzU88xRnJBpRS82h+gGQ8sTqXJE8UriTOALqyC4DnOYHN9Es2oJTy\nROHSyBOFK4n4g9bn5W25uPjP5hP8gR5sTTAa5yLtJgpJh0iaJ2m5pGWSvh7qB0iaK2mlpDlNj0wN\nr02VVC9phaRzY/VjJS0Nr90cq+8h6f5QP1/SsEKfqEtWPDX4/MSerWcwKzgSgJ5s4xP8IeGInOtY\nj2I78LdmdjTwceCrkkYCU4C5ZjaC6DnXUwAkjQIuA0YB44Fb1LyQ/q3AJDOrBqoljQ/1k4CNof4m\nYHpBzs6lwgH8ieNCeTtVeZeBut358JNLm3YThZltMLMlofwe8BIwBLgQmBGazYDcsywvAmaa2XYz\nWw2sAsZJOhjoY2ZNjwm+O7ZP/FgPAGd35qRcutRQmysv4iTeD0tqu9Z5onBps1dzFJIOA8YAC4CB\nZtYYXmoEBobyYKAhtlsDUWJpWb821BO+rwEwsx3AJkkD9iY2l15+Wezeid9PMY4F9GJLbtufeueS\n0OEn3Enan+jT/rVm9m78F9TMTFLRFwKdNm1arlxTU0NNTU2x39IVQPxGO5/Ibt+bHMQy4BigO9s5\ng6d4nKZRWn/qnWtbbW0ttbW1BT1mhxKFpG5ESeIeM3soVDdKGmRmG8Kw0huhfi1wSGz3oUQ9ibWh\n3LK+aZ9DgXWSqoB+ZvZWyzjiicJlw2H8kZGsAOADevpCgB00lyhRAJzPY7FE4VzbWn6IvuGGGzp9\nzI5c9STgDqDOzL4be+kRYGIoTwQeitVPkNRd0nCgGlhoZhuAzZLGhWNeCTzcyrEuIZocd2XgAh7N\nledxJlt7jcbwAAAPQElEQVTplWA02fForHxB3pZzpdeRHsWpwBeAFyUtDnVTgRuBWZImAauBSwHM\nrE7SLKAO2AFMNrOm/vJk4C6gF/Comc0O9XcA90iqBzYCEzp5Xi4lPs2vc+Vf8+kEI8mWp4D32I/9\neZ9qVlHNSuqTDspVLDX/DU83SZaVWCtd1GE0erGFjRxAr3DT2HBeZTXDW7am9XH3jtSVcv/Sv9cv\nuIiLQ6f7Om7iZv621f39/4VriyTMrFOTWX5ntiuaM5mXSxJ10EqScG15lAtyZR9+cknyROGKJv7H\n7ddttHOtiyeKT/Jb9kswFlfZPFG4fdbaNf3Nl01bi/kJt7fWMpQXGA1ADz70u1BdYjxRuE5q+YCi\naLx8FHUcxmsAbKIvTycVXsblDz85lwxPFK4o4r2JOZzLjgRjybL4lWKfBhSWaneulDxRuKL4HA/k\nyn5Z7L6bz8fZSLSazVD8qXcuGZ4oXMEdBowjWvtxO1U8woWJxpNlO6niQT6b276M+xOMxlUqTxSu\n4C6NledwLm/j6zt2xqzYT/Qv+JkPP7mS80ThCu6yWPn+vC23L+ZxJm9wIABDWMdp/D7hiFyl8UTh\nCuoI6jkhlLfRnYe5KNF4ysFOqniAz+W2L2VW3uttX6bsXOd5onAFFf8j9hjns5l+CUZTPuLDT5fw\nc7qwM/Zq65coO1conihcQcUnW33YqXCe4gzWh/IgGvkkv000HldZPFG4ghlJHaNZCkTPnvgVn0k4\novKxi678PLZ9OTMTi8VVHk8UrmAm5h57Ht078R59Eoym/Pw0Vp7AT9mP9xKLxVUWTxSuIKrYnpco\n7uKLyQVTpp4BXuIoAPrwnvcqXMl4onAd0t5VNRfwKINoBGAdMNsf3VkUt3FVrnwVtyUYiaskHXkU\n6p2SGiUtjdUNkDRX0kpJcyT1j702VVK9pBWSzo3Vj5W0NLx2c6y+h6T7Q/18ScMKeYKukPZ8Zc2X\nuT1Xvovokk5XeDOYyFZ6AHASzzIm4XhcZehIj+InsNvHwynAXDMbQfR86ykAkkYR3W81Kuxzi5o/\net4KTDKzaqBaUtMxJwEbQ/1NwPROnI9LwMd4JffsiV2IOxKOp5y9zQB+xl/ktq9OMBZXOdpNFGb2\nO+DtFtUXQm5AegZwcShfBMw0s+1mthpYBYyTdDDQx8wWhnZ3x/aJH+sB8GX3s+Yr3EqX0MN4jPN5\nNeF4yl18+OkKYH/eTS4YVxH2dY5ioJk1hnIjMDCUBwMNsXYNwJBW6teGesL3NQBmtgPYJMkXB8qI\n/XiPSbE+xPe5JsFoKsPvOY06RgLQB/hL7k42IFf2Oj2QbGYmqSS3gk6bNi1XrqmpoaamphRv69rw\nZW7nI7wDQD1H8DjnJRxRJRC3MJnv8zUA/o7/5kdcvdu8UGvLeJj5Xdvlrra2ltra2oIeUx35xZF0\nGPBLMzs2bK8AasxsQxhWmmdmR0maAmBmN4Z2s4HrgddCm5Gh/nLgDDP7SmgzzczmS6oC1pvZga3E\nYP5Lnpzoj07+z78b4hWGckjoLF7ND7mNq4Hd27Zet6f6NO6frvfqzfu8zqEcwFsAXMr9/Cxv3d7W\n9/f/Q5VHEmbWqcW/9nXo6RFgYihPBB6K1U+Q1F3ScKAaWGhmG4DNksaFye0rgYdbOdYlRJPjLgO+\nALkk0chB3M1fJhtQBdnCfnnDfP/MN335cVc0Hbk8dibRvT5HSloj6a+AG4FzJK0EzgrbmFkdMAuo\nAx4DJse6AZOBHwP1wCozmx3q7wAOkFQPXEe4gsqlWzc+5N9i29/lOrbSK7F4KtH3uYb3Q/k4XuSz\nPJhoPK58dWjoKQ186ClZLYeeruaH/JCvAPAnDmA4f4wt2ZHG4ZzyGnpq8m2U+2S1nFGM5kV20XWP\n+/v/ocqT5NCTq2D78y7TmJbbns43fF2nhPw/4F32B+Bo6vgrfpJsQK4seaJwu2lvuY4p3JhbrqOB\nIfyAryYRpgM2Av/FP+a2/y//4vdVuILzROH2oPXlOo6gnr/nO7ntqXybD+hd8uhcs+/w9zSE25IG\n0ZjX23OuEHyOwu2mtUthozHvXTzJ2ZzFPAAWAh9nJ7bb5400jvuX5xxFU90V3Mu9fAGAnXThRHax\npNW2rfP/W+XL5yhcSX2JH+eSxA66cjW0kiRcEu7jCp7kLAC6sos7ia5M250/NtXtPf9f7jqkmugS\n2Cbf5TqWJBeO2434G36YW1l2DPAf/GuyIbmy4UNPbjcth556sYWn2S+3pHUdIxnLc2ylN9kZzinv\noacmX+N7fI9rgWgl38/wKx7jgnbfy/9vlS8fenIlYNzGVbkk8SHd+Dz3+s11KfV9rmF2WG+rC8a9\nfJ4jqE84Kpd1nigqWGuXwba8FPYGrucL3Jvb/hr/wxJ/XE5qGV24knt4PWx/hHd4jPM5kDcSjctl\nmyeKirfnyc3ruIl/4z9y27fzpbxnIbh0+hMH8jngA3oCcASv8Djn0dba/e19YHCVzROFa9UUvs1N\n/F1u+zFgMrfQ1iWWLj2eBS5nJjvDf/ExLGEecEiur9GSXwnl9swThcvTkw+4Dfg2/5Sr+x2n8RfA\nDrolFpfbew9zMZO4g10huY8GFnESn+CZZANzmeOJwuUcy4ss4iS+HKv7DWdyPo/lVil12TKDL/JF\n7mJ7eKjRQN6glhqu4X/ows6Eo3NZ4YmiQrQ1Bt2b9/kX/oNFnMQxLM/V/y+fD0li/yRCdgVyD3/J\np3iCN8N2d7bzP3ydpzmVY3kx0dhcNvh9FBWitWU5eiK+yC38G//OwWzI1W8BruNH3M6XaZ6TSNf9\nAqXfP/vvNQzxMKM5LpYcdtCVmexkOktZzjEt9t+d/x/MnrK6j0LSeEkrJNVL+kbS8WRV+5e8Gsez\nmJv5OuuAW5mclySeZwwnALdzFT5xXV5eA8axgGlczza6A1DFTq4ElnEsjzGeL3APfdgc9vAJbheY\nWeJfQFdgFXAY0A1YAoxs0caybN68eSV5H8DAYl+77AhW2hVgP+LLtoYh1qKBGdgahthf82PryvZW\njtHacfdUV6y2ndl/XoZiba1+XsHf60hessc5p9XfhQ/oYbPBpvAtO4XfW2/ei+3f+lcafveLJevx\nh3+fTv2NripWAtpLJxM9HnU1gKSfAhcBLyUZVCHV1tZSU1NT0GP2kBgAHAC5718ChvNPHM4rHM4r\nVFNPv9wnxNt3O8YqDudWvsItTC7ju61rkw6gk2qBmoIe8WWO4jzmcBIL+Qbj+HNEl9Br6Mk2zgPO\nC1e+7UK8wuG8BLzOV1nDIbzOoaznYN6hP5sYA2+/DX37Qteuu0dfhN/9Usp6/IWQlkQxBFgT224A\nxhXt3S68EHaFB9GbNX+Plwtd9/rrMHfu3u+7axds2xZ9bd3a/LVtG9v2eILfbvP03+IjPMoF/IR7\nmcdKXwG2gi3iZC4BhvEql3E/lzOT43khr00XjGpWUQ3AD1o/0IDodr53gW3Ah8DgQw+F7t3hnXfg\noYeicrdu0KULSPlfsHtdW1+ltHIlPPts8/bnPw+XX17aGBKWlkRhJX23X/2q+Q9xKa1Z036bIniT\nj7KQP7GI65nLOSxgHDupAu4lRdNULkGvcRj/yTf4T77BMFbzSYZzJhM5mYUcyct0ZVeHjtMnfAHR\nh6Mmf/pToUMurfrYelknn5xcHAlJxVVPkj4OTDOz8WF7KrDLzKbH2iQfqHPOZZB18qqntCSKKuBl\n4GxgHdHD0y43s7KZo3DOuaxKxdCTme2QdA3wONEVUHd4knDOuXRIRY/COedceiU+k9nejXaS/kHS\n4vC1VNIOSf3Da6slvRheW1j66Dsdf39JP5f0kqS6MFeTifglHRmrXyxpk6SvZyX+8NpUSctD/X2S\nemQs/mtD3TJJ15Y69hBDe/F/VNJsSUtCnF/s6L7F1snY75TUKGlpSYPOj2+f4pd0iKR54Xd/WYf+\n33b2RozOfNGBG+1atP8M8ERs+4/AgAzHPwP461CuAvplKf5YfRdgPXBIVuIP+7wK9Ajb9wMTMxT/\nMcBSoGc4zlzg8LTFD0wDvh3KHwU2ht/1vTr3NMUetk8nejT50lL+zAv0sx8EHB/q9yeaH27zZ590\njyJ3o52ZbQeabrTbkyuAmS3qklxnYp/jl9QPON3M7oRonsbMNhU74BYK8fMH+BTwipmV+vrfzsS/\nGdgO9A4XU/QG1hYz2FZ0Jv6RwAIz22pmO4HfAp8tarS760j864G+odwX2GhmOzq4bzF1JnbM7HfA\n26UKthX7HL+ZbTCzJQBm9h7Rjc2D23qzpBNFazfaDWmtoaTewHnAA7FqA56Q9KykL7e2X5F1Jv7h\nwJuSfiLpeUm3hzal1Nmff5MJwH0Fj659+xy/mb0FfAd4nehKu3fM7ImiRru7zvz8lwKnSxoQXvs0\nMLSIsbamI/HfDhwtaR3wAnDtXuxbTJ2JPQ0KEr+kw4h6RgvaerOkE8XezKT/GfB7M3snVneqmY0B\nzge+Kun0gkbXvs7EXwWcANxiZicA7wNTChxfezr780dS9/DazwoZWAftc/ySDgeuI+q6Dwb2l/T5\ngkfYtn2O38xWANOBOUQPIFwMHbwrrnA6Ev8/AUvMbDBwPPADSX3a2acUshw7FCB+SfsDPweuDT2L\nPUo6UawFDoltH0KUGVszgRbDHma2Pnx/E/gFUXeslDoTfwPQYGaLwvbPiRJHKXXq5x+cDzwX/g1K\nrTPxnwg8Y2ZNwwkPAqcUJco96+zv/51mdqKZfRJ4h2isuZQ6Ev8phA8RZvYK0bzikaFdR8+9GDoT\nexp0Kn5J3Yh6p/9rZg+1+25JTMTEJluqgFeIPtV1Zw8TWkA/oomYXrG63kCfUN4PeBo4Nyvxh/qn\ngBHWPPE0PUvxh9d+SokngQv0+3McsAzoRTTPNQP4albiD/UHhe+HEo0z901b/MB/A9eH8kCiP2YD\nOnruaYw99vphJDeZ3ZmfvYC7gZs6/H5JnGSLkzmf6JPQKmBqqLsauDrWZiJwX4v9hocfzpLwH35q\nluIP9ccBi4jGDx+kxFc9FSD+/YA/ERJ2Bn/+/wdYTjTePwPolrH4nwrxLwHOTOPPn+hqm1+G3/Gl\nwBVt7Zuh2GcSzW1tI5or+KusxA+cRjRMuYRoyHIxML6t9/Ib7pxzzrUp6TkK55xzKeeJwjnnXJs8\nUTjnnGuTJwrnnHNt8kThnHOuTZ4onHPOtckThXPOuTZ5onDOOdem/w8vcEYwukDC8AAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xef9ab7ac>"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Part 2"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now repeat the experiment above with the following values:\n",
      "\n",
      "$$M_1=40\\times10^4\\pm2\\times10^4\\rm{kg}$$\n",
      "$$M_2=30\\times10^4\\pm10\\times10^4\\rm{kg}$$\n",
      "$$r=3.2\\pm1.0~\\rm{m}$$\n",
      "\n",
      "and as above, produce a plot.\n",
      "\n",
      "In this case, which method do you think is more accurate? Why? What do you think are the advantages of using a Monte-Carlo technique?"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "make_plot(4e5, 2e4, 3e5, 1e5, 3.2, 1, 1e90)\n",
      "# your solution here\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Est. error: 0.0021650635094610966\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAETCAYAAADNpUayAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEkxJREFUeJzt3X+sZPVZx/HPZ1moIlRAmi2FbTcpC4JSISWUgg1DQs2y\nJlQDFrANlaA06JbGxEqoJGxMjC0xpgVaShtAagkobSXbuIQiMIBpQSnLAoVd2QSSBetq+FVhNYHy\n+McclnGYO8+5e2fuzHP3/Uome2bme2cezpDPPPf7PedcR4QAAEvLsmkXAAAYP8IdAJYgwh0AliDC\nHQCWIMIdAJYgwh0AlqBFC3fb19veYfuxFmP/2vam5rbV9ovN46f2Pb7J9v/YPmPy1QNALV6s49xt\nf0TSK5K+GRHHzOPn1kk6NiJ+f+DxAyVtk3RoRPzvWIsFgOIWrXOPiPslvdj/mO33277d9kO277N9\n5JAf/V1JNw95/HckbSTYAeDtpj3n/nVJn4mI4yV9TtJX+5+0/T5JqyTdPeRnz9Hw0AeAPd7yab2x\n7f0kfVjSrbbffHifgWHnSLo1BuaObB8i6Vcl3THpOgGgoqmFu3q/NbwUEceNGHO2pD8c8vjHJX03\nIn42kcoAoLiR0zK2V9q+x/aPbT9u++I5xl1p+ynbm22PCutdIuKnkp62fVbzGrb9gb7X/GVJB0bE\nA0N+/FwxJQMAc8rm3F+T9McR8SuSTpT0R7aP6h9ge62kwyNitaQLJV0z7IVs3yzpB5KOtL3d9vmS\nPiHpAtuPSHpcUv9hjWdrSIDbXqXeETL35v95ALBnmtehkLZvk3RVRNzV99jXJN0TEX/X3N8i6ZSI\n2DHuYgEA7bQ+WqbpmI+T9ODAU4dK2t53/1lJhy20MADA7msV7s2RLd+W9NmIeGXYkIH7/AUQAJii\n9GgZ23tL+o6kb0XEbUOGPCdpZd/9w5rHBl+HwAeA3RARgw10KjtaxpKuk/RERHxpjmEbJJ3XjD9R\nvcMbh863R0TZ2+WXXz71GvbE2ql/+jfqn+5td2Wd+8mSPinpUdubmsc+L+m9TVhfGxEbba+1vU3S\nq5LO3+1qAABjMTLcI+Kf1WJePiLWja0iAMCCTfMM1Rp++EPp1VfV2W8/aedOad99p13RvHU6nWmX\nsCDUP13UX9NiXvI3Fuu9xmr1amnbtt721q3SEUdMtx4AexTbinEvqELSsr5dVPHLCcAeiXDPuO8L\n8403plcHAMwD4Z6hcwdQEOGeoXMHUBDhnukPdzp3AEUQ7hmmZQAURLhnmJYBUBDhnqFzB1AQ4Z6h\ncwdQEOGeYUEVQEGEe6Z/WobOHUARhHuGzh1AQYR7hgVVAAUR7hkWVAEURLhn6NwBFES4Z+jcARRE\nuGdYUAVQEOGeYVoGQEGEe4ZpGQAFEe4ZOncABRHuGTp3AAUR7hk6dwAFEe4ZOncABRHuGQ6FBFAQ\n4Z5hWgZAQYR7hmkZAAUR7hk6dwAFEe4ZOncABRHuGTp3AAUR7hk6dwAFEe4ZDoUEUBDhnmFaBkBB\nhHuGaRkABRHuGTp3AAUR7hk6dwAFEe4ZFlQBFES4Z/qnZejcARRBuGfo3AEURLhnWFAFUBDhnmFB\nFUBBhHuGzh1AQYR7hs4dQEGEe4YFVQAFEe4ZpmUAFES4Z5iWAVBQGu62r7e9w/Zjczzfsf2y7U3N\n7bLxlzlFdO4AClreYswNkq6S9M0RY+6NiDPGU9KMoXMHUFDauUfE/ZJeTIY5eb4uOncABY1jzj0k\nnWR7s+2Nto8ew2vODjp3AAW1mZbJPCxpZUTstH26pNskHTFs4Pr163dtdzoddTqdMbz9hHEoJIBF\n1O121e12F/w6jhaBZXuVpO9FxDEtxj4t6YMR8cLA49HmvWbOxRdLV13V2/7yl3v3AWCR2FZEzHvq\ne8HTMrZX2L321vYJ6n1hvJD8WB1MywAoKJ2WsX2zpFMkHWx7u6TLJe0tSRFxraSzJF1k+3VJOyWd\nM7lyp4AFVQAFpeEeEecmz39F0lfGVtGsoXMHUBBnqGZYUAVQEOGe4c/sASiIcM/QuQMoiHDPsKAK\noCDCPcOCKoCCCPcMnTuAggj3DJ07gIII9wwLqgAKItwzTMsAKIhwzzAtA6Agwj1D5w6gIMI9Q+cO\noCDCPUPnDqAgwj1D5w6gIMI9w6GQAAoi3DNMywAoiHDPMC0DoCDCPUPnDqAgwj1D5w6gIMI9w4Iq\ngIII9wx/Zg9AQYR7hs4dQEGEe4YFVQAFEe4ZFlQBFES4Z+jcARREuGfo3AEURLhnWFAFUBDhnuFQ\nSAAFEe4ZOncABRHuGRZUARREuGdYUAVQEOGeoXMHUBDhnqFzB1AQ4Z5hQRVAQYR7hmkZAAUR7hmm\nZQAURLhn6NwBFES4Z+jcARREuGdYUAVQEOGe4doyAAoi3DN07gAKItwzLKgCKIhwz7CgCqAgwj1D\n5w6gIMI9Q+cOoCDCPcOCKoCCCPcMh0ICKCgNd9vX295h+7ERY660/ZTtzbaPG2+JU0bnDqCgNp37\nDZLWzPWk7bWSDo+I1ZIulHTNmGqbDSyoAigoDfeIuF/SiyOGnCHpxmbsg5IOsL1iPOXNABZUARQ0\njjn3QyVt77v/rKTDxvC6s4HOHUBBy8f0Oh64PzQF169fv2u70+mo0+mM6e0niM4dwCLqdrvqdrsL\nfh1Hi27U9ipJ34uIY4Y89zVJ3Yi4pbm/RdIpEbFjYFy0ea+Zc9dd0mmn9bZPPVW6++7p1gNgj2Jb\nETHYQKfGMS2zQdJ5TREnSnppMNhLY1oGQEHptIztmyWdIulg29slXS5pb0mKiGsjYqPttba3SXpV\n0vmTLHjRMS0DoKA03CPi3BZj1o2nnBlE5w6gIM5QzdC5AyiIcM/QuQMoiHDP0LkDKIhwz3BtGQAF\nEe4ZpmUAFES4Z5iWAVAQ4Z6hcwdQEOGeoXMHUBDhnmFBFUBBhHuGP7MHoCDCPUPnDqAgwj3DgiqA\nggj3DAuqAAoi3DN07gAKItwzdO4ACiLcMyyoAiiIcM8wLQOgIMI9w7QMgIII9wydO4CCCPcMnTuA\nggj3DJ07gIII9wydO4CCCPcMh0ICKIhwzzAtA6Agwj3DtAyAggj3DJ07gIII9wydO4CCCPcMC6oA\nCiLcM/yZPQAFEe4ZOncABRHuGRZUARREuGdYUAVQEOGeoXMHUBDhnqFzB1AQ4Z5hQRVAQYR7hkMh\nARREuGfo3AEURLhnWFAFUBDhnmFBFUBBhHuGzh1AQYR7hs4dQEGEe4YFVQAFEe4ZpmUAFES4Z/o7\nd4mAB1AC4d4GUzMAiiHc22BRFUAxhHsbdO4AiknD3fYa21tsP2X7kiHPd2y/bHtTc7tsMqVOEdeX\nAVDM8lFP2t5L0tWSTpP0nKR/tb0hIp4cGHpvRJwxoRqnj84dQDFZ536CpG0R8UxEvCbpFkkfGzLO\nQx5bOjgcEkAxWbgfKml73/1nm8f6haSTbG+2vdH20eMscCawoAqgmJHTMuoFd+ZhSSsjYqft0yXd\nJumIYQPXr1+/a7vT6ajT6bSrctro3AEskm63q263u+DXcYwIK9snSlofEWua+5dKeiMivjjiZ56W\n9MGIeGHg8Rj1XjNt//2lV17pbb/8svTOd063HgB7DNuKiHlPfWfTMg9JWm17le19JJ0tacPAG6+w\ne/MWtk9Q7wvjhbe/VGEsqAIoZuS0TES8bnudpDsk7SXpuoh40vanm+evlXSWpItsvy5pp6RzJlzz\n4uNQSADFjJyWGesbVZ6WOfBA6aWXetvPPy8ddNB06wGwx5jUtAwkFlQBlEO4t8GhkACKIdzboHMH\nUAzh3gadO4BiCPc2OBQSQDGEextMywAohnBvg2kZAMUQ7m3QuQMohnBvg84dQDGEexssqAIohnBv\ng2vLACiGcG+Dzh1AMYR7GyyoAiiGcG+DBVUAxRDubdC5AyiGcG+Dzh1AMYR7GyyoAiiGcG+DQyEB\nFEO4t0HnDqAYwr0NFlQBFEO4t8GCKoBiCPc26NwBFEO4t0HnDqAYwr0NFlQBFEO4t8GhkACKIdzb\noHMHUAzh3gYLqgCKIdzbYEEVQDGEextMywAohnBvgwVVAMUQ7m3QuQMohnBvgwVVAMUQ7m2woAqg\nGMK9DTp3AMUQ7m3QuQMohnBvgwVVAMUQ7m1wKCSAYgj3NujcARRDuLfBgiqAYgj3NlhQBVAM4d4G\nnTuAYgj3NujcARRDuLfBgiqAYgj3NjgUEkAxhHsbdO4AiiHc22BBFUAxabjbXmN7i+2nbF8yx5gr\nm+c32z5u/GVOGQuqAIoZGe6295J0taQ1ko6WdK7towbGrJV0eESslnShpGsmVOv0LFum7pvbBTv3\nbrc77RIWhPqni/pryjr3EyRti4hnIuI1SbdI+tjAmDMk3ShJEfGgpANsrxh7pdNkvxXuBTv36v9z\nU/90UX9Ny5PnD5W0ve/+s5I+1GLMYZJ2LLi6WdE/LXPFFdJNN02vlt2xdav0ox9Nu4rdR/3TtRTq\nf8c7pEsvnXYliyoL97ZzEB64X2/uYoQdzz//1p1Nm3q3ap56atoVLAz1T1f1+it/Oe0mx4g5ZNsn\nSlofEWua+5dKeiMivtg35muSuhFxS3N/i6RTImLHwGstqcAHgMUSEYMNdCrr3B+StNr2Kkn/Luls\nSecOjNkgaZ2kW5ovg5cGg313iwMA7J6R4R4Rr9teJ+kOSXtJui4inrT96eb5ayNio+21trdJelXS\n+ROvGgAw0shpGQBATRM7Q9X2QbbvtP1vtr9v+4A5xj1j+1Hbm2z/y6TqaaP6CVtZ/bY7tl9u9vUm\n25dNo85hbF9ve4ftx0aMmeV9P7L+Wd73kmR7pe17bP/Y9uO2L55j3Ex+Bm3qn9XPwPbP2X7Q9iO2\nn7D9l3OMm9++j4iJ3CRdIelPm+1LJH1hjnFPSzpoUnXMo969JG2TtErS3pIekXTUwJi1kjY22x+S\n9MC0655n/R1JG6Zd6xz1f0TScZIem+P5md33Leuf2X3f1PduScc22/tJ2lrs//829c/sZyBp3+bf\n5ZIekPTrC933k7y2zK6Tm5p/f2vE2FlYbK1+wlab+qXZ2NdvExH3S3pxxJBZ3vdt6pdmdN9LUkT8\nR0Q80my/IulJSe8ZGDazn0HL+qUZ/QwiYmezuY96jdoLA0Pmve8nGe4r4q2jZnZImquQkPRPth+y\n/QcTrCcz7GSsQ1uMOWzCdbXVpv6QdFLza91G20cvWnULN8v7vo0y+745Ou44SQ8OPFXiMxhR/8x+\nBraX2X5Evay8JyKeGBgy732fHQqZFXSner8ODfqz/jsRESOOcz85In5i+12S7rS9pemCFlv1E7ba\n1PGwpJURsdP26ZJuk3TEZMsaq1nd922U2Pe295P0bUmfbTrgtw0ZuD9Tn0FS/8x+BhHxhqRjbf+i\npDtsdyKiOzBsXvt+QZ17RHw0Io4ZctsgaYftd0uS7UMk/eccr/GT5t//kvQP6k0vTMNzklb23V+p\n3rfjqDGHNY/NgrT+iPjvN3/9i4jbJe1t+6DFK3FBZnnfpyrse9t7S/qOpG9FxG1Dhsz0Z5DVX+Ez\niIiXJf2jpOMHnpr3vp/ktMwGSZ9qtj+l3rfk/2N7X9v7N9u/IOk3JM15tMSE7Tphy/Y+6p2wtWFg\nzAZJ50m7zt4desLWlKT1215h9y6UY/sE9Q6FHZzbm1WzvO9Ts77vm9quk/RERHxpjmEz+xm0qX9W\nPwPbB7s5mtD2z0v6qKTBa5zMe98vaFom8QVJf2/7AknPSPp4U9h7JH0jIn5TvSmd7zb7e7mkmyLi\n+xOsaU5R/IStNvVLOkvSRbZfl7RT0jlTK3iA7ZslnSLpYNvbJV2u3lE/M7/vpbx+zfC+b5ws6ZOS\nHrX9ZrB8XtJ7pRKfQVq/ZvczOETSjbaXqddw/21E3LXQ7OEkJgBYgvgzewCwBBHuALAEEe4AsAQR\n7gCwBBHuADAhbnFBvL6x77N9V3MG7T22B88wnxfCHQAm5wZJa1qO/StJfxMRvybpzyUNvTpkW4Q7\nAEzIsAvK2X6/7dub62ndZ/vI5qmjJN3dbHc1/MJ/rRHuALC4vi7pMxFxvKTPSfpq8/hmSWc2278t\naX/bB+7um0zyDFUAQJ/mwmYflnRrc2a+1LvMryT9iaSrbf+epPvUu3bMz3b3vQh3AFg8y9S7Lszb\n/pJScxHFM6VdXwJnRsRPF/JGAIBF0IT107bPknoXPLP9gWb7l5rry0jSpepdCG23Ee4AMCHNBeV+\nIOlI29ttny/pE5IuaP44x+Pq/ZUlSTpV0hbbWyW9S9JfLOi9uXAYACw9dO4AsAQR7gCwBBHuALAE\nEe4AsAQR7gCwBBHuALAEEe4AsAQR7gCwBP0fAATYYCItYOMAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0xeface2cc>"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}