mamacneil/GammaGLM.ipynb

## GammaGLM.ipynb
{
    "metadata": {
        "name": "",
        "signature": "sha256:295922906e09ca3fdcf9b7d9c33def73124780dec5e063873de3d0bca79b8727"
    },
    "nbformat": 3,
    "nbformat_minor": 0,
    "worksheets": [
        {
            "cells": [
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "# Example of Gamma GLM in [PyMC2](https://github.com/pymc-devs/pymc)\n",
                        "\n",
                        "Example of Gamma GLM, based on simulated data from [Sean Anderson's blog](http://seananderson.ca/2014/04/08/gamma-glms.html)"
                    ]
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "%matplotlib inline\n",
                        "import numpy as np\n",
                        "import pymc as pm\n",
                        "import matplotlib.pyplot as plt\n",
                        "\n",
                        "plt.style.use('fivethirtyeight')"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [],
                    "prompt_number": 1
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Simulate data under shape-rate configuration of the gamma used in PyMC"
                    ]
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "N = 100\n",
                        "x = pm.runiform(-1,1,N)\n",
                        "a_true = 0.5\n",
                        "b_true = 1.2\n",
                        "y_true = np.exp(a_true+b_true*x)\n",
                        "shape_true = 10\n",
                        "y = pm.rgamma(alpha=shape_true,beta=shape_true/y_true,size=N)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [],
                    "prompt_number": 2
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "plt.scatter(x,y)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "metadata": {},
                            "output_type": "pyout",
                            "prompt_number": 3,
                            "text": [
                                "<matplotlib.collections.PathCollection at 0x10e8833d0>"
                            ]
                        },
                        {
                            "metadata": {},
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEWCAYAAAC5XZqEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwVPXdx/HPJhAupkAikAgafZAgIgjocCkGIlZAGRWo\nF6Bar1Ex0MtYHMGhVK0aq+hUJAEloCiGGa0igtjHqcZIGB8vtTXWqlkfDOKFjRoicklCsvv8QTdP\nEnaTPWfPnj179v2a8Q9idvf78+D57O96PPX19QEBAGCjlHgXAABIPoQPAMB2hA8AwHaEDwDAdoQP\nAMB2hA8AwHaEDwDAdobDZ9SoUcrIyDjmn7lz58aiPgCAC3Uz+oKKigq1tLS0/vmbb77Rueeeqzlz\n5lhaGADAvQyHT2ZmZrs/b9iwQX369CF8AAARi2rOJxAI6Omnn9YVV1yhHj16WFUTAMDlogqf8vJy\nffHFF7rmmmusqgcAkASiCp8NGzbo7LPP1hlnnGFVPQCAJGA6fL799lu98soruvrqq62sBwCQBEyH\nT1lZmXr27KnLLrvMynoAAEnAVPgEAgE99dRT+vnPf67evXtbXVPC83q98S4hbmh78knWdkvJ3fZo\nmQqfHTt26PPPP2ehAQDAFMP7fCRpypQpqqurs7oWAECS4Gw3AIDtCB8AgO0IHwCA7QgfAIDtCB8A\ngO0IHwCA7QgfAIDtCB8AgO0IHwCA7QgfAIDtCB8AgO0IHwAJw+fzyOfzxLsMWIDwAZAQKitTlZ+f\nrvz8dFVWpsa7HESJ8AHgeD6fRwUFvbV3b4r27k1RQUFvekAJjvABANiO8AHgeFlZAZWWHlJ2tl/Z\n2X6Vlh5SVlYg3mUhCqYeJgcAdsvLa1FFxQFJInhcgPABkDAIHfdg2A1A0mMJt/0IHwBJjSXc8UH4\nAEhaLOGOH8IHAGA7wgdA0mIJd/yYCp+9e/dqwYIFGjp0qLKzszVx4kTt3LnT6toAIOaCS7grKg4o\nL68l3uUkDcNLrevr6zVjxgxNmjRJzz33nI4//njV1NRowIABsagPAGKO3o79DIfPypUrNWjQIK1e\nvbr1Zzk5OZYWBQBwN8PDbi+//LLOOussXXfddcrNzdXkyZO1du3aWNQGAHApw+FTU1OjdevWaciQ\nIXrhhRe0YMEC3XXXXQQQACBihofd/H6/zj77bP3+97+XJI0aNUq7du1SaWmpbrzxRssLBAC4j+Hw\nyc7O1mmnndbuZ7m5ufryyy/Dvsbr9RqvLMElY5uDaHvySdZ2S8nZ9tzc3Kjfw3D4TJw4UdXV1e1+\n9tlnn3W66MCKQhOJ1+tNujYH0fbka3uytltK7rZHy/CcT2Fhod577z099NBD2rVrl1588UU9/vjj\nKigoiEV9AAAXMhw+Y8eO1TPPPKPNmzdr0qRJuvfee7Vs2TLdcMMNsagPAGKCk6zjy9TzfKZPn67p\n06dbXQsA2KKyMlUFBb0lSaWlhzjZIA442w1AUuEka2cgfAAAtiN8ACQVTrJ2BlNzPgCQyIInWUsc\nKhovhA+ApEToxBfDbgAA2xE+AADbET4AANsRPgAA2xE+AADbET4AYILP51Fzc/94l5GwCB8AjuP0\nQz8rK1OVn5+uOXNyVFmZGu9yEhLhA8BRgjf2/Px0R97YORvOGoQPAMfgxp48CB8AMICz4azB8ToA\nHCN4Y2/7rB0n3tiDZ8PV1dXp9NMz4l1OQiJ8ADhKohz6mZUV0P7930kifMwgfAA4jpNDB9ZgzgcA\nYDvCBwBgO8IHAGA7wgcAYDvCBwBgO8PhU1RUpIyMjHb/DB8+PBa1AQBcytRS62HDhmnbtm2tf05N\ndd75SwAA5zIVPqmpqRowYIDVtQAAkoSpOZ+amhqdfvrpGj16tG644QbV1NRYXBYAwM0Mh8+4ceO0\nevVqPf/881q5cqV8Pp9mzJihffv2xaI+AIALGR52O//889v9edy4cRo9erTKysq0cOFCywoDALiX\np76+PupDlC6++GKddtppWrFiRch/7/V6o/0IAIBD5ObmRv0eUR8s2tDQoOrqak2ZMiXs71hRaCLx\ner1J1+Yg2p58bU/WdkvJ3fZoGQ6fZcuW6cILL9TgwYP13Xff6cEHH9Thw4c1f/78WNQHAHAhw+Hz\nzTffqKCgQN9//7369++vcePG6W9/+5tOPPHEWNQHAHAhw+Gzbt26WNQBAEginO0GwDCfzyOfzxPv\nMpDACB8AhlRWpio/P135+emqrExVSgq3ERjH3xoAEfP5PCoo6K29e1O0d2+KCgp6q6kpM95lIQER\nPgASFsN/iYvwARCxrKyASksPKTvbr+xsv0pLDyktrS4utXQc/kNiiXqTKYDkkpfXooqKA5KOhpHX\n67e9hrbDf5JUUNBbFRUHlJUV9YEtsAnhA8AwbvKIFsNuACLipPmVUMN/BGJioecDoEuVlakqKOgt\nSSotPaS8vJY4V3Ts8B8SCz0fAJ0KtbzaST0ggicxET4AANsRPgA6xfwKYoE5HwBdYn4FViN8AESE\n0IGVGHYDgE44aYm5mxA+ABAGR/jEDuEDACE4eYm5GxA+AADbET4AEoad8y8sMY8tVrsBSAjxOOKH\nJeaxQ88HgOPFc/6FI3xig/ABANiO8AHgeMy/uA9zPgASAvMv7hJVz+fhhx9WRkaGbrvtNqvqAYCw\nmH9xD9Ph8+6772rDhg0644wz5PGw8QoAEDlT4fPDDz/opptuUnFxsfr162d1TQBihHPK4BSmwue3\nv/2tZs+erby8PAUCdIGBRMA5ZXASw+GzYcMG1dTUaNmyZZLEkBuQADinDE5jaLWb1+vVH//4R/31\nr39VaurRb06BQKDL3o/X6zVfYYJKxjYH0XbnaW7uLym93c/q6uq0f/93lry/U9tth2Rse25ubtTv\n4amvr4943OyZZ57RokWLWoNHklpaWuTxeJSamqqvv/5a3bt3j7qoROf1ei25OImItju37bE6niZe\n7Q723OK5+s3p19zJDPV8LrroIp199tmtfw4EAlq4cKGGDh2qW2+9leABHMxN+2Ticc4brGUofPr2\n7au+ffu2+1mvXr3Ut29fDR8+3NLCAFgv0UNHaj9/JUkFBb1VUXHAFW1LJlEfr+PxeFh0AAAwJOrj\ndbZt22ZFHQAQkeA5b22H3aLt9Thh/ijZcLYbgIQTyfxVx0AJFzDMH8UHp1oDSEidnfPWdkPt//xP\nStgNtqH2P9XWchKEHQgfwMWS8SbaMVD++7+7R7zBtl8/v/79b06CsAPhA7gUx+l0reNzglatOqwF\nCzgJwg6ED+BQ0fRakvk4nY6BMmPGkU4fRBecP6qoOKCTTmLBgV1YcAA4kJMmwRNxJVioBQmdLVBo\n+zOrV9IhNHo+gMNY0Wux6rHTVg3dxWPuqeOChEgfRNe2J8TKt9ih5wO4VLTH6Vh1koDZXlw8e1z0\ndmKPng/gMFb1WoLvFc8bqdleHIsl3I/wARzICUM/kYZgc3N/S4fUknmxRDJh2A1wKCcM/XQ1dHd0\nSC1HUughtVgchdNWIi6GwFH0fAB0KtzQXaQ9FKO9uEh7XAzNJTZ6PgBizmjPpKseF49VSHz0fACY\nYuXCiHDvT5i4Fz0fAKb4fB7l5vr10kt71Ldvhq1BEeu5JMQe4QPAsLZ7d0pKmjVsmP03fjc9FjwZ\nMewGwJCOCw0KC/vFbSk0Q3OJi/ABANiO8AFgSMeFBiUl9fQ+YBhzPgAMazvfcuBAjaRT41oPEg/h\nA8CUYG9n/35/nCtBImLYDQBgO8IHAGA7w+Gzdu1anXPOOcrJyVFOTo6mT5+uV199NRa1AbBBPB70\nBhgOn8GDB+vuu+/Wm2++qTfeeENTpkzRlVdeqQ8//DAW9QGIIQ7nRLwYDp+ZM2fqZz/7mU455RQN\nGTJEy5YtU3p6ut5///1Y1AcgRnhuDuIpqtVuLS0tevHFF9XY2KhJkyZZVRMAwOVMhc9HH32k6dOn\nq7GxUb169dITTzyh3Nxcq2sDYEKw9+LxBBQIeMJuAM3KCmjjxoN6+eXueuWVblqxooHNorCNp76+\n3vDftiNHjuirr77SDz/8oC1btujxxx/X1q1bNXbs2JC/7/V6oy4UQOdSUlK0e/cpKizsJ0lasqRB\na9ak6b779uvkk2vk9/vD/m5x8Y8aMuRzNTc3x6FyJBorOhumwqej2bNna9CgQSopKYm6IDfwer1J\n2xOk7fFru8/nUX5+eusD1rKz/br44iPaurX7MQ9aC/W7Zh/GFu92x1Mytz1aluzzaWlpafetCgCA\nzhgOnzvvvFNvvfWWdu/erY8++kh33XWXdu7cqblz58aiPgAR6njg55IlDdqxIzXkg9Zi/RRSoCuG\nFxzU1tbqpptuUm1trfr06aORI0fq+eef19SpU2NRHwAD2h746fEEdOGFzWFDhYexIZ4Mhw/zOoCz\ntQ+SzkOF0EG8cLYbkAQ4QgdOQ/gALscROnAiwgdwMY7QgVMRPgAA2xE+gIuxpBpOxWO0AZdjSTWc\niPABkgChA6dh2A0AYDvCBwBgO8IHgGFsWkW0CB8AhnTctJqSwm0ExvG3BkhyRnoxoTatNjVlxrhC\nuBHhAyQxjt5BvBA+gM2MzpfEan7FzNE7oTatpqXVWV4b3I/wAWxktKfhxJ5JcNNqRcUB5eXxFGOY\nQ/gANjHa04j1oaDRHL2TlRVg4yqiwgkHQBLj6B3ECz0fwCZGexp2HQpKLwbxQM8HsJHRngY9E7gV\n4QPYzGiIOD10mpv7y+fzOL5OOAvDboBN3HgkTWVlqubMyXHUajwkBsIHsIETl0xHi0d0IxqEDxBj\nTrpJu7H3hcRkOHwefvhhTZ06VTk5ORo6dKjmzZunjz/+OBa1AbCQ1b0vHtGNaBgOn507d+rGG2/U\nq6++qpdeekndunXT7NmzVV9fH4v6gIRnx026qx5NrHpfeXkt2rz5i9bTDoBIGV7t9vzzz7f782OP\nPaacnBy9/fbbmjFjhmWFAW4SyyXTlZWpKijoLUkqLT1kewh06/adsrIyjvl5MNzoDSGUqOd8fvzx\nR/n9fvXr18+KegDXisVmzkh7NHYPkblxgQWsFfU+nyVLlujMM8/U+PHjragHQIzYtWG1bSBKUkFB\nb1VUHKAHhHai6vnccccdeuedd/TUU0/J42EFDZKHU1aNmTmyJysr4Jj6kbw89fX1pr6OLF26VC++\n+KK2bt2qoUOHdvq7Xq/XVHGAnVJSUlqfypmWVhfyUQEpKSnavfsUFRYeHWYuKanXySfXxPWxApHU\n3fZ3Y12/E/8bwVq5ublRv4ep8Ln99tu1ZcsWbd261ZIi3Mbr9Sbtfxentj2Sye9IJu59Po/y89Nb\nh5Sys/2tQ0pObXtbndVvVrh2J8OCg0S45k5leM5n8eLFevbZZ7Vx40b16dNHPp9PkpSenq7jjjvO\n8gKBzlgZKkbmKTIz/br00iNKTw/I43HvzTUabg4dRM/wnM+6det04MABzZo1S8OHD2/9Z9WqVbGo\nDwgrkhVVVu9vycoK6MknD2r58gZt3dpdmzalqbo6cVZzsTEUTmG457Nv375Y1AEYYvWKquBNuW0P\nKdx7/dd/BXTttT2P+exEwWMa4AQ8UgGu1jFU1qwJHyrxuClbMS9i5j0IHcQbB4siIRkZPsrLa9HG\njQc1f36Tlizp2emmx0g2glo1dGXFRkw2cyJR0fNBwoq0p+LzeXTVVcdZuukx2l6SFcOG8dzMmQwr\n2RBbhA8SWjxvfsl64+24evCEExhAgXH8rYHrOXGFlxU1xaNdoVYPBje4AkbQ80FScOIKLytqcmK7\ngEgQPkgaTrw5W1GTne0KtSQ9La1O0rGPVAA6Q/gAMKRjb8vr5cw2GEf4ABbw+Txqbu5vyftIzuyl\nteX0+uB8LDgAohTcazNnTk5Ue23Ys4NkQvgAEQr1DByrzo6z+gw6wOkIHyAC9Eqc8wA9uAPhA3Sh\ns16JVXttnLgXqS3CF1ZjwQEcywmT7z6fRwcPSv36+VuPsekouPqrrq5Op58efslxV+1x6p6deB7j\nA/ei5wNHsvKbttnhomANM2em6+67GzV8eHPYXklWVkDdun3X5Xt11Z5IDjYF3IDwgeO0/abd1CS9\n9lo37dplbq7BTIj5fB7t2tV+qO3Xv+6lsrJDqqg4EPJJqF29XyIvJnD6kCASE+EDx8rM9OuOOxq1\naVOaZs403gMyc9MPhtXTT6cd8++OO85Zw2F2Cg4JmglfIBTCB44T/KZ9zTVNWrGih209hrZhtWFD\nmpYsabDk275beg4MCcJKLDiAI+XltWjQIL82bTq2BxIpI4/G7qiuLkVr1qRp+/YDlvR4ulpM4ITF\nFYCd6PnAsYYMib7HYGS4qGMPZcWKBg0ZYt23/XA9B5YxIxnR84EjhPvmb8Xy43CvC/WZdi93Zhkz\nkhU9H8RdV9/8YzHX0NlnxmpugxMCgP9H+CCu4rEM2exnRhMe4cLOLYsRAKNMhc/OnTs1b948jRgx\nQhkZGSorK7O6LrhUon77j2ZepquwYxkzkpGp8Dl06JBGjhypoqIi9erVSx5P4t1MYL9QN/B4fPM3\n+pl29M5YxoxkY2rBwbRp0zRt2jRJ0sKFCy0tCO7U2cR6PM40i8VnpqSkhFzEEM2Sb8CtWO0GR4jH\nzTjSz4w0PHbvPkWFhemtv9N2CM2ph4YC8UL4wBZO/vYfyQbPSDaJFhb263TJtFPaCzgB4QPbOPHb\nf2VlartA7GzC3yk1A27gqa+vj+r/qBNPPFEPPvig5s+fH/Z3vF5vNB8BxERzc3/NmZPT2lvJzvZr\n8+YvOn00QjgpKSn/GXbrJ0kqKanXySfXyO/3W1ky4Ai5ublRv4ctPR8rCk0kXq836doc5IS2R3pO\nWqgVa5mZmcrKCv9AuM79ryoqhv7ns1MlnWr4HRLxjDcnXPN4Sea2R8vUUuuDBw+qqqpKVVVV8vv9\n2rNnj6qqqvTll19aXR9gSCT7cYJ7jaxe5u33+6NaMs0Zb0gmpsLn/fffV35+vvLz89XQ0KCioiLl\n5+erqKjI6vqSlpWbMWOxsdOJm0Uj2Y/T8QbvlA2eif7AOcAoU8NukydP1r59+6yuJWn4fB55PAEF\nAqGHWIxMgnelq/cyM8xjZX2RsmI4KlaHeAZrS0nhtCogUvzfYrPKylTNmtVb27d3b/0G/sYb/z/E\nYuU34K7ey+wjpqOpz0yPKdI643FaQtvadu8+xfT7cMYbkg3hY6PgjXvy5Bbdf3/P1hv4ggW99d57\nRy+FxxPQ/PlNuvHGRmVm+lt/ZvUwl9kQOXgwsp+FYkfYdTaMZvUNvmNthYX9orpGThkCBOzAPh+H\nePnl7srJaVR1dWrr0zuXL2/QyJHNqq42N8wVi42d6ekBLVnSoPvv7ylJWrKkQenpXb+nnc+tiWaz\naCRiORdDbwfJgp6PjYJhsGNHqoqKDrd+A1+8uFGvvNJNBw60/yZ9//09lZmpqIa5wn2bNtsLGDhQ\nGjOmWX/842Fde22jhg9v0cCBoX/XbG+t7etiMRxl1Yq0zz/3tKutpKSe8AAiRM/HZnl5Ldqy5ZBS\nUgLauPGgXn65u0pLu2vFigYdd5z1n9fZRL2ZXkDbxQZr1hzSxImhN1GGWpQQqhfWsb5Qr3PKyQgd\ne2/XXnuc3nzzx9baDhyokZm9PUAyInziIHgDHTDAr5NOatKCBU2tP+t4gx4yJKCVKw/r17/uJUla\nufKwqaAIN1wXfK9IVpN1vPkuWBB66CzcEFvHEOlYX26uXwUFvdXUJF166RG99lo3DRt2tGfl1B5F\nIOBprW3/fk4zACLFsFsMNDf3j3i4qeMQUMdhMp/Po+XLe+jii4/o4ouPaPnyHmHfu+1wlZGJ+uBQ\n0qxZRxc+xHpOI9jj6VjfwYNSv35+3XFHo7Zu7a5Nm9L07387Z7MlK9IA6xA+FqusTNWcOTlR7VLv\nGEj19Slau7aH1q7tofr60JfMyEqyUCHV1CQVFBzRVVcdF/Y9Ir35mr1Jp6cHtGrVYa1Y0aPdSkAn\nbbZkRRpgDcLHQrHYpR7JjTzU53o8oV/XMaQ8nqPvdemlR9rd9MPV3tXNN/hAtWHDOv+9UO0aOFA6\n6STn9yR46igQPeZ8EoCZCfdAwHPM68LNxZSWHtJrr0X+V6GzGjp7oFpHwfoOHlTrcm0nP/cHgHXo\n+VgoFpsY2y4ECPdenX1uJN/S8/JadMstjVqzJrra2z5QLdKen9ebopkz0zVlyk9ah/oY2gLcj56P\nxfLyWrR58xf/OZrffPAYPT8tkt5RZ72KgQOlgQPtXdLc2cZTejuAuxE+MdCt23dRPBPG/GkAkdyw\nuwqpaG76WVkBlZTUtz5QjSEzAOEQPkmg4x6eWAbCySfXtHmgWtdhyfwOkJyY84kBI/t8QrFy7sju\nB5QZfaAa8ztAcqLnY7GjczU5kqJ71o1VB2DadZhnNJxWD4DYo+djIav3+TDxDsCtCB8X4zgYAE7F\nsJuFnDiB7pQToQGgLcLHYlbt87GSU+oAgCCG3WLg6D4fbvgAEA7hAwCwHeEDALCd6fApLS3VmWee\nqezsbJ177rl66623rKwLAOBipsLnhRde0NKlS7V48WLt2LFD48eP1+WXX64vv/zS6voAAC5kKnyK\ni4t15ZVX6uqrr1Zubq4eeOABZWVlaf369VbXBwBwIcPh09TUpA8++EBTp05t9/PzzjtPb7/9tmWF\nAQDcy3D4fP/992ppadHAgQPb/bx///6qra21rDAAgHux2i0GcnNz411C3ND25JOs7ZaSu+3RMhw+\nxx9/vFJTU4/p5Xz77bfKysqyrDAAgHsZDp+0tDSNGTNG5eXl7X5eXl6uCRMmWFYYAMC9TJ3ttnDh\nQt18880666yzNGHCBK1fv161tbW67rrrrK4PAOBCpsJnzpw5qqur04oVK+Tz+TRixAg9++yzOvHE\nE62uDwDgQp76+npOwAQA2Mqy1W5PPvmkLrroIuXk5CgjI0N79uzp8jXPPPOMMjIy2v2TmZmppqYm\nq8qyhZm2S9KWLVs0YcIEZWVlaeLEidq2bVuMK7VWY2OjbrvtNp166qkaPHiw5s+fr6+//rrT1yTy\nNTd6pNRHH32kmTNn6oQTTtCIESP0wAMP2FSptYy0e/fu3cdc34yMDL3++us2Vhy9nTt3at68eRox\nYoQyMjJUVlbW5Wvccr2Ntt3sNbcsfA4fPqzzzz9fS5cuNfS63r17y+v1qrq6WtXV1fr000+VlpZm\nVVm2MNP2d955RzfccIPmzp2ryspKXX755br22mv197//PYaVWmvp0qXatm2b1q9fr+3bt+vHH3/U\n3Llz5ff7O31dIl5zo0dK7d+/X3PmzFF2drbKy8tVVFSkRx99VKtWrbK58uiYPUrrhRdeaL2+1dXV\nmjx5sk0VW+PQoUMaOXKkioqK1KtXL3k8nk5/3y3XWzLe9iCj19yyh8ndcsstkqR//OMfhl7n8XjU\nv39/q8qICzNtX716taZMmaJbb71VkvS73/1OO3bs0OrVq1VaWhqTOq30ww8/aOPGjSopKVF+fr4k\n6bHHHtOoUaP0xhtv6Lzzzgv72kS85m2PlJKkBx54QK+99prWr1+v5cuXH/P7zz33nBoaGrR69Wr1\n6NFDw4cPl9frVUlJiRYtWmR3+aYZbXdQRkaGBgwYYFeZlps2bZqmTZsm6egCq6645XpLxtseZPSa\nx32T6eHDhzVq1CidccYZmjt3rqqqquJdki3efffdhD6i6J///KeOHDnSLmQGDx6s0047rcs2JNo1\nN3Ok1DvvvKOf/vSn6tGjR7vf/+abb/TFF1/EtF6rRHOU1lVXXaXc3FxdcMEF2rJlSyzLdAQ3XO9o\nGb3mcQ2fYcOGqbi4WJs2bVJpaal69uypCy64QLt27YpnWbaora095oiiAQMGJMwRRbW1tUpNTVVm\nZma7nw8YMEDffvtt2Ncl4jU3c6RUuOsb/HeJwEy7f/KTn+iee+7Rhg0b9Nxzz2nKlCm6/vrr9eyz\nz9pRcty44XqbZfaadzrsds899+ihhx7q9A22bdumc845x3jFksaNG6dx48a1/nnChAmaPHmyHnvs\nMf3pT38y9Z5WiXXbnSrSdpvl5GtupUjHyd0mMzOz3VDNmDFjtG/fPj3yyCO64oor4lhZbCXr9ZbM\nX/NOw6ewsFDz5s3r9IMHDx5ssNTwUlJSNHr0aEd8C4512wcOHBjyiKKO357sFmm7m5ub1dLSorq6\nuna9n9raWk2aNCniz3PSNQ/HzJFS4a5v8N8lAquO0ho7dqw2btxodXmO4obrbaVIrnmn4ZOZmXnM\nsEosBQIB/etf/9Lo0aNt+8xwYt328ePHq7y8XL/61a9af1ZeXq6JEyfG7DMjEWm7x4wZo+7du+v1\n11/XZZddJkn66quvVF1dbeiYJSdd83DaHik1a9as1p+Xl5dr9uzZIV8zfvx43XnnnWpsbGydBygv\nL9egQYOUk5NjS93RMtPuUD788ENlZ2fHokTHcMP1tlIk19yyOR+fz6eqqip99tlnkqRPPvlEVVVV\nqq+vb/2dSy65RHfffXfrn++//369/vrrqqmpUVVVlRYtWqRPPvlE119/vVVl2cJM2xcsWKA333xT\nf/7zn1VdXa2HH35YlZWVrSvnnK5v37765S9/qT/84Q+qqKjQBx98oJtvvlkjR47Uueee2/p7brnm\nCxcuVFlZmZ566il9+umnuv3229sdKXXXXXe1u0Ffdtll6tWrlwoLC/Xxxx/rpZde0iOPPKLCwsJ4\nNcEUo+0uKyvTX/7yF3366afyer169NFHtW7dOt10003xaoIpBw8eVFVVlaqqquT3+7Vnzx5VVVW1\nLjF36/WWjLfd7DW3bKn1+vXrWzdVeTweXXHFFfJ4PCouLtb8+fMlSTU1NTrppJNaX7N//3795je/\nUW1trfr06aPRo0dr+/btGjt2rFVl2cJM28ePH69169bp3nvv1X333achQ4boiSee0FlnnRWXNphR\nVFSk1NRUXXfddWpoaFB+fr4ef/zxduPfbrnmXR0p5fP5VFNT0/r7ffr00ebNm7V48WJNnTpVGRkZ\nWrRokaFFFx0bAAAAgElEQVSlq05gtN0ej0crVqzQnj17lJqaqqFDh6q4uFiXX355nFpgzvvvv69L\nLrlE0tE2FRUVqaioSL/4xS9UXFzs2ustGW+72WvO8ToAANvFfZ8PACD5ED4AANsRPgAA2xE+AADb\nET4AANsRPgAA2xE+AADbET4AANsRPgAA2/0fLBkIMOkkeqkAAAAASUVORK5CYII=\n",
                            "text": [
                                "<matplotlib.figure.Figure at 0x10675f990>"
                            ]
                        }
                    ],
                    "prompt_number": 3
                },
                {
                    "cell_type": "markdown",
                    "metadata": {},
                    "source": [
                        "Build pymc model"
                    ]
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "# Priors\n",
                        "a = pm.Normal('a', mu=0.0, tau=0.0001)\n",
                        "b = pm.Normal('b', mu=0.0, tau=0.0001)\n",
                        "shape = pm.Uniform('shape', lower=0., upper=100.)\n",
                        "\n",
                        "# Model\n",
                        "mu = pm.Lambda('mu', lambda a=a,b=b,shape=shape: np.exp(a+b*x) )\n",
                        "\n",
                        "# Likelihood\n",
                        "Yi = pm.Gamma('Yi', alpha=shape, beta=shape/mu, value=y, observed=True)\n",
                        "\n",
                        "# Predictive model fit\n",
                        "xnew = np.linspace(min(x),max(x),N)\n",
                        "Ex_mu = pm.Lambda('Ex_mu', lambda a=a,b=b: np.exp(a+b*xnew))"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [],
                    "prompt_number": 4
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "M = pm.MCMC(locals())\n",
                        "M.sample(10000,9000)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "output_type": "stream",
                            "stream": "stdout",
                            "text": [
                                "\r",
                                " [--------         21%                  ] 2161 of 10000 complete in 0.5 sec"
                            ]
                        },
                        {
                            "output_type": "stream",
                            "stream": "stdout",
                            "text": [
                                "\r",
                                " [---------------  41%                  ] 4186 of 10000 complete in 1.0 sec"
                            ]
                        },
                        {
                            "output_type": "stream",
                            "stream": "stdout",
                            "text": [
                                "\r",
                                " [-----------------62%---               ] 6223 of 10000 complete in 1.5 sec"
                            ]
                        },
                        {
                            "output_type": "stream",
                            "stream": "stdout",
                            "text": [
                                "\r",
                                " [-----------------82%-----------       ] 8257 of 10000 complete in 2.0 sec"
                            ]
                        },
                        {
                            "output_type": "stream",
                            "stream": "stdout",
                            "text": [
                                "\r",
                                " [-----------------100%-----------------] 10000 of 10000 complete in 2.5 sec"
                            ]
                        }
                    ],
                    "prompt_number": 5
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "plt.hist(M.a.trace())\n",
                        "plt.plot(a_true,1,marker='.',markersize=30)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "metadata": {},
                            "output_type": "pyout",
                            "prompt_number": 6,
                            "text": [
                                "[<matplotlib.lines.Line2D at 0x10e91bd90>]"
                            ]
                        },
                        {
                            "metadata": {},
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEWCAYAAADoyannAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHYZJREFUeJzt3X9wVOXB9vFrWYiBoCQhuwkGoe8bthAkkCI/NBYTIkg7\n8tNBASsdoVqKTEU6BjaKpRZ8My5gqY6JIMRXi1gooCCiMk6ggyABq6ivA08yVFILJJsQgklIAiZ5\n/+Bxny6E/NyT7E2+n5mdkXPuPfd10ty5uns2J7aysrJ6AQBgoC4dHQAAgNaixAAAxqLEAADGosQA\nAMaixAAAxqLEAADGosQAAMZqtMReeOEFjR07Vv369dOAAQM0c+ZMHTt2zG/M/PnzFRER4fe45557\n/MbU1NQoLS1NcXFxio2N1axZs3T69OnAnw0AoFNptMQOHDigRx99VHv27NHOnTvVtWtXTZ06VWVl\nZb4xNptNY8eOVV5enu+xZcsWv+Okp6dr165dys7O1u7du1VeXq4ZM2aorq7OmrMCAHQKtpbcsaOy\nslL9+vXTpk2bNGHCBEmXX4mVlpZq8+bNDT7n/PnzcrlcyszM1PTp0yVJp06dUkJCgrZu3arU1NQA\nnAYAoDNq0TWx8vJy1dXVKTw83LfNZrPp0KFDcrlcGjFihBYuXKiSkhLf/qNHj+rSpUt+ZRUbG6uB\nAwcqNzc3AKcAAOisurZksNvt1tChQzVq1CjftnHjxmny5Mnq37+/CgoKtGLFCk2ePFn79u1TSEiI\nvF6v7Ha7IiMj/Y7lcDhUXFwcmLMAAHRKzS6xp556SocPH9b7778vm83m237ffff5/js+Pl6JiYlK\nSEjQhx9+qEmTJgU2LQAA/6FZbyemp6fr7bff1s6dO9W/f/9Gx8bExOjmm2/WN998I0lyOp2qra1V\naWmp3ziv1yun09nK2AAANKPElixZ4iuwAQMGNHnAkpISnTlzRtHR0ZKkxMREdevWTTk5Ob4xp06d\nUl5enkaPHt2G6ACAzq7REnvyySf11ltvad26dbrppptUVFSkoqIiVVZWSrr8acWlS5fqyJEjKigo\n0P79+zVr1iw5nU5NnDhRktSrVy/Nnj1by5Yt09///nd98cUXmjdvnoYMGaKUlBTLT7A95efnd3SE\nViN7xyB7xyD79aPRa2IbNmyQzWbTlClT/La73W4tWbJEdrtdx44d0+bNm3X+/HlFR0frrrvu0uuv\nv66wsDDf+IyMDNntds2ZM0fV1dVKTk7WunXr/K6tAQDQUo2W2Llz5xp9cmhoqLZt29bkJCEhIfJ4\nPPJ4PC1LBwBAI7h3IgDAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiU\nGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiN/j0xAOY6c6FWhRdqLTt+VahTFSUXG9wX08Ou\nPj3sls0N/IASA65ThRdqtehgmcWz1DS49U9J4ZQY2gVvJwIAjEWJAQCMRYkBAIxFiQEAjEWJAQCM\nRYkBAIxFiQEAjEWJAQCMRYkBAIxFiQEAjEWJAQCMRYkBAIxFiQEAjEWJAQCMRYkBAIxFiQEAjEWJ\nAQCMRYkBAIxFiQEAjEWJAQCM1WiJvfDCCxo7dqz69eunAQMGaObMmTp27NhV4zIyMhQfH68+ffpo\n4sSJOn78uN/+mpoapaWlKS4uTrGxsZo1a5ZOnz4d2DMBAHQ6jZbYgQMH9Oijj2rPnj3auXOnunbt\nqqlTp6qsrMw3Zs2aNcrMzJTH41FOTo4cDoemTZumiooK35j09HTt2rVL2dnZ2r17t8rLyzVjxgzV\n1dVZd2YAgOte18Z2btu2ze/fa9euVb9+/ZSbm6sJEyaovr5eWVlZWrRokSZNmiRJysrKksvl0tat\nW/Xwww/r/Pnz2rhxozIzM5WcnOw7TkJCgvbt26fU1FSLTg0AcL1r0TWx8vJy1dXVKTw8XJJUUFAg\nr9frV0ShoaFKSkpSbm6uJOno0aO6dOmS35jY2FgNHDjQNwYAgNZoUYm53W4NHTpUo0aNkiQVFRVJ\nkhwOh9+4qKgoeb1eSZLX65XdbldkZKTfGIfDoeLi4lYHBwCg0bcT/9NTTz2lw4cP6/3335fNZmty\nfHPGAADQFs0qsfT0dL3zzjt699131b9/f9/26OhoSVJxcbFiY2N924uLi+V0OiVJTqdTtbW1Ki0t\n9Xs15vV6lZSUdM058/PzW3YmQcLU3BLZO4pV2atCnZYct1lzV1UpP7+gw+ZvDr5n2pfL5bLkuE2W\n2JIlS7Rjxw69++67GjBggN++/v37Kzo6Wjk5OUpMTJQkVVdX69ChQ1q+fLkkKTExUd26dVNOTo6m\nT58uSTp16pTy8vI0evToa85r1QlbKT8/38jcEtk7ipXZK0ouSqqx5NhN6d69u1y3BO//JnzPXD8a\nLbEnn3xSW7Zs0caNG3XTTTf5roH17NlTYWFhstlsmj9/vlavXi2Xy6W4uDitWrVKPXv29BVWr169\nNHv2bC1btkwOh0Ph4eF6+umnNWTIEKWkpFh+ggCA61ejJbZhwwbZbDZNmTLFb7vb7daSJUskSQsX\nLlRVVZXS0tJUVlamESNGaPv27QoLC/ONz8jIkN1u15w5c1RdXa3k5GStW7eO62YAgDZptMTOnTvX\nrIO43W653e5r7g8JCZHH45HH42lZOgAAGtHsTycCpjpzoVaFF2o7ZO6YHnb16WHvkLmBzoASw3Wv\n8EKtFh0sa3qgBf6UFE6JARbiLvYAAGNRYgAAY1FiAABjUWIAAGNRYgAAY1FiAABjUWIAAGNRYgAA\nY1FiAABjUWIAAGNRYgAAY1FiAABjUWIAAGNRYgAAY1FiAABjUWIAAGNRYgAAY1FiAABjUWIAAGNR\nYgAAY1FiAABjUWIAAGN17egAwPUspIv0ecnFa+6vCnWqopH9bXGxtt6S4wLBhBIDLFRaU6dnjnzX\nxKgaS+ZePvImS44LBBPeTgQAGIsSAwAYixIDABiLEgMAGIsSAwAYixIDABiLEgMAGIsSAwAYixID\nABiLEgMAGIsSAwAYq8kSO3DggGbOnKnBgwcrIiJCmzZt8ts/f/58RURE+D3uuecevzE1NTVKS0tT\nXFycYmNjNWvWLJ0+fTqwZwIA6HSaLLELFy5oyJAhysjIUPfu3WWz2fz222w2jR07Vnl5eb7Hli1b\n/Makp6dr165dys7O1u7du1VeXq4ZM2aorq4usGcDAOhUmryL/fjx4zV+/HhJ0oIFC67aX19fr27d\nusnhcDT4/PPnz2vjxo3KzMxUcnKyJGnt2rVKSEjQvn37lJqa2pb8AIBOrM3XxGw2mw4dOiSXy6UR\nI0Zo4cKFKikp8e0/evSoLl265FdWsbGxGjhwoHJzc9s6PQCgE2vz3xMbN26cJk+erP79+6ugoEAr\nVqzQ5MmTtW/fPoWEhMjr9cputysyMtLveQ6HQ8XFxW2dHgDQibW5xO677z7ff8fHxysxMVEJCQn6\n8MMPNWnSpFYfNz8/v63ROoSpuaXrN3tVqLMdk/jryOu+HTl3VVWV8vMLOmz+5rhev9+DlcvlsuS4\nAf/LzjExMbr55pv1zTffSJKcTqdqa2tVWlrq92rM6/UqKSnpmsex6oStlJ+fb2Ru6frOXlFyUVb9\n9eSmdOnScb/F0pFzd+/eXa5bgvf76Xr+fu9sAv5dXlJSojNnzig6OlqSlJiYqG7duiknJ8c35tSp\nU8rLy9Po0aMDPT0AoBNp8pVYZWWlTpw4Ieny2xPffvutvvzyS0VGRioiIkIZGRmaMmWKnE6n/vWv\nf+mPf/yjnE6nJk6cKEnq1auXZs+erWXLlsnhcCg8PFxPP/20hgwZopSUFEtPDgBwfWuyxD777DNN\nnjxZ0uVPImZkZCgjI0MPPvigVq9erWPHjmnz5s06f/68oqOjddddd+n1119XWFiY7xgZGRmy2+2a\nM2eOqqurlZycrHXr1l31O2cAALREkyU2ZswYnTt37pr7t23b1uQkISEh8ng88ng8LUsHAEAjuHci\nAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWAG/7RRwpTMXalV4oday41eFOv/71lINu1hbb9nc\nADoWJQbLFV6o1aKDZRbPcu17Iy4feZPFcwPoKLydCAAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYA\nMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAW\nJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUG\nADBWkyV24MABzZw5U4MHD1ZERIQ2bdp01ZiMjAzFx8erT58+mjhxoo4fP+63v6amRmlpaYqLi1Ns\nbKxmzZql06dPB+4sAACdUpMlduHCBQ0ZMkQZGRnq3r27bDab3/41a9YoMzNTHo9HOTk5cjgcmjZt\nmioqKnxj0tPTtWvXLmVnZ2v37t0qLy/XjBkzVFdXF/gzAgB0Gk2W2Pjx47V06VJNmTJFXbr4D6+v\nr1dWVpYWLVqkSZMmKT4+XllZWaqoqNDWrVslSefPn9fGjRu1fPlyJScna9iwYVq7dq2+/vpr7du3\nz5KTAgB0Dm26JlZQUCCv16vU1FTfttDQUCUlJSk3N1eSdPToUV26dMlvTGxsrAYOHOgbAwBAa7Sp\nxIqKiiRJDofDb3tUVJS8Xq8kyev1ym63KzIy0m+Mw+FQcXFxW6YHAHRyXa068JXXzloqPz8/QEna\nl6m5JeuyV4U6LTluc3XktdfOOndVVZXy8ws6bP7mYK22L5fLZclx21Ri0dHRkqTi4mLFxsb6thcX\nF8vpvPyDy+l0qra2VqWlpX6vxrxer5KSkq55bKtO2Er5+flG5paszV5RclFSjSXHbo4rr+Uyt/W6\nd+8u1y3BuxZYq9ePNn2X9+/fX9HR0crJyfFtq66u1qFDhzR69GhJUmJiorp16+Y35tSpU8rLy/ON\nAQCgNZp8JVZZWakTJ05Iuvz2xLfffqsvv/xSkZGR6tu3r+bPn6/Vq1fL5XIpLi5Oq1atUs+ePTV9\n+nRJUq9evTR79mwtW7ZMDodD4eHhevrppzVkyBClpKRYenIAgOtbkyX22WefafLkyZIuX+fKyMhQ\nRkaGHnzwQb388stauHChqqqqlJaWprKyMo0YMULbt29XWFiY7xgZGRmy2+2aM2eOqqurlZycrHXr\n1rX5uhkAoHNrssTGjBmjc+fONTrG7XbL7XZfc39ISIg8Ho88Hk/LEwIAcA3cOxEAYCxKDABgLEoM\nAGAsSgwAYCxKDABgLEoMAGAsSgwAYCxKDABgLEoMAGAsSgwAYCxKDABgLEoMAGAsSgwAYCxKDABg\nLEoMAGAsSgwAYCxKDABgLEoMAGAsSgwAYKyuHR2gMzlzoVaFF2o7ZO6YHnb16WHvkLkBwCqUWDsq\nvFCrRQfLOmTuPyWFU2IArju8nQgAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWH7EHEHAhXaTP\nSy52yNz8TmTnQokBCLjSmjo9c+S7Dpmb34nsXHg7EQBgLEoMAGAs3k7sJJq6RlEV6lSFRdcwLtbW\nW3JcAKDEOonmXaOosWTu5SNvsuS4AMDbiQAAY1FiAABjUWIAAGO1ucQyMjIUERHh9xg0aNBVY+Lj\n49WnTx9NnDhRx48fb+u0AAAE5pXYj3/8Y+Xl5fkeBw8e9O1bs2aNMjMz5fF4lJOTI4fDoWnTpqmi\noiIQUwMAOrGAlJjdbpfD4fA9IiMjJUn19fXKysrSokWLNGnSJMXHxysrK0sVFRXaunVrIKYGAHRi\nASmxkydPKj4+XsOGDdOvfvUrnTx5UpJUUFAgr9er1NRU39jQ0FAlJSUpNzc3EFMDADqxNpfYyJEj\nlZWVpW3btunFF19UUVGRJkyYoHPnzqmoqEiS5HA4/J4TFRUlr9fb1qkBAJ1cm3/Zedy4cX7/Hjly\npIYNG6ZNmzZpxIgR13yezWZr69QAgE4u4Hfs6NGjhwYNGqRvvvlG9957rySpuLhYsbGxvjHFxcVy\nOp2NHic/Pz/Q0dpFY7mrQhs/ZyvV1dV1yrk7en7mbn9VVVXKzy9ocpypP2MkM7O7XC5LjhvwEquu\nrlZeXp7uuusu/ehHP1J0dLRycnKUmJjo23/o0CEtX7680eNYdcJWys/PbzT35XsTWnNrp6Z06dJx\nvxLYkXN39PzM3f66d+8u1y2N//xoaq0GM5OzW6HNJbZ06VL9/Oc/V2xsrEpKSrRy5UpVVVVp1qxZ\nkqT58+dr9erVcrlciouL06pVq9SzZ09Nnz69zeEBAJ1bm0vszJkzeuSRR3T27FlFRUVp5MiR+uij\nj9S3b19J0sKFC1VVVaW0tDSVlZVpxIgR2r59u8LCwtocHgDQubW5xDZs2NDkGLfbLbfb3dapAADw\nw70TAQDGosQAAMaixAAAxqLEAADGosQAAMaixAAAxqLEAADGCvhtpwCgI4V0kT4vudjomKpQ53/f\nBi6wYnrY1aeHPeDHxbVRYgCuK6U1dXrmyHfNGBn4+5j+KSmcEmtnvJ0IADAWJQYAMBYlBgAwFiUG\nADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAwFiUGADAWJQYAMBYlBgAw\nFiUGADAWJQYAMBYlBgAwFiUGADBW144O0N7+X+lFVX1vzbGrb3CqzHuxwX12m1RXX2/NxACCQkgX\n6fOShn8GBEpVqFMV15gjpoddfXrYLZ0/2HS6Evu//1WpT4svWThDTYNbe4XY9PTwmyycF0BHK62p\n0zNHvmuHmRr+OfOnpPBOV2K8nQgAMBYlBgAwFiUGADAWJQYAMBYlBgAwVqf7dCKA1ulV851mnPhA\ng879U2HfV/m2V3btruMR/1ub436m8zfwCVy0L0oMQJMGnTuhFYdf1I8qzjS4/67CzzT+24N6ZtRv\ndTwirp3ToTNr17cT169fr6FDhyomJkYpKSn65JNP2nN6AK0w6NwJ/eHIy9cssB/8r4rT+sORl5Vw\n9r/aKRnQjiW2fft2paen68knn9T+/fs1atQo3X///fr3v//dXhEAtFCvmu+04vCLGlDevHU6oPzf\n+j+5axR1/LDFyYDL2q3EXn75Zf3iF7/QL3/5S7lcLnk8HkVHRys7O7u9IgBooRknPmjyFdiV+lSd\n1aDd6xX6fbVFqYD/0S4ldvHiRX3xxRcaO3as3/bU1FTl5ua2RwQArTDo3D9b9bye3n9pVv57AU4D\nXK1dSuzs2bOqra2V0+n02x4VFSWv19seEQC0wn9+CrElbJKGleYFNgzQAFtZWZnlt1Y/c+aMBg8e\nrN27d+uOO+7wbX/++ee1detWHTlyxOoIAIDrULu8Euvdu7fsdvtVr7qKi4sVHR3dHhEAANehdimx\nkJAQJSYmau/evX7b9+7dq9GjR7dHBADAdajdftl5wYIFmjdvnoYPH67Ro0crOztbXq9Xc+bMaa8I\nAIDrTLuV2LRp01RaWqpVq1apqKhIgwcP1pYtW9S3b9/2igAAuM60ywc7AACwgqXXxFp7m6kTJ06o\nb9++Db5K+/jjj5WcnKyYmBglJibqtddeC3RsSYHPvnPnTk2bNk0DBgzQLbfconHjxun999+3Irol\nX/cffPLJJ+rdu7eSkpICFdePFdkvXryo5557TsOGDVN0dLSGDBmitWvXBjq6Jdn/+te/6s4779TN\nN9+sgQMH6te//rUlv5bSkuwFBQWKiIi46pGTk+M3LhjXanOyB+tabe7X/QfBtFabm701a9WyEmvt\nbaYuXryouXPn6s4775TNZvPbd/LkST3wwAO6/fbbtX//fv3ud7/T4sWLtXPnzqDPfvDgQaWkpOhv\nf/ub9u/fr/Hjx+uhhx4K+P0jrcj+g7KyMv3mN79RSkpKQDP/wKrsc+fO1d69e/Xiiy/q008/1Rtv\nvKFbb7016LN//PHHeuyxx/TQQw/p0KFDevPNN5WXl6dHH300KLJv375deXl5vseYMWN8+4J9rTaW\nPdjXamPZfxCsa7Wp7K1Zq5a9nXj33XcrISFBa9as8W277bbbNGXKFP3+97+/5vPS09NVXl6upKQk\nLV682O+LsmzZMr333nv69NNPfdsef/xxHT9+XHv27Anq7Nea54477tCKFSuMyP7QQw9p6NChqqur\n086dO3Xw4MGA5bYqe05Ojh5++GF98cUXioiICGheq7O/9NJLWrdunb766ivfto0bN8rtdgf0nqMt\nzV5QUOD7tHFiYmKDxwzWtdqc7Neap6PXakuyB9tabU721q5VS16JtfY2Ux9++KH27Nkjj8ej+vqr\nu/Xw4cMNHvPzzz9XbW1tUGdvSHl5eUB/sFqZff369Tp79qzS0tKafX4tYVX29957T8OHD9dLL72k\nW2+9VbfddpuWLFmiysrKoM+ekpKis2fP6oMPPlB9fb3Onj2r7du365577unw7NLlH5Qul0s/+9nP\ntGPHDr99wbxWm8rekGBZq1LT2YN1rUqNZ2/tWrXk04mtuc3UmTNn9MQTT+jNN99Ujx49GhxTXFx8\n1TEdDoe+//57nT179qp9wZT9Sq+++qoKCws1Y8aMNmf+gVXZv/76a3k8Hn300UfXfKuxrazKfvLk\nSR06dEg33HCD/vKXv6isrEyLFy9WYWGhXn/99aDOnpCQoLVr1+qRRx5RTU2Nvv/+e40dO1aZmZkB\nyd3a7DfeeKNWrFih22+/XXa7Xbt379bcuXOVlZWlBx54QFLwrtXmZL9SsKzV5mQP1rXanOytXatB\n80cx582bp7lz52r48OEdHaXFWpp9x44dWrZsmV577bUO/xWDprLX1NRo7ty5Wr58ufr169fO6RrX\nnK97XV2dunTpovXr1+vGG2+UJK1cuVL33XefSkpKFBUV1V5x/TQn+5EjR7RgwQItWbJEd999twoL\nC/XMM8/oiSee0CuvvNKOaf1FRkZqwYIFvn8nJibq3Llz+vOf/3zNIggWLc0eTGu1qezBvFab83Vv\n7Vq15O3E1txmav/+/Xr++ecVFRWlqKgoPf7446qsrFRUVJTeeOMNSZLT6WzwmF27dlXv3r2DOvsP\nduzYofnz5+uVV17RhAkTApLZyuxFRUXKy8vTggULfGNWrlypY8eOKSoqSvv27Qva7JIUHR2tmJgY\n36KQJJfLJUkBu65kVfbMzEwlJyfrt7/9rQYPHqzU1FStXr1amzdv1pkzLfvzKIHM3pCf/OQn+uc/\n/+eO98G6VhtyZfYfBNtabch/Zi8sLAzatdpUdqn1a9WSV2L/eZupKVOm+Lbv3btXU6dObfA5V37y\n57333tPq1auVk5OjmJgYSdKoUaO0a9cuv3F79+7V8OHDZbfbgzq7JL399tt67LHH9Morr2jy5MkB\nyWt19p49e141Zv369dq7d6/efPNN3XLLLUGbXZLuuOMO7dy5U5WVlQoLC5N0+SPtkoI+e319vbp0\n8f//mT/8u66ursOyN+Srr77y+14P1rXakCuzS8G5Vhvyn9ljY2ODdq02lV1q/Vq17O3Epm4z9eyz\nz+qzzz7zXdwbNGiQ3/P/8Y9/qEuXLn7b58yZo1dffVXp6el6+OGHlZubq7feeksbNmwI+uzbtm3T\nvHnz9Nxzz+n2229XUVGRpMvfEIG8YGxF9ivH9O7dWzfccMNV24Mx+/Tp07Vy5UotWLBAbrdbZWVl\ncrvdmjp1asBeEViV/d5779Vjjz2m7OxspaamqrCwUOnp6UpMTFRsbGyHZd+0aZNCQkKUkJCgLl26\n6IMPPtCGDRv07LPP+o4ZrGu1OdmDda02lb1r165Bu1ab83Vv7Vq1rMSaus1UUVGRTp482egxrrww\n2b9/f23ZskVPPfWUsrOz1adPH3k8Hk2aNCnos7/22muqq6uT2+2W2+32bf/pT3+qd999N6izt3R/\na1mRPSwsTO+8844WL16s1NRUhYeH695779Uf/vCHoM9+//3367vvvtOrr76qpUuXqlevXhozZozf\nwu+I7DabTatWrdK3334ru92uAQMG6OWXX9b999/vGxOsa7U52YN1rTYn+5WCZa02J3tr1yq3nQIA\nGKtd/hQLAABWoMQAAMaixAAAxqLEAADGosQAAMaixAAAxqLEAADGosQAAMaixAAAxvr/H0sQbAsw\n14MAAAAASUVORK5CYII=\n",
                            "text": [
                                "<matplotlib.figure.Figure at 0x10e91bb90>"
                            ]
                        }
                    ],
                    "prompt_number": 6
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "plt.hist(M.b.trace())\n",
                        "plt.plot(b_true,1,marker='.',markersize=30)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "metadata": {},
                            "output_type": "pyout",
                            "prompt_number": 7,
                            "text": [
                                "[<matplotlib.lines.Line2D at 0x10ecbafd0>]"
                            ]
                        },
                        {
                            "metadata": {},
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEWCAYAAADoyannAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHLdJREFUeJzt3X9wVPX97/HXshAS4AqJ2U1CkOgNKwkTIEV+jLEayIDe\ntiGITYXo0DG2lga9RSwpCVppjXN32IJSpiQVBcdKUSnQihTFaQNtGiTQIspYaLaokULMJoSEJOQH\nJrl/+GWnK5CQH5vNhzwfM5kh53z2nPfbT3ZfnnN2z1pqamraBQCAgQYFugAAALqLEAMAGIsQAwAY\nixADABiLEAMAGIsQAwAYixADABirwxB77rnnNGvWLI0dO1bjxo3TwoULdfz4cZ8xWVlZCg0N9fm5\n++67fcY0NzcrOztbsbGxio6OVkZGhs6cOdP73QAABpQOQ6y4uFiPPPKI3n33Xe3atUuDBw/Wvffe\nq5qaGu8Yi8WiWbNmqbS01Puzbds2n+3k5uZq9+7d2rx5s/bs2aO6ujotWLBAbW1t/ukKADAgWLpy\nx46GhgaNHTtWW7du1T333CPpyyOx6upqvfHGG1d8TG1trRwOh/Lz85Weni5JOn36tCZOnKjt27cr\nJSWlF9oAAAxEXbomVldXp7a2No0aNcq7zGKx6ODBg3I4HJo6daqWLl2qqqoq7/qjR4/q4sWLPmEV\nHR2t8ePHq6SkpBdaAAAMVIO7MjgnJ0eTJk3S9OnTvctmz56ttLQ0xcTEqKysTM8++6zS0tK0f/9+\nBQUFyePxyGq1KiwszGdbNptNlZWVvdMFAGBAuuYQW7lypQ4dOqS3335bFovFu/y+++7z/js+Pl6J\niYmaOHGi9u7dq7lz5/ZutQAA/JdrOp2Ym5ur3//+99q1a5diYmI6HBsZGanRo0frk08+kSTZ7Xa1\ntraqurraZ5zH45Hdbu9m2QAAXEOIrVixwhtg48aN63SDVVVVKi8vV0REhCQpMTFRQ4YMUWFhoXfM\n6dOnVVpaqhkzZvSgdADAQNdhiC1fvlyvvfaaNm7cqBtuuEEVFRWqqKhQQ0ODpC/frfjUU0/p8OHD\nKisrU1FRkTIyMmS325WamipJGjlypBYtWqRVq1bpL3/5iz744AMtXrxYCQkJmjlzpt8b7C/cbneg\nS+h112NPEn2Z5HrsSbp++/KHDq+Jbdq0SRaLRfPmzfNZnpOToxUrVshqter48eN64403VFtbq4iI\nCN1111165ZVXNHz4cO94p9Mpq9WqzMxMNTU1KTk5WRs3bvS5tgYAQFd1GGLnzp3r8MHBwcHasWNH\npzsJCgqSy+WSy+XqWnUAAHSAeycCAIxFiAEAjEWIAQCMRYgBAIxFiAEAjEWIAQCMRYgBAIzVpbvY\nA+i58gut+vxCa5/vtzHYrvqqFp9lkcOsihpm7fNagN5CiAF97PMLrVp2oKbzgX7R7PPb80mjCDEY\njRADBrCgQdL7Xzk6CwSOCNFdhBgwgFU3t+mnh88HugyOCNFtvLEDAGAsQgwAYCxCDABgLEIMAGAs\nQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIM\nAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYKwOQ+y5\n557TrFmzNHbsWI0bN04LFy7U8ePHLxvndDoVHx+vqKgopaam6sSJEz7rm5ublZ2drdjYWEVHRysj\nI0Nnzpzp3U4AAANOhyFWXFysRx55RO+++6527dqlwYMH695771VNTY13zLp165Sfny+Xy6XCwkLZ\nbDbNnz9f9fX13jG5ubnavXu3Nm/erD179qiurk4LFixQW1ub/zoDAFz3Bne0cseOHT6/v/DCCxo7\ndqxKSkp0zz33qL29XQUFBVq2bJnmzp0rSSooKJDD4dD27dv10EMPqba2Vlu2bFF+fr6Sk5O925k4\ncaL279+vlJQUP7UGALjedemaWF1dndra2jRq1ChJUllZmTwej08QBQcHKykpSSUlJZKko0eP6uLF\niz5joqOjNX78eO8YAAC6o0shlpOTo0mTJmn69OmSpIqKCkmSzWbzGRceHi6PxyNJ8ng8slqtCgsL\n8xljs9lUWVnZ7cIBAOjwdOJ/W7lypQ4dOqS3335bFoul0/HXMgYAgJ64phDLzc3VH/7wB7311luK\niYnxLo+IiJAkVVZWKjo62ru8srJSdrtdkmS329Xa2qrq6mqfozGPx6OkpKSr7tPtdnetEwPQkzn8\n2VdjsN1v2+6q/vLmqsbGRrndZd16LH+DZnA4HH7ZbqchtmLFCr355pt66623NG7cOJ91MTExioiI\nUGFhoRITEyVJTU1NOnjwoPLy8iRJiYmJGjJkiAoLC5Weni5JOn36tEpLSzVjxoyr7tdfDQeK2+2m\nJ0P4u6/6qhZJzX7bflcMGtQ/PioaEhIix01d/2/O3yA6DLHly5dr27Zt2rJli2644QbvNbARI0Zo\n+PDhslgsysrK0tq1a+VwOBQbG6s1a9ZoxIgR3sAaOXKkFi1apFWrVslms2nUqFF68sknlZCQoJkz\nZ/q9QQDA9avDENu0aZMsFovmzZvnszwnJ0crVqyQJC1dulSNjY3Kzs5WTU2Npk6dqp07d2r48OHe\n8U6nU1arVZmZmWpqalJycrI2btzIdTMAQI90GGLnzp27po3k5OQoJyfnquuDgoLkcrnkcrm6Vh0A\nAB3oHyfEAQDoBkIMAGAsQgwAYCxCDABgrGu+YwdgsvILrfr8Qus1jW0Mtv/PZ7n8o6W13W/bBgYa\nQgwDwucXWrXsQE3nA73892HkvGk3+G3bwEDD6UQAgLEIMQCAsQgxAICxCDEAgLEIMQCAsQgxAICx\nCDEAgLEIMQCAsfiwM4CACxokvd+Nu6T44+4qkcOsihpm7dVtwn8IMQABV93cpp8ePt/NR/fu3VWe\nTxpFiBmE04kAAGMRYgAAYxFiAABjEWIAAGMRYgAAYxFiAABjEWIAAGMRYgAAYxFiAABjEWIAAGMR\nYgAAYxFiAABjEWIAAGMRYgAAYxFiAABjEWIAAGMRYgAAYxFiAABjEWIAAGMRYgAAYxFiAABjdRpi\nxcXFWrhwoSZMmKDQ0FBt3brVZ31WVpZCQ0N9fu6++26fMc3NzcrOzlZsbKyio6OVkZGhM2fO9G4n\nAIABp9MQu3DhghISEuR0OhUSEiKLxeKz3mKxaNasWSotLfX+bNu2zWdMbm6udu/erc2bN2vPnj2q\nq6vTggUL1NbW1rvdAAAGlMGdDZgzZ47mzJkjSXr00UcvW9/e3q4hQ4bIZrNd8fG1tbXasmWL8vPz\nlZycLEl64YUXNHHiRO3fv18pKSk9qR8AMID1+JqYxWLRwYMH5XA4NHXqVC1dulRVVVXe9UePHtXF\nixd9wio6Olrjx49XSUlJT3cPABjAOj0S68zs2bOVlpammJgYlZWV6dlnn1VaWpr279+voKAgeTwe\nWa1WhYWF+TzOZrOpsrKyp7sHAAxgPQ6x++67z/vv+Ph4JSYmauLEidq7d6/mzp3b7e263e6eltbv\n0FPgNAbbA12CV3+6FtxfaukvdUhSY2Oj3O6yQJdhzHPrWjkcDr9st8ch9lWRkZEaPXq0PvnkE0mS\n3W5Xa2urqqurfY7GPB6PkpKSrrodfzUcKG63m54CqL6qRVJzoMuQJA0a1H8+2dJfaukvdUhSSEiI\nHDcF9u/apOdWoPX6X05VVZXKy8sVEREhSUpMTNSQIUNUWFjoHXP69GmVlpZqxowZvb17AMAA0umR\nWENDg06ePCnpy0P+U6dO6cMPP1RYWJhCQ0PldDo1b9482e12ffbZZ3rmmWdkt9uVmpoqSRo5cqQW\nLVqkVatWyWazadSoUXryySeVkJCgmTNn+rU5AMD1rdMQO3LkiNLS0iR9+U5Ep9Mpp9OpBx54QGvX\nrtXx48f1xhtvqLa2VhEREbrrrrv0yiuvaPjw4d5tOJ1OWa1WZWZmqqmpScnJydq4ceNlnzkDAKAr\nOg2xO++8U+fOnbvq+h07dnS6k6CgILlcLrlcrq5VBwBAB/rP1VQAALqIEAMAGIsQAwAYixADABiL\nEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixAD\nABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAY\nixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABiLEAMAGIsQAwAYixADABir0xArLi7W\nwoULNWHCBIWGhmrr1q2XjXE6nYqPj1dUVJRSU1N14sQJn/XNzc3Kzs5WbGysoqOjlZGRoTNnzvRe\nFwCAAanTELtw4YISEhLkdDoVEhIii8Xis37dunXKz8+Xy+VSYWGhbDab5s+fr/r6eu+Y3Nxc7d69\nW5s3b9aePXtUV1enBQsWqK2trfc7AgAMGJ2G2Jw5c/TUU09p3rx5GjTId3h7e7sKCgq0bNkyzZ07\nV/Hx8SooKFB9fb22b98uSaqtrdWWLVuUl5en5ORkTZ48WS+88II++ugj7d+/3y9NAQAGhh5dEysr\nK5PH41FKSop3WXBwsJKSklRSUiJJOnr0qC5evOgzJjo6WuPHj/eOAQCgO3oUYhUVFZIkm83mszw8\nPFwej0eS5PF4ZLVaFRYW5jPGZrOpsrKyJ7sHAAxwg/214a9eO+sqt9vdS5X0H/QUOI3B9kCX4NWf\nrgX3l1r6Sx2S1NjYKLe7LNBlGPPculYOh8Mv2+1RiEVEREiSKisrFR0d7V1eWVkpu/3LFw273a7W\n1lZVV1f7HI15PB4lJSVdddv+ajhQ3G43PQVQfVWLpOZAlyFJl11bDqT+Ukt/qUOSQkJC5LgpsH/X\nJj23Aq1HfzkxMTGKiIhQYWGhd1lTU5MOHjyoGTNmSJISExM1ZMgQnzGnT59WaWmpdwwAAN3R6ZFY\nQ0ODTp48KenLQ/5Tp07pww8/VFhYmMaMGaOsrCytXbtWDodDsbGxWrNmjUaMGKH09HRJ0siRI7Vo\n0SKtWrVKNptNo0aN0pNPPqmEhATNnDnTr80BAK5vnYbYkSNHlJaWJunL61xOp1NOp1MPPPCANmzY\noKVLl6qxsVHZ2dmqqanR1KlTtXPnTg0fPty7DafTKavVqszMTDU1NSk5OVkbN27s8XUzAMDA1mmI\n3XnnnTp37lyHY3JycpSTk3PV9UFBQXK5XHK5XF2vEACAq+g/V1MBAOgiQgwAYCxCDABgLEIMAGAs\nQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIM\nAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABg\nLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgLEIMAGAsQgwAYCxCDABgrB6HmNPpVGhoqM9P\nXFzcZWPi4+MVFRWl1NRUnThxoqe7BQCgd47Ebr31VpWWlnp/Dhw44F23bt065efny+VyqbCwUDab\nTfPnz1d9fX1v7BoAMID1SohZrVbZbDbvT1hYmCSpvb1dBQUFWrZsmebOnav4+HgVFBSovr5e27dv\n741dAwAGsF4JsU8//VTx8fGaPHmyvve97+nTTz+VJJWVlcnj8SglJcU7Njg4WElJSSopKemNXQMA\nBrAeh9i0adNUUFCgHTt2aP369aqoqNA999yjc+fOqaKiQpJks9l8HhMeHi6Px9PTXQMABrjBPd3A\n7NmzfX6fNm2aJk+erK1bt2rq1KlXfZzFYunprgEAA1yPQ+yrhg0bpri4OH3yySf61re+JUmqrKxU\ndHS0d0xlZaXsdnuH23G73b1dWsANxJ7a/1e4qi72g09yWIcEugKvtra2QJfg1V9q6S91SFJjY6Pc\n7rJAl3HdvV44HA6/bLfXQ6ypqUmlpaW66667dPPNNysiIkKFhYVKTEz0rj948KDy8vI63I6/Gg4U\nt9s9IHt6v6pFKw/V9FFFV5c3bWigS/AaNKgfhPr/6C+19Jc6JCkkJESOmwL7XL0eXy/8pcch9tRT\nT+kb3/iGoqOjVVVVpV/84hdqbGxURkaGJCkrK0tr166Vw+FQbGys1qxZoxEjRig9Pb3HxQMABrYe\nh1h5ebm+//3v6+zZswoPD9e0adP0pz/9SWPGjJEkLV26VI2NjcrOzlZNTY2mTp2qnTt3avjw4T0u\nHgAwsPU4xDZt2tTpmJycHOXk5PR0VwAA+Og/J6IBAOgiQgwAYCxCDABgLEIMAGAsQgwAYCxCDABg\nrF6/YwcAoHeUX2jV5xdaA12GIodZFTXMGugyrogQA4B+6vMLrVp2IPC3bXs+aVS/DTFOJwIAjEWI\nAQCMxenE61BfnUdvDLarvqqlwzEtre1+rwPAwEWIXYf69jx6c4dr86bd0Ed1ABiIOJ0IADAWIQYA\nMBYhBgAwFiEGADAWIQYAMBYhBgAwFiEGADAWIQYAMBYhBgAwFiEGADAWIQYAMBYhBgAwFiEGADAW\nIQYAMBYhBgAwFiEGADAWX4rZSzr7NuVr+Rbk3sK3KQMYKAixXnJt36bc8bcg9xa+TRnAQMHpRACA\nsQgxAICxjD6dWNfSpiNVLWq++qWoPjMyyBLoEgBgwDE6xC62tWv9sXqdbW4LdCl6lutQANDnOJ0I\nADCW0UdiQG8Y2XxeC06+o7hzH2v4F43e5Q2DQ3Qi9H/rjdj/o9qhHGkPFEGDpPf76OMwV3PpIzl8\nXKZzhBgGtLhzJ/XsofW6ub78iuvv+vyI5pw6oJ9O/786ERrbx9UhEKqb2/TTw+cDXYakZj4ucw36\n9HTiSy+9pEmTJikyMlIzZ87Ue++915e7B3zEnTupnx3ecNUAu+SW+jP62eENmnj2X31UGYBr1Wch\ntnPnTuXm5mr58uUqKirS9OnT9Z3vfEf/+c9/+qoEwGtk83k9e2i9xtVd29/fuLr/6P+VrFNS+RE/\nVwagK/osxDZs2KAHH3xQ3/3ud+VwOORyuRQREaHNmzf3VQmA14KT73R6BPZVUY1n9fixVxX8RZOf\nqgLQVX0SYi0tLfrggw80a9Ysn+UpKSkqKSnpixIAH3HnPu7W426uO60M9x97uRoA3dUnIXb27Fm1\ntrbKbrf7LA8PD5fH4+mLEgAf//0uxK4YJGlydWnvFgOg2yw1NTV+fw9neXm5JkyYoD179uj222/3\nLl+9erW2b9+uw4cP+7sEAMB1qE+OxG688UZZrdbLjroqKysVERHRFyUAAK5DfRJiQUFBSkxM1L59\n+3yW79u3TzNmzOiLEgAA16E++7Dzo48+qsWLF2vKlCmaMWOGNm/eLI/Ho8zMzL4qAQBwnemzEJs/\nf76qq6u1Zs0aVVRUaMKECdq2bZvGjBnTVyUAAK4zffLGDgAA/MFv18SKi4u1cOFCTZgwQaGhodq6\ndWunj/noo4/0zW9+U1FRUZowYYJcLpfP+qKiIoWGhl728+9//9tfbfjoak/Nzc3KysrSHXfcIZvN\nptTU1CuO+9vf/qbk5GRFRkYqMTFRL7/8sj/Kvyp/9BXouZK63ldRUZEyMjIUFxen0aNH64477tCW\nLVsuGxfI+fJHTybO1YkTJ5Samqpbb73VOw95eXm6ePGizzjTnlvX0leg56s7r+2XnDx5UmPGjLni\nGbjuzpXfTideuHBBCQkJysjIUFZWliyWjr808vz585o/f76+/vWva9++ffrXv/6lxx57TMOGDdNj\njz3mM7akpEShoaHe32+88Ua/9PBVXe2ptbVVISEhWrx4sfbu3avz5y+/qeinn36q+++/X4sWLdJL\nL72k9957Tz/+8Y914403Ki0tzV+t+PBHX5cEaq6krvd1+PBhJSQkaNmyZYqIiNCf//xnPf744woO\nDlZ6erqkwM+XP3q6xKS5Gjp0qB588EFNmjRJI0eO1LFjx7R06VK1tLQoLy9PUuDnSvJPX5eY8jp4\nSUtLix5++GHdcccdOnDggM+6nsyV30Jszpw5mjNnjqQv39TRmd/97ndqampSQUGBhg4dqri4OLnd\nbuXn518WYuHh4QoLC/NL3R3pak/Dhg3Tc889J0k6duyYamtrLxvz8ssva/To0Vq9erUkyeFw6O9/\n/7t+9atf9dkTzR99XRKouZK63tcTTzzh8/vDDz+soqIi7dq1y/uCH+j58kdPl5g0V7fccotuueUW\n7+9jxoxRUVGRzx2AAj1Xkn/6usSU18FLVq1apYkTJyopKUnFxcU+63oyV/3mSzEPHTqk22+/XUOH\nDvUuS0lJUXl5uT777DOfsTNnzlRcXJzmzZunoqKivi61Vx06dOiKt+N6//331draGqCqeo/pc3X+\n/Hmf/9u9Hubrqz1dYvJcffzxxyosLPSZm+thrq7U1yUmzdfevXv17rvvyuVyqb398rdh9GSu+k2I\neTyey25LZbPZvOskKSoqSs8//7xeffVVvfrqq3I4HJo3b57RX+lSWVl5xb6/+OILnT17NkBV9dz1\nMFfvvPOO/vrXv+qhhx7yLjN9vq7Uk8lzdffddysyMlK33XabbrvtNuXm5nrXmTxXHfVl2nyVl5fr\n8ccf14svvqhhw4ZdcUxP5qrffCnmtZxXHTdunMaNG+f9fdq0afrss8+0fv16n9tZIfBMn6uDBw/q\nBz/4gVwul772ta8FupxecbWeTJ6rl19+WQ0NDTp27JiefvppPf3003rmmWcCXVaPddSXafO1ePFi\nPfzww5oyZYpftt9vjsTsdvsVb0t1ad3VTJkyRR9/3L07kvcHV+t78ODBfXphvS+YMlfvvfee7r//\nfq1cufKyD+ObOl8d9XQlpsxVdHS0br31Vn3729/WqlWr9Otf/9p7+snUuZI67utK+vN8FRUVafXq\n1QoPD1d4eLh+9KMfqaGhQeHh4frNb34jqWdz1W+OxKZPn66f/exnam5u9l4X27dvn0aPHq2xY8de\n9XHHjh1TZGRkX5XZ66ZPn67du3f7LNu3b5+mTJkiq9UaoKr8w4S5uvT24dzcXP3whz+8bL2J89VZ\nT1diwlx9VWtrq9ra2tTW1iar1WrkXF3JV/u6kv48X189zfnHP/5Ra9euVWFhobfmnsyV30KsoaFB\nJ0+elCS1tbXp1KlT+vDDDxUWFqYxY8bo5z//uY4cOaI333xTkpSenq7Vq1dryZIlWr58udxut375\ny19qxYoV3m3m5+crJiZGcXFxamlp0bZt27Rnzx69+uqr/mqjRz1JX37uo6WlRWfPnvWeHmhvb9ek\nSZMkSZmZmXrxxReVm5urhx56SCUlJXrttde0adOmPunJX30Feq6601dRUZEWLFigRx55ROnp6aqo\nqJAkWa1WhYeHSwr8fPmjJxPn6vXXX1dISIji4+MVFBSk999/X3l5ebrvvvs0ZMgQSYGfK3/1Fej5\n6mpPcXFxPo//xz/+oUGDBvks78lc+e2OHUVFRd63RlosFu87Uh544AFt2LBBS5YsUXFxsT744APv\nY/75z39q+fLlOnLkiEJDQ5WZmamf/OQn3vXr16/XK6+8ojNnzig4OFjx8fF64oknNHv2bH+00Cs9\nTZo0SadOnfJ5jMViUXV1tXdMcXGxVq5cqRMnTigqKkqPP/64z4V3E/sK9Fx1p68lS5bo9ddfv+zd\nU2PHjvXpPZDz5Y+eTJyrHTt2aP369fr444/V3t6um266Sffff7+WLFni8w5n055b19JXoOerO68X\n/+23v/2tcnJyvK8fl3R3rrjtFADAWP3mjR0AAHQVIQYAMBYhBgAwFiEGADAWIQYAMBYhBgAwFiEG\nADAWIQYAMBYhBgAw1v8HTjHHHfAGCscAAAAASUVORK5CYII=\n",
                            "text": [
                                "<matplotlib.figure.Figure at 0x10ecbaf90>"
                            ]
                        }
                    ],
                    "prompt_number": 7
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "plt.hist(M.shape.trace())\n",
                        "plt.plot(shape_true,1,marker='.',markersize=30)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "metadata": {},
                            "output_type": "pyout",
                            "prompt_number": 8,
                            "text": [
                                "[<matplotlib.lines.Line2D at 0x10ef2b410>]"
                            ]
                        },
                        {
                            "metadata": {},
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEWCAYAAAA3h9P4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1UVHXiP/D3OMKCkDzEzECjwP5gAlwQtlBcOolyfFhP\nimHsonas2HzI3N+6emSF1MXU/U47QbHbiqlJuZmVq24+pOWpiVQU7KwZ+yuJyTU0AWcQh0AGhIHf\nH36bbTJnwBm4H8b365w5R+793Hvfdxx4M/feucjMZnMPiIiIBDVE6gBERESOsKiIiEhoLCoiIhIa\ni4qIiITGoiIiIqGxqIiISGgsKiIiEprDotq6dSseeOABhIeHIzw8HFOmTMGRI0fsxmi1WsTFxSEs\nLAzTp09HdXW13fyOjg7k5uYiKioKarUac+bMQV1dnfv3hIiIPJLDolKr1Vi3bh2OHj2KsrIyjB8/\nHo8++ij+/e9/AwCKi4tRUlICnU4HvV4PhUKBzMxMtLa22taRn5+PgwcPorS0FIcOHUJLSwuys7PR\n3d3dv3tGREQeQdbXO1P89Kc/xdq1a/HYY48hNjYWixYtwvLlywEA7e3t0Gg0WL9+PZ544gk0NzdD\no9GgpKQEWVlZAIBLly4hISEBu3fvRnp6uvv3iIiIPEqvz1FZrVbs2bMHHR0dSE1NRW1tLYxGo13Z\n+Pj4IDU1FZWVlQCAM2fOoLOz026MWq1GTEyMbQwREZEjQ50N+PzzzzFlyhR0dHTA19cXr776KjQa\nja1oFAqF3fiQkBA0NDQAAIxGI+RyOYKDg+3GKBQKmEwmd+0DERF5MKdFde+996K8vBzNzc3Yt28f\nnnzySRw4cMDhMjKZzG0BiYjozub00J+XlxciIyORmJiIP/7xj0hOTsbWrVuhUqkA4KZ3RiaTCUql\nEgCgVCphtVrR1NRkN8ZoNNrGEBEROdLnz1FZrVZ0d3cjMjISKpUKer3eNq+9vR0VFRVISUkBACQl\nJcHLy8tuzKVLl1BTU2MbQ0RE5IjDolq7di1OnjyJ2tpafP7553j22WdRXl6O7OxsAMDixYtRXFyM\nAwcO4IsvvsDTTz8Nf39/2xV+AQEBmDdvHgoKCvDxxx/js88+w6JFixAfH48JEyb0+84NNIPBIHWE\n28bs0mB2aTD74OLwHJXRaMTChQthNBoxfPhwxMfHY8+ePZg4cSIAYOnSpbBYLMjNzYXZbEZycjL2\n7t0LPz8/2zq0Wi3kcjlycnLQ3t6OtLQ0bNmyheexiIioV/r8OSq6NYPBAI1GI3WM28Ls0mB2aTD7\n4MJ7/RERkdBYVEREJDQWFRERCY1FRUREQmNRERGR0FhUREQkNBYVEREJjUVFRERCY1EREZHQWFRE\nRCQ0FhUREQmNRUVEREJjURERkdBYVEREJDQWFRERCY1FRUREQmNRERGR0FhUREQkNBYVEREJjUVF\nRERCGyp1ACJ3qm+zoqHNKmmG0GFyhA2TS5qByJOwqMijNLRZseyEWdIML6YGsqiI3IiH/oiISGgs\nKiIiEhqLioiIhMaiIiIiobGoiIhIaCwqIiISGouKiIiE5rCoXnjhBUycOBHh4eGIjo7G7Nmzcfbs\nWbsxixcvRlBQkN1jypQpdmM6OjqQm5uLqKgoqNVqzJkzB3V1de7fGyIi8jgOi6q8vBwLFizAkSNH\nsH//fgwdOhQPP/wwzOb/fqBSJpNh4sSJqKmpsT127dplt578/HwcPHgQpaWlOHToEFpaWpCdnY3u\n7u7+2SsiIvIYDu9MsWfPHruvN2/ejPDwcFRWVmLq1KkAgJ6eHnh5eUGhUPzoOpqbm7Fjxw6UlJQg\nLS3Ntp6EhASUlZUhPT3dHftBREQeqk/nqFpaWtDd3Y3AwEDbNJlMhoqKCmg0GiQnJ2Pp0qVobGy0\nzT9z5gw6OzvtCkmtViMmJgaVlZVu2AUiIvJkfbrXX15eHkaPHo2xY8fapk2aNAkZGRmIiIhAbW0t\nNmzYgIyMDJSVlcHb2xtGoxFyuRzBwcF261IoFDCZTO7ZCyIi8li9LqpnnnkGp06dwuHDhyGTyWzT\nZ82aZft3XFwckpKSkJCQgPfffx8zZsxwb1oiIrrj9Kqo8vPz8c477+DAgQOIiIhwODY0NBT33HMP\nzp8/DwBQKpWwWq1oamqye1dlNBqRmpp6y/UYDIbeRBPOYM0NeEZ2i49S4iSAxWKBwVDb6/Ge8LwP\nRsw+sDQazW0v67SoVq5ciX379uHAgQOIjo52usLGxkbU19dDpVIBAJKSkuDl5QW9Xo+srCwAwKVL\nl1BTU4OUlJRbrseVnZKKwWAYlLkBz8ne2ngdQIekeXx9faEZ2bvn0lOe98GG2QcXh0W1YsUK7Nq1\nCzt27MDw4cNx+fJlAIC/vz/8/Pxw7do1aLVazJw5E0qlEhcuXMC6deugVCoxffp0AEBAQADmzZuH\ngoICKBQKBAYGYtWqVYiPj8eECRP6fQeJiGhwc1hU27Ztg0wmw8yZM+2m5+XlYeXKlZDL5Th79ize\nfvttNDc3Q6VSYfz48di+fTv8/Pxs47VaLeRyOXJyctDe3o60tDRs2bLF7lwXERHRj3FYVFevXnW4\nsI+Pz02ftfox3t7e0Ol00Ol0fUtHRER3PN7rj4iIhMaiIiIiobGoiIhIaCwqIiISGouKiIiExqIi\nIiKhsaiIiEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGoiIhIaCwq\nIiISGouKiIiExqIiIiKhsaiIiEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMai\nIiIiobGoiIhIaCwqIiISGouKiIiE5rCoXnjhBUycOBHh4eGIjo7G7Nmzcfbs2ZvGabVaxMXFISws\nDNOnT0d1dbXd/I6ODuTm5iIqKgpqtRpz5sxBXV2de/eEiIg8ksOiKi8vx4IFC3DkyBHs378fQ4cO\nxcMPPwyz2WwbU1xcjJKSEuh0Ouj1eigUCmRmZqK1tdU2Jj8/HwcPHkRpaSkOHTqElpYWZGdno7u7\nu//2jIiIPMJQRzP37Nlj9/XmzZsRHh6OyspKTJ06FT09Pdi0aROWLVuGGTNmAAA2bdoEjUaD3bt3\n44knnkBzczN27NiBkpISpKWl2daTkJCAsrIypKen99OuERGRJ+jTOaqWlhZ0d3cjMDAQAFBbWwuj\n0WhXNj4+PkhNTUVlZSUA4MyZM+js7LQbo1arERMTYxtDRER0K30qqry8PIwePRpjx44FAFy+fBkA\noFAo7MaFhITAaDQCAIxGI+RyOYKDg+3GKBQKmEym2w5ORER3BoeH/r7vmWeewalTp3D48GHIZDKn\n43szxhGDweDS8lIZrLkBz8hu8VFKnASwWCwwGGp7Pd4TnvfBiNkHlkajue1le1VU+fn5eOedd3Dg\nwAFERETYpqtUKgCAyWSCWq22TTeZTFAqb/zAUCqVsFqtaGpqsntXZTQakZqaestturJTUjEYDIMy\nN+A52VsbrwPokDSPr68vNCN791x6yvM+2DD74OL00N/KlSvxz3/+E/v370d0dLTdvIiICKhUKuj1\netu09vZ2VFRUICUlBQCQlJQELy8vuzGXLl1CTU2NbQwREdGtOHxHtWLFCuzatQs7duzA8OHDbeek\n/P394efnB5lMhsWLF6OoqAgajQZRUVEoLCyEv78/srKyAAABAQGYN28eCgoKoFAoEBgYiFWrViE+\nPh4TJkzo9x0kIqLBzWFRbdu2DTKZDDNnzrSbnpeXh5UrVwIAli5dCovFgtzcXJjNZiQnJ2Pv3r3w\n8/OzjddqtZDL5cjJyUF7ezvS0tKwZcsWl89jERGR53NYVFevXu3VSvLy8pCXl3fL+d7e3tDpdNDp\ndH1LR0REdzze64+IiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGoiIhIaCwqIiISGouK\niIiExqIiIiKhsaiIiEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGo\niIhIaCwqIiISGouKiIiExqIiIiKhDZU6AJGn8R4CfNp4vVdjLT5KtPZybG+FDpMjbJjcreskkhKL\nisjNmjq6seaTb/uwRIdbt/9iaiCLijwKD/0REZHQWFRERCQ0FhUREQmNRUVERELjxRREHqYvVx26\nwtEVi7zykNzJaVGVl5fjpZdeQlVVFerr67Fx40bMnTvXNn/x4sV466237JYZM2YMjhw5Yvu6o6MD\nq1evxt69e9He3o7x48ejqKgI99xzjxt3hYiA27nq0BU/fsUirzwkd3J66K+trQ3x8fHQarXw9fWF\nTCazmy+TyTBx4kTU1NTYHrt27bIbk5+fj4MHD6K0tBSHDh1CS0sLsrOz0d3d7d69ISIij+P0HdXk\nyZMxefJkAMCSJUtumt/T0wMvLy8oFIofXb65uRk7duxASUkJ0tLSAACbN29GQkICysrKkJ6e7kp+\nIiLycC5fTCGTyVBRUQGNRoPk5GQsXboUjY2NtvlnzpxBZ2enXSGp1WrExMSgsrLS1c0TEZGHc/li\nikmTJiEjIwMRERGora3Fhg0bkJGRgbKyMnh7e8NoNEIulyM4ONhuOYVCAZPJ5OrmiYjIw7lcVLNm\nzbL9Oy4uDklJSUhISMD777+PGTNm3PZ6DQaDq9EkMVhzA56R3eKjlDgJJD/3KvX2AcBiscBgqJU6\nhkOe8HofTDQazW0v6/bL00NDQ3HPPffg/PnzAAClUgmr1Yqmpia7d1VGoxGpqam3XI8rOyUVg8Ew\nKHMDnpP9xuXS7r13Xl8NGSLtxxOl3j4A+Pr6QjNS3NeTp7ze7xRuf0U3Njaivr4eKpUKAJCUlAQv\nLy/o9XrbmEuXLqGmpgYpKSnu3jwREXkYp++orl27hnPnzgG4cUjh4sWLqKqqQnBwMIKCgqDVajFz\n5kwolUpcuHAB69atg1KpxPTp0wEAAQEBmDdvHgoKCqBQKBAYGIhVq1YhPj4eEyZM6NedIyKiwc9p\nUZ0+fRoZGRkAblzhp9VqodVqMXfuXBQVFeHs2bN4++230dzcDJVKhfHjx2P79u3w8/OzrUOr1UIu\nlyMnJwft7e1IS0vDli1bbvpMFhER0Q85LaoHH3wQV69eveX8PXv2ON2It7c3dDoddDpd39LRoFLf\nZkVDm3XAt/v9W/lct/YM+PaJqH/xXn/kNg1tViw7YZZo6zcuoFg/ZrhE2yei/iL95UFEREQOsKiI\niEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGoiIhIaCwqIiISGouK\niIiExqIiIiKhsaiIiEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGo\niIhIaCwqIiISGouKiIiExqIiIiKhsaiIiEhoLCoiIhIai4qIiITmtKjKy8sxe/ZsjBo1CkFBQdi5\nc+dNY7RaLeLi4hAWFobp06ejurrabn5HRwdyc3MRFRUFtVqNOXPmoK6uzn17QUREHstpUbW1tSE+\nPh5arRa+vr6QyWR284uLi1FSUgKdTge9Xg+FQoHMzEy0trbaxuTn5+PgwYMoLS3FoUOH0NLSguzs\nbHR3d7t/j4iIyKM4LarJkydj9erVmDlzJoYMsR/e09ODTZs2YdmyZZgxYwbi4uKwadMmtLa2Yvfu\n3QCA5uZm7NixA+vXr0daWhoSExOxefNmfP755ygrK+uXnSIiIs/h0jmq2tpaGI1GpKen26b5+Pgg\nNTUVlZWVAIAzZ86gs7PTboxarUZMTIxtDBER0a0MdWXhy5cvAwAUCoXd9JCQEDQ0NAAAjEYj5HI5\ngoOD7cYoFAqYTCZXNk9EgvIeAnzaeF2y7YcOkyNsmFyy7ZN7uVRUjvzwXFZfGQwGNyUZWIM1N+B6\ndouP0k1Jbp8I5z2lziD19gGg0dKFgn+1Oh/YT/7n5z9B6yWjwzF38veqFDQazW0v61JRqVQqAIDJ\nZIJarbZNN5lMUCpv/NBSKpWwWq1oamqye1dlNBqRmpp6y3W7slNSMRgMgzI34J7srY3XAXS4J9Bt\n+uF51Dsxg9TbFyGDr68vNCNv/Xq+079XBxuXXk0RERFQqVTQ6/W2ae3t7aioqEBKSgoAICkpCV5e\nXnZjLl26hJqaGtsYIiKiW3H6juratWs4d+4cgBuHFC5evIiqqioEBwdjxIgRWLx4MYqKiqDRaBAV\nFYXCwkL4+/sjKysLABAQEIB58+ahoKAACoUCgYGBWLVqFeLj4zFhwoR+3TkiIhr8nBbV6dOnkZGR\nAeDGeSetVgutVou5c+di48aNWLp0KSwWC3Jzc2E2m5GcnIy9e/fCz8/Ptg6tVgu5XI6cnBy0t7cj\nLS0NW7Zscfk8FhEReT6nRfXggw/i6tWrDsfk5eUhLy/vlvO9vb2h0+mg0+n6npCIiO5o0p91JSIi\ncoBFRUREQmNRERGR0FhUREQkNBYVEREJjUVFRERCY1EREZHQWFRERCQ0FhUREQmNRUVEREJjURER\nkdBYVEREJLR++wu/NPDq26xoaLPe1rIWH+X//uHD23fd2uPS8kREP4ZF5UEa2qxYdsLswhpc++u8\n68cMd2l5IqIfw0N/REQkNBYVEREJjUVFRERCY1EREZHQWFRERCQ0FhUREQmNRUVEREJjURERkdBY\nVEREJDQWFRERCY1FRUREQmNRERGR0FhUREQkNBYVEREJzeWi0mq1CAoKsnvExsbeNCYuLg5hYWGY\nPn06qqurXd0sERHdIdzyjuree+9FTU2N7XHixAnbvOLiYpSUlECn00Gv10OhUCAzMxOtra3u2DQR\nEXk4txSVXC6HQqGwPYKDgwEAPT092LRpE5YtW4YZM2YgLi4OmzZtQmtrK3bv3u2OTRMRkYdzS1F9\n/fXXiIuLQ2JiIp588kl8/fXXAIDa2loYjUakp6fbxvr4+CA1NRWVlZXu2DQREXk4l4tqzJgx2LRp\nE/bs2YO//vWvuHz5MqZOnYqrV6/i8uXLAACFQmG3TEhICIxGo6ubJiKiO8BQV1cwadIku6/HjBmD\nxMRE7Ny5E8nJybdcTiaTubppIiK6A7hcVD80bNgwxMbG4vz583jooYcAACaTCWq12jbGZDJBqVQ6\nXI/BYHB3tAEhZW6Lj+PntL91d3dLun1mEGP7ImSwWCwwGGodjhmsP2OAwZldo9Hc9rJuL6r29nbU\n1NRg/PjxiIyMhEqlgl6vR1JSkm1+RUUF1q9f73A9ruyUVAwGg6S5WxuvA+iQbPtDhkj/sTxmkH77\nImTw9fWFZuStvxel/l51xWDOfrtcLqrVq1dj2rRpUKvVaGxsxPPPPw+LxYI5c+YAABYvXoyioiJo\nNBpERUWhsLAQ/v7+yMrKcjk8ERF5PpeLqr6+HvPnz8eVK1cQEhKCMWPG4IMPPsCIESMAAEuXLoXF\nYkFubi7MZjOSk5Oxd+9e+Pn5uRyeiIg8n8tFtW3bNqdj8vLykJeX5+qmiIjoDiT9wWwiIiIHWFRE\nRCQ0FhUREQmNRUVEREJjURERkdDc/oFfIiIC6tusaGizun29Fh/l/36437nQYXKEDZO7PcNAY1ER\nEfWDhjYrlp0w99Pae3cHmhdTAz2iqHjoj4iIhMaiIiIiobGoiIhIaDxHRUQex3sI8KmDCw76ckHC\n7bpu7enX9d9JWFRE5HGaOrqx5pNvnYzq3z+Js37M8H5d/52Eh/6IiEhoLCoiIhIai4qIiITGoiIi\nIqGxqIiISGgsKiIiEhqLioiIhMaiIiIiobGoiIhIaLwzhZt809qFankozl2wSLL96AD+VxKRZ+JP\nNzc5c6UTRf9PmpICgCdihiHxbm/Jtk9E1F946I+IiITGoiIiIqGxqIiISGgsKiIiEhovpiAaQAEd\n3yL73HuIvfof+HX99+Kba0N9UR30f/B21C/R/BP+HSOi72NREQ2Q2KvnsOHUXxHZWv+j88c3nMbk\niyewZuz/RXVQ1ACnIxLXgB76e+WVVzB69GiEhoZiwoQJOHny5EBunkgysVfPYe0nG29ZUt/5aWsd\n1n6yEQlXvhygZETiG7Ci2rt3L/Lz87FixQocO3YMY8eOxa9+9St88803AxWBSBIBHd9iw6m/Irql\nd6/16JZv8D+VxUitP93PyYgGhwErqo0bN+LRRx/FY489Bo1GA51OB5VKhdLS0oGKQCSJ7HPvOX0n\n9UNhliv4/b9fh09Xez+lIho8BqSorl+/js8++wwTJ060m56eno7KysqBiEAkmdir/7mt5SJbLmGO\n4V03pyEafAbkYoorV67AarVCqVTaTQ8JCYHRaByICESS+f7VfX0xBEBiU417w9AdxXsI8Gnjdcm2\nHzpMjrBhcpfXIzObzT1uyONQfX09Ro0ahUOHDuEXv/iFbfqf//xn7N69G5988kl/RyAiokFqQA79\n3X333ZDL5Te9ezKZTFCpVAMRgYiIBqkBKSpvb28kJSXho48+spv+0UcfISUlZSAiEBHRIDVgH/hd\nsmQJFi1ahPvuuw8pKSkoLS2F0WhETk7OQEUgIqJBaMCKKjMzE01NTSgsLMTly5cxatQo7Nq1CyNG\njBioCERENAgNyMUUREREt0uIu6c3NDTgqaeeQnR0NEJDQzFu3DiUl5dLHcuphIQEBAUF3fTIzs6W\nOppTXV1dWLduHRITExEaGorExERs2LABVqtV6mi90tLSgry8PCQkJCAsLAxTp07Fp59+KnWsm5SX\nl2P27NkYNWoUgoKCsHPnzpvGaLVaxMXFISwsDNOnT0d1dbUESW/mLPv+/fsxa9YsREdHIygoCMeP\nH5co6c0cZe/q6kJBQQEeeOABqNVqxMbGYsGCBcLcJcfZ875hwwaMHTsWarUakZGRmDlzJk6dOiVR\nWnu9eb1/5/e//z2CgoLw0ksvOV2v5EVlNpsxdepUyGQy/OMf/8CpU6eg0+mgUCikjubUxx9/jJqa\nGtvj448/hkwmQ2ZmptTRnCoqKsKrr74KnU6HTz75BM899xy2bduGF154QepovfK73/0OZWVlePnl\nl3Hy5ElMnDgRM2fORH193+4A0d/a2toQHx8PrVYLX19fyGQyu/nFxcUoKSmBTqeDXq+HQqFAZmYm\nWltbJUr8X86yWywWjBs3Dn/6058A4Kb5UnKU/dq1a6iqqkJubi6OHj2KnTt34ptvvkFWVpYQv6g5\ne97vvfdeFBYW4sSJE3jvvfcQERGBRx55RIjPpDrL/p19+/bh9OnTCAsL69XrRvJDf+vWrcPJkydx\n+PBhKWO4RWFhIf72t7/hyy+/xE9+8hOp4ziUnZ2Nu+++GyUlJbZpTz31FMxmM9566y0JkzlnsVgw\ncuRIvP7665g2bZpt+oQJEzBp0iSsXr1awnS3NmLECDz//POYM2cOAKCnpwexsbFYtGgRli9fDgBo\nb2+HRqPB+vXr8cQTT0iY1t4Ps3/flStXEB0djYMHD+KBBx6QIJ1jjrJ/58svv8S4ceNw4sQJxMXF\nDWA6x3qT/dtvv0VERAT27t17091/pHSr7BcuXMAvf/lL7Nu3D4888ggWLlyI3/72tw7XJfk7qnff\nfRf33XcfcnJyoNFo8OCDD2Lr1q1Sx+qznp4evP766/j1r38tfEkBwOTJk3H06FEYDAYAQHV1NY4f\nP44pU6ZInMy5rq4uWK3Wm55nHx8fVFRUSJSq72pra2E0GpGenm6b5uPjg9TUVN5abIB9++23AIDA\nwECJk/TN9evXsX37dgQHByMpKUnqOE51dXVh/vz5yM3NhUaj6fVykv89qq+//hrbtm3DkiVLsHz5\nclRVVWHlypUAgAULFkicrvc++ugjXLhwAY8//rjUUXpl/vz5qKurw9ixYzF06FB0dXVhxYoV+M1v\nfiN1NKfuuusujB07Fs8//zzi4uKgVCptdziJiho8f8fp8uXLAHDTYe6QkBA0NDRIEemOdP36daxe\nvRrTpk1DWFiY1HF65b333sP8+fPR1taGkJAQ7Nq1C0FBQVLHckqr1SIkJKTPH0uSvKi6u7tx//33\nY82aNQBuXKDwn//8B6+88sqgKqrt27fj/vvvx89+9jOpo/TKyy+/jDfeeAOlpaWIjY1FVVUV8vLy\nEB4ejnnz5kkdz6nNmzdjyZIlGDVqFORyOZKSkvDII4/gs88+kzqaW4h0vseTdXV1YeHChWhpacHb\nb78tdZxeGz9+PI4fP44rV67gtddew+zZs/Hhhx8iPDxc6mi3dOzYMbz55ps4duyY3fSeHudnnyQ/\n9BcaGoqYmBi7aRqNRpgrcHrDZDLh8OHDeOyxx6SO0mtFRUVYvnw5MjMzERcXh+zsbCxZsgQvvvii\n1NF6JTIyEu+++y7q6urwxRdf4IMPPkBnZyciIyOljtZr390+zGQy2U03mUw33cCZ3K+rqwtPPvkk\nzp49i3379g2qw37Dhg1DZGQk7r//frz00ksYPny4wyvsRFBeXo6GhgbExMQgJCQEISEhuHjxItau\nXYv4+HiHy0peVOPGjUNNjf0dor/66iuhfzP4oZ07d8LHxwdZWVlSR+m1np4eDBli/98/ZMiQXv12\nIxJfX18olUqYzWbo9Xo89NBDUkfqtYiICKhUKuj1etu09vZ2VFRU8NZi/ayzsxM5OTk4e/YsDhw4\nMCiuMnbEarWiu7tb6hgOzZ8/HydOnMDx48dx/PhxHDt2DGFhYViyZAn27dvncFnJD/09/fTTmDJl\nCoqKipCZmYmqqips2bIFBQUFUkfrlZ6eHvz973/HrFmzMGzYMKnj9NpDDz2E4uJiREREICYmBlVV\nVSgpKXF4dZFI9Ho9rFYrNBoNzp8/jzVr1iAmJgaPPvqo1NHsXLt2DefOnQNw4zD3xYsXUVVVheDg\nYIwYMQKLFy9GUVERNBoNoqKiUFhYCH9/fyF+6XGW3Ww248KFC2hubgYAnDt3DnfddRdCQ0Mlf0fo\nKHtYWBgef/xxnDlzBm+++SZ6enps5wsDAgLg4+MjZXSH2QMCAvCXv/wF06ZNg1KpxJUrV7B161Y0\nNDQI8bEBDDKHAAABQUlEQVQYZ6+ZkJAQu/FDhw6FUql0em5Z8svTAeDIkSNYt24dvvrqK4wcORIL\nFizAwoULpY7VK0ePHsXDDz+MDz/8ED//+c+ljtNr165dg1arxf79+213sc/KysIf/vAHeHt7Sx3P\nqXfeeQfPPvss6urqEBQUhIyMDKxZswZ33XWX1NHsHDt2DBkZGQBunHf67h3r3LlzsXHjRgDAc889\nh9deew1msxnJyckoLCxEbGysZJm/4yz7G2+8Ybus+Pvz8/LybBdEScVR9pUrVyIxMdFu+ndE+GXN\nUfbCwkIsWLAA//rXv9DU1ITg4GDcd999WLFihRA/f3rzev++0aNH9+rydCGKioiI6FYkP0dFRETk\nCIuKiIiExqIiIiKhsaiIiEhoLCoiIhIai4qIiITGoiIiIqGxqIiISGgsKiIiEtr/B1DNwPm16NcL\nAAAAAElFTkSuQmCC\n",
                            "text": [
                                "<matplotlib.figure.Figure at 0x10ef2b3d0>"
                            ]
                        }
                    ],
                    "prompt_number": 8
                },
                {
                    "cell_type": "code",
                    "collapsed": false,
                    "input": [
                        "# Plot data and model\n",
                        "plt.scatter(x,y,c=\"dodgerblue\",s=30)\n",
                        "ymu = np.array([np.median(i) for i in M.Ex_mu.trace().T])\n",
                        "yl95 = np.array([np.percentile(i,2.5) for i in M.Ex_mu.trace().T])\n",
                        "yu95 = np.array([np.percentile(i,97.5) for i in M.Ex_mu.trace().T])\n",
                        "plt.plot(xnew,ymu,c=\"red\",linewidth=1.5)\n",
                        "plt.plot(xnew,yl95,c=\"grey\",linewidth=1)\n",
                        "plt.plot(xnew,yu95,c=\"grey\",linewidth=1)"
                    ],
                    "language": "python",
                    "metadata": {},
                    "outputs": [
                        {
                            "metadata": {},
                            "output_type": "pyout",
                            "prompt_number": 14,
                            "text": [
                                "[<matplotlib.lines.Line2D at 0x10fc27950>]"
                            ]
                        },
                        {
                            "metadata": {},
                            "output_type": "display_data",
                            "png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEWCAYAAAC5XZqEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9+PHXvTfrZkHCymSHERISRthbpqIgGqZ7K20d\n1Vr7VatWpcVqHUV/uCpVaQsCIlu2EEaEAIEEkkASSMhOyLoZN3f8/oiJhAxyb25ubnLfz8ejf+Te\nc899f3rwvO9nvY+isLDQiBBCCGFFyrYOQAghhP2R5COEEMLqJPkIIYSwOkk+QgghrE6SjxBCCKuT\n5COEEMLqJPkIIYSwOpOTT2hoKF5eXvX+t2jRotaITwghRAfkYOoHDh48iF6vr/07MzOTKVOmcOed\nd1o0MCGEEB2XycnH29u7zt9r1qzB09NTko8QQohma9Gcj9Fo5Ouvv2bhwoU4OztbKiYhhBAdXIuS\nz/79+7ly5Qr333+/peIRQghhB1qUfNasWcOIESMYMmSIpeIRQghhB8xOPrm5uezYsYP77rvPkvEI\nIYSwA2Ynn7Vr1+Li4sLdd99tyXiEEELYAbOSj9Fo5N///jcLFizA1dXV0jG1e0lJSW0dQpuRttsf\ne2032HfbW8qs5HPo0CFSUlJkoYEQQgizmLzPB2DSpEkUFBRYOhYhhBB2Qmq7CSGEsDpJPkIIIaxO\nko8QQgirk+QjhBDC6iT5CCGEsDpJPkIIIaxOko8QQgirk+QjhBDC6iT5CCGEsDpJPkIIIaxOko8Q\nQgirM6u2mxBCWFtasYKEAhUAA731BHoa2zgi0RKSfIQQNi+tWMHjP7qRXlqdfAI89KyeoZEE1I7J\nsJsQwuYlFKhqEw9Aeomqthck2idJPkIIIaxOko8QwuYN9NYT4KGv/TvAQ8/ALvomPiFsncz5CCFs\nXqCnkdUzNL8uOOiiJ9BD5nvaM0k+Qoh2IdDTSKCnrq3DEBYiyUcIYddkCXfbkOQjhLBbsoS77ciC\nAyGE3WrJEm6DwdBaYdkFST5CCGGGDRs2cO3atbYOo92S5COEsFvmLuHOyMigtLSUzp07t2Z4HZpZ\ncz5ZWVm89tpr7Nmzh9LSUnr37s27777L+PHjLR2fEEK0GnOXcJ88eZLhw4ejUChaO8QOy+TkU1hY\nyKxZsxg3bhzr16+nS5cupKam0q1bt9aITwghWpWpS7ivXbtGRkYGs2fPJjU1tfUC6+BMTj4ffvgh\nfn5+fPLJJ7Wv9ezZ06JBCSGErYqJiWHo0KE4Ojq2dSjtmslzPtu2bWP48OE8+OCDBAUFMXHiRD77\n7LPWiE0IIWxKWVkZiYmJhIWFtXUo7Z7JySc1NZUvvviCvn37snHjRp544glef/11SUBCiA7vzJkz\nDBgwAFdX17YOpd0zedjNYDAwYsQIXnnlFQBCQ0NJTk7m888/59FHH7V4gEIIYQuqqqqIjY1l4cKF\nbR1Kh2By8vHx8WHgwIF1XgsKCiI9Pb3RzyQlJZkeWTtnj22uIW23P/bQ7pSUFDp16kReXh55eXm1\nr9tD228UFBTU4nOYnHzGjBlDYmJindcuXrzY5KIDSwTaniQlJdldm2tI2+2v7fbQboPBwKFDh5gz\nZw6+vr61r9tD21uLyXM+Tz31FCdOnODdd98lOTmZ77//nk8//ZRHHnmkNeITQog2l5SUhIeHR53E\nI1rG5J7PsGHD+Pbbb3njjTd45513CAwM5OWXX+bhhx9ujfiEEMLiTKlkbTQaOXHiBGPHjrVWeHbB\nrAoHM2fOZObMmZaORQghWp2plazT0tLQ6/X06dPHmmF2eFLbTQhhV0ytZH3ixAlGjBghpXQsTJKP\nEEI0Iicnh4KCAgYNGtTWoXQ4knyEEHbFlErWJ0+eZNiwYahUzXvGj2g+eZKpEMKuNLeSdVFREZcv\nX2batGnWDtEuSPIRQtid5lSyPnHiBKGhoTg7O1spKvsiw25CCHEDjUZDYmIiw4YNa+tQOixJPkII\ncYOYmBgGDx4sBURbkSQfIYS4TkVFBXFxcQwfPrzpA403f+KpaJwkHyGEuM6ZM2fo27cvnp6ejR9U\nWIj7+PF4RkVZL7AORpKPEEL8QqvVcvr0aUaOHNnkcc4ffYQqPp6qrl2tFFnHI8lHCCF+ERsbS2Bg\nIN7e3o0eo8jKwunjT7gyO5LYgbeSViyVD8whS62FEDbFlKKfllRVVUVMTAwLFixo+ri3/45Bq+We\nEW9zOcrrprXhRMMk+QghbIapRT8t6dy5c/j6+tK1iaE0ZUoKXmu/4tsxj3K5W3/g19pwN9s3JOqS\nYTchhM0wteinpeh0Ok6ePMmoUaOaPM75L3/B4ODEhzNfafWYOjpJPkIIuxcXF0fXrl3p0aNHo8eo\nYmJw2riRgseW4xTQvfb1pmrDicbJsJsQwmbUFP1ML/l12K21b+w6nY6ff/6ZuXPnNn6Q0YjLq69i\n6NoV5+d/y2qqa8Npq7SE+qgarA0nmibJRwhhM5pb9NOS4uLi6NatGz4+Po0e47B7Nw6HD1O+ciV4\nehJIdW24pKRLBHoEtWp8HZUkHyGETWlO0U9Lqen13H777Y0fpNfj8tpr6Pv0QfvAA1aJyx5I8hFC\n2K1z587RrVu3Jud6HL/5BlV8PJo1a8DJyYrRdWyy4EAIYZdqej1jxoxp/KCSElzefBPd2LHo7rjD\nesHZAen5CCHsUmxsLD4+Pk32epw/+ABlbi5l//0vKKSSgSVJz0cIYXe0Wi0nTpxg7NixjR6jSEvD\n+Z//RLtwIfoRI6wYnX2Q5COEsDunT58mMDCwyWoGLm+8AUDFK7KhtDWYnHxWrFiBl5dXnf8NGjSo\nNWITQgiLq6io4NSpU4wePbrRY1TR0TitX0/lb36DMTDQitHZD7PmfAYMGMDWrVtr/1apWr/8hRBC\nWEJMTAx9+vRpvHK1wYDLn/6EwceHymeesW5wdsSs5KNSqejWrZulYxFCiFal0WiIjY1l6dKljR7j\n+N13OJw4QdmqVeDubsXo7ItZcz6pqakMHjyYsLAwHn74YVJTUy0clhBCWF50dDSDBw9u/CmlGg0u\nr7+OPiyMqiVLrBucnTE5+URERPDJJ5+wYcMGPvzwQ7Kzs5k1axbXrl1rjfiEEMIiioqKSEhIICIi\notFjnD/4AOXVq5SvWAFKWY/Vmkwedps+fXqdvyMiIggLC2Pt2rUsX77cYoEJIYQlHTt2jPDwcFxd\nXRt8X5GWhvOHH6JdsAD9uHFWjs7+KAoLC1tcte/2229n4MCB/P3vf2/w/aSkpJZ+hRBCmK2oqIif\nf/6ZKVOm4ODQ8G/uvn/6E51/+olz332HtokiowKCglpeTLXFFQ4qKipITExk0qRJjR5jiUDbk6Sk\nJLtrcw1pu/21vT20e+PGjYwbN47Bgwc3+L7q6FHcd++m4sUX6TVxYrPP2x7abqtMTj4vv/wyc+bM\nwd/fn7y8PN555x3Ky8tZIpNzQggbdOXKFYqLiwkJCWn4AJ0O9fPPYwgIoPLpp60bnB0zOflkZmby\nyCOPkJ+fT9euXYmIiGDPnj0EBAS0RnxCCGE2o9HIoUOHGDduXKP7EZ2++AJVXFx11epG5oOE5Zmc\nfL744ovWiEMIISwuISEBpVLZ6NCYIicHl7feomraNJOqVhuNRvbu3dtkeR7RNKlqLYQwSVqx4tcn\njXrrUdhotWedTkdUVBSzZ89uNEaXP/8ZysupWLnSpKrV8fHxZGdnEyild8wmyUcI0WxpxQoe/9GN\n9NLq5BPgoeed0Y0/kqAtxcTE0KNHD/z9/Rt8XxUVhdN//kPFs89i6N+/2ectKSnh8OHDLFiwgMLC\nQkuFa3dkF5UQotkSClS1iQcgvUTF5XK3NoyoYRqNhpiYGCZMmNDwAVot6t//HkPPnlS+8EKzz1sz\n3BYWFiYlxlpIej5CiHbpxuG/QM9ftywePXqU4OBgOnfu3OBnnVetQnXhApr//c+kRQbx8fGUlpY2\nWSVBNI8kHyFEsw301hPgoSe95Ndht15qDeBh1TgaGv5bPUNDoKeRnJwckpOTue+++xr8rCI1FeeV\nK6maOxfdrFnN/s7i4uLa4Tap5N9yknyEEM0W6Glk9QzNrz2OLnoqs7OxdvJpaPgvoUBFgEcVBw8e\nZMyYMbi4uNT/oNGI+g9/AKWS8r/+tdnfZzQa2b17N8OGDZPhNguR5COEMEmgp5FAT13t30lZLa7Q\nZTFJSUlUVlY2uqHU4YcfcPzxR8rfegujCXsTz5w5Q1VVFSNHjrRUqHZPFhwIIW4qrVjBnlQH9qQ6\nkFbc9kura4b/agR46OnXqYJDhw4xZcoUlA1VpC4uRv3HP6IPDUX7+OPN/q6CggKOHTvGrFmzGj6v\nMIv0fIQQTWpqfqWtNDT8l3Y2Gl9f30arrbi89RaKrCzKvvkGGikueiO9Xs+uXbsYO3YsXl5eFotf\nSPIRQtxEY/Mr1w+9tYXrh//y8/M5e/Ysy5Yta/BY1YkTOH32GdpHHkE/YkSzvyM6OhoXFxeGDh1q\nkZjFr6QPKYRo12r23owZMwb3hh57rdWi/u1vMfr6UvHKK80+b0ZGBmfPnmXmzJk2W8WhPZOejxCi\nSQ0trx7YRX+TT1lPXFwcer2e0NDQBt93fu89VOfPV+/paezx2TeorKxk165dTJs2DTc329tE2xFI\n8hFCNKmh+ZVAD9tY4abRaDhy5Ajz589vcDGAMj4e53ffRRsZadKengMHDhAYGEh/E8ruCNNI8hFC\n3NSNy6ttxcGDBxk8eDDdu3ev/2ZVFeqnnsLo6UnFihXNPmdCQgJZWVksXbrUgpGKG0nyEUK0S8nJ\nyWRnZzNjxowG33f+xz9wOH0azZo1GJv56IOioiIOHDjA/PnzcXR0bLKEj2gZST5CiHZHq9Wyf/9+\nZs6ciaOjY733lbGxOK9cifauu9DNm9esc+r1enbs2EFERAQ9evSwySXmHYmsdhNCtDuHDx8mMDCw\n4efpVFbi+uSTGLt0oeKdd5p9zmPHjuHi4sKwYcOAxpeYC8uQ5COEaFfS09NJTk5m0qRJDb7v/Ne/\nooqLo/z99zF6ezfrnFeuXCE+Pl6WVVuRDLsJIdqFtGIF53MMXNy/hxFjpzVYOFR19CjO77+P9t57\n0c2Z06zzajQadu3axezZs3G97vEKtr7EvL2T5COEsHk18y9+uXtQGwLYmjCU1T1vmH8pKcH1iScw\n9uxJ+dtvN+u8BoOBHTt2EBISUm8Iz5aXmHcEknyEEDYvoUCFtvAqfbXn+N7jCSoaKPGj/tOfUFy5\ngmb7dvBo3iMeoqOjARg9enSD79vqEvOOQOZ8hBA2z6DTMrFsM0fVt1KhrF9xwGHzZpy+/prKZ55B\nP3Zss8555coVzp49y+zZs6VadRuQ/8eFEDavMvknylx8uew0GKg7/6K4ehX100+jGzaMypdeatb5\nSktL2blzJ7Nnz264HpxodTLsJoSwaWlpaWRdTuKxO+9jhkYDXDf/YjDg+tRTKLRayj/7DJycbno+\nvV7P9u3bCQ8Pb3iptrCKFvV83nvvPby8vHjhhRcsFY8QQtSqrKzkxx9/ZPr06fTr5sz03jqm99bV\nTvw7ffghDgcPUv7Xv2JoZh22w4cP4+TkRERERGuGLm7C7OTz888/s2bNGoYMGSLr4oUQreLgwYP0\n6tWLPn361HtPdeIELm++iXb+fKruvbdZ50tISODSpUvMnj1b7lttzKzkU1RUxGOPPcaqVavo3Lmz\npWMSQrQCW3sU9s1cvHiRq1evNryZtKgI14cfxujrS/n770MzEkl+fj4HDhxg7ty5De4REtZl1pzP\nM888w/z585kwYQJGo6x7F8LWtbc6ZRqNhn379jF37lycbpzHMRpRP/ssivT06mXVzfgBXFlZydat\nW5kwYULDFbCF1ZmcfNasWUNqaiqff/45gHRdhWgHbPVR2A0xGo38+OOPhISE4OfnV+99p6++wmnj\nRipeeQV9I/tzbjzfzp076dmzJ0OGDGmNkIUZTEo+SUlJ/OUvf2Hnzp2oVNX/kI1G4017P0lJSeZH\n2E7ZY5trSNttj5Z+QN39MdoqLUlJlyxyfku2OyUlhaKiIoKDg+udV52YyOAXX6RozBiS5s6FZnxv\nYmIihYWFDBo0qFWuj61e89YUFBTU4nOYlHyio6PJz89nzJgxta/p9XqOHj3KV199RUZGRoPlzS0R\naHuSlJRkd22uIW23zba7FCvq1SkL9VER6NHyeC3Z7vz8fPbt28eiRYvqzyeXlOC+ZAl4e5Px8Tdc\nVvoATT9n5+LFi2RmZrJkyZJWeRy2LV9zW2dS8pk7dy4jRoyo/dtoNLJ8+XL69+/Pc88912DiEUK0\nvfZQp0yn07F9+3YmTJhQP/EYjah/9zuUycmk/HcLj8T0vun8VV5eHnv37mXevHmtknhEy5iUfDp1\n6kSnTp3qvKZWq+nUqRODBg2yaGBCCMuy9Tplhw4dwtvbm+Dg4HrvOX32GU6bNlHx5z9zasBk0g80\nPX9VXl7ODz/8wOTJk/Hx8bFK/MI0LS6vo1AoZNGBEKJFkpOTSU5O5pZbbql3P1GdOIHL//0fVbNm\nUfn00zc9l16vZ9u2bQQFBcmPYhvW4vI6W7dutUQcQgg7VVpayp49e7jtttvq7b9R5OXh+sADGH18\nKP9//w+Uyiafs2M0Gtm/fz+Ojo6MHz++Wd+fVqz4dTiyifkjYVlS200I0WZqnqcTFhaGv79/3Td1\nOlwfeghFbi6lu3Zh9PICmp6/On36NFlZWYybs4h9V6r3Bw301gENJ5j2tv+pI5HkI4RoM8ePH0ep\nVDZYZ83lL3/B4aefKFu1CkN4eJ33Gpq/SklJ4cSJE0y8bTG/2e9Vm1CeGVHOdwlODSaYG/c/ZWuU\nZJQqSSioHvqTnlDrkeQjRAdl68NJaWlpnD17lmXLltV7no7D5s04f/ABlQ89RNWyZTc9V25uLj/+\n+CO333478VqvOglFb1A0e4PtAyGVvHZEzVXpCbU6ST5CdEC2Ppyk0WjYuXMns2bNqrcMWhkfj+tT\nT6EbNYqKFStueq7S0lJ++OEHpk6dip+fH/GpzY/jxvkjd0djbeIB264E0d5J8hHCBrW012LL5XRq\n5nlCQkLo1atX3TcLC3Fdtgyjuztla9aAs3OT59Jqtfzwww+EhoYyYMAAoH5CUSmNBLjr6yTimgUK\nN84fVRks2VLRFEk+QtgYW+q1tMbQ3fHjx1EoFIy+sS6bXo/rI4+gTE9Hs2ULRl/fJs9jMBjYvn07\n3bt3rzNnVH9Bgo4Zvaoa3WB7/fxRWgOVIGoSlbAsST5C2BhL9FqaWo7cXJZKgtcnME/NReLi4liy\nZEm9eR6XV1/Fcc8eyt5/H/11JbwaYjQa2bdvH0ajkalTp9bbG1R/QULzNti2h0oQHYUkHyE6IEvc\nRC2RBK9PYO76a8zT7Gbm7NvqzfM4fvstzqtWUfnoo1Q98MBNe1zR0dFkZ2cTGRlZW+TYUmy9EkRH\nIclHCBtjiV4L2MZNtCaBqYxVTCtbxynnCYxw6gn8Gpfq6FHUzz6LbvJkKlasuGmP69y5c8TFxbFo\n0aL6z/oR7YYkHyFsjK0M/TQnCaYVK0iiH5dTHRqfEzIaGVu+nWJlV+KdRgNltW8pUlNxveceDIGB\nlH31FTg4NNnjSk5O5siRI0RGRkqx0HZOko8QNsgWei03S4LNmRMa6K1nrOI4XXWZbPV4iABPw68J\nrLAQt4ULQa+nbN262goGjcnIyGD37t3ccccdeHl52fw+JtE0ST5CiEY1lQSbMyekKE4jvPIQvW5Z\nyiS3KgZ2qahOYFVVuD7wAMqUFDSbNmHo16/2Mw31uHwU2WzdupVZs2bh6+trUysChXkk+QghWkVJ\nSQk7duzg1tmz6NXLndp5HqMR9XPP4XjgAGWrVqGfMKHO527scQU4FXB4+yYmTZpE7969AdvexySa\nR5KPEMIsTc0J6XQ6tm7dSnh4eL2NpM7/+AdOX39NxfPPN1o6p6bHpdFoWLduAyNHjpTHI3QwknyE\nECarmW/57bAKOjvpKdMpGdDFUDsndODAATw9PRk5cmSdzzlu2IDLG2+gjYyk8v/+r8nvqKioYOPG\njQwZMoTwGwqLWmpFoGg7knyEECapN9/irmf1TE1t4jl79iwZGRksXry4zuZPVVQU6iefRDd2LOUf\nfQRNPIRSq9WyadMmevfu3WDFa1tZESjMJ8lHCGGSevMtpb/Ot2RkZHDkyBEWLlxYZw+OMjER12XL\nMPTqRdnatXDDQ+OuV1VVxebNm+nevTsTJkxo9EnJtrAiUJivxY/RFkIIqK4uvW3bNmbOnInXdcum\nFRkZuN11Fzg6olm/vskl1Tqdji1btuDp6cm0adMaTTyi/ZPkI4QwSc18S40Adz39O1WydetWwsLC\n6NOnT+17ivx83O68E0VhYXXi+WW1WkNqFik4OzszY8YMSTwdnAy7CSFMcuN8i4+qgIToo3h4eNSd\nnykpwTUyEmVqKpoNG+o9jfR6Op2OH374AWdnZ2bPnl2v6KjoeCT5CCFMdv18y86dRygoKCAyMvLX\n3kpFBW5Ll6I6c4ayb76pt5fnelVVVWzZsgW1Ws2sWbMk8dgJucpCCLNdunSJlJQUbr/9dhwdHatf\n/KV6gerwYco/+QTdnDmNfr6yspLvv/8eNzc3STx2Rno+QgizZGRksGfPHkaMGIGHh0f1iwYD6qee\nwnHnTsr//neqFi5s9PMVFRVs2rSJ7t27y+ICO2Tyz4zPPvuM8ePH07NnT3r27MnMmTP58ccfWyM2\nIUQrSytWsCfVgT2pDqQVN//mn5+fz9atW+uubDMacXnhBZzWr6fi1VfRPvJIo5/XaDR89913+Pv7\nS+KxUyb3fPz9/XnjjTfo168fBoOBtWvXsmzZMvbt20doaGhrxCiEaAXmFucsLi5m06bqWmt9+vQh\nKSkJAOc338T5iy+ofPppKp99ttHPFxUVsWnTJgYPHsyoUaMk8dgpk3s+t956K7fccgu9e/emb9++\nvPzyy7i7uxMTE9Ma8QkhWkljxTmbotFo2LhxIyNGjKhTa83pgw9wefddKh94gIrXXmu0ekF+fj7f\nffcd4eHhjB49WhKPHWvRnI9er+f777+nsrKScePGWSomIYQNqpmjGTx4MMOGDat9vdv69ahXrkS7\nYAEV777baOLJyMhg69atTJw4kcGDB1srbGGjzEo+cXFxzJw5k8rKStRqNf/6178ICgqydGxCCBPV\nFPysMoC3i5FSLQzwNjQ4lBbgYeDZEeXoDAoq9LD7skOjxTm1Wi3ff/89gYGBjBo1qvZ1x//8h14r\nV1I1ezblq1eDquGeU2pqKrt27WLmzJl1NqEK+6UoLCw0uRpfVVUVV69epaioiM2bN/Ppp5+yZcuW\nOr+GrlczJiyEaD1atQ8vHPOtHUrzdzdwa18tO5IdeWdMJk7lWY0eG+Cu528TClGXpmI01r0l6PV6\noqOjcXd3JyQkpHaozHvXLvq8+iolw4eT9P77GJ2dG4wrIyODuLg4Ro4cWafsjmi/LNHZMKvn4+jo\nWPtQp7CwMGJiYvjss8/4+OOPGzze3npFSUlJdtfmGtL2tmv7nlSHOnM4V0uVuKiqC39m6ryYHuTR\n6LHppSqyK1yY3r9/nXPWVB7o0aMHM2fOrE08jhs3on71VfRjx3JxxQr6h4Q0GFNsbCyJiYksXLiQ\nrl27WrK5NqGtr3l7ZpF9Pnq9HoPBYIlTCSFshMFgYOfOnfVqrTlu3Ij60UfRjxmDZt06DBkZ9T5r\nNBqJjo4mLi6OyMhIOnfubO3whY0zOfm89tprzJo1Cz8/P0pLS/nuu++Iiopiw4YNrRGfEKKZbnzA\nmr+7gQp9ww9au9nD2IxGI3v37kWr1XLHHXfUVh5w/M9/UC9fXpt4cHOrF4fRaOTAgQNcvXqVRYsW\n4dbAMUKYnHxycnJ47LHHyMnJwdPTk5CQEDZs2MDUqVNbIz4hRDNdX/Dz+gUH84O09R601tTD2IxG\nI4cOHSI/P58FCxbg4FB9m3Bcswb1M8+gnzQJzdq1DSYenU7Hrl27KCsrIzIyEudG5oGEMDn5NDav\nI4Roe6Y8YK2xY6Ojo7ly5Qp333137QPhnD76CPUrr1A1fTplX38NanW9z1VWVrJlyxZcXFy48847\na5OWEA2Rfx1CdHA1y6+heritqQoGMTExnD9/nsjISFxcXMBoxPnNN3F591208+dXL6duoDdTWlrK\n999/j7+/P5MnT5YCoeKmJPkI0YGZUkLn7NmznDp1isjIyOp5Gr0elxdewPnLL9Hefz/l773X4D6e\n0tJS1q1bR0hICBEREVK1QDSL/DwRogNrbgmduLg4jh8/zl133YWnpydUVuL60EM4f/klFc88Q/n7\n7zeYeDIzMzl27BijR4+WOm3CJNLzEcLOxcfHc/ToUe66667qJdFFRbjdey8OP/1E+VtvoV2+vMHP\nxcXFcfjwYUJDQxkyZIiVoxbtnSQfITqwmy2pPnPmDD///DMLFizAy8sLRUYGbpGRKBMSKFu9mqpF\ni+qdU6fTcfDgQdLT04mMjCQ/P99q7REdhyQfITqwxpZU12wCjY+PJzIykk6dOqE8fx63yEgUhYWU\nrV+ProHtE0VFRWzfvh0PDw8WL16Ms7OzJB9hFkk+QnRwNy6pNhgMHDhwgIyMDBYuXIibmxuqAwdw\nu+8+jGo1pdu2YQgLq3eelJQUdu/ezciRIxk2bJjM74gWkeQjhB3RarVs374dg8FQuwnU8d//Rv3c\ncxgGDEDzv/9hDAys8xmDwcDRo0c5f/48t912G/7+/m0UvehIJPkIYSdKSkpqi4ROnToVlVKJ89/+\nhsuKFVTdcgtl//oXeHrW+YxGo2H79u2oVCqWLl2Kq6trG0UvOhpJPkLYgezsbLZs2UJ4eDgjRoxA\nYTBU7+H5/HO0S5ZQ/uGH4OhY5zN5eXls3ryZ4OBgRo8eXbtx9MZNqzL8JswhyUeIDu7ixYvs3buX\nW265hf79+0NpKa6PP47jtm1UPv10g4+9TklJ4ccff2Ty5Ml1Hpfd0KbVd0b3sGZzRAchyUeIDuzM\nmTNER0fIKUy5AAAgAElEQVQzf/58evTogSItDbclS1DGx1P+t79xackTJFz+tRcT4GHg1KlTnDhx\ngttvvx0/P78652to0+rlcjdkl48wlSQfITogo9FIVFQUly5dYuHChXTq1AnV8eO43nsviooKytat\nI2XUjDq9mEB3LY+qf6D0WjaLFy+urnQgRCuR5COEFZlS5NOc46F6ddqePXsoKChg4cKFqNXq6hVt\nv/89hsBANFu2YBg4kITUX3sxakMJQzPXU9DdhXsWLqytZn2jhjat9lJrAI8GjxeiMZJ8hLASU4p8\nmnM8VFcf2L59O3q9nrvuugtHhQKXP/wB508/pWraNMq+/BJueKpoV91VpmnWkeg8nGmjRuDkpG/k\n7A1vWq3MzkaSjzCVJB8hrKSxIp+NPX/H1ONLS0vZsmULXl5ezJgxA1VpKa4PPYTj3r1ULl9Oxeuv\nw3XP2BnorWe04gR9NfuJUt+OoUsQg7pqbtqOGzetJmXdvDcmxI0k+QjRAWRmZrJ161bCw8MZOXIk\nqpQUXJcsQXnpEmUffkjVfffVOV6n05EQfYAxhgx8pi0iwqMLA7to6j3xVIjWIslHCCu5WZFPc443\nGo2cO3eOI0eOMGPGDPr27YvD/v2oH3wQlEo0mzahnzixzmeur89277LFv8zvNO/pp0JYiiQfIayk\nsSKf5h6v0+nYv38/mZmZLFy4EK/OnXH65z9xefVVDIMGoVm7FmPv3nXOKfXZhK2Q5COEFd04X2Lu\n8UVFRWzbto3OnTuzePFinLRa1A89hNOmTVTdfjtln3wC7u61xxsMBo4cOcKFCxcsWp8trVhBEv24\nnOrQ7NV4QoAkHyGswpwl042p6b1EREQQHh6O6vx5XB96CGViIuWvv472d7+rU7GgpKSEHTt24Ojo\naNH6bOasxhOihiQfIVqZpW7SBoOBY8eOER8fX9178fPD6YsvcHn5ZYweHtXzO5Mn1/nMpUuX2Lt3\nL+Hh4URERFh0mM3U1XhCXE+SjxCtzBI3aY1Gw44dO1AqldW9F60W13vuwXHbNqqmT6f8448xdu9e\ne7xOp+PQoUMkJyczd+5c/Pz8LNr7EqKlTE4+7733Hlu2bOHSpUs4OTkxcuRI/vznPzN48ODWiE8I\nu5eamsru3bsJCQlh9OjROMbE4Prggyiysih/8020Tz0Fv1Schupq1Dt27MDb25tly5bh4uLSKkNk\npq7eE+J6JiefqKgoHn30UYYPH47BYODtt99m/vz5HD9+nM437JwWQph/k9br9Rw9epQLFy4we/Zs\nAv39cVq1CpfXX8fo54dm5070I0bU9miMRiMOWSdJjD3OhAkTCA4Orh1ma40hsprVeGez9Tg5Ot10\n9Z4Q1zM5+WzYsKHO36tXr6Znz54cP36cWbNmWSwwIToKU5dYAxQWFrJjxw7UajVLly7FraIC9ZIl\nOO7aRdXcuZT985/QuXNtj6agWMPEss24KStZdMdiBgd0skbTCPQ0UpF9iaDeQXVelyE+cTMtnvMp\nKSnBYDBIr0eIJjR3ibXRaCQ+Pp7Dhw8zatQowsPDcdy5E/Wzz6IoKKB85Uq0jz5au5rtQr4S54I4\n5pXvIt55NLHOE5igK2fwDZtGrTlEJqvgRHO0OPn88Y9/ZOjQoYwaNcoS8Qhht8rKyti7dy9FRUUs\nWLCAbo6OqJ94Aqf//Q/9kCFo1q3DMHRonePTovcztKKQH93uId/Bt9Fzm9P7MpesghPN0aLk86c/\n/Yno6Gh27NghO6WF3bD0kJLRaCQpKYkDBw4QHBzMnDlzcD5xAtdHH0WRkUHFCy9Q+cIL8MtjDmp6\nR1FRUQT0C2aD6i7yS52Bpns0Nb2vtGIFCfkqEvJlSEy0HUVhYaFZ//Jeeuklvv/+e7Zs2VL9aN4m\nJCUlmRWcENaiUCiodOnB5XI3AHqpNThXZGM01v3PQ6v24YVjvr8OKbnreWdMJk7lWWZ9b0VFBefO\nnUOj0TB06FC83dzw/fJLfP/1Lyp9fUl58000ISG1x5eUlHDu3Dn0ej2hoaF07ty5WXG3Vvxt9R2i\nbQUFBd38oJswq+fz4osvsnnz5mYlHrBMoO1JUlKS3bW5hi22vTk9lbRiBb+rM0/hyeoZnvWO3ZPq\nUHdIqVRFps6L6UEeJrXdaDRy9uxZjh49ytChQ4mIiMD5yBHUTz2F6uJFtIsWUfnOO/j98jRRrVbL\n8ePHiY+PZ9SoUYSFhaG8bnn1r4+x9qCpZ+s0Fb+5Gmr36pk3DvE1HVd7ZYv/3tsLk5PP888/z7p1\n6/jmm2/w9PQkOzsbAHd3d9zc3CweoBCNaW5Sac7ktynzFI5KIw+EVOLyy+EeTqYNHly7do09e/bU\nPvCtq4MD6qefxunbb9H37o1m40Z006YB1Unq/PnzREVF0bNnT+6555528d+ZqTXshP0xOfl88cUX\nKBQK5s2bV+f1P/7xj7z44osWC0yIprRGUmmOgd56lg+rYF2CMxml1T2PAHc9q2fe/CFsBoOB06dP\nEx0dXbuSzWnLFtQvvIAiP5+KZ5+l8g9/ALUagKysLA4cOIDRaGTu3Ln4+ja+oMCU+GVjqLAFJief\na9eutUYcQpikNZJKc27KgZ5GfN0MtYkHqoeuEgpU9Gri/AUFBezevRulUsnixYvxqqpC/eCDOG7e\njH7oUDTr12MICwOqn0gaFRXFlStXGDduXJ3Noi1lzVVvQjRFaruJDq1eUnHXE+BhqHecKTdlR2WD\nLzdIr9cTExPDyZMnGTNmDGGhoTj997+4vPIKlJZy7revkXj/MwzorqCHtpKTJ09y5swZQkNDuf/+\n+3950FvjzFl5J0NiwhZI8hHtkik9lfemlHE0wwGdQUGFHv7wk5pVt5TVu1E396bc2HdX3LCY6+rV\nq+zbtw93d3eWLFmCV0oK6jlzcIiOpmzEKB6d9zmHPUNRHDIwSnmSEZUH6RUYwNKlS/H8ZaFBU2Qz\np2jPJPmIdsmUnkp6iZJ/nFTXea0lQ3SNfXfSL8mnuLiYI0eOkJ6ezuTJkwnq1An1K6/g+PXXGLt2\npezjj9kx5h4OH3QnsCqBkeV7qFC44TdhPnOGdWt2HG21mfPG3pbs8RPmkOQj2q22HD5q6Lurqqo4\ndOgQcXFxhIWFMW3CBNy/+gqXlSuhrAztk09S8Yc/QOfOaE5kcFtpFI7GSk6oZ5DmEMQMrzLAtofD\nGuptvTO6RxtHJdojST6iw2vtFV46nY4zZ85w/PhxBgwYwL333EOnAwdwmTQJVXIyVdOnU/HWWxgG\nDiQ/P5+oH34gOyePnE5TOKEPw6hQmhVTW6xca6i3dbnc7bp9RkI0jyQf0eG11govo9FIQkICUVFR\ndOvWjTFjxhBhMOASGYnD0aPoBw1C89136KZPp6ioiOjdu0lOTiYiIoJbb72VOWWOJBSUmx2TrFwT\n7ZkkH2EXLD1El56ezqFDhwCYNWsWPcvLqXrpJdz37MHQrRvl//gH2nvvpUij4ec9e7h48SJDhw7l\n/vvvx8XFxWIxWXvosaHeVi+1ho5YvUC0Lkk+QpggLy+PqKgo8vPzGTduHIM8PVH/7W84fvMNBien\n6iKgv/sdJUD0wYMkJSURGhrKAw88UJt02rOGeluV2dlI8hGmkuQjRDOUlpZy9OjR2mGz2yZPxm3V\nKpxXrcJYVcWlhY8S+8Cz9OnvTWZsNGfPniUkJIT7778ftVp98y+g/TyA7cbeVlKWbcYpbJskHyGa\nUFVVxcmTJzl9+jRDhgzh/nvuwfN//8N50SKUubkU3r6Ax8auIFrdnwGnTzHi+BYG9Alk2bJleHg0\nvzcge3aEvZHkI0QDjEYjFy5cICoqCj8/P6bMXYJh12HUT09DnZyAbuxYyv7zH/Z0GU3W7jTuLPmE\nMqU7P7ouZmRoZzw8TJuHkQewCXsjyUeIG1y5coVDhw6hUqm4dc4c1DEXKJ97F4NTT3Cx+yBWLP+O\npS/cgjbvEpcOfE1YpZJo9UzSHfr/8njrmxcZbW/ay5CgaD8k+Qjxi4yMDI4ePUpxcTETxo9nYEYG\n6ocfxuHYMdK8evH7Jf/i+xFLCTQk0mnjt3RyUTIsYgw/XAglvbT6PyVz99rYcrVpGRIUrUGSj7BJ\n1vqlbTQauXr1KtHR0Vy7do3Ro0YRmpOD6+9+h8OxY5R19+PEHz9g75T7OHM2jTvKPkWrUOMzdBLz\nRweiUChY3bOMhAIV2iotoT6qBvfa3Kw9trxnR4YERWuQ5CNsjiV/aTd20zcajaSkpHDixAk0Gg0j\nR45kaHo6rs88g8PJk1T5+vPu0g9ZHf4IviQz6sS3TFQ7sVM9B6VXL54cUoZCUX2umtVfSUmXCPSo\n/1TL5rZHqk0LeyLJR9icml/a1z8xNCbHAdCZlIAauul/cksxZZkJnDhxAoVCwchhwxhy7hzq5ctR\nxcaiDexJzMsfEDfnHjZHZTOn4iv0CgeOuMzksUl+PO8IA7uUmdQrae89B1seEhTtlyQfYbMeCKlk\nW7LTr08MNbEHdP1N39FYQefcGLb87xg+XToxKTycoD17cI6MRJmRgT4oiPR3P2aZ6304FicyKepb\nhla6cFJ9C2kOQaBQ4O6oYXrv9pEwLMmWhwRF+yXJR9icml/aLirqPjHUjB6Dh/4awZXH6V91hnSH\n/gzuFcZte7/F6be/RVFejm7yZMrff5/SiRPZeiie0edXc03VHUXwTGLy+3NV0/Jf+x2h5yBDgsLS\nJPkIm1PzS7t6qM10BoOB1NRU8mLOMk+TxQWHYWRnDmT5oY8ZHbcbo7MzVYsWUfn44+T7+XH69Gku\nfPUVLt36sdt9KQUqHxzTq4f8enpW4uZgbNGv/Zv1HGQZs7BHknyETaq+Aeua3WMwGo3k5eVx/vx5\nEhIS8PDwYGj//sxJz8Lr/z2Fx5WLVPXwpeKVV6i8/35SS0s5ffo02YcPExISwr333ss1gwfrd7tR\nUAJVBgU7UhxZPVNjkSGmxnoOsoxZ2CtJPqLNNfbLvzlzDcXFxSQkJHDhwgW0Wi2DBg0icuBAfP/3\nPxyffbZ6aC0igrJXP0czZw7nL17k9LZtGBQOuPcaTp/QefTspsTd3Yg71p/baO+LEYQwlyQf0aZu\n9su/oR5DWVkZSUlJXLhwgWvXrhEUFMS0sWPpFRWFy5/+hOrMGYxublQtXIj2gQfI7dWLs2fPcv7r\nrwkICCB0zHRejh1AeooDpNT9ztaa25ChNSHqkuQj2lRzf/nr9XqSk5M5f/486enp9OnTh1GjRtG7\nrAz1mjU4/fe/KIqL0QcHU/7uu5QvWMCl3FzOnj1L/s8/ExwczJIlS+jUqRN7Uh1I1zjc9Duv15Lk\n0VSC7QiLEYQwh1nJJyoqio8++ojY2FgyMzNZtWoVS5cutXRsooMx5waem5tLfHw8Fy5coEuXLgQH\nBzN70iTcduzA6ckncTh2DKOTE1Xz5qF98EHyBg3iXFwc59evx9vbm5CQEPr374+Dg/m/s1o6L9NU\ngpVlzMJemfVfZFlZGSEhISxZsoQnn3wShUJh6bhEB9PYDbzhJ2MWExNzgfPnz1NRUcHgwYNZtHAh\nXVJScPz8c5y++666l9O/P+V/+Quau+8msaCA+Ph4rp07x+DBg7n77rvx9vZuMBZTexutPS8jy5iF\nPTIr+cyYMYMZM2YAsHz5cosGJDqmxm7g03vrWD1DQ1xmJcWZF9HlXGDPxmz69u3LxIkT6Wkw4Lxh\nA47PPIMqMRGjiwtV8+ZRec89XOndm/jz57n0/ff4+/szbNgw+vTpg0qlaiKS1ts0qVX7sCe1+j+p\n63t2MrQmRH0y5yPaTFV5CadPJ3Dp0iWys7Pp1asXQeEh9HUbjeu2bTi+9x4Ox48DVD8/54MPyJky\nhfNXr3LhwgUcU1MJDg5m/PjxuLm5mfTdpvQ2mpM80ooVvHDMt8GhORlaE6I+ST7CKgZ66wlw11Fa\ndI2eVRcIMlwg5UAByr59CAsLo7enJ+odO3BavRrV4cMoDAb0wcFUvPoq+XPmkFheTkJCAiW7djFw\n4EDmzp1Lt27dWjzk25x5qOYkj5sNzcnQmhB1SfIRrUqv15ORkUFycjJ3FiWj1Rnw8OvPkIFjGebt\ngvOOHTh+9hkOBw6g0OnQ9+1L5e9/T8Gtt5JgNJKUlMS1gwfp168f48aNIzAwEKVSefMvbgZTFhJI\n8hDCsqySfJKSkqzxNTbFHttcIz4+ntzcXLKzs8nNzcXV1ZUePXowPHwo3mVleP30E14f/gWPU6dQ\n6PVU+PuTt3QpaVOnkqJWk5WVhebIEXr06EFAQADh4eEolUq0Wi2XLl2yWJxJ9KvXWzmbraci2/Tv\n8FX7EODu8msic9fj63CNpKQsi8Vry+z537s9tj0oqP6jQ0xlleRjiUDbk6SkJLtrc2lpKcnJyZw9\ne5aioiJ8fX0ZPHgwt912Gx5ZWTjUzOGcOAGAfuBAKp95hvyZM0l0cCAxKYminBz69evH1KlTCQgI\nuOnCgYaYspz7cmr9f/5Ojk4E9Tbv2r0zJpNMnVf1d3fRE+jhAXg0+/PtdSOqPf57r2HPbW8ps5KP\nRqOp/QVqMBhIS0sjNjYWb29vAgICLBqgsE1Go5GcnBySk5NJSUmhqKiI3r17ExAQwN133onruXM4\nbNqE4/btqC5cAEAXHk7FK6+QPWUKF41GLl68SOGpU/Tr14+xY8eanXBqNGcY7fobfICHwaKr0JzK\ns5ge1PxkY2rsQnQkZiWfmJgY7rjjDgAUCgUrVqxgxYoVLF26lFWrVlk0QHtkyV/AljxXZWUlV65c\n4VxiKlevpKJ0dKZXr95MnDgRf1dXnA8coOybb/B++GGUeXkYVSr048ZRdv/9ZIwfT1JpKUlJSWhj\nYujbty9jxowhMDCwRQnnejeb9L/xBt+7k473ppSRXlI9h9SWq9CkxpuwN2Yln4kTJ3Lt2jVLx2IX\n0ooVJBYocXeCggoFjsq6ScHSj5Bu6lw3S0w1vZvLly9z+fJlcnNz8ermx8HSAZx1nErPrCt4H95K\nn/S3cY2JRmEw4NSpE7oZM9DOmkVqSAiXsrO5dOkSypMn6d+/PzNnzsTHx8ekVWqWSqA33uBTixxI\nL1G26AFx18fmq/Yx+zxC2BtZ7WZFNclgTt8qtiU7knHd5PR7U8sY1MVg0V/ATZ2rscTUiSKuXLnC\n5cuXSUtLQ61W07t3byIiIghwdSVhQxTemz7lgws76V5cPZleMHgYqueeQzNtGlGVlZRXVpKamorn\nL0Nq8+bNo0uXLqSXKIkrUBF3uflJxJRkbO3NnPVic3epfgSDGclRNqIKeyPJx4pqkoGLqqo28QCk\nl6o4muGAm6MWDycjvx1WAUCFHr6Nd8LDydjgzvkWx1KixMNwjR76K/TIvsyW/6WgMlTSs2dPevfu\nzaTRo+kcH4/DgQM47N+P6vRpxhmNXHP15qeBMzkweDanBo7irkHZuJZcIuf0aTp16kRISAjjx4/H\nw+PX+Q9ze3SmJOOb7cex9A2+Xmyl5v9QkI2owt5I8rEROoOCjFIVrx9R197Q/N0NvD2xvM5rpgzD\n3XizDXTX0k2fTkxMBlcuZrG4OAOAbIeeZDn0otfwYG7Xp+B4+DCq1atxOHYMRUVF9dzNqFFU/vGP\nnBgylffyfXDWpBJQdYkxyl04anszfPhwAgMDSU1NbXD1j7XmNJraj2OJG/z1w2xVhpbF2lB8Mscj\n7IUkHyuqSQYV+uqhtuuTTIUessqUdW7QV0uVFGsVZt+0uztX8PrgVBJTM9Hkp1OZlc3Zo53w8/Nj\nQP8+7KyagltiMqMuHeKJ1M+JeOMIKk0pAPrgYLQPPkjlxImk9+/Plfx8rly5Qk7qz4z39kfVtw/O\n3e6iwqkLE3pXEXjDIq8b52kUNO8mf+PnLN1backN/sbe2zMjyutcxwB3GSoTorkk+VhRzS/vxAIl\n4/10nMlVoTMoqNDD7ssOPBlWWe8zzZ2XNxqNlJSUkJGRQWZmJlevXuVaYREunX1w9fZn+MhRhHZ3\nw+3MGVRHjuBw5AhzTsagrPxliG/AIHSLF1E2YQJZwcGklZaSlpZGRmoqnQsLCQwMpN/QsXx9bgBX\nNM6QBQGl1XNVDSWeG4fY3ptSVi+JBHgY6gwnAg0MzZXazHDUjb23Vadc+Hi6hhJt9UXydbj2y94e\nIcTNSPKxsupf3tU32m6uhtqb6vwgLRipd4Pu10mPv7ueq7W9pOqbdlVVFTk5OWRlZZGZmUlmZiZG\noxFfX1/8/PzwChjEh3tV+B0/wciULQRejqLb1VgUBgNGpRL90KFUPfQg2jFjOBs4iLMF5Wjy0qjI\nTMe9uJjAwECGDBnCrFmzUKvVAOxJdahOPL9IL62eNxrUpe74U0NDbOklyjpJJMDDwB8Oqkktdqht\n62+HVZBdpuTRoRW4/PLxjFIVo/30NjkcVWVQUKJV1K6Wq65mIMlHiOaQ5GNhacUKkujH5VSHmy4O\naGgI6MZf+Yn5Sub0qcS5Mh99cSZVRVkc+OEqezQFdOnSBR8fH/r3709Q8GhKz6TS9dTP9Fq7HqcT\n0WzIrV6NpnFy41TvMege+QN+M0dytWdPMgoLuXr1KhlXM7iWVsJVVS8G9QnGvf8sunu7EtTV8jvs\nr2/vnlSH2sQD1QlKa1DwQEgl25KdyCit3nsT4K43ewWZpcmKNCEsR5KPBVlij46vaxWOZZnk5uaS\nlJzLxbRcKq/lcU3pQZ7Kn1yVH0uHBzHXMRvns2dRbd2K/ucYHOPO4aivAuBq1z5URUziS9eJxPYd\nTlEXD7oaM4hwvYz2wgWcrubh6u1PUP+huAbfygtHu/Po0Irqm36WstHYm3vzNfcm7eNqoKBcUZt4\noGUryCxNVqQJYTmSfCzIlBVdNXM0eXl55OXlkZubS15eHiUlJXh5edGtWze6du3KqO59+GZXKV3O\nn2dm2o8MzzzJoIyzKCur54eMnp7kDgznuynPcarPKK4E9sbRuZxZXS5TlZXFIEMUuboAKtz8CY2Y\nwNvn+nKlzAXyIUCr589jywFwUVH3pt9A7M25+WrVPlwpUPLnseW/bqJt4LiGEpSfh56CCtt+Kq6s\nSBPCMiT5tDKj0YhGo6GgoID8/Hzy8vLIz88nPz8fR0dHunbtSteuXavLzQwaRJf0dJzOnUN16BCq\n2FiUiYlM0VXf7LTundCGDqXqtkepGDqUrD59yHJw4HRyHilXc+hhuICT8hq5Bn/cfIOYN3YC6VXe\nKBSK2iG8K2UutbGll6goqFAQ4NH8oaOmbr6NPlCtgd5BTSI7l1c93ObjagAjhHSVoS0h7IEkHwu6\n/tf86LKd+CkySNqZy0WMdOnShS5dutCtWzcG9ulD94IC3C9eRHXqFMr4eFTx8SgzMmrPVdXdh9yg\noRTefyuq4cFUDQgkW6kkNy+PnJwcinJy8Nbp6N69O30CfdihHcuFCh+MChUBHnpGhVXf9Afz6407\nIb9+r8JRWT3PlFGqrLts2Iybvjl7eT6KcbHZ1W1CiNYjyceCan7Nn83WU5bVlb6dujOosgDPlBRU\nCQmozp9Hef48ytRUFMbqG6rRyQnDwIHoJkygasgQCgYM4Lx7d75JNOJQnou3Pgt1XhbdFVUE+HTD\n39+fYcOG0aVLlzoFOfsOUpBQUL1surEbdmNzMYEe1SvwVs+07k2/4WTlwPTeOhnaEqKDk+RjYYGe\nRiqyLjLsD/ejvHDh1yTj4IChXz/0YWFoFy2iMCiI3B49yHNyIu/aNfLz8ykoKMA1MxOjq44SrR8F\nTqH8rJpBsdKb98aVMa2JApjNmYu42ZxNS+czqh+V3bLekxDCPkjyaQ0KBbqZM9HNm8e1vn3J7d6d\nfBcXCoqKKCgooKCgAKeCArooFHh7e+Pv709oaChdu3bFycmJPakOfHDAzSKhNFQRurV6FYGexgYe\nqNZ470mWLgthvyT5WFjNPp99/YZQWZSNe2kp3k5OdZKMt7c3Li4ujZ7DUjfltnhAmSkPVJOly0LY\nL0k+FnT9zd5Dfxdevm58MrPS5Ju9pW7K7eEBZbJ0WQj7JMnHgq6/2ZeovCkphYQC8ybP5aYshOjI\nlDc/RLRXNcN3NWRORQhhK6TnY0G2NoEucypCCFslyceCrt/n4+ToZBM3exm+E0LYIkk+FhboaaQi\n+xJBves/zVMIIUQ1mfMRQghhdZJ8hBBCWJ3Zyefzzz9n6NCh+Pj4MGXKFI4ePWrJuIQQQnRgZiWf\njRs38tJLL/H8889z6NAhRo0aRWRkJOnp6ZaOTwghRAdkVvJZtWoVy5Yt47777iMoKIiVK1fSo0cP\nvvzyS0vHJ4QQogMyOflotVrOnDnD1KlT67w+bdo0jh8/brHAhBBCdFwmJ5/8/Hz0ej3du3ev83rX\nrl3JycmxWGBCCCE6Llnt1gqCgux3j4+03f7Ya7vBvtveUiYnn5onaN7Yy8nNzaVHjx4WC0wIIUTH\nZXLycXJyIjw8nP3799d5ff/+/YwePdpigQkhhOi4zCqvs3z5ch5//HGGDx/O6NGj+fLLL8nJyeHB\nBx+0dHxCCCE6ILOSz5133klBQQF///vfyc7OJjg4mHXr1hEQEGDp+IQQQnRAisLCQqmxL4QQwqos\nttrtq6++Yu7cufTs2RMvLy/S0tJu+plvv/0WLy+vOv/z9vZGq9VaKiyrMKftAJs3b2b06NH06NGD\nMWPGsHXr1laO1LIqKyt54YUX6NevH/7+/ixZsoSMjIwmP9Oer7mpJaXi4uK49dZb8fX1JTg4mJUr\nV1opUssypd2XL1+ud329vLzYt2+fFSNuuaioKBYvXkxwcDBeXl6sXbv2pp/pKNfb1Labe80tlnzK\ny8uZPn06L730kkmfc3V1JSkpicTERBITE0lISMDJyclSYVmFOW2Pjo7m4YcfZtGiRRw+fJjIyEge\neOABTp482YqRWtZLL73E1q1b+fLLL9m+fTslJSUsWrQIg8HQ5Ofa4zU3taRUcXExd955Jz4+Puzf\nv+82b6oAAAV9SURBVJ8VK1bw0Ucf8c9//tPKkbeMuaW0Nm7cWHt9ExMTmThxopUitoyysjJCQkJY\nsWIFarUahULR5PEd5XqD6W2vYeo1t9jzfJ588kkATp06ZdLnFAoFXbt2tVQYbcKctn/yySdMmjSJ\n5557DoDf//73HDp0iE8++YTPP/+8VeK0pKKiIr755hs+/vhjJk+eDMDq1asJDQ3lwIEDTJs2rdHP\ntsdrfn1JKYCVK1eyd+9evvzyS1599dV6x69fv56Kigo++eQTnJ2dGTRoEElJSXz88cf85je/sXb4\nZjO13TW8vLzo1q2btcK0uBkzZjBjxgygeoHVzXSU6w2mt72Gqde8zTeZlpeXExoaypAhQ1i0aBGx\nsbFtHZJV/Pzzz+26RNHp06epqqqqk2T8/f0ZOHDgTdvQ3q65OSWloqOjGTt2LM7OznWOz8zM5MqV\nK60ar6W0pJTWPffcQ1BQELNnz2bz5s2tGaZN6AjXu6VMveZtmnwGDBjAqlWr+M9//sPnn3+Oi4sL\ns2fPJjk5uS3DsoqcnJx6JYq6devWbkoU5eTkoFKp8Pb2rvN6t27dyM3NbfRz7fGam1NSqrHrW/Ne\ne2BOuz08PHjzzTdZs2YN69evZ9KkSTz00EOsW7fOGiG3mY5wvc1l7jVvctjtzTff5N13323yBFu3\nbmX8+PGmRwxEREQQERFR+/fo0aOZOHEiq1ev5m9/+5tZ57SU1m67rWpuu81ly9fckpo7Tt7ReHt7\n1xmqCQ8P59q1a3zwwQcsXLiwDSNrXfZ6vcH8a95k8nnqqadYvHhxk1/s7+9vYqiNUyqVhIWF2cSv\n4NZue/fu3RssUXTjrydra267dToder2egoKCOr2fnJwcxo0b1+zvs6Vr3hhzSko1dn1r3msPLFVK\na9iwYXzzzTeWDs+mdITrbUnNueZNJh9vb+96wyqtyWg0cu7cOcLCwqz2nY1p7baPGjWK/fv389vf\n/rb2tf379zNmzJhW+87maG67w8PDcXR0ZN++fdx9990AXL16lcTERJPKLNnSNW/M9SWl5s2bV/v6\n/v37mT9/foOfGTVqFK+99hqVlZW18wD79+/Hz8+Pnj17WiXuljKn3Q05e/YsPj4+rRGizegI19uS\nmnPNLTbnk52dTWxsLBcvXgTgwoULxMbGUlhYWHvMHXfcwRtvvFH791//+lf27dtHamoqsbGx/OY3\nv+HChQs89NBDlgrLKsxp+xNPPMFPP/3E+++/T2JiIu+99x6HDx+uXTln6zp16sS9997Ln//8Zw4e\nPMiZM2d4/PHHCQkJYcqUKbXHdZRrvnz5ctauXcu///1vEhISePHFF+uUlHr99dfr3KDvvvtu1Go1\nTz31FOfPn+eHH37ggw8+4KmnnmqrJpjF1HavXbuW7777joSEBJKSkvjoo4/44osveOyxx9qqCWbR\naDTExsYSGxuLwWAgLS2N2NjY2iXmHfV6g+ltN/eaW2yp9Zdfflm7qUqhULBw4UIUCgWrVq1iyZIl\nAKSmphIYGFj7meLiYp5++mlycnLw9PQkLCyM7du3M2zYMEuFZRXmtH3UqFF88cUXvPXWW7z99tv0\n7duXf/3rXwwfPrxN2mCOFStWoFKpePDBB6moqGDy5Ml8+umndca/O8o1v1lJqezsbFJTU2uP9/T0\nZNOmTTz//PNMnToVLy8vfvOb35i0dNUWmNpuhULB3//+d9LS0lCpVPTv359Vq1YRGRnZRi0wT0xM\nDHfccQdQ3aYVK1awYsUKli5dyqpVqzrs9QbT227uNZfyOkIIIayuzff5CCGEsD+SfIQQQlidJB8h\nhBBWJ8lHCCGE1UnyEUIIYXWSfIQQQlidJB8hhBBWJ8lHCCGE1UnyEUIIYXX/H2BUE0rgrkFMAAAA\nAElFTkSuQmCC\n",
                            "text": [
                                "<matplotlib.figure.Figure at 0x10f8862d0>"
                            ]
                        }
                    ],
                    "prompt_number": 14
                }
            ],
            "metadata": {}
        }
    ]
}