sampling_trunc_expon.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First lets make our true distribution, a truncated exponential from 0 .. b=1 with beta=1.0\n",
    "\n",
    "We Make a histogram of samples from it to ensure it looks right.\n",
    "\n",
    "Our goal is to recover the value of b=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAERxJREFUeJzt3X+MZWV9x/H3x12stUUx7pgYlnVpuiRuUaMdAWNTMWAD\n1EAaq2GttTToRltMU00jjQ0Y/EdKtLURxa2lKyaCaI3d1FVMLXQTZS1LpciPQFakMmi6KyJJaxQJ\n3/5xL3aczuw9M3Pm3rnPvF/JhHvueXLP93B3P/Pd5zz33FQVkqS2PG3SBUiS+me4S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhq0eVIH3rJlS23fvn1Sh5ekqXT77bd/v6pmRo2bWLhv\n376dQ4cOTerwkjSVkvxnl3FOy0hSg0aGe5JrkxxJctcxxpyZ5I4kdyf5135LlCQtV5fOfS9wzlI7\nk5wAfAQ4v6p+DXh9P6VJklZqZLhX1QHgB8cY8kbgc1X1neH4Iz3VJklaoT7m3E8BnpPkliS3J3lz\nD68pSVqFPlbLbAZ+HTgL+EXg1iQHq+r+hQOT7AZ2A2zbtq2HQ0uSFtNH5z4H3FRV/1NV3wcOAC9Z\nbGBV7amq2aqanZkZuUxTkrRCfYT7PwK/kWRzkmcCpwP39vC6kqQVGjktk+R64ExgS5I54HLgOICq\nuqaq7k3yJeBO4Eng41W15LJJSdLaGxnuVbWrw5irgKt6qaiD7Zd+4WePH3z/b4/rsJI0NfyEqiQ1\naGL3lunL/C4e7OQlCRoI94WcspEkp2UkqUmGuyQ1yHCXpAY1N+fehfPyklrXdLgb4pI2qqbDfb6F\nSyYlqWUbJty7sNOX1IoNH+529JJa5GoZSWqQ4S5JDTLcJalBhrskNWjDX1BdiitnJE0zO3dJalCX\nr9m7FngtcKSqTj3GuJcDtwIXVtVn+ytx8uziJU2bLtMye4EPA9ctNSDJJuBK4Mv9lLV+GfSSpkGX\n71A9kGT7iGHvAP4BeHkPNU0Ng17SerXqC6pJTgR+B3g1Gyzc5zPoJa0nfVxQ/Wvg3VX15KiBSXYn\nOZTk0NGjR3s4tCRpMX0shZwFbkgCsAU4L8kTVfX5hQOrag+wB2B2drZ6OLYkaRGrDveqOvmpx0n2\nAv+0WLBvVE7XSJqELkshrwfOBLYkmQMuB44DqKpr1rS6KeWdJiVNWpfVMru6vlhVXbSqaiRJvfD2\nA2PkFI2kcfH2A5LUIDv3CbGLl7SW7NwlqUF27uuMHb2kPhju64BLJyX1zWkZSWqQ4S5JDXJaZko4\nFy9pOQz3dcy5eEkr5bSMJDXIzn0KLdXRO10j6Sl27pLUIMNdkhrktExDXFEj6Sl27pLUIDv3Ri28\n6GonL20sdu6S1KCR4Z7k2iRHkty1xP7fS3Jnkm8m+VqSl/RfpiRpObp07nuBc46x/9vAq6rqRcD7\ngD091CVJWoUuX5B9IMn2Y+z/2rzNg8DW1ZelvrmSRtpY+r6gejHwxaV2JtkN7AbYtm1bz4dWVwa9\n1L7eLqgmeTWDcH/3UmOqak9VzVbV7MzMTF+HliQt0EvnnuTFwMeBc6vqkT5eU+PhfWqkNq26c0+y\nDfgc8PtVdf/qS5IkrdbIzj3J9cCZwJYkc8DlwHEAVXUNcBnwXOAjSQCeqKrZtSpYkjRal9Uyu0bs\nfwvwlt4q0rrgdI003fyEqiQ1yHvLaFns6KXpYOcuSQ0y3CWpQU7LqBdO10jri527JDXIzl1j4z1t\npPEx3LWmlpqukbS2DHdNnB291D/DXRNhRy+tLS+oSlKD7Ny1bjldI62cnbskNchwl6QGOS2jqeAU\njbQ8hrvWlS6raAx6aTSnZSSpQV2+Zu9a4LXAkao6dZH9AT4EnAf8CLioqv6970KlUezopf/TZVpm\nL/Bh4Lol9p8L7Bj+nA58dPhfac35YShpcSOnZarqAPCDYwy5ALiuBg4CJyR5fl8FSpKWr48LqicC\nD83bnhs+970eXltaNadrtBGN9YJqkt1JDiU5dPTo0XEeWpI2lD7C/WHgpHnbW4fP/T9VtaeqZqtq\ndmZmpodDS5IW00e47wPenIEzgMeqyikZSZqgLkshrwfOBLYkmQMuB44DqKprgP0MlkEeZrAU8g/X\nqlipq6VW0RxrdY3z8WrJyHCvql0j9hfwx71VJElaNW8/IC3CFTaadoa7NOQHotQSw10awS5e08gb\nh0lSg+zcpWWwi9e0MNylFTLotZ45LSNJDbJzl3pmR6/1wHCXeuAySq03hru0hpYKfTt6rTXDXZoA\np2601gx3acLs7rUWXC0jSQ2yc5emgNM4Wi47d0lqkJ27tE65vFKrYecuSQ0y3CWpQZ2mZZKcA3wI\n2AR8vKrev2D/NuATwAnDMZdW1f6ea5W0gBdatZQuX5C9CbgaeA0wB9yWZF9V3TNv2F8AN1bVR5Ps\nZPCl2dvXoF5pw1vul38b+htTl879NOBwVT0AkOQG4AJgfrgX8Kzh42cD3+2zSEkrZ3e/MXUJ9xOB\nh+ZtzwGnLxjzXuDLSd4B/BJw9mIvlGQ3sBtg27Zty61VUo8M/bb1tRRyF7C3qj6Q5BXAJ5OcWlVP\nzh9UVXuAPQCzs7PV07EldeTyyo2jy2qZh4GT5m1vHT4338XAjQBVdSvwDGBLHwVKkpavS+d+G7Aj\nyckMQv1C4I0LxnwHOAvYm+SFDML9aJ+FSlo7Czt6p2mm38hwr6onklwC3MRgmeO1VXV3kiuAQ1W1\nD3gX8LdJ/pTBxdWLqsppF2lKOR8//TrNuQ/XrO9f8Nxl8x7fA7yy39IkrQcG/XTy3jKSVsTQX98M\nd0mdudpmenhvGUlqkJ27pFXz1gfrj527JDXIzl3SmvGi6+QY7pLGwqAfL8Nd0tg5R7/2DHdJ64ah\n3x/DXdK6Z+gvn+EuaWo5j780l0JKUoMMd0lqkNMykprgFM3PM9wlNafLDc5a/wVguEvakFr/BeCc\nuyQ1qFO4JzknyX1JDie5dIkxb0hyT5K7k3yq3zIlScsxclomySbgauA1wBxwW5J9w6/We2rMDuDP\ngVdW1aNJnrdWBUvSuBxr6ma9T9l06dxPAw5X1QNV9ThwA3DBgjFvBa6uqkcBqupIv2VKkpajywXV\nE4GH5m3PAacvGHMKQJKvApuA91bVl3qpUJLWufW4DLOv1TKbgR3AmcBW4ECSF1XVD+cPSrIb2A2w\nbdu2ng4tSeO31JTNegn6LtMyDwMnzdveOnxuvjlgX1X9tKq+DdzPIOx/TlXtqarZqpqdmZlZac2S\npBG6dO63ATuSnMwg1C8E3rhgzOeBXcDfJ9nCYJrmgT4LlaRpM8m7WY7s3KvqCeAS4CbgXuDGqro7\nyRVJzh8Ouwl4JMk9wM3An1XVI2tVtCTp2DrNuVfVfmD/gucum/e4gHcOfyRJE+YnVCWpQYa7JDXI\ncJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3\nSWqQ4S5JDTLcJalBncI9yTlJ7ktyOMmlxxj3uiSVZLa/EiVJyzUy3JNsAq4GzgV2AruS7Fxk3PHA\nnwBf77tISdLydOncTwMOV9UDVfU4cANwwSLj3gdcCfy4x/okSSvQJdxPBB6atz03fO5nkrwMOKmq\nvtBjbZKkFVr1BdUkTwM+CLyrw9jdSQ4lOXT06NHVHlqStIQu4f4wcNK87a3D555yPHAqcEuSB4Ez\ngH2LXVStqj1VNVtVszMzMyuvWpJ0TF3C/TZgR5KTkzwduBDY99TOqnqsqrZU1faq2g4cBM6vqkNr\nUrEkaaSR4V5VTwCXADcB9wI3VtXdSa5Icv5aFyhJWr7NXQZV1X5g/4LnLlti7JmrL0uStBp+QlWS\nGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB\nhrskNchwl6QGGe6S1CDDXZIaZLhLUoM6hXuSc5Lcl+RwkksX2f/OJPckuTPJV5K8oP9SJUldjQz3\nJJuAq4FzgZ3AriQ7Fwz7BjBbVS8GPgv8Zd+FSpK669K5nwYcrqoHqupx4AbggvkDqurmqvrRcPMg\nsLXfMiVJy9El3E8EHpq3PTd8bikXA19cbEeS3UkOJTl09OjR7lVKkpal1wuqSd4EzAJXLba/qvZU\n1WxVzc7MzPR5aEnSPJs7jHkYOGne9tbhcz8nydnAe4BXVdVP+ilPkrQSXTr324AdSU5O8nTgQmDf\n/AFJXgp8DDi/qo70X6YkaTlGhntVPQFcAtwE3AvcWFV3J7kiyfnDYVcBvwx8JskdSfYt8XKSpDHo\nMi1DVe0H9i947rJ5j8/uuS5J0ir4CVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtS\ngwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqUKdwT3JOkvuSHE5y6SL7\nfyHJp4f7v55ke9+FSpK6GxnuSTYBVwPnAjuBXUl2Lhh2MfBoVf0q8FfAlX0XKknqrkvnfhpwuKoe\nqKrHgRuACxaMuQD4xPDxZ4GzkqS/MiVJy9El3E8EHpq3PTd8btExVfUE8Bjw3D4KlCQt3+ZxHizJ\nbmD3cPO/k9y3wpfaAny/n6qmhue8MXjOG0CuXNU5v6DLoC7h/jBw0rztrcPnFhszl2Qz8GzgkYUv\nVFV7gD1dCjuWJIeqana1rzNNPOeNwXPeGMZxzl2mZW4DdiQ5OcnTgQuBfQvG7AP+YPj4d4F/qarq\nr0xJ0nKM7Nyr6okklwA3AZuAa6vq7iRXAIeqah/wd8AnkxwGfsDgF4AkaUI6zblX1X5g/4LnLpv3\n+MfA6/st7ZhWPbUzhTznjcFz3hjW/Jzj7IkktcfbD0hSg9Z1uG/E2x50OOd3JrknyZ1JvpKk07Ko\n9WzUOc8b97oklWTqV1Z0Oeckbxi+13cn+dS4a+xbhz/b25LcnOQbwz/f502izr4kuTbJkSR3LbE/\nSf5m+P/jziQv67WAqlqXPwwu3n4L+BXg6cB/ADsXjPkj4Jrh4wuBT0+67jGc86uBZw4fv30jnPNw\n3PHAAeAgMDvpusfwPu8AvgE8Z7j9vEnXPYZz3gO8ffh4J/DgpOte5Tn/JvAy4K4l9p8HfBEIcAbw\n9T6Pv547941424OR51xVN1fVj4abBxl87mCadXmfAd7H4J5FPx5ncWukyzm/Fbi6qh4FqKojY66x\nb13OuYBnDR8/G/juGOvrXVUdYLB6cCkXANfVwEHghCTP7+v46zncN+JtD7qc83wXM/jNP81GnvPw\nn6snVdUXxlnYGuryPp8CnJLkq0kOJjlnbNWtjS7n/F7gTUnmGKzOe8d4SpuY5f59X5ax3n5A/Uny\nJmAWeNWka1lLSZ4GfBC4aMKljNtmBlMzZzL419mBJC+qqh9OtKq1tQvYW1UfSPIKBp+dObWqnpx0\nYdNoPXfuy7ntAce67cEU6XLOJDkbeA9wflX9ZEy1rZVR53w8cCpwS5IHGcxN7pvyi6pd3uc5YF9V\n/bSqvg3czyDsp1WXc74YuBGgqm4FnsHgvjOt6vT3faXWc7hvxNsejDznJC8FPsYg2Kd9HhZGnHNV\nPVZVW6pqe1VtZ3Cd4fyqOjSZcnvR5c/25xl07STZwmCa5oFxFtmzLuf8HeAsgCQvZBDuR8da5Xjt\nA948XDVzBvBYVX2vt1ef9BXlEVebz2PQsXwLeM/wuSsY/OWGwZv/GeAw8G/Ar0y65jGc8z8D/wXc\nMfzZN+ma1/qcF4y9hSlfLdPxfQ6D6ah7gG8CF0665jGc807gqwxW0twB/Naka17l+V4PfA/4KYN/\niV0MvA1427z3+Orh/49v9v3n2k+oSlKD1vO0jCRphQx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIa9L9HPX4TPL5RTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107d100d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Truncated distribution between 1 and 2 (stupid numpy form)\n",
    "true_dist=scipy.stats.truncexpon(loc=0,b=1,scale=1)\n",
    "_=hist(true_dist.rvs(1000000), normed=True, bins=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do an example where we have 20 samples from this distribution as our data sets.\n",
    "\n",
    "Our log-PDF for any given sample is just the sum of the log-pdfs for each sample, as they are independent.\n",
    "\n",
    "Note that I'm using a stupid uniform prior here.  You should use something else, maybe a Jeffreys prior or something.  Depends on the actual problem.\n",
    "\n",
    "It also turns out that the likelihood flattens out basically completely for b>>sample_max, so we need to add in an upper limit or the emcee samples wander off to infinity when they're trying to explore the tails if you're unlucky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sample = true_dist.rvs(20)\n",
    "\n",
    "def log_prob(b_beta, sample):\n",
    "    b,beta=b_beta\n",
    "    if beta<0:\n",
    "        return -np.inf\n",
    "    if b<=sample.max():\n",
    "        return -np.inf\n",
    "    if b>10.0:\n",
    "        return -np.inf\n",
    "    return scipy.stats.truncexpon.logpdf(sample, loc=0,b=b,scale=beta).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can just plot the log-pdf for this example, but Steve doesn't want to do that for some reason.  He loves emcee.\n",
    "Or maybe in the extended version there are more parameter.  Even for three parameters it's worth doing a grid over emcee.\n",
    "\n",
    "But first let's plot it in two slices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHZtJREFUeJzt3Xt4XHW97/H3d2Zyv9/aNE3a0BulLbSUFkuhogKCiNbL\nUbcKylE36vaCbpWz3Xjc6vZ41OP96KOgG1RQ1A0oiigFFNhCQVIuhbZAL/RO25S0TZo0ze27/5hJ\nm4akSTszWZk1n9fzzDOz7t+Zp/3ML7/1m7XM3RERkfCIBF2AiIikloJdRCRkFOwiIiGjYBcRCRkF\nu4hIyCjYRURCRsEugTGzK83sbwOmD5rZtDQfs9HM3Mxi6TzOaCVqmTHKdb9gZjenuybJfAp2SSsz\nO8/MHjazA2bWYmYPmdniodZ192J335Ti4282swtTuU+R8W5ctFoknMysFLgT+DDwGyAXWAYcDrKu\nE2FmMXfvCboOkROhFruk0ywAd7/F3Xvd/ZC7r3D31UOtPLBbwswKzOybZrYl0dr/m5kVJJYtSfwV\nsN/MnjKzVw2zv5uAKcAfEt081wxY/G4z22pme83s2gHbfMHMbjWzm82sFbjSzPLM7DtmtjPx+I6Z\n5SXWP6Y7aYj3UWVmfzCzVjN7zMy+PHh94EIzW594Pz8wMzvOZ5pvZr82szYze9zM5h9nXclSCnZJ\np+eBXjP7mZm9zswqTmDbbwBnAUuBSuAaoM/MJgN/BL6cmP9p4DYzqxm8A3e/AtgKvCHRzfP1AYvP\nA04FLgA+b2anDVi2HLgVKAd+AVwLLAEWAPOBs4HPjfJ9/ABoB2qB9yYeg10GLAbOAN4OXHyc/S0H\n/pP4e/8l8DszyxllLZIlFOySNu7eSjxAHfgx0GxmvzezicfbzswiwPuAq919R6K1/7C7HwYuB+5y\n97vcvc/d7wGagEtPsLwvJv6CeAp4inhg91vp7r9L7P8Q8G7gS+6+x92bgS8CV4x0ADOLAm8F/s3d\nO9x9LfCzIVb9qrvvd/etwF+Jf4EMZ5W73+ru3cC3gHziXzoiRyjYJa3cfZ27X+nu9cA8oA74zgib\nVRMPrI1DLJsKvC3RbbHfzPYT//KYdIKl7RrwugMoHjC9bdC6dcCWAdNbEvNGUkP8PNbA/Q3e90i1\nDHZke3fvA7aPshbJIgp2GTPu/izwU+IBfzx7gU5g+hDLtgE3uXv5gEeRu391uMOeTKmDpncS/0Lp\nNyUxD+LdLIX9C8ysdsB6zUAPUD9gXsNJ1DPQke0Tf9nUD6hFBFCwSxqZ2Wwz+5SZ1SemG4B3Ao8c\nb7tES/QG4FtmVmdmUTM7J3HC8mbgDWZ2cWJ+vpm9qv8YQ9gNJDs2/hbgc2ZWY2bVwOcTdUC8G2eu\nmS0ws3zgCwPeRy9wO/AFMys0s9nAe5Ks5Swze0tiHP4niI8wOu7nKdlHwS7p1Aa8AnjUzNqJB9Az\nwKdGse2ngaeBx4AW4GtAxN23ET+B+K/EW8TbgM8w/L/l/0s8lPeb2adP8n18mXg//upETY8n5uHu\nzwNfAu4F1gODR7x8FCgj3t1yE/EviWSGe94BvAPYR7yf/y2J/naRI0w32hAZO2b2NaDW3YcaHSOS\nEmqxi6RRojvqDIs7G3g/8Nug65Jw0y9PRdKrhHj3Sx3x/v5vEu9OEUkbdcWIiISMumJEREImkK6Y\n6upqb2xsDOLQIiIZa9WqVXvd/WWXzxgskGBvbGykqakpiEOLiGQsM9sy8lrqihERCR0Fu4hIyCjY\nRURCRsEuIhIyCnYRkZBRsIuIhIyCXUQkZDIq2P/y7G5+eP9QN9UREZF+GRXsD214ie/dtx5d30ZE\nZHgZFeyNVYUc6u6luS2Z+xSIiIRbRgX71KoiADa/1BFwJSIi41dGBXvjkWBvD7gSEZHxK6OCva48\nn1jE2KJgFxEZVkYFeywaoaGyUF0xIiLHkVHBDjC1qlAtdhGR48i4YG+sKmLL3g4NeRQRGUZKgt3M\nLjGz58xsg5n9Syr2OZypVYW0He6hpb0rnYcREclYSQe7mUWBHwCvA+YA7zSzOcnudziNGvIoInJc\nqWixnw1scPdN7t4F/ApYnoL9DmlqVSGA+tlFRIaRimCfDGwbML09Me8YZnaVmTWZWVNzc/NJH6y+\nopCIwea9CnYRkaGM2clTd7/e3Re5+6KamhFvsj2s3FiEyRUF6ooRERlGKoJ9B9AwYLo+MS9tGquK\n1BUjIjKMVAT7Y8BMMzvFzHKBfwB+n4L9DmtqlX6kJCIynKSD3d17gI8CdwPrgN+4+5pk93s8jVVF\nHDjUzT4NeRQReZlYKnbi7ncBd6ViX6PRP+TxhZfaqSjKHavDiohkhIz75SnA9AnFAGzcczDgSkRE\nxp+MDPaGigJyosYmDXkUEXmZjAz2WDRCY1WRWuwiIkPIyGAHmFZTxMZmBbuIyGAZG+zTa4rZ8lIH\n3b19QZciIjKuZHSw9/Q5W1s0nl1EZKDMDXaNjBERGVLGBvu0mvhY9o3NGhkjIjJQxgZ7aX4OE0ry\ndAJVRGSQjA12iPezb1Kwi4gcI7ODfUIRG5vbdf9TEZEBMjrYp1UXc+BQNy/pYmAiIkdkdLBrZIyI\nyMtldrAnRsZsUD+7iMgRGR3sdWUFFOVGWb9bwS4i0i+jgz0SMWbVlvDcrragSxERGTcyOtgBTp1Y\nwnO72zQyRkQkIeODfdbEElrau9h7UCNjREQgBMF+am0JAM/vVneMiAiEINhnTYwHu/rZRUTiMj7Y\nq4tzqSzKVYtdRCQh44PdzJg1sZjnFOwiIkAIgh3iI2Oe36WRMSIiEJZgry2lvauX7fsOBV2KiEjg\nQhLs8WvGqJ9dRCQkwT6zf2SMgl1EJBzBXpqfQ11ZvoY8iogQkmAHmD2plHUvtgZdhohI4EIT7HPr\nStnY3E5nd2/QpYiIBCo0wT5nUim9fa7uGBHJeqEJ9rl1ZQCs2anuGBHJbqEJ9obKAkryY6zZeSDo\nUkREAhWaYDcz5kwqVYtdRLJeUsFuZv/PzJ41s9Vm9lszK09VYSdjTl0pz+5qpbdPlxYQkeyVbIv9\nHmCeu58BPA98NvmSTt7cujI6u/t4Ya/ugSoi2SupYHf3Fe7ek5h8BKhPvqSTN7euFNAJVBHJbqns\nY38f8KfhFprZVWbWZGZNzc3NKTzsUTMmFJMbjSjYRSSrxUZawczuBWqHWHStu9+RWOdaoAf4xXD7\ncffrgesBFi1alJZO8JxohFm1xaxVsItIFhsx2N39wuMtN7MrgcuAC3wcXBB97qQyVqzdhbtjZkGX\nIyIy5pIdFXMJcA3wRnfvSE1JyZk3uZR9Hd3s2K9rs4tIdkq2j/37QAlwj5k9aWY/SkFNSTmjPj7i\ncvV2/VBJRLLTiF0xx+PuM1JVSKrMnlRCbjTCU9v2c+npk4IuR0RkzIXml6f98mJRTqsr5clt+4Mu\nRUQkEKELdoAF9WU8s+OAfoEqIlkplMF+Rn057V29bGzWL1BFJPuEMtjnN8RPoKo7RkSyUSiDfVp1\nESV5MZ5SsItIFgplsEcixhkNZRryKCJZKZTBDjC/vpx1L7bqHqgiknVCG+xn1JfT0+esfVHXjRGR\n7BLaYF/QfwJ1q/rZRSS7hDbYa8vyqSvLZ9XWfUGXIiIypkIb7ACLGitp2tzCOLjopIjImAl5sFew\nu/Uw2/fpSo8ikj1CHexnTa0AYNUWdceISPYIdbDPri2lOC9G05aWoEsRERkzoQ72aMQ4c0o5TZvV\nYheR7BHqYAdYNLWS53a3ceBQd9CliIiMidAH++LGCtzhCQ17FJEsEfpgXzClnGjE1B0jIlkj9MFe\nmBtjzqRSnUAVkawR+mAHWNxYyRNb93O4RxcEE5Hwy4pgP2d6FYd7+nhC140RkSyQFcF+9imVRAxW\nbnwp6FJERNIuK4K9rCCHuXVlrNykYBeR8MuKYAdYOr2KJ7bu41CX+tlFJNyyJtiXTK+iu9d13RgR\nCb2sCfbFjZVEI8bKTXuDLkVEJK2yJtiL82LMry/jYZ1AFZGQy5pgh/iwx9XbD3DwcE/QpYiIpE12\nBfu0anr7nMde0K9QRSS8sirYFzVWkBeL8OD65qBLERFJm6wK9vycKEumVfHA8wp2EQmvrAp2gPNn\n1bCpuZ1tLR1BlyIikhbZF+yn1gCo1S4ioZWSYDezT5mZm1l1KvaXTtOqi6ivKFCwi0hoJR3sZtYA\nvBbYmnw56WdmnD+rhoc37KWrpy/ockREUi4VLfZvA9cAnoJ9jYnzZ9XQ3tWrywuISCglFexmthzY\n4e5PjWLdq8ysycyampuD7QZZOqOaWMTUHSMioTRisJvZvWb2zBCP5cC/Ap8fzYHc/Xp3X+Tui2pq\napKtOynFeTEWNVZw/3N7Aq1DRCQdRgx2d7/Q3ecNfgCbgFOAp8xsM1APPG5mtektOTVeM3sCz+5q\n07BHEQmdk+6Kcfen3X2Cuze6eyOwHVjo7rtSVl0aXTQn/v1z77rdAVciIpJaWTeOvd8p1UXMmFDM\nPWsV7CISLikL9kTLPaMudv7aORN59IUWDnR0B12KiEjKZG2LHeCiORPp7XP+qpOoIhIiWR3s8+vL\nmVCSp+4YEQmVrA72SMS44LSJ3P/cHg736CbXIhIOWR3sEO9nb+/q1S3zRCQ0sj7Yl86ooiQvxl2r\nXwy6FBGRlMj6YM+LRblozkTuXrNLFwUTkVDI+mAHuGz+JFo7e/jbBl07RkQyn4IdOG9GDWUFOdz5\nlLpjRCTzKdiB3FiEi+dOZMXa3XR2a3SMiGQ2BXvCZWfUcfBwjy7lKyIZT8GesHR6FRWFOdyp0TEi\nkuEU7AmxaITXnT6Je9fu5uDhnqDLERE5aQr2Ad5y5mQOdffyp6fVaheRzKVgH+CsqRU0VhVy2+Pb\ngy5FROSkKdgHMDPeurCeRza16M5KIpKxFOyDvHnhZABuf3xHwJWIiJwcBfsg9RWFLJ1exW2Pb8fd\ngy5HROSEKdiH8NaF9Wxt6eCxzfuCLkVE5IQp2IfwutNrKcqN8uvHtgVdiojICVOwD6EwN8abzpzM\nnat3sr+jK+hyREROiIJ9GO9+xVQO9/Rx6yoNfRSRzKJgH8aculIWTinnl49u1UlUEckoCvbjuHzJ\nVDbtbWelbpsnIhlEwX4cl54+ifLCHG5+dEvQpYiIjJqC/Tjyc6K87ax6VqzZze7WzqDLEREZFQX7\nCK5Y0kifOz97eHPQpYiIjIqCfQRTqgq5eG4tNz+yhXZdzldEMoCCfRQ+sGwarZ09/GeTfrAkIuOf\ngn0UzppawcIp5dzw0GZ6+zT0UUTGNwX7KH1g2TS2tnSwYs2uoEsRETkuBfsoXTy3lobKAq57cJN+\nsCQi45qCfZSiEeOqV07nyW37eWiDfrAkIuNX0sFuZh8zs2fNbI2ZfT0VRY1Xb19UT21pPt+973m1\n2kVk3Eoq2M3s1cByYL67zwW+kZKqxqm8WJQPv2o6j23ex8pNarWLyPiUbIv9w8BX3f0wgLvvSb6k\n8e0dixuYUJLH9+5bH3QpIiJDSjbYZwHLzOxRM3vAzBanoqjxLD8nygfPn84jm1p4VK12ERmHRgx2\nM7vXzJ4Z4rEciAGVwBLgM8BvzMyG2c9VZtZkZk3Nzc0pfRNj7V1nT6GmJI9vrHhOfe0iMu6MGOzu\nfqG7zxvicQewHbjd4/4O9AHVw+znendf5O6LampqUvsuxlhBbpSrL5jJY5v3cd+60Pc+iUiGSbYr\n5nfAqwHMbBaQC+xNtqhM8I7FDUyrLuJrf36Wnt6+oMsRETki2WC/AZhmZs8AvwLe61nSN5ETjfCZ\ni09l/Z6D3P74jqDLERE5IpbMxu7eBVyeoloyziXzalnQUM637nmeN8yvoyA3GnRJIiL65WkyzIzP\nvm42u1o7ue7BjUGXIyICKNiT9oppVbz+jEn88P6NbGvpCLocEREFeyp87vWnETHj3+9cG3QpIiIK\n9lSYVFbAxy6YwYq1u7n/OQ1/FJFgKdhT5P3nncIp1UV88Q9r6ezuDbocEcliCvYUyYtF+dLyubyw\nt53v/2VD0OWISBZTsKfQspk1vHVhPT98YCPP7DgQdDkikqUU7Cn2vy87jcqiXK65dTXd+kWqiARA\nwZ5i5YW5/Pvyeax9sZXrHtDYdhEZewr2NLhkXi2vP2MS371vvbpkRGTMKdjT5MvL51FZlMvHf/UE\nHV09QZcjIllEwZ4mFUW5fPvtC3hhb7t+uCQiY0rBnkZLZ1TzofOnc8vft/Gnp18MuhwRyRIK9jT7\n54tmMb++jGtuW80Le9uDLkdEsoCCPc1yohG+/66FxCLGB29qov2w+ttFJL0U7GOgobKQ773zTDbs\nOcg1t63WfVJFJK0U7GNk2cwaPnPxbP64+kWue3BT0OWISIgp2MfQh86fxutPn8TX/vwsf35GJ1NF\nJD0U7GPIzPjG2+azoKGcq3/1JKu27Au6JBEJIQX7GCvIjfKT9yyitiyff/x5E5s1UkZEUkzBHoCq\n4jxuvHIx7s57b/w7u1s7gy5JREJEwR6QaTXF3HDlYva2HeZdP36EvQcPB12SiISEgj1AZ06p4IYr\nF7Nj/yEu/8mj7O/oCrokEQkBBXvAXjGtih+/ZxGb9rZzxX/8nX3tCncRSY6CfRxYNrOG6y4/i+d2\nt/GO61eqz11EkqJgHydePXsCP/2fi9mx7xBv+9FKtrV0BF2SiGQoBfs4snR6Nb/4xyW0dnbzP370\nsG7SISInRcE+zixoKOc3HzyHqBlvv24l96zdHXRJIpJhFOzj0KyJJfzuI+cyc0IxV93UxPUPbtSF\nw0Rk1BTs49SE0nx+/cFzuHTeJL5y17P882+e0i32RGRUFOzjWH5OlP//zjP55IWz+N2TO3jTDx5i\nw562oMsSkXFOwT7ORSLG1RfO5Ob3v4KW9i7e+P2H+O0T29U1IyLDUrBniHNnVPPHjy9jXl0Zn/z1\nU3z0l0/Qoh8zicgQkgp2M1tgZo+Y2ZNm1mRmZ6eqMHm5iaX53HLVEq655FRWrN3Fa7/9IPdq1IyI\nDJJsi/3rwBfdfQHw+cS0pFE0YvzTq2bw+4+eR01JHh/4eRMfu+UJ9ujXqiKSkGywO1CaeF0G7Exy\nfzJKp00q5Y6PnMsnLpzJ3Wt28ZpvPsCND71AT29f0KWJSMAsmZNwZnYacDdgxL8klrr7lmHWvQq4\nCmDKlClnbdky5GpyEl7Y286//X4NDz7fzGmTSrn20tM4b2Z10GWJSIqZ2Sp3XzTieiMFu5ndC9QO\nseha4ALgAXe/zczeDlzl7heOdNBFixZ5U1PTSKvJCXB3/vTMLv7PH9exY/8hls2s5n9dMpt5k8uC\nLk1EUiRlwT7CQQ4A5e7uZmbAAXcvHWk7BXv6dHb3cvMjW/j+Xzewv6Ob5Qvq+OirZzBzYknQpYlI\nkkYb7Mn2se8Ezk+8fg2wPsn9SZLyc6J8YNk0HvjMq/nQ+dNZsWY3F337QT500yqe3q6Liolkg2Rb\n7OcB3wViQCfwT+6+aqTt1GIfOy3tXdz40Av89OHNtHX2sGxmNe879xTOn1VDJGJBlyciJ2BMumJO\nloJ97LV1dnPTI1u48aHNNLcdZmpVIVcsmcrbzmqgrDAn6PJEZBQU7DKkrp4+7l6zi5+v3Mxjm/eR\nnxPh0nmTeMvCes6ZXkVUrXiRcUvBLiNas/MANz+ylTtX76Sts4fa0nzedOZk3nzmZGZNLCZ+PlxE\nxgsFu4xaZ3cv963bw+2Pb+f+55vp7XOmVRfx2rm1XDKvlvn1ZQp5kXFAwS4npbntMHev2cXda3ax\ncuNL9PQ5taX5XHDaBF45q4al06soyVefvEgQFOyStP0dXdy3bg9/XrOLhzbspaOrl2jEWDilnGUz\nazh3RjWnTy4jN6aLhIqMBQW7pFRXTx+Pb93Hf61v5r/W7+XpHQdwh7xYhPkN5SyaWsHixkoWTqnQ\nKBuRNFGwS1q1tHfx6KaXaNqyj6bNLazZ2UpPX/zf0rSaIubVlTFvcinz6sqYW1emsBdJAQW7jKmO\nrh6e3LafVZv3sXrHAdbsOMDOA0cvJdxQWcDs2lJmTChmRk0xMyYUM31CMcV5sQCrFsksow12/a+S\nlCjMjbF0ejVLpx+9quRLBw+zZmcra3a28syOAzy3u42/PrvnSMseYFJZPjMmFNNYVURDZQENFYU0\nVMYfZQVq5YucDAW7pE1VcR6vnFXDK2fVHJnX3dvHlpc62LDnIBubD7JhT/xxx7YdtHb2HLN9aX6M\nhspC6isKqC3NZ0JpPhNL85lYmndkujQ/pqGYIoMo2GVM5UQj8e6YCcUvW3bgUDfbWjrYvq+DbS2H\n2NrSwbZ9HWxqbueRTS0cONT9sm3ycyJMLM2nsiiXysJcKopyqSzKpaIwl6qi/ukcKgpzKS/MpTgv\nplE8EnoKdhk3ygpyKJtcNuw15A919bKnrZPdrYfZ3do54HGYfR1d7GrtZN2LrbzU3sXhnuHvJJUX\ni1CSn0NpfoyS/BjF+TFK8nIoyY9Rkh9/Ls6LUZAbpSAnSmFulPwBrwtyouT3v86Nkh+L6oJqMq4o\n2CVjFORGmVpVxNSqohHXPdTVS0tHFy0Hu2jp6GJfexf7O7po6+zh4OEeWjt7aOvsPjLd3HaQts6e\nI9MnKj8nQn5OlNxohNxY5OhzLEJO9Oh0TjRCXixCTtSOXR6LkBeNEI1EiEWNaMSIRQY+R4hGiC+P\nDFoeNSJmxCKRI9MDl0es/wE26DlihiWe+18PnI4YGIZFODJ9ZD0GTatLbNxQsEsoFeRGmZxbwOTy\nghPetrfP6ejq4VB3L4e6eo8+97/u7qWjq5fOxPOhxOtD3b109fTFH71Hn7sTrzsO9S/vpbvXjy7v\n6eNwYp1MFg/7o18KJKb7896wI+vFp49+GRz5SrCjr48sO876NmDDgesdrefl+xi472OWjbD+kF9b\nQ8wcPGvwF95X3nw6Z59SOdTeUkbBLjJINGKJLpmxHZXj7vR5/Iult8/p6etLPPvR516n153evj56\n+pyeXj9mncHbdvc67o4DfYn9x4/j9PVxZH7/seOvOWa6fxsfMN2/TXzesftwjq5DYgCUD3iP8eeB\n8/rXcQaPvj6y/qD1jt2OAdsN2P8o1vcBBb6sxgH1DTUofKih4i+bM8SGRXnRIfaWWgp2kXHCzIga\nAy6dnP4AkHDS8AARkZBRsIuIhIyCXUQkZBTsIiIho2AXEQkZBbuISMgo2EVEQkbBLiISMoHcaMPM\nmoEtY37g1KsG9gZdxDihz+JY+jyO0mdxrGQ+j6nuXjPSSoEEe1iYWdNo7maSDfRZHEufx1H6LI41\nFp+HumJEREJGwS4iEjIK9uRcH3QB44g+i2Pp8zhKn8Wx0v55qI9dRCRk1GIXEQkZBbuISMgo2E+Q\nmTWY2V/NbK2ZrTGzq4OuKWhmFjWzJ8zszqBrCZqZlZvZrWb2rJmtM7Nzgq4pSGb2ycT/k2fM7BYz\nyw+6prFiZjeY2R4ze2bAvEozu8fM1ieeK9JxbAX7iesBPuXuc4AlwEfMbE7ANQXtamBd0EWME98F\n/uzus4H5ZPHnYmaTgY8Di9x9HvFbQv1DsFWNqZ8Clwya9y/Afe4+E7gvMZ1yCvYT5O4vuvvjiddt\nxP/jTg62quCYWT3weuAnQdcSNDMrA14J/AeAu3e5+/5gqwpcDCgwsxhQCOwMuJ4x4+4PAi2DZi8H\nfpZ4/TPgTek4toI9CWbWCJwJPBpsJYH6DnAN0Bd0IePAKUAzcGOia+onZlYUdFFBcfcdwDeArcCL\nwAF3XxFsVYGb6O4vJl7vAiam4yAK9pNkZsXAbcAn3L016HqCYGaXAXvcfVXQtYwTMWAh8EN3PxNo\nJ01/ameCRP/xcuJfeHVAkZldHmxV44fHx5qnZby5gv0kmFkO8VD/hbvfHnQ9AToXeKOZbQZ+BbzG\nzG4OtqRAbQe2u3v/X3C3Eg/6bHUh8IK7N7t7N3A7sDTgmoK228wmASSe96TjIAr2E2RmRrwPdZ27\nfyvoeoLk7p9193p3byR+Uuwv7p61LTJ33wVsM7NTE7MuANYGWFLQtgJLzKww8f/mArL4ZHLC74H3\nJl6/F7gjHQdRsJ+4c4EriLdOn0w8Lg26KBk3Pgb8wsxWAwuArwRcT2ASf7ncCjwOPE08b7Lm8gJm\ndguwEjjVzLab2fuBrwIXmdl64n/RfDUtx9YlBUREwkUtdhGRkFGwi4iEjIJdRCRkFOwiIiGjYBcR\nCRkFu4hIyCjYRURC5r8B3okY2NtVGwgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d07bf90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYFFW+xvHvbxLDkHPOOSsMSA4SRV0UM64Kq6KAgBFX\n3XV117t71zULoqAYUNKqBAURUXJ0QHKQnMMgOTPMuX904x1xhhmY6a6e6ffzPPPQXVXd59dl+Xb1\nqVNV5pxDRERyvgivCxARkeBQ4IuIhAkFvohImFDgi4iECQW+iEiYUOCLiIQJBb4EjZn1NLO5KZ4f\nN7PKAW6zopk5M4sKZDsZ5a+lagaXfcHMPg10TRI+FPiSpcyspZnNN7MjZnbQzOaZWePUlnXO5XXO\nbc7i9reaWYesfM/sSF8WkpqQ2OuRnMHM8gNfA32AcUAM0Ao442Vdl8PMopxzSV7XIRII2sOXrFQd\nwDk32jl33jl3yjk3zTm3IrWFU3ZvmFluM3vVzLb5fx3MNbPc/nlN/b8aDpvZcjNrm8b7jQTKA1/5\nu4sGpZh9t5ltN7MDZvZcite8YGafm9mnZnYU6GlmuczsDTPb7f97w8xy+Zf/TbdUKp+jiJl9ZWZH\nzexHM3vp4uWBDma2wf95hpiZXWKdxprZWDM7ZmZLzaxBinZLm9kXZpZoZlvMbIB/ehfgWeAO/3pY\n7p/ey8zW+t9rs5k9dIl2JQdS4EtW+hk4b2Yfm9l1ZlboMl77CtAIaA4UBgYByWZWBpgMvOSf/iTw\nhZkVu/gNnHP3ANuBG/3dRS+nmN0SqAG0B543s1op5nUDPgcKAp8BzwFNgauABkAT4C8Z/BxDgBNA\nSeA+/9/FbgAaA/WB24HOl3i/bsB/8X32UcAEM4s2swjgK2A5UMb/uR41s87OuanAP4Gx/vVw4Uti\nv7/t/EAv4HUza5jBzyU5gAJfsoxz7ii+YHXAcCDRzCaZWYlLvc4fXn8CBjrndvl/Hcx3zp0B/ghM\ncc5Ncc4lO+e+AxKArpdZ3ov+XxzL8YVkgxTzFjjnJvjf/xRwN/B359x+51wi8CJwT3oNmFkkcAvw\nN+fcSefcGuDjVBb9X+fcYefcdmAGvi+WtCxxzn3unDsHvAbE4vsyagwUc8793Tl31n8sZDhwZ1pv\n5Jyb7Jzb5HxmAdPwdblJmFDgS5Zyzq11zvV0zpUF6gKlgTfSeVlRfEG2KZV5FYDb/N0fh83sML4v\nlVKXWdreFI9PAnlTPN9x0bKlgW0pnm/zT0tPMXzHxVK+38XvnV4tF/v19c65ZGCnv5YKQOmL1suz\nQJpfrv5fXQv9B9MP4/vSLJrOZ5IcRAdtJWCcc+vM7CMgvb7iA8BpoAq+ve+UdgAjnXMPZrTZyyoy\n9dfsxheoq/3Py/unga+7Ju7CgmZWMsXrEoEkoCy+7i2AcldQT0q/vt7/S6isv5YkYItzrloar/vN\nZ/Ifg/gCuBeY6Jw7Z2YTgEsdP5AcRnv4kmXMrKaZPWFmZf3PywF3AQsv9Tr/nusI4DX/gchIM2vm\nD6lPgRvNrLN/eqyZtb3QRir2AZkd2z8a+IuZFTOzosDz/jrA94VUx8yuMrNY4IUUn+M88CXwgpnF\nmVlNfAGbGY3MrLv5ziN4FN+Ip4XAYuCYmT3tP+AdaWZ1UwyB3QdU9H9JgG/EVC78X0pmdh3QKZO1\nSTajwJesdAy4BlhkZifwBdMq4IkMvPZJYCXwI3AQ+DcQ4Zzbge/A5bP4wmoH8BRpb7v/whfWh83s\nySv8HC/hO06wwl/TUv80nHM/A38HpgMbgItH4DwCFMDXbTMS35dHZoalTgTuAA7hO47Q3Tl3zv/l\ncgO+/v8t+H4lve9vG3wHegF+MbOlzrljwAB8w2UPAT2ASZmoS7Ih0w1QRALHzP4NlHTOpTZaRySo\ntIcvkoX83Vr1zacJcD8w3uu6REAHbUWyWj583Til8fWjv4qvW0bEc+rSEREJE+rSEREJEyHVpVO0\naFFXsWJFr8sQEclWlixZcsA597vLjVwspAK/YsWKJCQkeF2GiEi2Ymbb0l9KXToiImFDgS8iEiYU\n+CIiYUKBLyISJhT4IiJhQoEvIhImFPgiImEiRwT+6XPneWHSag6dOOt1KSIiIStHBP6KnUcYtWg7\nt7w7nx0HT3pdjohISMoRgd+kUmE+ub8JB46dofvQ+azadcTrkkREQk6OCHyAppWL8Hmf5kRHGHe8\nt4A5GxK9LklEJKQEPPDNrIuZrTezjWb250C2Vb1EPr7s24JyhePo9eGPfLl0ZyCbExHJVgIa+GYW\nCQwBrgNqA3eZWe1AtlmyQCzjHm5Gk0qFeXzccobM2Iiu+S8iEvg9/CbARufcZufcWWAMvhtSB1T+\n2Gg+6tWEbleV5j/fruevE1dxPlmhLyLhLdCXRy4D7EjxfCdwTcoFzKw30BugfPnyWdZwTFQEr99+\nFSULxPLerM3sO3qGt+68mtwxkVnWhohIduL5QVvn3DDnXLxzLr5YsXSv339ZIiKMZ66rxQs31mb6\n2n30eH8hBzVWX0TCVKADfxdQLsXzsv5pQdWzRSXe6dGQ1buPcutQjdUXkfAU6MD/EahmZpXMLAa4\nE5gU4DZTdV29Unz2wDX8cuIsN78zn5U7NVZfRMJLQAPfOZcEPAJ8C6wFxjnnVgeyzUtpXLEwX/Rp\nRq6oCO4YtoCZ6/d7VYqISNAFvA/fOTfFOVfdOVfFOfc/gW4vPVWL5+PLvs2pUCQP93+cwH8TdqT/\nIhGRHMDzg7ZeKJE/lnEPNaVZ5SI89fkK3v5+g8bqi0iOF5aBD5AvNpoRPRvT/eoyvPrdzzw7fhVJ\n55O9LktEJGACPQ4/pMVERfDq7Q0oWSCWd2ZuYs+RUwzu0ZC8ucJ6tYhIDhW2e/gXmBmDutTkX93r\nMWfDAW57dwF7j5z2uiwRkSwX9oF/wV1NyjOiZ2N2HDzJTUPmsWb3Ua9LEhHJUgr8FNpUL8Z/H26G\nGdz27nxmaNimiOQgCvyL1CqVn/F9W1ChSB4e+DiBUYu2e12SiEiWUOCn4sIllltXK8qz41fyr2/W\nkqyrbYpINqfAT0PeXFEMvzeeu68pz3uzNtN/9E+cPnfe67JERK6Yxh9eQlRkBC/dVJcKReL455R1\n7DlyiuH3xlMkby6vSxMRuWzaw0+HmdG7dRXeudt3tc3uQ+ezOfG412WJiFw2BX4Gda1XitG9m3L8\ndBLdh85n8ZaDXpckInJZFPiXoWH5Qozv24LCeWL44/uLmLgs6Jf2FxG5Ygr8y1S+SBxf9mnOVeUL\nMnDMMt7ShddEJJtQ4F+BgnExjLy/Cd0bluG1737m0bHLNIJHREJewALfzP5jZuvMbIWZjTezgoFq\nywu5oiJ59bYGPNW5BhOX7abH8IUkHjvjdVkiImkK5B7+d0Bd51x94GfgmQC25Qkzo1+7qrz7x4as\n2XOUm4bMY91eXYNHREJTwALfOTfNf4tDgIX4bmCeI3WpW4r/PtScpORkbnlnPj+s2+d1SSIivxOs\nPvw/Ad8EqS1P1CtbgIn9WlKpmO8aPO/P2ayDuSISUjIV+GY23cxWpfLXLcUyzwFJwGdpvEdvM0sw\ns4TExMTMlOO5kgViGfdQMzrVLslLk9fy7PhVnNNdtEQkRFgg90LNrCfwENDeOXcyveXj4+NdQkJC\nwOoJluRkx6vfrWfIjE00r1KEoXc3okBctNdliUgOZWZLnHPx6S0XyFE6XYBBwB8yEvY5SUSE8VTn\nmrx2ewMSth7i5nfmseXACa/LEpEwF8g+/MFAPuA7M1tmZu8GsK2Q1L1hWT578BoOnzrHTUPmMX/T\nAa9LEpEwFshROlWdc+Wcc1f5/x4OVFuhrHHFwkzo24Li+XJx7weLdUMVEfGMzrQNgvJF4viib3Na\nVPXdUOX5iTqYKyLBp8APkvyx0Yzo2ZjerSvzyYJt3PvBYg6eOOt1WSISRhT4QRQZYTzbtRav3d6A\nJdsP0W3IXJ2ZKyJBo8D3QPeGZRn3UDPOnEum+zvz+Xb1Xq9LEpEwoMD3yFXlCvJV/5ZUK5GPh0Yu\n0WWWRSTgFPgeKpE/lrG9m9L9at9llvuNWsrJs0npv1BE5Aoo8D0WGx3Jq7c34LmutZi6ai+3DF3A\nzkNhdZ6aiASJAj8EmBkPtq7MiJ6N2XnoJN0Gz9M9c0UkyynwQ0jbGsWZ0K8FBXJH02P4Qp2kJSJZ\nSoEfYqoUy8v4fi1+PUnrLxNWcjZJJ2mJSOYp8ENQgdy+k7Qeal2ZTxdu567hC9l/9LTXZYlINqfA\nD1GREcYzXWsxuMfVrNl9lBvensuSberXF5Erp8APcTfUL82Efi3IHRPJncMWMnLhNo3XF5ErosDP\nBmqUzMekfi1pWbUof52wiqe/WMHpc+e9LktEshkFfjZRIC6aD+5rzIBrqzIuYSd3vLeA3YdPeV2W\niGQjCvxsJCLCeLxTDYbd04hNiSe48e25LNj0i9dliUg2EfDAN7MnzMyZWdFAtxUuOtUpyYR+LSgY\nF80fP1jEB3O3qF9fRNIV0MA3s3JAJ0BnEGWxqsXzMqFfCzrUKs4/vl7Do2OXceqs+vVFJG2B3sN/\nHd+NzLX7GQD5YqMZencjnupcg0nLd9N96Hy2/6Lr8IhI6gIW+GbWDdjlnFueznK9zSzBzBISExMD\nVU6OFRFh9GtXlQ97Nmb34VPc8PYcfli3z+uyRCQEWWb6fs1sOlAylVnPAc8CnZxzR8xsKxDvnDtw\nqfeLj493CQkJV1xPuNv+y0n6fLaE1buP0q9dFR7vWIPICPO6LBEJMDNb4pyLT2+5qMw04pzrkEbj\n9YBKwHIzAygLLDWzJs453d4pQMoXieOLPs158avVDJmxiZ+2H+atu66maN5cXpcmIiEgIF06zrmV\nzrnizrmKzrmKwE6gocI+8GKjI/lX9/r859b6LNl2iOvfmkPCVl2SQUQ0Dj/Hui2+HOP7tiB3tO+S\nDO/P2ayhmyJhLiiB79/Tv2T/vWS92qXzM6l/S66tWZyXJq+l36ilHDt9zuuyRMQj2sPP4fLHRvPe\nPY14tmtNvl29j26D57F+7zGvyxIRDyjww4CZ0bt1FUY9cA3HziRx05B5jP9pp9dliUiQKfDDyDWV\nizB5QEvqlS3AY2OX8+z4lbrqpkgYUeCHmeL5Yhn1wDU81KYyoxZt59Z357PtlxNelyUiQaDAD0NR\nkRE8c10tht8bz/ZfTnLDW3OZvGKP12WJSIAp8MNYx9olmDKwFVWK56XfqKX8dcIqdfGI5GAK/DBX\ntlAc4x5qxoOtKjFy4TZuGTqfLQfUxSOSEynwhZioCJ67vjYf3BfPrsOnuPHtuXy1fLfXZYlIFlPg\ny6/a1yrB5AGtqFEyH/1H/6RRPCI5jAJffqNMwdyM6d3011E8N78zn82Jx70uS0SygAJffifaP4rn\nw56N2XvE18Uzcdkur8sSkUxS4Eua2tUszpSBrahdOj8DxyzjmS9XqItHJBtT4MsllSqQm9EPNqVv\n2yqMXrxD1+IRycYU+JKuqMgIBnWpySd/asIvJ87yh8FzGblwmy63LJLNKPAlw1pXL8Y3A1vRtHIR\n/jphFQ+NXMKhE2e9LktEMiiggW9m/c1snZmtNrOXA9mWBEexfLn4sGdj/nJ9LWas3891b85h4eZf\nvC5LRDIgYIFvZu2AbkAD51wd4JVAtSXBFRFhPNCqMuP7tiAuJpK7hi/k1WnrSTqf7HVpInIJgdzD\n7wP8r3PuDIBzbn8A2xIP1C1TgK/6t+TWhmV5+4eN3P7eAnYcPOl1WSKShkAGfnWglZktMrNZZtY4\ntYXMrLeZJZhZQmJiYgDLkUDIkyuK/9zWgLfvupoN+47T9a05uiyDSIjKVOCb2XQzW5XKXzcgCigM\nNAWeAsaZmV38Hs65Yc65eOdcfLFixTJTjnjoxgalmTKwFdWK56X/6J8Y9PlyTp5N8rosEUkhKjMv\nds51SGuemfUBvnS+sXuLzSwZKApoNz6HKlfYd+XNN7/fwOAZG0nYeoi37rqaumUKeF2aiBDYLp0J\nQDsAM6sOxAAHAtiehICoyAie6FSDUQ805eTZ89z8zjyGztzE+WSN2RfxWiADfwRQ2cxWAWOA+5zO\n1AkbzaoU4ZuBrehYuwT/nrqOu4YvZOchHdAV8ZKFUgbHx8e7hIQEr8uQLOSc48ulu/jbpNUY8I+b\n6tLtqtKkcjhHRK6QmS1xzsWnt5zOtJWAMjNuaVSWbwb6rrP/6Nhl9B/9E0dOnvO6NJGwo8CXoChX\nOI6xDzXjqc41mLpqL13enM38jTqkIxJMCnwJmsgIo1+7qozv24LcMZH0eH8R/zN5DWeSdMllkWBQ\n4EvQ1StbgMn9W3FvswoMn7OFboPnsW7vUa/LEsnxFPjiidwxkfy9W10+7NmYA8fP8oe35/H+nM0k\na/imSMAo8MVT7WoW59tHW9GmRjFemryWe0YsYvfhU16XJZIjKfDFc0Xy5mLYPY349y31+Gn7YTq/\nPpvPl+zUDVZEspgCX0KCmXFH4/JMHdiaWqXy8+R/l9N75BISj53xujSRHEOBLyGlfJE4xvRuyl+u\nr8WsnxPp9Pospqzc43VZIjmCAl9CzoUbrEwZ0JJyhePo+9lSBoz+icMndTtFkcxQ4EvIqlo8H1/2\nac4THaszZeUeOr0+mxnrdB8dkSulwJeQFhUZQf/21ZjQrwWF88TQ66Mf+fMXKzh2WpdmELlcCnzJ\nFuqWKcDER1rQt20VxiXsoMsbc5i/SZdmELkcCnzJNnJFRTKoS03++3BzYqIi6DF8ES9MWs2ps7o0\ng0hGKPAl22lUoRBTBrSiZ/OKfDR/K13fmsOPWw96XZZIyAtY4JvZVWa20MyW+W9S3iRQbUn4yR0T\nyQt/qMOoB67h3Plkbn9vAS9MWq376IpcQiD38F8GXnTOXQU8738ukqWaVy3Kt4+25t6mFfho/lY6\nvzFbffsiaQhk4Dsgv/9xAWB3ANuSMJYnVxQvdqvL2N5NiTCjx/BFPDd+JcfPaG9fJKWA3eLQzGoB\n3wKG74uluXNuWyrL9QZ6A5QvX77Rtm2/W0Qkw06dPc8r09YzYt4WShfIzb+616N19WJelyUSUBm9\nxWGmAt/MpgMlU5n1HNAemOWc+8LMbgd6O+c6XOr9dE9bySpLth1i0OfL2ZR4gjviy/HcDbXIHxvt\ndVkiARGUwE+ngCNAQeecM98dq4845/Jf6jUKfMlKp8+d543pGxg2exPF88Xyr+71aFezuNdliWS5\nULiJ+W6gjf/xtcCGALYl8jux0ZH8+bqajO/bgvy5o+j10Y88Pm6ZrskjYSuQgf8g8KqZLQf+ib+f\nXiTYGpQryFf9W9L/2qpMXLabjq/PZuqqvV6XJRJ0AevSuRLq0pFAW7XrCIM+X8GaPUfpUqckL3ar\nQ4n8sV6XJZIpodClIxJyLlyTZ1CXGsxYv58Or81i9OLtupeuhAUFvoSd6MgI+ratytRHW1OndH6e\n+XIldw5fyKbE416XJhJQCnwJW5WK5mH0g0359y31WLfnKNe9OYfBP2zgbFKy16WJBIQCX8LahXvp\nTn+iDR1rleCVaT/zh8FzWbbjsNeliWQ5Bb4IUDxfLEPubsiwexpx+OQ5bn5nHi9+tZoTujyD5CAK\nfJEUOtUpyXePt+aP11Tgw3lbfbdVXK/bKkrOoMAXuUi+2Gj+cVNdPn+4GbljIun14Y8MHPMTvxw/\n43VpIpmiwBdJQ3zFwkwe0JKB7asxZeUe2r82i7E/aginZF8KfJFLyBUVyWMdqzNlQCuqF8/H01+s\n5I5hC/h53zGvSxO5bAp8kQyoViIfY3o35eVb67Nx/3G6vjmHf09dp/vpSraiwBfJoIgI4/b4cnz/\nRFtuuroMQ2duouPrs5ixTgd1JXtQ4ItcpsJ5YnjltgaM6d2U2OhIen30I30/W8LeI6e9Lk3kkhT4\nIleoaeUiTBnQiic7Vef7tftp/+pMRszdwnkd1JUQpcAXyYSYqAgeubYa0x5rTaOKhfn712voNmQu\nK3bqTF0JPQp8kSxQoUgePu7VmME9rmbf0TN0GzKPv01cxdHT57wuTeRXCnyRLGJm3FC/NN8/0YZ7\nm1bgk4Xb6PDqLCYu20Uo3XdCwlemAt/MbjOz1WaWbGbxF817xsw2mtl6M+ucuTJFso/8sdG82K0u\nE/q2oET+WAaOWcZdwxdq7L54LrN7+KuA7sDslBPNrDZwJ1AH6AK8Y2aRmWxLJFtpUK4gE/q14KWb\n6rJ2zzG6vjmH/5m8huO6IJt4JFOB75xb65xbn8qsbsAY59wZ59wWYCPQJDNtiWRHkRHGH5tWYMaT\nbbm1UVmGz9lC+1dnMmn5bnXzSNAFqg+/DLAjxfOd/mm/Y2a9zSzBzBISExMDVI6ItwrnieF/b6nP\n+L7NKZ4vlgGjf6LH8EXq5pGgSjfwzWy6ma1K5a9bVhTgnBvmnIt3zsUXK1YsK95SJGRdXb7Qr908\na/Ycpeubc/jnlLXq5pGgiEpvAedchyt4311AuRTPy/qniYS9C908XeuV4uWp6xg2ezMTl+3iuetr\nc2P9UpiZ1yVKDhWoLp1JwJ1mlsvMKgHVgMUBakskW0qrm2eDunkkQDI7LPNmM9sJNAMmm9m3AM65\n1cA4YA0wFejnnNNlBUVScXE3z3X+0TzHdNKWZDELpZEC8fHxLiEhwesyRDxz8MRZXp66jrEJOyiS\nJ4ZBnWtya6OyRESom0fSZmZLnHPx6S2nM21FQsiFbp6J/VpQvnAcg75YwU3vzGPJtkNelyY5gAJf\nJATVL1uQL/o05407rmLf0dPcMnQ+j41dxr6jugSzXDkFvkiIMjNuuroMPzzRln7tqjB5xR7avTKT\nITM2cvqcDonJ5VPgi4S4PLmieKpzTaY/3oaWVYvyn2/X0+n12UxbvVdn68plUeCLZBPli8Qx7N54\nPr3/GnJFRdB75BLuHbFYwzglwxT4ItlMy2pFmTKwFX+7sTbLdxymy5tzePGr1Rw5pWGccmkKfJFs\nKDoygl4tKjHjybbc0bgcH83fSrtXZvLZom0knU/2ujwJUQp8kWysSN5c/PPmenzdvyVVi+flufGr\nuP6tuczZoAsRyu8p8EVygDqlCzC2d1OG3t2QU+fOc88Hi+n14WI27lf/vvw/Bb5IDmFmXFevFN89\n3ppnu9YkYeshOr8xh+cnruLgibNelychQIEvksPkioqkd+sqzHyqLT2alOezRdtp858ZDJu9iTNJ\nGr8fzhT4IjlUkby5+MdNdZk6sBXxFQrxzynr6PjabL5ZuUfj98OUAl8kh6tWIh8f9mrCJ39qQmx0\nBH0+W8od7y1kxc7DXpcmQabAFwkTrasXY8qAVvzz5npsPnCcPwyex+Njl7HnyCmvS5MgUeCLhJGo\nyAh6XFOeGU+2pU/bKny90nd9ntemrddtFsNAZm+AcpuZrTazZDOLTzG9o5ktMbOV/n+vzXypIpJV\n8sVG83SXmnz/eBs61CrBWz9spO1/ZvLpwm2c04lbOVZm9/BXAd2B2RdNPwDc6JyrB9wHjMxkOyIS\nAOUKxzG4R0O+7NucykXz8JcJq+j8+mymrtKF2XKiTAW+c26tc259KtN/cs7t9j9dDeQ2s1yZaUtE\nAqdh+UKMfagpw++NJyLCePjTJdz67gKWbDvodWmShYLRh38LsNQ5dya1mWbW28wSzCwhMVGng4t4\nxczoWLsEUwe24l/d67H94EluGbqAh0cuYXPica/LkyyQ7j1tzWw6UDKVWc855yb6l5kJPOmcS7jo\ntXWASUAn59ym9IrRPW1FQsfJs0m8P2cL783axOmkZHo0Kc+A9tUolk8/1kNNRu9pG5XeAs65DldY\nQFlgPHBvRsJeREJLXEwUA9pX464m5Xnr+w2MXrydL5fupHfrKjzQqhJ5cqUbHxJiAtKlY2YFgcnA\nn51z8wLRhogER7F8vjN2pz3WmtbVi/H69J9p+8pMRi3arksxZzOZHZZ5s5ntBJoBk83sW/+sR4Cq\nwPNmtsz/VzyTtYqIhyoXy8vQPzbiiz7NqFA4jmfHr6TzG7rVYnaSbh9+MKkPXyR7cM4xbc0+/j11\nHZsTT9CwfEEGdalJ08pFvC4tLGW0D19n2orIZTMzOtcpybRHW/O/3eux+/Bp7hy2kHtHLGbVriNe\nlydp0B6+iGTa6XPn+WTBVt6ZuYnDJ89xQ/1SPNGpBpWK5vG6tLCQ0T18Bb6IZJmjp88xfPZmPpi7\nhTNJydweX46B7atRskCs16XlaAp8EfFM4rEzDJmxkc8WbSPCjJ7NK/JwmyoUyhPjdWk5kgJfRDy3\n4+BJXp/+M+N/2kXemCgealOZXi00hj+rKfBFJGSs33uMV6at57s1+yiaN4b+1/pO6IqJ0riRrKDA\nF5GQs2TbIV6euo5FWw5StlBuHutQnZuuLkNkhHldWramYZkiEnIaVSjEmN5N+fhPTSiQO5on/ruc\nzm/M5usVu0lODp2dz5xKgS8iQWVmtKlejK8eacnQuxtiwCOjfqLrW3N01m6AKfBFxBMREcZ19Uox\n9dHWvHHHVZw+d57eI5fQbcg8Zqzfr+APAPXhi0hISDqfzJc/7eLN6RvYdfgUjSoU4olO1WlepajX\npYU8HbQVkWzpbFIy4xJ2MPiHjew9eppmlYvwRKfqxFcs7HVpIUuBLyLZ2ulz5xm1aDvvzNzEgeNn\naFO9GI93rE6DcgW9Li3kKPBFJEc4eTaJkQu28e6sTRw6eY4OtUrweMfq1C6d3+vSQoYCX0RylGOn\nz/HRvK0Mm7OZY6eTuL5eKR7tUI1qJfJ5XZrngjIO38xuM7PVZpZsZr9rzMzKm9lxM3syM+2IiOSL\njaZ/+2rMHXQt/a+tysz1++n0xmz6jVrK+r3HvC4vW8jssMxVQHdgdhrzXwO+yWQbIiK/KhAXzROd\najD36Wvp27YKM9ftp/Mbs+n32VLW7T3qdXkhLVNXMHLOrQXfiRQXM7ObgC3Aicy0ISKSmkJ5Yniq\nc00eaFmZD+Zu4aP5W5m8cg/X1S3JgPbVqFVKffwXC9RNzPMCTwMvZmDZ3maWYGYJiYmJgShHRHKw\nQnlieLLwGgR7AAAJ2UlEQVRzDeY+3Y7+11ZlzoYDXPfmHB4euYQ1u7XHn1K6e/hmNh0omcqs55xz\nE9N42QvA686546nt/afknBsGDAPfQdv06hERSU3BuBie6FSD+1tWYsTcLXw4bytTV++lc50SDGhf\njTqlC3hdoufSDXznXIcreN9rgFvN7GWgIJBsZqedc4Ov4L1ERDKsYFwMj3eqwf0tKzNi3hZGzNvC\nt6v30am2L/jrlgnf4A/IXQicc60uPDazF4DjCnsRCaYCcdE81rE6f2pZiQ/nbeGDuVuYtmYfHWqV\n4NEO4Rn8mR2WebOZ7QSaAZPN7NusKUtEJGsUyB3Nox2qM/fpa3msQ3UWb/mFG96eywMf/8iyHYe9\nLi+odOKViISVo/4TuD6Yu4Ujp87RqlpRHmlXlWsqF/G6tCumM21FRC7h+JkkPl24jffnbObA8bM0\nrliIR66tRutqRVMdah7KFPgiIhlw6ux5xv64nfdmb2bPkdPUK1OAR66tSsdaJYjIJrdeVOCLiFyG\ns0nJfLl0J0NnbWLbLyepXiIv/dpV5Yb6pUP+nrsKfBGRK5B0PpmvV+xhyIyNbNh/nIpF4ujTtgo3\nX12WmKjQvEmgAl9EJBOSkx3T1uxl8IyNrNp1lNIFYnm4bRVujy9HbHSk1+X9hgJfRCQLOOeY+XMi\nQ37YSMK2QxTNm4verStx9zUVyJMrIKcyXTYFvohIFnLOsWjLQQb/sJG5Gw9QMC6ans0rcl+zihTK\nE+NpbQp8EZEA+Wn7IYbM2Mj0tfvJHR3JnU3K8UCrypQpmNuTehT4IiIB9vO+Y7w7axOTlu0GoNtV\nZXi4TeWg34VLgS8iEiS7Dp/i/TmbGbN4B6fOnadDrRL0aVuZRhUKB6V9Bb6ISJAdPHGWj+dv5eMF\nWzl88hxNKhbm4baVaVejeEDP3lXgi4h45OTZJMYs3sH7czaz+8hpapbMx8NtqnBD/VJERWb9WH4F\nvoiIx86dT2bSst28O2sTG/Yfp2yh3DzYqjK3x5cjd0zWjeVX4IuIhIjkZMf36/YzdOZGlm4/TOE8\nMfRqXpF7mlWgYFzmh3Qq8EVEQoxzjh+3HuLdWZv4Yd1+4mIiubNxee5vVSlTQzoV+CIiIWztnqO8\nN2sTX63YA0Cv5hX5yw21r+i9Mhr4mb3j1W1mttrMks0s/qJ59c1sgX/+SjOLzUxbIiI5Sa1S+Xnj\nzquZPagdPZtXpFzhuIC3mdkLQawCugPvpZxoZlHAp8A9zrnlZlYEOJfJtkREcpwyBXPz1yvcs79c\nmQp859xaILXxpZ2AFc655f7lfslMOyIiknmBurhzdcCZ2bdmttTMBqW1oJn1NrMEM0tITEwMUDki\nIpLuHr6ZTQdKpjLrOefcxEu8b0ugMXAS+N5/UOH7ixd0zg0DhoHvoG1GCxcRkcuTbuA75zpcwfvu\nBGY75w4AmNkUoCHwu8AXEZHgCFSXzrdAPTOL8x/AbQOsCVBbIiKSAZkdlnmzme0EmgGTzexbAOfc\nIeA14EdgGbDUOTc5s8WKiMiVy+wonfHA+DTmfYpvaKaIiISA0LwFu4iIZLmQurSCmSUC27Lo7YoC\nB7LovbJSKNYVijWB6rocoVgThGZdoVgTZK6uCs65YuktFFKBn5XMLCEj15YItlCsKxRrAtV1OUKx\nJgjNukKxJghOXerSEREJEwp8EZEwkZMDf5jXBaQhFOsKxZpAdV2OUKwJQrOuUKwJglBXju3DFxGR\n38rJe/giIpKCAl9EJExku8A3sxFmtt/MVqUx38zsLTPbaGYrzKxhinn3mdkG/999Qa7rbn89K81s\nvpk1SDFvq3/6MjPLsns8ZqCmtmZ2xN/uMjN7PsW8Lma23r8e/5xVNWWwrqdS1LTKzM6bWWH/vECt\nq3JmNsPM1vjv0jYwlWWCvm1lsK6gblsZrCno21YG6/Ji24o1s8Vmttxf14upLJPLzMb618kiM6uY\nYt4z/unrzaxzpopxzmWrP6A1vitvrkpjflfgG8CApsAi//TCwGb/v4X8jwsFsa7mF9oDrrtQl//5\nVqCoB+uqLfB1KtMjgU1AZSAGWA7UDlZdFy17I/BDENZVKaCh/3E+4OeLP7MX21YG6wrqtpXBmoK+\nbWWkLo+2LQPy+h9HA4uAphct0xd41//4TmCs/3Ft/zrKBVTyr7vIK60l2+3hO+dmAwcvsUg34BPn\nsxAoaGalgM7Ad865g853cbfvgC7Bqss5N9/fLsBCoGxWtX2lNV1CE2Cjc26zc+4sMAbfevWirruA\n0VnVdlqcc3ucc0v9j48Ba4EyFy0W9G0rI3UFe9vK4LpKS8C2rSuoK1jblnPOHfc/jfb/XTxaphvw\nsf/x50B7MzP/9DHOuTPOuS3ARnzr8Ipku8DPgDLAjhTPd/qnpTXdC/fj21O8wAHTzGyJmfUOci3N\n/D81vzGzOv5pIbGuzCwOX3B+kWJywNeV/+f01fj2xFLydNu6RF0pBXXbSqcmz7at9NZVsLctM4s0\ns2XAfnw7B2luW865JOAIUIQsXl+ZvYm5XCYza4fvf8qWKSa3dM7tMrPiwHdmts6/FxxoS/Fdg+O4\nmXUFJgDVgtBuRt0IzHPOpfw1ENB1ZWZ58YXAo865o1n1vpmVkbqCvW2lU5Nn21YG/xsGddtyzp0H\nrjKzgsB4M6vrnEv1GFYg5cQ9/F1AuRTPy/qnpTU9aMysPvA+0M2luLG7c26X/9/9+C43fcU/2S6H\nc+7ohZ+azrkpQLSZFSUE1pXfnVz0kzuQ68rMovEFxWfOuS9TWcSTbSsDdQV920qvJq+2rYysK7+g\nblsp2jgMzOD3XX6/rhfz3TSqAPALWb2+svoARTD+gIqkfSDyen57YG2xf3phYAu+g2qF/I8LB7Gu\n8vj635pfND0PkC/F4/lAlyDVVJL/P/muCbDdv96i8B14rMT/H1irE6x15Z9fAF8/f55grCv/5/4E\neOMSywR928pgXUHdtjJYU9C3rYzU5dG2VQwo6H+cG5gD3HDRMv347UHbcf7HdfjtQdvNZOKgbbbr\n0jGz0fhGABQ13922/obvIAjOuXeBKfhGU2zEdwP1Xv55B83sH/juwgXwd/fbn3OBrut5fH1y7/iO\nxZDkfFfGK4HvJx74/mcY5ZybGqSabgX6mFkScAq40/m2siQzewTfrSojgRHOudVZUVMG6wK4GZjm\nnDuR4qUBW1dAC+AeYKW/rxXgWXxh6uW2lZG6gr1tZaQmL7atjNQFwd+2SgEfm1kkvl6Vcc65r83s\n70CCc24S8AEw0sw24vsyutNf82ozG4fvFrFJQD/n6x66Irq0gohImMiJffgiIpIKBb6ISJhQ4IuI\nhAkFvohImFDgi4iECQW+iEiYUOCLiISJ/wMjQL0yAp4IKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d07bcd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b_plot = np.arange(0.1, 15.0, 0.001)\n",
    "y = [log_prob([bi,1.0],sample) for bi in b_plot]\n",
    "plot(b_plot,y)\n",
    "title(\"Slice through b\")\n",
    "figure()\n",
    "beta_plot = np.arange(0.1, 3.0, 0.001)\n",
    "y_beta = [log_prob([1.0,beta],sample) for beta in beta_plot]\n",
    "plot(beta_plot,y_beta)\n",
    "_=title(\"Slice through beta\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, now let's run emcee on our log-posterior.  We'll use 20 walkers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import emcee\n",
    "walkers=20\n",
    "sampler=emcee.EnsembleSampler(nwalkers=walkers, dim=2, lnpostfn=log_prob, args=[sample])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Need to make sure that none of our initial points has bad posteriors, so just sample around 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p0=emcee.utils.sample_ball([2.0,1.0], [0.01,0.01], size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the sampler and pull out the parameters of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nsample=20000\n",
    "b_samples=np.zeros((nsample,walkers))\n",
    "beta_samples=np.zeros((nsample,walkers))\n",
    "\n",
    "for i,(params,logp,rstate) in enumerate(sampler.sample(p0, iterations=nsample)):\n",
    "    b_samples[i] = params[:,0]\n",
    "    beta_samples[i] = params[:,1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, that worked.  So lets histogram our results and plot against the true posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEN1JREFUeJzt3X+MXWldx/H3h64V+SGgOxJsu7TRAjaILI4FJAECS9IF\n0pIApk0wkKw0JhQQiFoiaUhNDKCC/lENFRGCQFlXoqOMFgIYowHSWVgX2lqYlB+dAu7w2x+BUvn6\nx9zFy2Wmc2bmTO/06fuVTHp+PPfc755tP/eZ55zz3FQVkqR23WfcBUiS1pdBL0mNM+glqXEGvSQ1\nzqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcDeN64xtvvLG2b98+rreXpGvSnXfe+ZWqmljJa8YW9Nu3\nb2dmZmZcby9J16Qkn1/paxy6kaTGGfSS1DiDXpIaZ9BLUuMMeklqXKegT7Inybkks0kOL7L/piQf\nTvKJJHcneWb/pUqSVmPZoE+yCTgG3ArsAg4k2TXS7DXA7VV1M7Af+JO+C5UkrU6XHv1uYLaqzlfV\nJeAEsG+kTQE/Plh+EPDF/kqUJK1FlwemtgAXhtbngMePtHkt8P4kLwXuD9zSS3WSpDXr68nYA8Db\nquoPkzwReEeSR1fV94YbJTkIHAS46aab1vym2w+/7/vLn3vds9Z8PElqUZehm4vAtqH1rYNtw24D\nbgeoqo8A9wVuHD1QVR2vqsmqmpyYWNFUDZKkVeoS9KeAnUl2JNnMwsXWqZE2XwCeDpDk51gI+vk+\nC5Ukrc6yQV9Vl4FDwEngLAt315xOcjTJ3kGzVwEvTvJvwLuBF1VVrVfRkqTuOo3RV9U0MD2y7cjQ\n8hngSf2WJknqg0/GSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4/qavXLs\nnMlSkhZnj16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1rpnbK4d5q6Uk/b9OPfoke5KcSzKb5PAi\n+9+U5K7Bz6eTfKP/UiVJq7Fsjz7JJuAY8AxgDjiVZGrw9YEAVNUrhtq/FLh5HWqVJK1Clx79bmC2\nqs5X1SXgBLDvCu0PsPAF4ZKkDaBL0G8BLgytzw22/ZAkDwd2AB9aYv/BJDNJZubn51daqyRpFfq+\n62Y/cEdV/e9iO6vqeFVNVtXkxMREz28tSVpMl6C/CGwbWt862LaY/ThsI0kbSpegPwXsTLIjyWYW\nwnxqtFGSRwEPAT7Sb4mSpLVYNuir6jJwCDgJnAVur6rTSY4m2TvUdD9woqpqfUqVJK1Gpwemqmoa\nmB7ZdmRk/bX9lSVJ6otTIEhS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGNTkf/bDhuenB+eklXX/s\n0UtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqXKegT7Inybkks0kO\nL9HmV5KcSXI6ybv6LVOStFrLznWTZBNwDHgGMAecSjJVVWeG2uwEXg08qaq+nuSn1qtgSdLKdOnR\n7wZmq+p8VV0CTgD7Rtq8GDhWVV8HqKp7+i1TkrRaXYJ+C3BhaH1usG3YI4BHJPnXJB9NsmexAyU5\nmGQmycz8/PzqKpYkrUhfF2NvAHYCTwUOAH+W5MGjjarqeFVNVtXkxMRET28tSbqSLkF/Edg2tL51\nsG3YHDBVVd+tqs8Cn2Yh+CVJY9Yl6E8BO5PsSLIZ2A9MjbT5GxZ68yS5kYWhnPM91ilJWqVlg76q\nLgOHgJPAWeD2qjqd5GiSvYNmJ4GvJjkDfBj4zar66noVLUnqrtNXCVbVNDA9su3I0HIBrxz8SJI2\nkOa/M3bU8HfI+v2xkq4HToEgSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN\nM+glqXEGvSQ1zqCXpMYZ9JLUuOtu9sphzmQp6Xpgj16SGmfQS1LjOgV9kj1JziWZTXJ4kf0vSjKf\n5K7Bz6/1X6okaTWWHaNPsgk4BjwDmANOJZmqqjMjTd9TVYfWoUZJ0hp06dHvBmar6nxVXQJOAPvW\ntyxJUl+6BP0W4MLQ+txg26jnJrk7yR1Jti12oCQHk8wkmZmfn19FuZKklerrYuzfAdur6jHAB4C3\nL9aoqo5X1WRVTU5MTPT01pKkK+kS9BeB4R761sG276uqr1bVdwarbwF+sZ/yJElr1SXoTwE7k+xI\nshnYD0wNN0jysKHVvcDZ/kqUJK3FsnfdVNXlJIeAk8Am4K1VdTrJUWCmqqaAlyXZC1wGvga8aB1r\nliStQKcpEKpqGpge2XZkaPnVwKv7LU2S1Ifreq6bYc57I6lVToEgSY0z6CWpcQa9JDXOoJekxhn0\nktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOGevXIQzWUpqiT16SWqcQS9J\njesU9En2JDmXZDbJ4Su0e26SSjLZX4mSpLVYNuiTbAKOAbcCu4ADSXYt0u6BwMuBj/VdpCRp9br0\n6HcDs1V1vqouASeAfYu0+13g9cC3e6xPkrRGXYJ+C3BhaH1usO37kjwO2FZV7+MKkhxMMpNkZn5+\nfsXFSpJWbs0XY5PcB3gj8Krl2lbV8aqarKrJiYmJtb61JKmDLkF/Edg2tL51sO1eDwQeDfxTks8B\nTwCmvCArSRtDl6A/BexMsiPJZmA/MHXvzqr6ZlXdWFXbq2o78FFgb1XNrEvFkqQVWTboq+oycAg4\nCZwFbq+q00mOJtm73gVKktam0xQIVTUNTI9sO7JE26euvayNY3g6BHBKBEnXHp+MlaTGGfSS1DiD\nXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXO74xdIb9PVtK1xh69JDXOoJekxhn0ktQ4g16SGmfQ\nS1LjvOtmDbwDR9K1wB69JDWuU9An2ZPkXJLZJIcX2f/rST6Z5K4k/5JkV/+lSpJWY9mgT7IJOAbc\nCuwCDiwS5O+qqp+vqscCbwDe2HulkqRV6dKj3w3MVtX5qroEnAD2DTeoqm8Nrd4fqP5KlCStRZeL\nsVuAC0Prc8DjRxsleQnwSmAz8LReqpMkrVlvF2Or6lhV/Qzw28BrFmuT5GCSmSQz8/Pzfb21JOkK\nugT9RWDb0PrWwbalnACes9iOqjpeVZNVNTkxMdG9SknSqnUJ+lPAziQ7kmwG9gNTww2S7BxafRbw\nmf5KlCStxbJj9FV1Ockh4CSwCXhrVZ1OchSYqaop4FCSW4DvAl8HXrieRUuSuuv0ZGxVTQPTI9uO\nDC2/vOe6JEk98clYSWqcc930xHlvJG1U9uglqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS47y9ch14\nq6WkjcQevSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqct1euM2+1lDRu9uglqXGdgj7JniTnkswm\nObzI/lcmOZPk7iQfTPLw/kuVJK3GskGfZBNwDLgV2AUcSLJrpNkngMmqegxwB/CGvguVJK1Olx79\nbmC2qs5X1SXgBLBvuEFVfbiq/mew+lFga79lSpJWq0vQbwEuDK3PDbYt5TbgH9ZSlCSpP73edZPk\nBcAk8JQl9h8EDgLcdNNNfb61JGkJXYL+IrBtaH3rYNsPSHIL8DvAU6rqO4sdqKqOA8cBJicna8XV\nXuOGb7UEb7eUdHV0Gbo5BexMsiPJZmA/MDXcIMnNwJuBvVV1T/9lSpJWa9mgr6rLwCHgJHAWuL2q\nTic5mmTvoNnvAw8A/irJXUmmljicJOkq6zRGX1XTwPTItiNDy7f0XJckqSc+GStJjTPoJalxBr0k\nNc6gl6TGOU3xGDmFsaSrwR69JDXOoJekxhn0ktQ4x+g3CMfrJa0Xe/SS1DiDXpIa59DNBuQwjqQ+\n2aOXpMYZ9JLUOINekhrnGP0G53i9pLWyRy9JjesU9En2JDmXZDbJ4UX2PznJx5NcTvK8/ssULPTu\n7/2RpK6WDfokm4BjwK3ALuBAkl0jzb4AvAh4V98FSpLWpssY/W5gtqrOAyQ5AewDztzboKo+N9j3\nvXWoUZK0Bl2CfgtwYWh9Dnj8+pSjrkaHb7xQK2kpV/VibJKDSWaSzMzPz1/Nt5ak61aXHv1FYNvQ\n+tbBthWrquPAcYDJyclazTG0OG/DlLSULj36U8DOJDuSbAb2A1PrW5YkqS/LBn1VXQYOASeBs8Dt\nVXU6ydEkewGS/FKSOeD5wJuTnF7PoiVJ3XV6MraqpoHpkW1HhpZPsTCkow3AYRxJw3wyVpIa51w3\njVvqKVp7+tL1wx69JDXOHv116krz5djbl9pi0OuHeDFXaotBrysy9KVrn0Gvzgx96dpk0GtVHOOX\nrh0GvXrX5YtR/DCQrh6DXmNxpWEg7/2X+mXQa+y6fjXiar5CcfjDYa0fIF6jUB+6/j3q8++bQa+m\ndflwWM0HSJcPjbV+t++VjtXlH37X35qWep+1/rd0fX2f52ytVnou1vr/qKu1npdUjWda+MnJyZqZ\nmVnTMcb9l0KSrrbPv/7Zd1bV5Epe4xQIktQ4g16SGmfQS1LjDHpJalynoE+yJ8m5JLNJDi+y/0eT\nvGew/2NJtvddqCRpdZYN+iSbgGPArcAu4ECSXSPNbgO+XlU/C7wJeH3fhUqSVqdLj343MFtV56vq\nEnAC2DfSZh/w9sHyHcDTk6S/MiVJq9Ul6LcAF4bW5wbbFm1TVZeBbwI/2UeBkqS1uapPxiY5CBwc\nrP5XknNX8/2vghuBr4y7iA3I87I0z83iPC9Le+RKX9Al6C8C24bWtw62LdZmLskNwIOAr44eqKqO\nA8dXWuS1IsnMSp9Yux54XpbmuVmc52VpSVY8pUCXoZtTwM4kO5JsBvYDUyNtpoAXDpafB3yoxjW3\ngiTpByzbo6+qy0kOASeBTcBbq+p0kqPATFVNAX8OvCPJLPA1Fj4MJEkbQKcx+qqaBqZHth0ZWv42\n8Px+S7smNTsstUael6V5bhbneVnais/N2GavlCRdHU6BIEmNM+h7kGRbkg8nOZPkdJKXj7umjSTJ\npiSfSPL3465lo0jy4CR3JPn3JGeTPHHcNW0USV4x+Hf0qSTvTnLfcdc0LknemuSeJJ8a2vYTST6Q\n5DODPx+y3HEM+n5cBl5VVbuAJwAvWWSaiOvZy4Gz4y5ig/lj4B+r6lHAL+D5ASDJFuBlwGRVPZqF\nG0Cu55s73gbsGdl2GPhgVe0EPjhYvyKDvgdV9aWq+vhg+T9Z+Ec7+vTwdSnJVuBZwFvGXctGkeRB\nwJNZuFuNqrpUVd8Yb1Ubyg3Ajw2eybkf8MUx1zM2VfXPLNzJOGx4ypm3A89Z7jgGfc8GM3feDHxs\nvJVsGH8E/BbwvXEXsoHsAOaBvxgMab0lyf3HXdRGUFUXgT8AvgB8CfhmVb1/vFVtOA+tqi8Nlr8M\nPHS5Fxj0PUryAOCvgd+oqm+Nu55xS/Js4J6qunPctWwwNwCPA/60qm4G/psOv35fDwbjzftY+DD8\naeD+SV4w3qo2rsGDqcveOmnQ9yTJj7AQ8u+sqveOu54N4knA3iSfY2HW06cl+cvxlrQhzAFzVXXv\nb313sBD8gluAz1bVfFV9F3gv8Mtjrmmj+Y8kDwMY/HnPci8w6HswmJL5z4GzVfXGcdezUVTVq6tq\na1VtZ+GC2oeq6rrvnVXVl4ELSe6dnOrpwJkxlrSRfAF4QpL7Df5dPR0vVI8annLmhcDfLvcCg74f\nTwJ+lYUe612Dn2eOuyhtaC8F3pnkbuCxwO+NuZ4NYfBbzh3Ax4FPspBR1+1TskneDXwEeGSSuSS3\nAa8DnpHkMyz8BvS6ZY/jk7GS1DZ79JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG\n/R/l7Zvt6jNHOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1084d13d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSBJREFUeJzt3X+s3fV93/Hnaw4kSoISU9+0zNgYVKb8aBOgV5AmaDFT\nQhy64FTtVLO0NRmR1yhk7SpVIotEIqJqaSutWwQbsVKLUrUmGwmdu5oSt5CylTr1debwKyEYhxVb\naHYwJWFkYabv/XG+l3y53ON77r3nnnN9v8+HdHS/38/n8z3nfb7++v39nM/3V6oKSVJ3/INxByBJ\nGi0TvyR1jIlfkjrGxC9JHWPil6SOMfFLUseY+CWpY+ZM/EnWJbknycNJHkryq7O0SZLPJjmY5P4k\nF7XqtiZ5tHltHfYXkCTNT+a6gCvJWcBZVfW1JGcA+4EPVNXDrTZXAB8DrgAuAf5DVV2S5ExgCpgE\nqln2p6rq6SX5NpKkOb1irgZV9STwZDP9vSTfANYCD7eabQZurd5eZG+S1zc7jI3Anqo6DpBkD7AJ\n2Hmyz1yzZk1t2LBh/t9Gkjpq//7936mqiUHazpn425JsAC4Evjqjai3wRGv+cFPWr/ykNmzYwNTU\n1HxCk6ROS/K/Bm078MHdJK8Fvgj8WlV9dyGBzfH+25JMJZk6duzYsN9ektQYKPEnOY1e0v/DqvrS\nLE2OAOta82c3Zf3KX6aqtlfVZFVNTkwM9GtFkrQAg5zVE+D3gG9U1b/r02wX8MvN2T1vB55pjg3c\nBVyeZHWS1cDlTZkkaUwGGeN/J/BLwANJDjRl/wZYD1BVNwO76Z3RcxB4DvhQU3c8yaeBfc1yN0wf\n6JUkjccgZ/X8DyBztCngo33qdgA7FhSdJGnovHJXkjrGxC9JHWPil6SOMfFLUsfM68pdLb0N1/3p\nS+Yf/8zPzFrXLpek+bDHL0kdY+KXpI4x8UtSx5j4JaljTPyS1DEmfknqGBO/JHWM5/EvAzPP3Zek\npWSPX5I6xh7/Mtfv14BX8UpaKHv8ktQx9vhXAHv/kubDHr8kdcycPf4kO4B/Chytqp+Ypf43gA+2\n3u9NwETzvN3Hge8BLwAnqmpyWIFLkhZmkB7/LcCmfpVV9TtVdUFVXQB8HPjLGQ9Uv6ypN+lL0jIw\nZ+KvqnuB43O1a1wF7FxURJKkJTW0Mf4kr6b3y+CLreICvpxkf5Jtw/osSdLCDfOsnvcDfzVjmOfS\nqjqS5A3AniTfbH5BvEyzY9gGsH79+iGGtTx5ta6kcRnmWT1bmDHMU1VHmr9HgTuAi/stXFXbq2qy\nqiYnJiaGGJYkqW0oiT/J64B3Af+1VfaaJGdMTwOXAw8O4/MkSQs3yOmcO4GNwJokh4FPAqcBVNXN\nTbOfBb5cVf+nteiPAnckmf6cP6qqPxte6JKkhZgz8VfVVQO0uYXeaZ/tskPA2xYamCRpaXjlriR1\njIlfkjrGm7SN0ChO4fSGbZLmYo9fkjrGxC9JHeNQzwrmsI+k2Zj4O8KdgKRpDvVIUseY+CWpY0z8\nktQxJn5J6hgTvyR1jIlfkjrGxC9JHWPil6SO8QKuJbYcn63rxVxSt9njl6SOMfFLUsfMmfiT7Ehy\nNMmsD0pPsjHJM0kONK/rW3WbkjyS5GCS64YZuCRpYQbp8d8CbJqjzX+vqgua1w0ASVYBNwHvA94M\nXJXkzYsJVpK0eHMm/qq6Fzi+gPe+GDhYVYeq6nngNmDzAt5HkjREwxrj/+kkX09yZ5K3NGVrgSda\nbQ43ZZKkMRrG6ZxfA86pqmeTXAH8MXD+fN8kyTZgG8D69euHEJYkaTaL7vFX1Xer6tlmejdwWpI1\nwBFgXavp2U1Zv/fZXlWTVTU5MTGx2LAkSX0sOvEn+bEkaaYvbt7zKWAfcH6Sc5OcDmwBdi328yRJ\nizPnUE+SncBGYE2Sw8AngdMAqupm4OeBjyQ5AXwf2FJVBZxIci1wF7AK2FFVDy3Jt9CCeRWv1D1z\nJv6qumqO+huBG/vU7QZ2Lyy0U9dyvE2DJE3zyl1J6hgTvyR1jHfn1Isc75e6wR6/JHWMiV+SOsbE\nL0kdY+KXpI7x4K5m5YFeaeWyxy9JHWPil6SOcahHc3LYR1pZTPyaF3cC0qnPoR5J6hgTvyR1jEM9\nWrB+t592CEha3uzxS1LHmPglqWNM/JLUMXMm/iQ7khxN8mCf+g8muT/JA0nuS/K2Vt3jTfmBJFPD\nDFyStDCDHNy9hd4zdW/tU/9t4F1V9XSS9wHbgUta9ZdV1XcWFeUpwOfsSjpVDPKw9XuTbDhJ/X2t\n2b3A2YsPS5K0VIY9xn8NcGdrvoAvJ9mfZNuQP0uStABDO48/yWX0Ev+lreJLq+pIkjcAe5J8s6ru\n7bP8NmAbwPr164cVliRphqH0+JO8Ffg8sLmqnpour6ojzd+jwB3Axf3eo6q2V9VkVU1OTEwMIyxJ\n0iwWnfiTrAe+BPxSVX2rVf6aJGdMTwOXA7OeGSRJGp05h3qS7AQ2AmuSHAY+CZwGUFU3A9cDPwL8\nxyQAJ6pqEvhR4I6m7BXAH1XVny3Bd9Ay4x08peVtkLN6rpqj/sPAh2cpPwS87eVLSJLGySt3Jalj\nvDunlpR38JSWH3v8ktQxJn5J6hiHejQWnvkjjY89fknqGBO/JHWMiV+SOsYxfo2d4/3SaNnjl6SO\nMfFLUseY+CWpY0z8ktQxJn5J6hgTvyR1jIlfkjrGxC9JHWPil6SOGSjxJ9mR5GiSWR+Wnp7PJjmY\n5P4kF7XqtiZ5tHltHVbg6rYN1/3piy9J8zPoLRtuAW4Ebu1T/z7g/OZ1CfCfgEuSnEnv4eyTQAH7\nk+yqqqcXE/RyYdJZWt7KQVoaAyX+qro3yYaTNNkM3FpVBexN8vokZwEbgT1VdRwgyR5gE7BzMUFr\n5XJnKi29Yd2kbS3wRGv+cFPWr1yaF3cI0vAsm4O7SbYlmUoydezYsXGHI0kr1rAS/xFgXWv+7Kas\nX/nLVNX2qpqsqsmJiYkhhSVJmmlYiX8X8MvN2T1vB56pqieBu4DLk6xOshq4vCmTJI3JQGP8SXbS\nO1C7JslhemfqnAZQVTcDu4ErgIPAc8CHmrrjST4N7Gve6obpA72SpPEY9Kyeq+aoL+Cjfep2ADvm\nH5okaSksm4O7kqTRMPFLUseY+CWpY4Z1AZc0Nv0u7vI2D9Ls7PFLUseY+CWpY0z8ktQxJn5J6hgT\nvyR1jIlfkjrGxC9JHWPil6SOMfFLUseY+CWpY0z8ktQxJn5J6hgTvyR1jHfn1IrVvmund+qUfmig\nHn+STUkeSXIwyXWz1P9ukgPN61tJ/q5V90Krbtcwg5ckzd+cPf4kq4CbgPcAh4F9SXZV1cPTbarq\nX7fafwy4sPUW36+qC4YX8nj1u/e7JJ0qBunxXwwcrKpDVfU8cBuw+STtrwJ2DiM4SdLwDZL41wJP\ntOYPN2Uvk+Qc4Fzg7lbxq5JMJdmb5AMLjlSSNBTDPri7Bbi9ql5olZ1TVUeSnAfcneSBqnps5oJJ\ntgHbANavXz/ksCRJ0wbp8R8B1rXmz27KZrOFGcM8VXWk+XsI+AovHf9vt9teVZNVNTkxMTFAWJKk\nhRgk8e8Dzk9ybpLT6SX3l52dk+SNwGrgr1tlq5O8spleA7wTeHjmspKk0ZlzqKeqTiS5FrgLWAXs\nqKqHktwATFXV9E5gC3BbVVVr8TcBn0vy9/R2Mp9pnw0kSRq9gcb4q2o3sHtG2fUz5j81y3L3AT+5\niPgkSUPmLRskqWNM/JLUMSZ+SeoYE78kdYyJX5I6xsQvSR1j4pekjjHxS1LH+AQudYJP45J+yB6/\nJHWMiV+SOsbEL0kdY+KXpI7x4O4AfMC6pJXEHr8kdYyJX5I6xsQvSR1j4pekjhko8SfZlOSRJAeT\nXDdL/dVJjiU50Lw+3KrbmuTR5rV1mMFLkuZvzrN6kqwCbgLeAxwG9iXZNctD079QVdfOWPZM4JPA\nJFDA/mbZp4cSvSRp3gbp8V8MHKyqQ1X1PHAbsHnA938vsKeqjjfJfg+waWGhSpKGYZDEvxZ4ojV/\nuCmb6eeS3J/k9iTr5rmsJGlEhnVw90+ADVX1Vnq9+t+f7xsk2ZZkKsnUsWPHhhSWJGmmQa7cPQKs\na82f3ZS9qKqeas1+Hvjt1rIbZyz7ldk+pKq2A9sBJicna4C4pAXxFs3qukF6/PuA85Ocm+R0YAuw\nq90gyVmt2SuBbzTTdwGXJ1mdZDVweVMmSRqTOXv8VXUiybX0EvYqYEdVPZTkBmCqqnYB/yrJlcAJ\n4DhwdbPs8SSfprfzALihqo4vwfeQJA1ooJu0VdVuYPeMsutb0x8HPt5n2R3AjkXEKEkaIq/claSO\nMfFLUseY+CWpY0z8ktQxPoFLneY5/eoie/yS1DH2+PvwObuSVip7/JLUMSZ+SeoYE78kdYyJX5I6\nxsQvSR1j4pekjvF0TqnhxVzqCnv8ktQxJn5J6hgTvyR1jGP80iwc79dKNlCPP8mmJI8kOZjkulnq\nfz3Jw0nuT/IXSc5p1b2Q5EDz2jVzWUnSaM3Z40+yCrgJeA9wGNiXZFdVPdxq9j+Byap6LslHgN8G\nfqGp+35VXTDkuCVJCzRIj/9i4GBVHaqq54HbgM3tBlV1T1U918zuBc4ebpiSpGEZJPGvBZ5ozR9u\nyvq5BrizNf+qJFNJ9ib5wAJilCQN0VAP7ib5RWASeFer+JyqOpLkPODuJA9U1WOzLLsN2Aawfv36\nYYY1MO/BL6kLBunxHwHWtebPbspeIsm7gU8AV1bVD6bLq+pI8/cQ8BXgwtk+pKq2V9VkVU1OTEwM\n/AUkSfMzSI9/H3B+knPpJfwtwD9vN0hyIfA5YFNVHW2Vrwaeq6ofJFkDvJPegV/plOGpnVpp5kz8\nVXUiybXAXcAqYEdVPZTkBmCqqnYBvwO8FvgvSQD+tqquBN4EfC7J39P7dfGZGWcDSZJGbKAx/qra\nDeyeUXZ9a/rdfZa7D/jJxQQoSRoub9kgSR3jLRukeXC8XyuBPX5J6hgTvyR1jEM90gI57KNTlT1+\nSeoYE78kdYyJX5I6xsQvSR3T+YO73pFTUtfY45ekjul8j18ahn6/HD3NU8uRPX5J6hh7/NIS8peA\nliN7/JLUMfb4pTHwdg8ap04mfk/h1HLlDkGj0MnELy0ndkQ0agON8SfZlOSRJAeTXDdL/SuTfKGp\n/2qSDa26jzfljyR57/BClyQtxJw9/iSrgJuA9wCHgX1Jds14aPo1wNNV9eNJtgC/BfxCkjcDW4C3\nAP8Q+PMk/6iqXhj2F5FWGs8I0lIZZKjnYuBgVR0CSHIbsBloJ/7NwKea6duBG5OkKb+tqn4AfDvJ\nweb9/no44Q/On9NaKRayLffbWXhMoZsGSfxrgSda84eBS/q1qaoTSZ4BfqQp3ztj2bULjnaeTPZS\nzyD/F4b1/6W9A1nMjmW57JSWSxzDtGwO7ibZBmxrZp9N8siQP2IN8J0hv+dSOBXiNMbhORXinFeM\n+a35lS/mPVtGsh4X8x0aSxnnOYM2HCTxHwHWtebPbspma3M4ySuA1wFPDbgsAFW1Hdg+WNjzl2Sq\nqiaX6v2H5VSI0xiH51SI0xiHZ7nEOchZPfuA85Ocm+R0egdrd81oswvY2kz/PHB3VVVTvqU56+dc\n4Hzgb4YTuiRpIebs8Tdj9tcCdwGrgB1V9VCSG4CpqtoF/B7wB83B2+P0dg407f4zvQPBJ4CPekaP\nJI3XQGP8VbUb2D2j7PrW9P8F/lmfZX8T+M1FxDgsSzaMNGSnQpzGODynQpzGODzLIs70RmQkSV3h\n3TklqWNWROIf4JYSv57k4ST3J/mLJOe06l5IcqB5zTxoPcoYr05yrBXLh1t1W5M82ry2zlx2xHH+\nbivGbyX5u1bdqNbljiRHkzzYpz5JPtt8h/uTXNSqG8m6HCDGDzaxPZDkviRva9U93pQfSDI1xhg3\nJnmm9W96favupNvJCGP8jVZ8Dzbb4JlN3UjWY/NZ65Lc0+SZh5L86ixtxr5dvqiqTukXvQPOjwHn\nAacDXwfePKPNZcCrm+mPAF9o1T27TGK8GrhxlmXPBA41f1c306vHFeeM9h+jd7B/ZOuy+Zx/DFwE\nPNin/grgTiDA24GvjmFdzhXjO6Y/G3jfdIzN/OPAmmWwHjcC/22x28lSxjij7fvpnVE40vXYfNZZ\nwEXN9BnAt2b5Pz727XL6tRJ6/C/eUqKqngembynxoqq6p6qea2b30rueYFnFeBLvBfZU1fGqehrY\nA2xaJnFeBexcolj6qqp76Z091s9m4Nbq2Qu8PslZjHBdzhVjVd3XxADj2SYHWY/9LGZ7npd5xjiW\n7RGgqp6sqq81098DvsHL71Iw9u1y2kpI/LPdUuJkt4W4ht5ed9qrkkwl2ZvkA0sRIIPH+HPNT8Db\nk0xf+Dbf77cYA39WM1x2LnB3q3gU63IQ/b7HKNflfMzcJgv4cpL96V3RPk4/neTrSe5M8pambNmt\nxySvppcsv9gqHst6TO/uxBcCX51RtWy2y2Vzy4ZRSPKLwCTwrlbxOVV1JMl5wN1JHqiqx8YQ3p8A\nO6vqB0n+JfD7wD8ZQxyD2gLcXi+9LmO5rMtTRpLL6CX+S1vFlzbr8Q3AniTfbHq+o/Y1ev+mzya5\nAvhjehdhLkfvB/6qqtq/Dka+HpO8lt7O59eq6rtL+VmLsRJ6/APdFiLJu4FPAFdW726hAFTVkebv\nIeAr9PbUI4+xqp5qxfV54KcGXXaUcbZsYcbP6hGty0H0+x6jXJdzSvJWev/Wm6vqqeny1no8CtxB\nb2hl5Krqu1X1bDO9GzgtyRqW2XpsnGx7HMl6THIavaT/h1X1pVmaLJ/tchQHPpbyRe9XyyF6ww7T\nB5reMqPNhfQORp0/o3w18Mpmeg3wKEtwkGrAGM9qTf8ssLd+eODn202sq5vpM8e1Lpt2b6R34Cyj\nXpetz9tA/4OSP8NLD6L9zajX5QAxrgcOAu+YUf4a4IzW9H3ApjHF+GPT/8b0kubfNut0oO1kFDE2\n9a+jdxzgNWNcjwFuBf79Sdosi+2yqk79xN+suCvoHUV/DPhEU3YDvd49wJ8D/xs40Lx2NeXvAB5o\nNtwHgGvGGOO/BR5qYrkHeGNr2X/RJImDwIfGuS6b+U8Bn5mx3CjX5U7gSeD/0RsPvQb4FeBXmvrQ\ne3jQY00sk6NelwPE+Hng6dY2OdWUn9esw68328Mnxhjjta1tci+tndRs28k4YmzaXE3vuR/t5Ua2\nHpvPu5TeMYX7W/+mVyy37XL65ZW7ktQxK2GMX5I0DyZ+SeoYE78kdYyJX5I6xsQvSR1j4pekjjHx\nS1LHmPglqWP+P3fg10dGSl1gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108a9e7d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist(b_samples[200:].flatten(), bins=100, normed=True)\n",
    "xlim(0.5, 10)\n",
    "figure()\n",
    "_=hist(beta_samples[200:].flatten(), bins=100, normed=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a log-scale version of the same thing.  You can see the flattening in this one, which is why we need the upper limit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can easily free the beta and lower limit parameters if you need to - you might need to check there are reasonable limits on them but it should behave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10adbd210>"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAo8AAAHjCAYAAABCTvXHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UpnlZH/jvVdVv0zMMCI2IMyCYRY5koqIjaPSsIBoH\nw8LxiCsYdWHR2bzgy2pMNFkhIefsWc3G1Wx8SS8S8CUgDuhOWAiii4dVgTCgwWEQdxYVZkSHeWHe\np7ur6to/qoa0TVf374Z65rmf6s/nnOd0PVXX+d2/eu56qq+67vt3/aq7AwAAI9aWPQEAAFaH5BEA\ngGGSRwAAhkkeAQAYJnkEAGCY5BEAgGGSRwAAhkkeAQAYJnkEAGDYgWVPYKpDdbiP5OJlTwMA2Efu\nzh23dvejlzmHb3jmxX3b7Zt7Pu5733/ird191V6Nt3LJ45FcnKfXs5Y9DUZUjcdO2SZzdFxbbwIw\n6Df7mj9b9hxuu30z/+mtj9/zcdcf+/8e28vxVi55BADYjzrJVraWPY3zcs8jAADDVB4BAGahs9kq\njwAA7CMqjwAAM7B9z+P8F3tKHgEAZsKCGQAA9hWVR6ZZW1/MuL33TVEBYJV0Opsr0KNY5REAgGEq\njwAAM2HBDAAAQzrJ5gokjy5bAwAwTOURAGAmVuGytcojAADDVB4BAGagk5Vo1SN5ZJopG7ZPeQNM\n6R+5pSckAPvT/PeXcdkaAIAJVB4BAGag01r1AACwv6g8AgDMQSeb8y88qjwCADBO5REAYAY6q7Ha\nWvIIADALlc3UsidxXi5bAwAwTOWRaaY0/q4Jfz1NaT4+Zdzh46/AHcoA7GudZGsF/jtSeQQAYJjK\nIwDATKzCPY+SRwCAGeisRvK4sMvWVfWqqrqlqq4/T9yXV9VGVT1/UXMBAGBvLPKex1cnuepcAVW1\nnuTHkvzGAucBALAStrr2/LHXFpY8dvc7ktx+nrDvSfKGJLcsah4AAOydpd3zWFWXJfmmJM9M8uXn\nib06ydVJciRHFz85AICH2Krc87jMBTM/meQfd/dWnadvX3cfT3I8SS6tR65AByQAgGk6lc0V6KK4\nzOTxyiSv20kcjyX5xqra6O5fX+KcOJ+19YUMW2vjf2n15uZC5gAAnN/SksfufuKDH1fVq5O8SeII\nAFzIFrHAZa8tLHmsqtcmeUaSY1V1U5KXJzmYJN39c4s6LgAAi7Ow5LG7Xzgh9kWLmgcAwCqwYAYA\ngAkqmz3/BTPznyEAALOh8ggAMAOdZGsF6nrznyEAALOh8ggAMBMWzLD/9NZw6Nrhw+PDbmx8OrMB\nAPZAVf2PSb4r21fP/zDJi7v7gbPFumwNADAD3durrff6cT5VdVmS701yZXdfkWQ9yQt2i1d5BACY\nia3lXbY+kOSiqjqV5GiSP98tUOURAGB/O1ZV1532uPr0L3b3zUn+1yQfSfKxJHd292/sNpjKIwDA\nDGzvMLOQut6t3X3lbl+sqs9K8rwkT0zyiSS/WlXf3t2/dLZ4lUcAgAvb1yX5k+7+eHefSvLGJH9z\nt2CVRwCAWVja9oQfSfIVVXU0yf1JnpXkut2CJY8AADOwrB1muvvdVXVNkvcl2Ujy+0mO7xYveWSS\nOnRoOLa7xwdeXx+P3dwcjwUAzqu7X57k5SOxkkcAgJnY7PnvMGPBDAAAw1QeAQBmoFOLatWzpySP\nAAAzsbWc1daTzH+GAADMhsojAMAMLHCHmT01/xkCADAbKo8AADPQqZVo1SN5ZJoJDbrr8OHh2D5x\nYnwONVgw763xMQGAIZJHAICZWMb2hFNJHgEAZqA72dSqBwCA/UTlEQBgFipbmf+CGZVHAACGqTwC\nAMxAZzXueZQ8AgDMhB1mAADYV1QemaQOHVpIbLrHY0+eHBtyvJ95UhNuUB6d6yLGBGDf6lS2VmCH\nGZVHAACGqTwCAMzEKtzzKHkEAJiBTrK1Aqut5z9DAABmQ+URAGAWKpt2mAEAYD9ReQQAmIFVuedR\n8sg06+vjsb01HjulJ+LoHLbGeyfW2vjxe3OwgeSU3o2L6gk5Oq4+kwAMkjwCAMzEKtzzKHkEAJiB\n7lqJy9bznyEAALOh8ggAMBObKo8AAOwnKo8AADPQSbYsmAEAYEy5bA0AwP6i8sgkdejghODx0nvV\nhL9jRhtajzbzTtJTGooPNimfMuakhupTaP4NsDK2d5iZ/2VrlUcAAIapPAIAzMTmCtT1JI8AADPQ\nKZetAQDYX1QeAQBmYmsJdb2qenKSXzntU5+f5GXd/ZNni5c8AgBcwLr7Q0m+JEmqaj3JzUl+bbd4\nySMAwAx0J5vLv+fxWUn+v+7+s90CJI8AAPvbsaq67rTnx7v7+C6xL0jy2nMNJnlkkjo43iS8N6c0\nvl5Qk+wF6NHm41Man0+xNtakPEmyNd4oHYDlW9Bq61u7+8rzBVXVoSTPTfIj54pb2F2ZVfWqqrql\nqq7f5et/p6reX1V/WFW/V1VfvKi5AADM3XarnrU9f0zw7CTv6+6/PFfQIpf0vDrJVef4+p8k+Zru\n/htJ/kWS3cqnAAAs3gtznkvWyQIvW3f3O6rqCef4+u+d9vRdSS5f1FwAAFbBZpazYKaqLk7y9Un+\nh/PFzuWex5ckectuX6yqq5NcnSRHcvShmhMAwAWhu+9N8qiR2KUnj1X1zGwnj1+9W8zOiqDjSXJp\nPbIfoqkBADxkOgtbMLOnlpo8VtUXJXllkmd3923LnAsAwHLV1AUuS7G0GVbV45O8Mcl3dPcfL2se\nAACMW1jlsapem+QZ2W5MeVOSlyc5mCTd/XNJXpbta+s/U1VJsjHSg4glWx/vMViHD42Pe+LkeOxg\nn8XhfoxJakLrxD61MR48PIEJf8dN6d1Yg5c/2t0gAHOwtaQFM1MscrX1C8/z9e9K8l2LOj4AAHtv\n6QtmAACYzd7W5yV5BACYCQtmAADYV1QeAQBmYHtv6/lftlZ5BABgmMojAMBMrEKrHpVHAACGqTwy\nydaxh48HjzaoTrL2l7ePj3tw7Md2yl9GfXJCk/K1se+rtxbUeHvC66r5N8DqsLc1AACTaNUDAMC+\novIIADAHrVUPAAD7jMojAMAMdFajVY/kEQBgJly2BgBgX1F5BACYAX0e2Zf6wHix+oHHHB2OPXrP\n/eOT2Ngcizs4PmQ2NsZj19eHwiqD88y0huI1ePwk6c3BOUxpJq5JOcAFTfIIADATKo8AAAzp6PMI\nAMA+o/IIADATq9DnUeURAIBhKo8AAHPQq7FgRuURAIBhKo8AADOgSTj70slHXTQcu3Vo/A1w8vLP\nGo49PNgkvO+8e3jMOrD3b4Up7bGnNBRfiCmNvwFYmFVIHl22BgBgmMojAMAMaBIOAMC+o/IIADAT\nvQKVR8kjAMBM2GEGAIB9ReURAGAGekV2mJE8MskDjxz/kbnv0eOF7c++5cRw7NbDLx6KW3tgfMze\n3BqOzcmT47HLVgu4uNATXqsp/SN7SmfMFXGhf//AviR5BACYCQtmAAAYtLw+j1X1iCSvTHJFtjdJ\n+++7+51ni5U8AgDwU0n+Y3c/v6oOJTm6W6DkEQBgJpZx2bqqHp7kv07you059Mkku97gr1UPAMD+\ndqyqrjvtcfUZX39iko8n+XdV9ftV9cqq2nV1qsojAMAMdBbWqufW7r7yHF8/kORLk3xPd7+7qn4q\nyQ8n+dGzBas8AgBc2G5KclN3v3vn+TXZTibPSvIIADAHvd3yda8f5z1s918k+WhVPXnnU89KcsNu\n8S5bM8kDjxgvp9972XjT43v/4shw7NFbxv7mOXjvrgvFPtXd9w6H1pHDY4GnNobH7K0JDaLXxs9B\nZXPvjz+l8fiiGooPH3/JjbeXfXxg5Sxxb+vvSfLLOyutP5zkxbsFSh4BAC5w3f0HSc51X+QnSR4B\nAGagsxo7zLjnEQCAYSqPAACzsLztCaeQPAIAzMQqrLNz2RoAgGEqjwAAM2HBDAAA+4rKI5N84gvH\nb8ZYf/QDw7F3TGjoffDesR/bA6PNvJOsnTg1HJsDY8fvu+4eH3NC4+9J1teHwkabiSdJb47HTmoo\nvjVh3EUYbVI+5YakKY3PV+FGJ2ChtneEmX/lUfIIADATq7Da2mVrAACGqTwCAMzEKtzBovIIAMAw\nlUcAgJmwYAYAgCGdWonkcWGXravqVVV1S1Vdv8vXq6r+dVXdWFXvr6ovXdRcAADYG4usPL46yb9J\n8gu7fP3ZSZ6083h6kp/d+ZcZO3LL+N8b/+gb3jwc+4qN5wzH3vOxsf6NaxuXDI955OR4n8cajK0j\nR8bHnNAPcEqfxd7YGAsc7AeZJFP+Ju6t8Tu/a7R/5ir1mVyFO9+BWVmF3xoLqzx29zuS3H6OkOcl\n+YXe9q4kj6iqxy5qPgAAfOaWec/jZUk+etrzm3Y+97EzA6vq6iRXJ8mRjO9EAgCwMlZkh5mVaNXT\n3ce7+8ruvvJgxrecAwBgby2z8nhzksed9vzync8BAFyYVuCmx2VWHq9N8p07q66/Ismd3f0pl6wB\nAC4U3bXnj722sMpjVb02yTOSHKuqm5K8PMnBJOnun0vy5iTfmOTGJPclefGi5gIAwN5YWPLY3S88\nz9c7yT9Y1PEBAFbNKnT4WokFMwAAzIPtCZnk8T9+3Xjwd4yHPvazPzEce9vjHzMUt3bq4PCYW+uP\nHI49euNtQ3F1+NDwmD2hQXWtTbh/ZUpD7UE9qaH4hIbmow3FpzT+7q3x2NFG7YsqC0xoFL8SpQlg\nss5qtOqRPAIAzEEnWYHk0WVrAACGqTwCAMzEKtyVovIIAMAwlUcAgLlYgcqj5BEAYBYWsyPMXnPZ\nGgCAYSqPAABz4bI1+02fOjkc+xP/x/OHY9//gz8zHPvkO79zKO7EbZcMj3nxXw6HZuOzLx2KO3DL\nXcNjTrlI0adOjQdvDjbJnjLmaDPvZFJD7xrsPd4b43OtCQ3Ne7Sh+tr4mJNMaWgOsESSRwCAOejV\n2GHGPY8AAAxTeQQAmAv3PAIAMM5lawAA9hGVRwCAuVjSZeuq+tMkdyfZTLLR3VfuFit5BAAgSZ7Z\n3beeL0jyyMI89l/93nDsS77lq4djv/BzxpoyXv+4o8NjHrh//K1w6Z+N/Vm4dvLi4THXavwel7pv\nODQZfQlOTLiD5YETEyYwrgd7UtahQ+NjntoYjh3tCdlT+lxe6L0bJ/xcp1dglQA8FFbgreCeRwCA\nOegkXXv/SI5V1XWnPa7e5ei/UVXv3eXrn6TyCACwv916rnsYd3x1d99cVZ+d5G1V9Ufd/Y6zBao8\nAgDMRPfeP8aO2zfv/HtLkl9L8rTdYiWPAAAXsKq6uKoe9uDHSf5Wkut3i3fZGgBgLpazYOYxSX6t\nthe5HUjy77v7P+4WLHkEAJiLfuh3mOnuDyf54tF4l60BABim8ggAMBO1An0eJY/Mws0vfcJw7H/4\nP18zFHfVqW8aHvMvP375cOz99429bQ7ctzk8Zvqi4dD1KY2vRw+/OT7XunismXaSZMJcc/LkUFhP\naCZda3t/+WfKmL21qIs7g83Ha8LxV6mhuebjsFSSRwCAOejYYQYAgP1F5REAYBZqKautp5I8AgDM\nhcvWAADsJyqPAABzofIIAMB+ovIIADAXK1B5HE4eq+pvJ/nrSY48+LnufsUiJsWFp6+7fjj2ile9\ndCjuQy/52eExv+DW7xyOvf/EJUNxBx44NDzmkVvHV9et3Xd4ODYHx97iNaXp8qlTw6GTfgcONnOu\nrfFm1r213N/CNaGfeqY0Hx9sqL5vTWn8vTZ4EhbVJF2TcqborMRq66HL1lX1c0m+Ncn3ZHvTim9J\n8nkLnBcAADM0es/j3+zu70xyR3f/8yRfmeQLFjctAIALT/XeP/baaPJ4/86/91XV5yY5leSxez8d\nAADmbPSexzdV1SOS/Msk78v2VflXLmxWAAAXohW4TXY0efzx7j6R5A1V9aZsL5p5YHHTAgBgjkYv\nW7/zwQ+6+0R333n65wAAuDCcs/JYVZ+T5LIkF1XVU7O90jpJLk1ydMFzAwC4oCxigcteO99l629I\n8qIklyf5idM+f1eSf7KgOQEAMFPnTB67+zVJXlNV39zdb3iI5gTn9IQfHbtj4h8/50uGx/zuK353\nOPZn7vzaobi1jYPDY9ak/sRjTcqT5NCf3zkU14enzHXCZDc2xsc9NDaHvn/8dusabJKeJBn8vnpK\n0+cpTconxNb6lO7jo8bH7M3NscCasgPugpp0j1pUM+/RBvyr1Ex8yqYCq/R9zcV+aRKe5Her6uer\n6i1JUlVPqaqXLHBeAADM0Gjy+O+SvDXJ5+48/+Mk37+QGQEAXIh6QY89Npo8Huvu12fnukJ3byQZ\nvG4BAMCQfZQ83ltVj3pwClX1FUnGbqYCAGDfGL2b/AeSXJvk86vqd5M8OsnzFzYrAIAL0H5o1fOg\nG5L8WpL7ktyd5Nezfd8jAAAXkNHk8Rey3dvxf955/m1JfjHJtyxiUgAAF6R9VHm8orufctrzt1fV\nDYuYEOyV67/p8cOxP/RbbxqOfeoX/ulQ3Pvv+WvDY9bGeD+89fvH+/HVYy4dijtw94nhMaeY1K1s\ntHfgFKfG+0z2yZN7fvgpfSaHeycmqRobtzfHeydOOf5o/8ZaG/8J6K0JPSF7Qk/I0di1Cb0zpxx/\nlUzp38jirEDyOPpufd/OIpkkSVU9Pcl1i5kSAABzdb69rf8w2znwwSS/V1Uf2Xn+eUn+aPHTAwC4\nMFTvjwUzz3lIZgEAwEo4397Wf/ZQTQQA4IK3j/a2/rRU1VVV9aGqurGqfvgsX398Vb29qn6/qt5f\nVd+4yPkAAMzaPtphZrKqWk/y00meneQpSV5YVU85I+x/SvL67n5qkhck+ZlFzQcAgM/ceB+J6Z6W\n5Mbu/nCSVNXrkjwv2w3HH9RJHuwl8vAkf77A+QAAzNp+WDDzmbgsyUdPe35TkqefEfPPkvxGVX1P\nkouTfN3ZBqqqq5NcnSRHcnTPJwoAwJhFJo8jXpjk1d39r6rqK5P8YlVd0f1XO7B29/Ekx5Pk0nrk\nCuTkzMHGn35kOPafvfQlw7Gv/bf/21Dc8+5+8fCYt28dG46tHn/bXnLz2I3XtTHe9Hhtffxm7rUJ\nTYfr5KmxuPXxZs69NjbmpHE3xhuPT2l8XhNeq+69/zU46XVdQEP3Scef8PM62tB8kkWMmQV9Txqa\nj5tLdjGXeZzDIhfM3Jzkcac9v3znc6d7SZLXJ0l3vzPJkSTj/4sCAPCQWmTy+J4kT6qqJ1bVoWwv\niLn2jJiPJHlWklTVF2Y7efz4AucEADBP/V8ahe/lY68tLHns7o0kL03y1iQfzPaq6g9U1Suq6rk7\nYT+Y5Lur6j8neW2SF/UirscAAKyCJbbqqar1nfaJbzpX3ELveezuNyd58xmfe9lpH9+Q5KsWOQcA\nAIZ8X7YLfpeeK2ihTcIBAJhgSZXHqro8yd9O8srzxUoeAQD2t2NVdd1pj6vPEvOTSf5RBpb+L7tV\nDwAAOxbUJPzW7r5y12NWPSfJLd393qp6xvkGU3kEALiwfVWS51bVnyZ5XZKvrapf2i1Y5RGSHH7L\ne4Zj/5tX/NBQ3Cv/6U8Oj/n8T/zd4dj777poOHbz8Nhb/NLx/sw5dOd48MHNBTSzntJ4fGu8QXIP\nNlOuKb8218f/Pu/BJunbcxgNnPL9j5+ryoQfmOEJLKZJdq2NvVq9tUKNPlap8fcimnmz57r7R5L8\nSJLsVB7/YXd/+27xkkcAgLlYgb9jJI8AACRJuvu3k/z2uWIkjwAAc7CgHWH2muQRAGAuViB5tNoa\nAIBhKo8AAHOh8ggAwH6i8ggTPeqV7xyKe8na9w+P+U1/913Dsdfc/2XDsf3RQ0Nx954Y79vX6+N9\n29ZObA7Hrg/2Gay7xsfMgfFfcbU2+Lf0lN6Rpyb0blyfcA5Ge0JO6Yk5oc9jb20Mx45PYELvxgW0\nmVzEmEkm9GRczAR6c8L7ZcI5WLrR13WVvqds93BdhQUzq/WqAgCwVCqPAABzsQKVR8kjAMAcrEif\nR5etAQAYpvIIADAXKo8AAOwnKo8AAHOxApVHySMAwEyswoIZySMsyLHjY83Ek+Ttm185HPvlL/7j\n4dj39OcPRh4cHrM2JjRz3jw8HHvwnrFfRwdHm3knqVPjDZLX7ntgLHBCM+1Jjb9PnBwf9/Bg8+/N\n8YbmCzGhSXmmNLOeYvDnZVIz7a0F/O8+4WdlyvHrwPj7ZfQ1qLXx89pT5jpp3LHva8qYWdCP4H4k\neQQAmIsVqDxaMAMAwDCVRwCAOeisROVR8ggAMBOrsGDGZWsAAIapPAIAzIXKIwAA+4nKIwDATKzC\nPY+SR5iBR/38eEPxj93x9OHYQ996/1DcA4+d0Ei3JvzaqPHGx0cH43r9ouExD9029v0nyVaODMXV\nlGbeE5pk14EJTaIfODEWN+FU9QODTdKT1IHBgac0aJ7SUH1jY3zcrbFG6VMaumdC6Ojxe0Lz+fSE\n5u81oan+6PmaNOZiGtUPn68prxXDJI8AAHOh8ggAwJAV6fNowQwAAMNUHgEAZqB2HnOn8ggAwDCV\nRwCAuViBex4ljwAAM6HPI7Dnjr7x3cOxl9/5ZUNxn/j+TwyPecfRS4ZjT118eDi2B/u2XfwXw0Nm\n7dLx49epsX5w6/eON/mrzc3h2LW77huOzeFDY3FbU3oHToidMu7wmOOv1ZR7wob7Jw72Y0ySrO39\nHV+T7nOb0D90Uv/ISQ0sB42f1kn9I5c6JpJHAIDZWIHKo5QcAIBhKo8AAHOxApVHySMAwBz0aiyY\ncdkaAIBhKo8AAHOh8ggAwH6i8ggAMBPLuOexqo4keUeSw9nODa/p7pfvFi95hH3swG+9dyjucz50\n2fCYt//ow4djtw6PN16++wljF0K21sd/bT18Y/y3cG2ONUjeuHj8+IdvPzEcu/WI8ebrdf/JsbiT\np8bHXD8yHDvsgfHvf9L/lxeNn4Ph12BCQ/cpRpvfL0pNaCi+kNfgwPi5mtbQnD12IsnXdvc9VXUw\nye9U1Vu6+11nC5Y8AgDMxRJy6N7O3O/ZeXpw57HrTNzzCAAwE9V7/0hyrKquO+1x9acct2q9qv4g\nyS1J3tbdu+6Fq/IIALC/3drdV54roLs3k3xJVT0iya9V1RXdff3ZYlUeAQDmoBf0mDKF7k8keXuS\nq3aLkTwCAFzAqurROxXHVNVFSb4+yR/tFu+yNQDAXCxn0fljk7ymqtazXVh8fXe/abdgySMAwAxU\nltPnsbvfn+Spo/EuWwMAMEzlEcjGTTcPxz75e28djr3t73zpcOytT98Yijtxz8HhMe/eHP8Vd+ju\nsT/3D9w/XhY4cN/48dfWx5s5r22MNV/fuuTw+Jh33T8cO9x4+6LxxuM1oUF03zs+1xwa/XkZ/7nK\nlObra4PndUKD7kUVpurg2M9rb443/8/WeOykn4EFNXWfhRXola7yCADAsIUmj1V1VVV9qKpurKof\n3iXmv62qG6rqA1X17xc5HwCAOavuPX/stYVdtt5ZsfPT2V7ufVOS91TVtd19w2kxT0ryI0m+qrvv\nqKrPXtR8AABm7dPoy7gMi6w8Pi3Jjd394e4+meR1SZ53Rsx3J/np7r4jSbr7lgXOBwCAz9Aik8fL\nknz0tOc37XzudF+Q5Auq6ner6l1VddZu5lV19YP7MZ7KiQVNFwBguRa0t/WeWvZq6wNJnpTkGUku\nT/KOqvobO1vjfFJ3H09yPEkurUeuQEEXAGB/WmTl8eYkjzvt+eU7nzvdTUmu7e5T3f0nSf4428kk\nAMCFZ8l7W49YZPL4niRPqqonVtWhJC9Icu0ZMb+e7apjqupYti9jf3iBcwIAmK0L+rJ1d29U1UuT\nvDXJepJXdfcHquoVSa7r7mt3vva3quqGJJtJfqi7b1vUnIDPXJ8Yv+/4ka9653Dso95/xVDch/7+\n+G/Ck48Yb/x8yUfH/pY+fMf48bcOjjfpPnz7eOPp9QNjc12/f3zM8WbaydZFh4bi6v6T48ef0E6k\nHnbx+LD33jcWuDHWpD7JpNdquKH4+vrwkLU1IRsYbVKepE8NvgZrE+pOU2InqBr8vhbQpoYF3/PY\n3W9O8uYzPvey0z7uJD+w8wAAuLCtQL5rhxkAAIYte7U1AABJsqB7FPeayiMAAMNUHgEA5mIFKo+S\nRwCAGai4bA0AwD6j8gjMQl93/VDck//eeO/Eu5/31OHYW75sLG7jovG+eQfuH49dPzmhz99gZaIP\njtcHDkzpxzfYO68Pjf8Xs3b3A+PH3xjvH1kHB3syXnRk/Pij/RCT5NBYT8wpfSZ7bcLPyvqE81qD\nr+uE3pHZ3ByPnWLw++rNrcUcf5FWoDelyiMAAMNUHgEAZmIV7nmUPAIAzEFnJVZbu2wNAMAwlUcA\ngJmoFVjjo/IIAMAwlUcAgLlYgXseJY8AADNhtTXAHusTJ4ZjL3n9u4ZjH/62zxqKu/uZXzA85j2f\nO97M+b5j47FH1seaNK9tjP8vVBNi106NNX6ujfGbt/rwgv47GuwpXydPTRhzsPF3ktw/2Px8QY2h\ne8K4dfHRsTFPTXitakJD8ZpwJ92psYbmk5qkM0zyCAAwBx07zAAAsL+oPAIAzMQq3POo8ggAwDCV\nRwCAuViByqPkEQBgBiouWwMAsM+oPAIAzEH3SrTqkTwCJNm8446huKNvfPfwmJc+9nOGYz/x1Z83\nHHvq6NhFo4P3jTfpfuDYeOPrA/eNNQlfPzF+/AN3L/c/zDq1MR583/3jsUcGu5RPSRjuG2w8nqSm\nNOneHDuvk2yO/wwkE2LXB5vqL+J7QvIIADAXq3DPo+QRAGAuViB5tGAGAIBhkkcAgJmo3vvHeY9Z\n9biqentV3VBVH6iq7ztXvMvWAAAXto0kP9jd76uqhyV5b1W9rbtvOFuw5BEAYA46ydZDf9Njd38s\nycd2Pr67qj6Y5LIkkkcAgFlbTO54rKquO+358e4+frbAqnpCkqcm2bUvmeQRYEE2PvYXw7GX/Op4\n7PqjHjkUt/nXLhse877Ljw7HHnhgrHfe1vr4bfUbDxvsh5jk4O33DcfWibH+jT1hrnnkw8djT54a\nCpvSZ7JRyknfAAALW0lEQVSOHhk//saEPoejPREn9E6sw+P9Q/vEyeHYYQekOTtu7e4rzxdUVZck\neUOS7+/uu3aL86oCAMzEsvo8VtXBbCeOv9zdbzxXrNXWAAAXsNreiujnk3ywu3/ifPGSRwCAuXhw\nf+u9fJzfVyX5jiRfW1V/sPP4xt2CXbYGALiAdffvJBneCF3yCAAwE/a2BgBgTMfe1gAA7C8qjwAA\nM1BJamyBy1JJHgFWzOZtt48FjsYlOXrd+nDsgcc8eiju1BMfMzxmTdmSbWs8NGuDawA2x49fd94z\nfvzRJtmTvv8pL8AENfZa1ZEpTcrHm59r6L06nCkAgLlY0N8Ge0nyCAAwE6tw2dqCGQAAhqk8AgDM\ngVY9AADsNyqPAACzMLwX9VJJHgEAZmIVtid02RoAgGEqjwAkW5vDoRsf+4uhuBqMS5K1iy8ejq2j\nR4dj8+jPGgrrI4fHj39qvPF1r4/VaGprQbWcKU26By+X9gMPDA9ZUxp/jzZ0T5KTp/Z+zLlYgcvW\nKo8AAAxTeQQAmINOagV2mFF5BABgmMojAMBcrMA9j5JHAIC5mH/u6LI1AADjVB4BAGaiLvTL1lV1\nVZKfSrKe5JXd/b/sEvfNSa5J8uXdfd0i5wTA/Gzde+948JTYj398KKwOHhoesi++aPz4g+P2hD6b\ndeTI+PHX14dD+9Rg78Qav2jZU/pMTokdNf6yMsHCkseqWk/y00m+PslNSd5TVdd29w1nxD0syfcl\nefei5gIAsBJWoPK4yHsen5bkxu7+cHefTPK6JM87S9y/SPJjScZb1gMA7DedZGsBjz22yOTxsiQf\nPe35TTuf+6Sq+tIkj+vu/+tcA1XV1VV1XVVddyon9n6mAAAMWdqCmapaS/ITSV50vtjuPp7keJJc\nWo+cfz0XAGCiSq/EgplFVh5vTvK4055fvvO5Bz0syRVJfruq/jTJVyS5tqquXOCcAAD4DCyy8vie\nJE+qqidmO2l8QZJve/CL3X1nkmMPPq+q307yD622BgAuWCtQeVxY8tjdG1X10iRvzXarnld19weq\n6hVJruvuaxd1bACAlXQhJ49J0t1vTvLmMz73sl1in7HIuQAA8JmzwwwAF7w+dXI4dvMT47HDqoZD\n1w4fHh93QpPwrI0tg6gJc53UpHxzQk+ZrcHYwe9pNh5s1TNzK/aqAgCwTCqPAAAzcaG36gEAYJ9R\neQQAmIsVqDxKHgEAZqFXInl02RoAgGEqjwAAc9BReQQAYH9ReQSAZZtQbdp64IEFTmTAlCbhNV6j\nWjt0cDi2R1+vrflX8T7FEpqEV9WrkjwnyS3dfcX54lUeAQBmorr3/DHg1UmuGp2j5BEA4ALW3e9I\ncvtovMvWAABzsZgFM8eq6rrTnh/v7uOf7mCSRwCA/e3W7r5yrwaTPAIAzEFnJRb5SB4BAGbBDjMA\nAMxcVb02yTuTPLmqbqqql5wrXuURAGAullB57O4XTomXPAIA46YkN705HLr1wHgsyyV5BACYC/c8\nAgCwn6g8AgDMgVY9AACM66S3lj2J83LZGgCAYSqPAABzYcEMAAD7icojAMAcWDADAMAkLlsDALCf\nqDwCAMyFyiMAAPuJyiMAwCz0SlQeJY8AAHPQSbbsMAMAwD6i8ggAMBcrcNla5REAgGEqjwAAc6Hy\nCADAfqLyCAAwC21vawAABnXSrVUPAAD7iMojAMBcrMBla5VHAACGqTwCAMzFCrTqkTwCAMxBt72t\nAQDYX1QeAQDmYgUuW6s8AgAwTOURAGAmegXueZQ8AgDMQrtsDQDA/qLyCAAwBx07zAAAsL+oPAIA\nzEXPf8GMyiMAAMNUHgEAZqCT9Arc8yh5BACYg26XrQEA2F8WmjxW1VVV9aGqurGqfvgsX/+Bqrqh\nqt5fVb9VVZ+3yPkAAMxZb/WeP0acL2c73cKSx6paT/LTSZ6d5ClJXlhVTzkj7PeTXNndX5TkmiQ/\nvqj5AADwqQZztk9aZOXxaUlu7O4Pd/fJJK9L8rzTA7r77d19387TdyW5fIHzAQCYt97a+8f5nTdn\nO90iF8xcluSjpz2/KcnTzxH/kiRvOdsXqurqJFfvPD3xm33N9XsyQxbtWJJblz0JhjhXq8X5Wh3O\n1ep48rIncHfueOtv9jXHFjD0kaq67rTnx7v7+GnPJ+Vss1htXVXfnuTKJF9ztq/vfIPHd2Kv6+4r\nH8Lp8WlyrlaHc7VanK/V4VytjjOSq6Xo7quWPYcRi0web07yuNOeX77zub+iqr4uyT9N8jXdfWKB\n8wEA4FMN5WwPWuQ9j+9J8qSqemJVHUrygiTXnh5QVU9N8m+TPLe7b1ngXAAAOLvz5mynW1jlsbs3\nquqlSd6aZD3Jq7r7A1X1iiTXdfe1Sf5lkkuS/GpVJclHuvu55xn6+Hm+znw4V6vDuVotztfqcK5W\nxwV7rnbL2XaLr+75b4MDAMA82GEGAIBhkkcAAIbNNnkc2NrwRVX18ar6g53Hdy1jniRV9aqquqWq\nztp/s7b9651z+f6q+tKHeo5sGzhXz6iqO097X73soZ4jSVU9rqrevrN96weq6vvOEuN9NQOD58r7\naiaq6khV/aeq+s875+ufnyXmcFX9ys57691V9YSHfqbzNos+j2c6bZucr892o8r3VNW13X3DGaG/\n0t0vfcgnyJleneTfJPmFXb7+7CRP2nk8PcnP5twN41mcV+fc5ypJ/p/ufs5DMx12sZHkB7v7fVX1\nsCTvraq3nfE70PtqHkbOVeJ9NRcnknxtd99TVQeT/E5VvaW733VazEuS3NHd/1VVvSDJjyX51mVM\ndq7mWnmctE0Oy9Xd70hy+zlCnpfkF3rbu5I8oqoe+9DMjtMNnCtmoLs/1t3v2/n47iQfzPYOEKfz\nvpqBwXPFTOy8X+7ZeXpw53HmyuHnJXnNzsfXJHlW7bSEYdtck8ezbZNztjfjN+9crrmmqh53lq8z\nD6Pnk3n4yp1LOm+pqr++7Mlc6HYumT01ybvP+JL31cyc41wl3lezUVXrVfUHSW5J8rbu3vW91d0b\nSe5M8qiHdpbzNtfkccR/SPKE7v6iJG/Lf/krAfj0vS/J53X3Fyf535P8+pLnc0GrqkuSvCHJ93f3\nXcueD7s7z7nyvpqR7t7s7i/J9i4qT6uqK5Y9p1Uz1+TxvNvkdPdtp21n+MokX/YQzY3pJm17xPJ0\n910PXtLp7jcnOVhVx5Y8rQvSzv1Yb0jyy939xrOEeF/NxPnOlffVPHX3J5K8PcmZ+0l/8r1VVQeS\nPDzJbQ/t7OZtrsnjyNaGp9/b89xs32fCPF2b5Dt3Vod+RZI7u/tjy54Un6qqPufBe3uq6mnZ/h3h\nl+ZDbOcc/HySD3b3T+wS5n01AyPnyvtqPqrq0VX1iJ2PL8r2wtw/OiPs2iT/3c7Hz0/yf7cdVf6K\nWa62Htza8Hur6rnZXul2e5IXLW3CF7iqem2SZyQ5VlU3JXl5tm9CTnf/XJI3J/nGJDcmuS/Ji5cz\nUwbO1fOT/L2q2khyf5IX+KW5FF+V5DuS/OHOvVlJ8k+SPD7xvpqZkXPlfTUfj03ymp2uLmtJXt/d\nbzojv/j5JL9YVTdmO794wfKmO0+2JwQAYNhcL1sDADBDkkcAAIZJHgEAGCZ5BABgmOQRAIBhkkdg\n36mqJ1TV9cueB8B+JHkEAGCY5BHYrw5U1S9X1Qer6pqqOrrsCQHsB5JHYL96cpKf6e4vTHJXkr+/\n5PkA7AuSR2C/+mh3/+7Ox7+U5KuXORmA/ULyCOxXZ+69ai9WgD0geQT2q8dX1VfufPxtSX5nmZMB\n2C8kj8B+9aEk/6CqPpjks5L87JLnA7AvVLcrOQAAjFF5BABgmOQRAIBhkkcAAIZJHgEAGCZ5BABg\nmOQRAIBhkkcAAIb9/zjem8h/atqkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ad95190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(12,8))\n",
    "_=hist2d(b_samples[200:].flatten(), beta_samples[200:].flatten(), bins=50, range=[(0.5,3),(0.2,1.5)], normed=True)\n",
    "colorbar()\n",
    "xlabel(\"b\")\n",
    "ylabel(\"beta\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}