Created
April 26, 2017 16:53
-
-
Save jasonmhite/e60f7cd893e9fe4e3d66240a3c368561 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "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": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import pymc3\n", | |
| "\n", | |
| "from theano.tensor import dscalar, dvector\n", | |
| "from theano.compile.ops import as_op\n", | |
| "\n", | |
| "import scipy.optimize as opt" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The main difference is that you need to make your response model using Theano, a la https://pymc-devs.github.io/pymc3/notebooks/getting_started.html#Arbitrary-deterministics.\n", | |
| "\n", | |
| "This requires specifying the types of your function inputs and outputs.\n", | |
| "http://deeplearning.net/software/theano/library/tensor/basic.html\n", | |
| "\n", | |
| "I will be using Dr. Smith's (in)famous spring model test problem, see his book example 8.2. This is a single parameter, $C$, with a nominal value of 1.5." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# Some parameters and the problem grid\n", | |
| "K = 20.5 # fixed\n", | |
| "C0 = np.float64(1.5) # Just used to generate some data\n", | |
| "\n", | |
| "CMIN, CMAX = 1.4, 1.6\n", | |
| "SIGMA = 0.1 # amount of noise used to generate data\n", | |
| "\n", | |
| "NT = 500 # Number of time steps\n", | |
| "TMIN, TMAX = 0, 5\n", | |
| "T = np.linspace(TMIN, TMAX, NT) # grid of time domain points" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def response_model(C):\n", | |
| " return 2. * np.exp(-C * T / 2.) * np.cos(np.sqrt(K - (C ** 2) / 4.) * T)\n", | |
| "\n", | |
| "# Turn this into a theano op\n", | |
| "response_model_op = as_op(itypes=[dscalar], otypes=[dvector])(response_model)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Generate some made up data" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "Z_nominal = response_model(C0)\n", | |
| "Z_data = Z_nominal + np.random.normal(loc=0.0, scale=SIGMA, size=len(T))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl4VOXZuO/nzCQk7BD2EEjCjspmDJtRVLQIVGq/2qq1\ngrWltdpqtbbS9uv2a6td7GettW6folWg9qsLyqKiIkGWEPYdQhIIe9i3bDPn/f1xZjKTYcKWSWaS\nee7rypXMmTPzvudk5nne91nFGIOiKIoSf1jRnoCiKIoSHVQBKIqixCmqABRFUeIUVQCKoihxiioA\nRVGUOEUVgKIoSpyiCkBRFCVOUQWgKIoSp6gCUBRFiVPc0Z7AuejQoYNJT0+P9jQURVEaDStXrjxk\njOl4IefGtAJIT08nPz8/2tNQFEVpNIjIzgs9V01AiqIocYoqAEVRlDhFFYCiKEqcogpAURQlTlEF\noCiKEqeoAlAURYlTmqYCKMmD3Ced34qiKEpYYjoP4JIoyYNXbwFvJbgSYfJsSMuO9qwURVFijqa3\nAyjOdYS/8YKnAuY/Bu//UHcDiqIoITQ9BZCe46z8sQAb9qyE/Jdh+kRVAoqiKEHUWQGISJqIfCoi\nm0Rko4g8GOYcEZGnRaRARNaJyLC6jlsradmO2afXGEACx70VsPBxVQKKoig+IrED8ACPGGMGAiOA\n+0VkYMg5NwN9fD9TgX9EYNzaScuGMdPAlVDzeOFCxz+gSkBRFKXuCsAYs88Ys8r390lgM5Aactok\n4DXjsAxoKyJd6zp2bWxds4SiZe9An5ug/0RIvRLEAmM7/oHi3PoaWlEUpdEQUR+AiKQDQ4HlIU+l\nAiVBj3dztpKICPM/y2XCrFJ+sqo9ZvP7sP0jGHo3uJqBuBz/QHpOfQytKIrSqIhYGKiItAT+Azxk\njDlRh/eZimMmokePHhf9+hGefNrQmTwzgNn2SCZ5l0HZYccvUJzrCH8NC1UURYnMDkBEEnCE/xvG\nmLfCnLIHSAt63N137CyMMS8YY7KMMVkdO15QT4MatO07mh8lOlP4h2cSRtxwfLfzZM4jKvwVRVF8\nRCIKSID/BTYbY/5Sy2mzgbt90UAjgOPGmH11HTssadl8+ZuP0jmxnC2mBwvtIbDyVXX+KoqihBCJ\nHcBo4BvA9SKyxvczXkS+KyLf9Z0zFygECoAXge9FYNxaaZY+nHvHDgHg1arrnKQwdf4qiqLUoM4+\nAGPMYmoE3Ic9xwD313Wsi+G/upTyB7zk2ldwyLSmg6tCnb+KoihBNL1MYB8pBz7nGms9XlzM8Y6A\noXeq/V9RFCWIJqsASM/hSwmOzX+2PRoG3xnlCSmKosQWTVcBpGUz9u6fkmjZrLL7cKT9kGjPSFEU\nJaZougoAaNFrOMN7dcIAi7aVRns6iqIoMUWTVgAAY/p1AuDTrQejPBNFUZTYoskrgOv7Owrgs40l\neHeGVqhQFEWJX5q8Asgo20ialHKsysXm6d/XZDBFURQfTV4BUJzLCGsTAMuqemsymKIoio+mrwDS\ncxjh3g7AMjMAklOiPCFFUZTYoOkrgLRshl9/KwB53n54501TM5CiKArxoACA7lYpqVLKCVqwuaoT\nrJ0R7SkpiqJEnbhQAKTnkG05ZqDVdm9YPUN3AYqixD3xoQDSshnSsz0Aa+xeYHvUGawoStwTHwoA\nGDxsBABrTW9tC6koikIcKYABQ0aQaMEO040Tt7+jlUEVRYl74kYBNHO7GJDaFoOwnt7Rno6iKErU\niRsFADA0rS0AazZsgNwn1RGsKEpcU+eOYI2JwWltAFizYjEk/sXxBUyereYgRVHikrjaAQxJawfA\nGm86xtY+wYqixDdxpQDSz2ygDacopR37aQ+WW6OBFEWJWyKiAETkZRE5KCIbanl+jIgcF5E1vp9f\nRGLci0XWzeIyqxiATXZP6HOjmn8URYlbIrUDmA6MO885ucaYIb6f30Ro3IvEMEB2AbDZ9ICWHaMz\nDUVRlBggIgrAGLMIOBKJ96pXBt/JQNduADaZTG0UryhKXNOQPoBRIrJOROaJyGUNOG6AtGwG3PoT\nADa3HKHmH0VR4pqGUgCrgB7GmEHA34B3ajtRRKaKSL6I5JeWRr6Re+8rRpDgEopPeDld4Yn4+yuK\nojQWGkQBGGNOGGNO+f6eCySISIdazn3BGJNljMnq2DHyNvpEt0XvTq0wBrbsPxnx91cURWksNIgC\nEJEuIiK+v7N94x5uiLHDMaBrKwA27zsRrSkoiqJEnYhkAovITGAM0EFEdgO/BBIAjDHPAV8B7hMR\nD1AG3G6MMZEY+1IY2LU1b7GHTaoAFEWJYyKiAIwxd5zn+WeAZyIxViQY2LU1oDsARVHim7jKBPYz\nwKcAtuw9hr1Ii8IpihKfxKUCaNcikY7NLco8sOfj5+HVW1QJKIoSd8SlAgDom+xEAG23u2pROEVR\n4pL4VAAlefSREgC2mzRtEakoSlwSfwqgJA+mT6DPscUAbDPdYdwTmhWsKErcEX8KoDgXvFX0sfYA\nsN3uBmVRS0lQFEWJGnHVEQxwTD2uBPoYRwEUmFTsnt3jUBMqihLvxJ/cS8uGKXNod9VtdEio4AxJ\n7Gl5RbRnpSiK0uDEnwIARwlMfIq+nVoCULBxZZQnpCiK0vDEpwIAJxLo4HwAtn30kuYBKIoSd8Sv\nAijOpQ9Od7Dt3i6aB6AoStwRvwogPYc+7oMAbDepkJwS5QkpiqI0LPGrANKy6Tv2XgC226mYuY/C\n+z9UU5CiKHFD/CoAoJ19iA4cdyKBvK0h/xWtC6QoStwQ1wqA9JyghLDugNG6QIqixA3xrQCA3tY+\nAHaYbiCW1gVSFCVuiG8FUJxLpuwFfAogcwxMnq11gRRFiQviWwGk55DpKgWg0HSDMdNU+CuKEjfE\ntwJIyybztt8CUOjKiPJkFEVRGpb4VgBAaiuLZlRSWtWMk9O/qhFAiqLEDRFRACLysogcFJENtTwv\nIvK0iBSIyDoRGRaJcSOBtXMxGbIfgEJPikYAKYoSN0RqBzAdGHeO528G+vh+pgL/iNC4dSc9h0zL\nyQguJE0jgBRFiRsiogCMMYuAI+c4ZRLwmnFYBrQVka6RGLvOpGWTOeQaAAqveFCdwIqixA0N5QNI\nBUqCHu/2HTsLEZkqIvkikl9aWtogk8vs1QeAwoo2DTKeoihKLBBzTmBjzAvGmCxjTFbHjh0bZMyM\nDi0A2LH3IOQ+qY5gRVHigoZSAHuAtKDH3X3HYoLMjk5jmOLDZ7A//r3WA1IUJS5oKAUwG7jbFw00\nAjhujNnXQGOflzbJCXRI9FBOM/aatloPSFGUuCAiTeFFZCYwBuggIruBXwIJAMaY54C5wHigADgD\n3BOJcSNJZkoyh/ZVUWhS6e46rdFAiqI0eSKiAIwxd5zneQPcH4mx6ovM7l3I21dCUfdJXNNtd7Sn\noyiKUu/EnBM4WmR2dBzBhSV7YOWr6gdQFKXJowrAR2YHxxFc6O0Exqt+AEVRmjyqAHxU7wBMNxCX\n9gVQFKXJExEfQFMgrX1z3Jawx06h7Nqfk9zras0Kro2SPGd3lJ6j90hRGjGqAHwkuCx6pDSnsPQ0\nRX2/xcBuraM9pdjCL/STU2D+Y46JzJWoDXQUpRGjCsBPSR6Z1n4KaUXhoVPxqQBCV/b502Hzu9Bl\nECx/3hH6ImBs58fvJ1EFoCiNElUA4Ai+V2+hV/mXWcAECgu2wqBu0Z5Vw+K7B9Ur++Hfgc+fcp7b\n8QkggAFjgWU5j9VPoiiNGlUA4KxivZVk+KpTFO+JmSTlhsN3D6ojoDbPrvm8CGCB5YI+N0HLTjD4\nDl39K0ojRqOAwFnFuhLJsA4AUFjWPMoTigLJKb4/BCw3DLil5vOjfgBX3u08v3UerJnZ0DNUFCXC\nqAIAZxU77gkyXL7GMEcqMbuWR3lSDUhJHsx71Fn9Yxz7fv8JMPGv0Ot65/eNv4Y2aWB7auZJlORp\nBVVFaaSoCchP2WE6mmO05AwnaMGRbUtJ6TE82rNqGIpzwVsVeGx7qCrMJbfzXazsNIKqg4b05bsY\n23EUnVyJAT9BcgpMnxh4POV9NQkpSiNCFYCf9BzEnUiG7Ge9yaTo4DFSSvLiQ6Cl54ArwRHkwCIz\nlF/m9qboVH6N035pweQrZvBwlzU0b9UeVr8G3grnSW8FrJ0ZH/dLUZoIqgD8pGXD5NlkvrmO9aVQ\nuHk1WYXPNe049+CwzylzYO0MXtrTk98V9cFUCOmyn/HWclrKGVba/fjEHsJLa8v4fPdAXix7hO52\naNE8E5XLUBTl0lAFEExaNhltN0EpFNmdm3ace2jY5+TZvN7hh/x28QYAHnL/h/td75AgXt8L3mOD\nnc73+QmbD7fhdnmUfyf+mq7iawUtLugyJDrXoijKJaFO4BAy0jMBKDLdmm6ce0keLHzcMdv4HLrL\nV+bzi3cd4f944is85P5PkPB3uNwq5h3rUQZbO9htOjK58jFOSUtH+GOcDGF1BitKo0EVQAiZ/QYD\nUNTiiqZp/vGv/HcsdKJ9xOK41ZaHNmRgG7iv12HucH0SOF/8HxEBhDac4tXEP9Er+RTbTHd+3Pz/\nYQw1M4MVRWkUqAIIIb2DkwNQdCYJuyi36a1o/Qlf2IAFmWP4U8aL7DtlM6RtGQ8Ptp2dj7jAnQwT\n/gdu+AWMftBJAhOLtu4qXrg1jVbN3Mw93Jm3zLVaQVVRGiHqAwihVVICHZtblJ6x2fvxc3RP+FPT\n2gn4kt7wVIBlsaH5CN5YUYYbL38s+zUJH5bCuCeg7HCgJpB/12DbThmIcU/Qa9Aoflm1mx/9ey2/\n5tvkjLqWTv1HNZ37pChxgO4AwpCRdAqAIrtT0zNr+JLesCywvfxhpY1BmOL6gL6yy7ne/WtqviZ4\n12CMoxyA/xqWynX9OnKiwvDzfVo+W1EaG6oAwpDZRgAoJLVpmjXKDoMxrLJ7kWsPoiVneMD9Dk4Z\nCBesngGf/C7QFtO/a/CbeZJTIPdJZPcKHv/yIFo2c/PhpgMs3n4o2lemKMpFEBEFICLjRGSriBSI\nyGNhnh8jIsdFZI3v5xeRGLdeKMkjc897ABSZrs5quamtbNNzwHLxtOfLAEx2fUhbqxyy7nEKvQUX\nhfOHwU6eDdf/zLkf8x+rVhBdTqzje9f1AuD/vfkZnp1xVEJDURo5dVYAIuIC/g7cDAwE7hCRgWFO\nzTXGDPH9/Kau49Ybxblk4CQ4FdqdzzaHNAXSslnX+z4W2kNoTjn3uuc5x7sMhu0fUp3QZbkDu5+0\nbMh5xNk9hCiIb/Y8THcpZevJJP718l+cXYPWCFKUmCcSO4BsoMAYU2iMqQRmAZMi8L7RIT2HDFcp\n4NsBrJ7RJIXYK2XXAHCXawHt5SRgnOYvtj/2X2DonWfvfkLNQek5JO1ezGPuWQA8XTGR8tWzHPNR\nsBlJUZSYIxIKIBUoCXq823cslFEisk5E5onIZbW9mYhMFZF8EckvLS2NwPQukrRsegwdi4XNbtOB\nCi9NywkMHDldyZwdVQiGbyR84sT6u5rBgElBIaBJMPjOmi/0l44Y94RTGnrIHc7x9BzGJ65lgOzk\nAO351/7Us81IiqLEHA0VBroK6GGMOSUi44F3gD7hTjTGvAC8AJCVlRWV4jKJQ28nbekGdprO7LJS\n6dOUnMAlebz58UYqvZ24rl8n0sb+s2YbyM4Dwzd8Dy4dYbkAcUpDr5kJk2dj3fw4Dy5fy3dLevKP\n3enc7k6imSlzGslU9xpQFCWWiMQOYA+QFvS4u+9YNcaYE8aYU76/5wIJItIhAmPXD2nZZKQ5l1R4\n7dNNxwlckoc9fRIztjp69a6UgoBt33+NoY/91OgYVlVzhb92Bsx/jJsOvkp/2cX+yiTerLo60D9Y\nS0QoSkwSCQWwAugjIhkikgjcDtToJygiXUREfH9n+8Y9HIGx642MNKcncFFRQdMRXsW5LKnsxS7T\nmVRKGbP6BxfusK1h+09w/sZyhPypUvBWYomXH7jfAuDZqolU2JaWiAB1iCsxS51NQMYYj4g8AHwA\nuICXjTEbReS7vuefA74C3CciHqAMuN0YE9O1gzPdTkx70fZNsHta08gGTs/hbdvxq/yXaxEuvM7q\nfc2sGlVBw16nPxTUbx46sAnmPuII+O0fOhFDXptx1gr6SQlbTRrvmRy+Igvj2wwUpurqJX2Ogkt3\nN/bPoRIzRCQPwBgz1xjT1xjTyxjzO9+x53zCH2PMM8aYy4wxg40xI4wxSyIxbn2SWb4JCCkL3cgp\n73IlH8hoAL7kXuo4fpELd9gGm4d8yWQY24kc6jMWLBeWwLcTnLDSl5KnYMSKbzNQDdPZJX6O/EpE\no6qUCKOZwLWQ0c+pbV/YhMpCL9h8gFMei0Etj5N51ThnNTr4jrPCOi+I0HDQlp3xlQXli9ZSOjbz\nsOW4m8+9A+PbDBQmbPaiiYQSUZQwaDG4WujSfwRJ7rkc8rTh+O3v0Kaxb7tL8nhn/kagE5Mq3oM1\nnzrCP9S0c6HXGfo6cCKCvJU0swyTOxTz5529eck7nqtdG5uMEr1oznd/g007EP48vxLxm5Hi8T4q\n9YLEsik+KyvL5Ofnn//EemLcU4vYsv8k794/msFpbaM2jzpTksfxV24j68zTeLFY1ux+OlmnnNIO\nOY9EdBzWzoBVr3PUm8TIiqcppxkf9fo/+gy9tvYQ03gh1I7vN+14Khw/iTgF+rAsGP8kZE2p/bWK\nUgsistIYk3Uh5+oO4BxkdmzBlv0nKTp0uvEqgJI8mP8Yn1ZdRhVuRlgb6SQnwJVUPyvJfevArqKd\nVPEV1yJe997I/xan8MT+R6nOHaiLM7SxEs4ZXJzrCH9sp/qG8WVh27bjYO88sGZ4bjzdL6VBUB/A\nOcjo0AKAwkOnozyTS8QvdPas4kOvsyD4gpUPqcMiL4CDxvLzTdc8BJu3vKM55EmObzt2ODt+eo6z\n8g+HsePvHikNjiqAc5DRoSUARY1VAfiETrlxs9B2Wl3emLC+fiqcVvcMCJgUM92HucG1lkoSmWmP\nrbsztDFTmzM4pXf4862E+LtHSoOjJqBzkNnR2QEUHToV5ZlcIj6hs8Q7gDMkcVmLE3Sf8mr9mBKC\nHZWWyykr3bIT90gGC3Lhnwm38d1vfIeEksXxaccOdgYnpzi+ktUznKzqcPQZG/4exYsvIF6uM8qo\nAjgHmT4TUFHpaYwxSG3b9VjFJ3Q+nL0dSuCmkVmQFrYEU8TGqhZw8x8DbyWjrJn0bT+dbUe8zD3W\nnUmRdDpHG7+QSk6p2UKzNnzP2dMnUe6xceElERsRC1KHVvtPANi+APKn+8qRixOxFZx852pW04zX\nlARmpJLnLmScpnLPLhFVAOegbfNE2iUJR8u9lG5dRqf+I6M9pYvGm3oVC44eByr5wuWd63cwv6My\n98lqe7fYlUxJ3c1Pj3Rl+vylTErp0TS+bMERPNiBiqphhPKZ1Kv56ER3PttWypqtu9l9+lkqSQSg\nJWVkWAcY2rw/16UvZdSOp2gmvlpL7z9EtUlt1T8dwe93FHvKYeHjMGaa87ghBGZDEc5fEuns6YZS\nMjGOKoBzUZJHRuV2jtKbHbMeo9O9f2p0H5LVu45y6FQlPdo3p1/nVg0zaEjc+q293Pxh/WlWH2vB\nmpcfZMg3/9ro7uNZBPdJBkc4eyoCQvnAJo7M+SXPV41npieNE/jbZTYDIIkKvLg4RTLr7XTWry/n\nNYbSnme4zbWQe93z6CTHAuPZoaYiA4ULYedSpyx3JARmrBCJvIdQAT/uCWeXVn4C9q+DhOZN655d\nIqoAzkVxLhkcZBW9KfK0Z6T/y92IPigfbT4AwI0DOzecCSsk+Sl57Uxud5XxvPeLTK+4jqcWPu70\nHrgQs0ms4hdS/h0A4vwuXIhd9Dn/rBrDk1V/4gSOGXFIq2OMvyKVUat/TIa9ixauKhj/JEcG3Mm2\nAydZsuMwH27cz5b98Lz3i7zqvYkprg/4rvs92sppwHJ8K7bXGcpfhsPveG9KiWKXmpwYTPAuwlMG\ncx7yZaoHYSU4vy/mnjUxs5EqgHORnEKmtRZsKDJdYMdMZ8XVWLaLJXl8tqoESOL6/p0a9sPrNweV\n5MHq1/mGuzUveicwxx7OTwveoNOOT8KbTRoLoT6Pze9C4UIO2q15uOo+FttXAHC1tZ4fud9kiHcX\nlPQCsw3EgC2wfw3tOw9kxJ5cRrRL4YfDDrM6eQTPrz7DB4XwnPcW3vSO4ZcJr3GLtQQx4vgK2veC\njW/5fAGJTuOewXc2DcEU/BkN9ReFS6QLvWb/sfITNQV+uITXroOg/4Twr48Ts5EqgHNRdpgM2Q/4\n2kNinBVfY9guluRx4JW72XLmSZKpIOvw+/DRYw3/4S3OBdtDdznETVY+8+1sXvfcwMMJ/6lZIyjW\n72cw4YRU54GsK9zLvRU/oJS2tOckv094iS9YK5xQfxso3RL0JgbyX4VVr/vacDrO4GGuZjw/eTZr\nN2/j90tOsbw8jQerHuAdazR/SHiRTntWwp6VgDg7An9Ib1MoEHcuARvOpDM/5PMMMH0ieCsubLyh\nd5+dbX0uAR8p30QMoXkA5yI9hwzXQQAKTVffQbtxlDYuzmVRVT8ARlqbaLbt3egkYiWnOIIeuMc9\nH4AZ3rFUGHfNcxoLtVTmXHCyJ1+r+gWltGV4twTmT6hkXA8bcSXg2GzCYfts+0F+BF+DncF5P2KW\neYw/uF+gFaf51B7K+IrHWeS9wvda46xqyw5Ht1popHodlOQ5/hNvRfjPaKjw3Rz0efb7Xj7/6wUK\nf4HRD9UU/uHGCP2OnKuwXyPt+aA7gHORlk3GzQ/B27DLdMJjLNzi+9LFOskpfGYPAuBa1zroMsgx\nXzW0nbjsMM46wyZbtjJAdrLZ9GSOPYIvuxYHSkX7yx7Euo117UwnAgdTLSQ+PNGD772xEo8NX+lZ\nzu/bzCTxk7nOOZYbeo5y7r1f0NeGWM7/5lQpeMoRMXzNvZAxrjX8sOp7LLEv5+6qaXzPfpeHE/6D\n2+VLFovWyjQSJhF//ajqnAg7cB+CP6PBPheRwOfZU+68Zscn1L6eFV93OuPsmkLrLIWOUdt3xG/2\nWzvDec8Dm84Ke25spiFVAOchafgUUj+YzZ4zbnbTmXT3sdh3spXk4Z37Y3K9fwPgWlkDyxbCzX9s\neMdreg64m4G3EnEnck+mlx9vhlc847jVWoyIqbnaimUbq8+fUR2aabn5+GhX7p+7Ao+xmOp6n2n7\nZyAHgl7jrYRdS53XiAv63QxlR2Hn54FzUq90zBFlhx1hMu/RoDES6Nz3av7ZYj/PHsvgfzY251nv\nJNa1G8szk3rQ1n9/ouEEvljFU2sxPJ9CBcCCzDGBYIvg14x7IpAHsfx5GP4dZ9VfjV/B+sxjfcdB\ny05ODoV/vufyH1yo83nNrJrhv/6eF43QpKkK4ALI6NaZPQWHKBr0EOnZw2L/n1ucy1pPD47Tkp6y\nn3TrAHjFETANnYgV4iy95dQRnthhs74yk1VyGVfKloDQinUba3Guz14PICztcif3LWlOFRb3uuYy\nzT0jfGkf4xdMxqnDlJ4D0yc4K15XQs3SHLlP1hiDYXfBxKdwAd8Hrio8zAMzVrH4EEz691FeHPI3\n+l4x/PyCqz52VhcTrhm8W7BcMPQu53iN8iHiLBb8uQ3vP+TsDPwFBIfcUTP6af+6oNcGY2DkA3Dj\nr2se9l93/vTaE+qCi+7lT3dMTQMmBXYM4cJ/jQlUc70QBRxDu1xVABdARocWLC44RGHnL3BdWka0\np3N+0nP4zDgOx2usdc6xaIYH+j/kr95CkreSO8zX+DsTeaX7b7iyy/tU28hjse598Jc1aH4F0oPv\nFI+hkgTudn3Iz92vhxH+Qg0BJVbgSz9lzoXV/h98Z413HJGZwrsPXM3Ulxax8ZCHWxd15X+W/5yb\nrros0N8h3DXUx87qYsI1ayh3L+S/4ig/y+3IUr9S8K/WQ3cGnnI4dQDbasYuTxu2ksG2snHsrcrk\nsGnFYdOacl9yHUDSkua025tP+xYJdOYIGVUFpPceSO92LlrPfcRRKuD4DPy7z+C+DJ8/BVvmOH/v\n+MT5nTWllvBf40xTLOh9w7nvWYxFEqkCuAD8VUEbU02gRUljoAqu7dEMUr9Zu3BoKIIEwF2uD3nO\nM4F5OyrYt/cDupoDTjOZybPrHv8dScJ9WSfP5vDWJdyT148TZ2xucq3ilwmvI5Y7sDoVyzH19L7R\nsQ17KgI1/s9X3vkCGsikFufyf3328uiRZN63RzK14kEeWfZvHlh9CzIljECpERNf7gi31Csjc4/P\nV6Y6uFyGKzFIqBtnp3Pl3dAmreZc/JnkGIyB7SaVxfYVLFk/kOVyByc9Plt/IUD38OOWAZuDbXHd\nYK2TWJcpTzBEdjDY2kGWu4ABSSlY0yc6Y4rlzM2E+GtWv+YogFrCf6uztLfMgYKPA4I9+Pr3r3HK\nfXgrYsZcFBEFICLjgL/iNIV/yRjzRMjz4nt+PHAGmGKMWXXWG8UoGb6icIWljaAqaEkeJ6Z/jbWn\nn8aNh5Fj/wt6DY/2rGqsbLu6TjOuexJzCip4vTKHR93/CnwZ/D2HY4EagtOJNPFe8xjf2zGckhNH\nGNSmjKeG2rhKr3HMBOEa3lxKE5zahGqQQkoW4W8JHgZ4d/Jnz1d50nMbm+2e/KlgMS1CX1sjyso4\nQmrrPMf84c+QPZdt/EIJrY0U6hwd94QjBIPNOoPvPHuc9BwKpCfvVQ3jfe9wdpjUGk93bt2MfqaI\n/mdWkSYH6SAnaW+dooU5AyKYQbdTNuzbHDldydG177N383IK7c4Umm7sII1C041C04237BzwQPt3\nKxlpT2W0tZHR1gZ6WgfPvrZ9a53rC55r54HO76JF1e1QnVDxcp+jmJrlQoK5UHNRPVNnBSAiLuDv\nwI3AbmCFiMw2xmwKOu1moI/vZzjwD9/vRkF1UbjGUBa6OJcVVRnYWAy1CmixtwASJfqr6pCV7Tft\n3swpWMoMz/V83/U2SS6i/mU4i9DtfuFC/ra9I8srJ9FRjvFS+c9ovuyo82X2JwjmPBIICfTf70jd\n82CFZAQrL7bXAAAgAElEQVQRuN89m/5SwkNV9zPXzqYw38ULHZfS4+iSwKpz9YyzV7T+0hVzH3GE\nV22x9Rc693C1kUKdo2WHYeJTtSatlRctZ/7ydbx+oAf5Z35bfbwdJxnjWsso9zZG3TCJ1L3zHSXm\nS+TFSoDxfw4f4NBmEBT9vPqaKu96ly2HPazdvIXVVT1Ztl/YexzmMJI5tlPrq4ccIMdazzXWOkZa\nG2ktZY7Cmv+Y46z33yPLhdPkyG8KgmolsNoXKRTsLwimXTqMejDqi51I7ACygQJjTCGAiMwCJgHB\nCmAS8Jpx+k8uE5G2ItLVGLMvAuPXO93bNSfBJew7Xs6ZSg/NE2PYcpaew1KzE4CRrm2QPDx2bI5B\nwnCYMVzR0c360lbM9o7iq+6lgbC6WDD/QEBpLXwcCheyzNuXpyu/iGB4yv0sneSoc15o0/v6ut+h\nJbdtLxgvN7hW87b8gqny32w53IYvvl7CMwnvkONax1l+CD/hBPTmMLkiFzr3sM5R24l8Ok8PiNKT\nFbwy73NmrTrAEboBHlomCOMHpzJxUDdGJhaSULIHkvvB/B87pR2C6ToofFgnnLXwSEzLZlA6DOqU\nwDfWzsB0EIpaDuPzT99jiac/S+zL2GU684a3M294x+LCyzDZzjWudeSUrOeKPY/gEv89C3Luh2J7\nqC7REW4HcLTYifbavzaq5tlISLJUoCTo8W7OXt2HOycVaBQKwGUJPVNaUHDwFMWHzjCwW+toT6l2\n0rJZ1mIvHIMRI66Gsv0xGVkjIkzptotHSrvxivcmbvMsQoJXo9EOAQ02hYyZxtGi1TxUeT82Fg+4\n3ma0a0PQxQRt5+szkinUP3BgE8x5GIyX3tZe3rGm8VDrn/HJiVQmV/2YaWYG33LN9TmnxXG69rkp\n8H5lR2HXssD8B0yqmSuSnFJzJxPu3viPn+UcDbo3w74RcO4GRfbslq682OtpZm0qp8IL0JqBUsxd\nrgVMysqkxaQnfW/SEdKHB/kGQhh699nHQucYOv/pE5zQZCDT9QaZt/yRb5QdxptUxrqCjeQe7cCi\n0z1YfdhihenPCk9/nuSrtOUko10budZaT457E12tY2c1QnLudUiJjmAfwN7VAaWb/4oTWux3gDdw\nLkzMLWVFZCowFaBHjx5Rnk2AjA6OAig6dDqmFcDxJa+y8Vh7EvBw5appcPNvYi+yxsfEkYN4fN0O\nNpt0ljOAEWbz2VU1o6EEQpy/5u53edQ9jf2kkCVbecj9H9+JAv3HOw7V5JSazs76ut/Bwiwt21lB\n5r8CGFrbJ3ipchp/cd3KM95b+Z3nLtbYvfl94iu06X+t45Se9+jZQlR8JSWypgR8FrUlN50rimXI\n7c49OXXQF0FjfOWrTVCyVyUFdlee836Rd7yj8awvB2Cslc997vcYJtsdheVKP1sQBu+ARKDr4LPL\nOYT5/4Uv6RBUXdVb4ex+xkzDlZbN0KtgKPAD4HjBMpa+9zKLDrVgkT2I3aYTc7wjmOMdAVWQ2dbF\nkBYlDDn2MYOkgAGuPTQbdnvNVX2o8qkR4WQCimDNzLqZ4S6BSCiAPUBa0OPuvmMXew4AxpgXgBcA\nsrKywgX5RoXMxhAJVJJH3vx/YvghQ2QbyfZJxy4aS5E1QTRLH86dV57h6fwzTO/wI0Yc/14Ne3vU\nCu+FrOKnz1/MghP9ac1pnkr8O27x2XzdSU5JAQhfergh7vfgOxzB4RPqlvHyo4R/c5lVzCNV9zHH\nHsGa8l78detLZPFR+BW0sQPZ7X4FE9TTocZOprZyCcHXP/w7VK+Ije2rd+Rhg92TZz2TmGdfhcHC\nwuZLrbZyX8X/0s/aHZiPuKD480APBH+IaJfBASVzLrPJ+XZh6TnObij4Xuz41Blzyvs1zm3TewTj\nvmIxbvoEjKeSYtOFXDOIRfZglnIFhcegkG68xTcAcFvQY4OLzO2FZHY9RY/0XnSo2EXHY+vp0Gso\nLdKHkXj7OyTmPUPCtnlUGYsKEqjAzenK5hxZuZSjlZdxxG7BRLOC5vW8Y4+EAlgB9BGRDByhfjtw\nZ8g5s4EHfP6B4cDxxmL/99MoGsQX57LU2x9w6v/UiDuPIcEPVK/u7rp8FP9YXcaHe5ux+6636b7q\nj4GwumiZrIJWmhtMJo8X9gLgjwkv0F0OO7Hrw+4OCKFQYdmQCXfB5QlWvV7dMOZm1wr6SwkPVt3P\nOtOLr5b9hO8V5/GASSBJQnoLhNuphJZeSE5x/mfHSwKx++HMXp4yWP8mft+DMcKyqj48672FXF9p\nkkSq+IrrU77jnkNPTylY3sC4bXrA8V1waGvgmNcL+S87f/sryPpNSuFITnHmTC2RNv48jLUznEVG\n6RaclXiFs/oO1zN76NeRfevI2LuaDPMhd7s+pNK42ZzYh7VX/YE1p9qwbvdxdhw8ReExL4W0ggNl\nsMZvKuwJS48AC3yP7/D9hFAU+DMrcSeZ9bxjr7MCMMZ4ROQB4AOcMNCXjTEbReS7vuefA+bihIAW\n4ISB3lPXcRuajFiPBPJ9OZcZJzRthHtrzbjzWCJoi97JlciE1L/zzq4k/rmhjGljpkWnZlEwPqF6\nquBzvv9ZKpUk8A3Xh4xzrYRe1wUyVf2r32gksIWaR4pzz4r0ybD283+Jv+J/PF/hOe8XeebYCOZI\nOr91v8zohK3Q9wuBUgnhbPzDvwNLn3Hed96jOBEvHmdFfuXkmq+zXI6gBjixF9sIC0wW//Dcwmrb\naXzfnHLudH3CtztvofPRVT5HseWs+P2lrVt2dBRAbZxvYVCS5whx23ZyL8IJcwgsit5/qGaV1j2r\nnM+mfxcXbAqzXE7Ekc/mnygeBrONwa3yuXu8o/DPfPoXihe+RqHdiSLTlT2mA4dMG0pNGw7ThrKE\ndlTYQqXHpgoXiXhIpIpmVNFCymnHSdrLKdq1aUPixOcahw/AGDMXR8gHH3su6G8D3B+JsaKFPxeg\naN/hs+OBo41PoB7zJLDZfpZEsRk25UnIiNFI25D4+ikHnuAdfsXM1Yd4IKs5raJhsgrjNPzFkkSK\nyvfQ3yrhZwmzapYpCJMg1mBzDmfjDrWP+zJdE8XLTxL+xQ2u1Uyr+hbbTXe+XvVTru+WyI/6tWDg\nqaU1rz9Y4IkERQn5dw3GF/1iatq4h94F+S9zwiTztjeHf3rHUmCcJK127kqmMJvJrg+c5jZHfG8l\nlrOb6HNTQBEd2OQrd+2jyxXOMeMb83zx88HRSEbOX7hx8J2w+o0gc5CpGR4bfA9snMQ1CBSv8++O\nfDTvfTUDFz/OQG9RyEA+k6G/bHVt+QF+KtzQet655x4BYs4JHKt0PLqWlpRxrCqZo698jXb3/Ct2\nlIDvQ7/cexkGiyHtzpAUq8IfzhJWQ0wB2bKZPDOAVz9czgOXexpe+Pu/kJYFIx/gP4e689babiQn\nuHjmK4NIOvHjwJzC2ccbMoEtnI0755GaGar+DGQBjE2WtY05idN40TuRZ81X+GRnJZ/uLGesdZrJ\niY8y2tqMGE9NgWcs5374q2naPiHsj3P3lanwFOayXLKY7UngPc9wzpAEQNeWFt8e05/bux2k+Yx5\nPsenHwu6DYX9G3xJaYmOAvA7dDe/61T8XP58oIrnyAcgqXXA4Q5hk8guajcWbA7yJ6jVdg/8cf+D\n74AuQwL1hIKr2QJhm8+kDnN2Ff7/37gnamYRh2J7nTk1hh1APCDrZpEpGawzvSj0pHDl2pmxowB8\nH/qllZcBMLJXhyhP6DyEptPPf4wHE97h65UDeGlnF6Yc+CEt3X9q2KY1/tWYbbMj91/8d+XvAPh1\nTnN6Dx4FjAqcH+2aRbWN779XfgETXF3UW0mieLm/2Xy+9tXv8fdPt/NGcSs+sq/ko/Ir6Sn7ucnK\n5wbXKoZYhSSJJ5Ap7E8mI2Cr3+9txfL332Hp7io+9A7jCK2BawEY2eYId41I56ZrRpPgsoCMgJ8i\nOAu462DYu+ZsZ23WFOenOuzTt5pPau1c67kifC6mPlHwa9Kya4Zs1pbBvPJVx+keWpgu2EnurzPk\nxx9lBTUL4vW5qYZJqSZBSjbGncBxgiFD9rPO9GKH3Y0rw1YhjBK+D/2yVwrhDIxMTTj/a6JNsGO6\n80BGFeWS9dlB8k934lXPDdwvc8LbeeurqqVlgW1TYdx8v+r7nCGJW1xLuC2pPXDN2XOPZmRVbePX\nFv7YeWCghv3gO+iQls0vm1vc98pd/KtyNG94r2en6cKL3om86J2IhU2m7CO9laHNhu60KndRXu7m\ntGnGXpNCkenqCPydgSllyl7Gdy/nS7dNoXenVuHnHCxk/UrLH8F0Lke0X2Ae3+30YzhfnsWlBj2E\nfCZr3N/cJx3B7h+3tj7MoRFG4oIJfwmzc/TVDbJc0KEvHN5+9k7A9tR7EIQqgAtl8J30Xv4c2FBA\nd2cLGEMcKbPZcqYViVQy5KPbIfWt2NmhnI+0bOTAJh6seplv8BNe8kxgSuJCWoQmIkW6kmKwMhn/\nJMz5IY9X3ckmk04POcDvEqcjzX8XPhkq2pFV4cavLfwx3Llp2XS653W+X5zLfWcKWPX5fBZ4h/GZ\nPYgCk+r8nABOVADtgetqvLwVZ7jS2ka2tZkx1loGyC7kSDOoGI5THCCE2nr9nkuRBkc4+Vffluvs\nKKT6IPSehavSGq6kRbBJCXFCV/1d26ojq4IL4nng0LaAeelCfR0RQhXAhZKWTZ+sAlgK2+zUs+1+\nUSZvzVqgO1da20myz8RMxu8FUZIHcx/havEwTLaxyvTl5db38f3QSpplhyOXZRtGmXw09O9MX9KG\nBDw8k/gsrUZ9q3F1eroU+zfgfuVmsi0P2dYWfsoMyk0C20x39poUTpiWnOpwBUndLqe5fZJOnn30\n2v4ynTh6dvnr2las51Lc51OkwWYV4w04YkMriNY359r1hfoj/NdUSzXZ4KS4akVgTOC6/MX04jET\nOJbpk3gYaMN2OzWmyioALK3qDZQzwtoSHbt0XSjOBdtGBB51/4s7qv6b50ov43Z3MzpKmeN8nPuI\nowQiZXsPWS3vWf4Wj6517Pw/GXCYQdc/G/sNakK5FNOU794HkyRVDJIiBvmD0o8vglO+xikMhqK/\nBbsDfEj4PrnFuY7ppi73MdzqOxr/h1BlFa5RfbDgrs1Z7zeFhfpEonBdqgAugp4Ds0n8bB976Mhp\nqxUtYkjILtvjhOmNHNATronxlWooySk+G7xhZMJ2buhQwcf7mvEUX+Z3Cb4EIH+2anBPVj8X6xcI\nSWiqlETuX9WNY7aL61zr+OZ11zvvc2DTuROKYpGLNU35W3YGRUCx/PmgCCKfo9Nf4njiU46J4/On\nYOt88Pc9Dq5lA2d3AKuL2aY+fS518SmFlgsPrWV1rh1ZOJ9IFL6zqgAuAnfP4WSmzGfLYS8FN89g\ncIwI2cPblrH1iJdmVDK48EW45uZoT+nCqZG444LxTzItoQsLZx5ilvc67nZ96JQJsBKCHIeznC/V\npdROCRZMItBtCI8fn8CaM71IpZS/uJ/F+izPKYx2IQlFjZ207EBIor/1Yf8JgWiY6tpBIVEpt884\nt/AMFo6RMNvUh8+lrj6l0NyL0L7AwaG5tV13lH1JqgAukj7dO7Pl8F620YPB0Z6Mj+Vr1gGpZFnb\naOY9Hd1CahdLmMSd3vv/zdddNq95b+KnVffy78TfYA29M3wM/sWWMK5RUx/m7XLxSlV/n93/GdrJ\nCSc2u2iRzyF3gQlFjRW/AvZWOhnYfr+W/x4GFZs7y8Zfm/AK3mF5fclUXYbULNoWC31x62riCxPO\nHDY0N4a/h1a0J9DY6NOpJQAFB2OnKNzSSqdWzQhrM9WF1F69xfmSxTr+VVSNmvGGR9xv0pGjrDT9\nmGXfEOiNG3r+gElhXn8B4yFst1P5cdVUAKZlFDC0T1qgRr7tb2riqlkauTHc04uhtuJufgbf4WSw\niisQihnuHvib4ORPdz57K19z3tO/Mp7/WOB1/pX3x7+FV252XhMNwn72LpK0bGelnzXFUQbX/yz8\nTsJ/f2Ls86M7gIukb2dHAWzbuRty58ZEhc1lB51/48i0JNhvRbeQ2sVSi323zeo3+KX9Gg9UPcjj\nTOHaVleQGu58CFSI7DI4fIZo6Gpz8mwO5/+bb+YN5CTNmeBewT033wySU7MO0fDvwP51TkZqY4oG\nuhjOFznkv99+m//K6YH+zeHs/TWyaIO6ZIUmSwUl3jH3kehE1EXat3CuHVGsNGUKQRXARdKns5Pk\nsm3XXjj4u6j/Q0tPVrD94CmSE1wM+sI98Po70S2kdimEfnF8sdQTinJ5Z1MiC4oreWjWamZ+ewRu\nlxU+zE6E6r6srmaB/0n+9EDKvu94eZcrmbrfQ4k5yuCWx/lzZiGybpaz2vU7mU+VwrJ/OGYPvzmo\nMSnWcxFGIZ5XCG77oLrSKN6KmveghlktTOkEf5RLcLKUL/EOcO5rtO5pQ5hoYjiaTBXARdKzfXMS\nLZs9dgdO225aEN1/6PIixzadld6OxPThMVv7/6JJy0bSsvnDsApu/msuK4qP8vTbn/Fw51Xhw+yC\nE7P9Agoc4e9PzfdW4Fkzk4d2FLJyfyu6tbR40fsLkreVOs+vfgNu/qPjZK5O1qGmUGtMijUcta1G\nz+c3CQ4V9ZcZ9xO6iwgOh/S/PjRZavyTNRVzY76n5yPapUPOgSqAi8TtsshsZbPluEWB6c5gqyR6\n/9CSPJYu3ga0Y0SmryJhjDudLpaUls146mtD+PpLy3k6/wy9Ej5hUrM/OUKm2tHoL1Lmwy+gQgSX\n11g8vKwl872taMUZXuq9jk5bDgVeF+xUrn4/ccIkG7LJS31yKavR0FDR0DLj59tFhHv/4O5jjf2e\nno/6DGOtI6oALpaSPHqfWskWhrPddGOwOUft8nqeB6/ewrLTvwHaMaLlAaB3dOZSz4xqVsjPOi/j\ntwdG8KOqqaTwR672x1xbLidscftHTnneUAHlE1xlNONBeZQPvf1pyRlea/ZHBtKHGorDctfsi+vv\nRBXFpt0R51JWoxciwC5l4dHEFivnJNhcFvw4yqgCuFiKc+kru4DhbLe7N0jBptrmcdCTzA6TSnPK\nGXRmGTC6YefQEPgU3bc85ex1HeFl73jurfoRz/A3bnStDKz+h97JWbVXfIJrz/K3uH9VN9bYvWjN\naV5O/DNDXYVA30DUDwLD7mr6K9OLWY2eq7G6cnHEqCNYFcDFkp5DH+sjALaZ7tGz6aXnsMx8BkCW\nazsJmdc3/Bwaguo8AcPP3W9QTiIzvGOZWvVD7rPf4wfut0jaMifQXIQ3qp2O9hee4O2CKn61YQQn\n7QRSKeXVxD/Q27XfWfhv/9BJMAtOxYemL+wu5PpiVGA1WmLUEawK4BLo63L62W8zaY7TMErRC0t7\nPQibyhlx1YiY+DDVC+k51e0GLTH8LmE6qa7jPFl5K896J/GePZJ7XXO5yZVPN3MUgIOmDZ9UDeO1\nt86wyfQEYKyVzx8SXiRFTgasPrY3OoXFGgMxKrAaLTHqCFYFcLEU55Ju9pJEBXtMB46f2EabKE1l\n+UEX4Kv/H65kcZPBH09uIf3Hc//o73LV3kp+/u4mttmp/MozhV95ptAcp+OUvyMVQGeO8GjCm/yX\ntejs6pVwdoaq4hCjAqvREqOOYFUAF0t6Di53Av1kN2tNLzY3v4oRUZjGgRPlFB46TYsE4fL5XwW7\nrGlu1Ytzg/rR2rBtPox+kOzhVzO3SwIfLl/Lf/Z1IO+gxUmvI/hbWJVc2fYMt5z+DxPlc5JcNtiu\nQBy7n3Dt/BSHGBVYjZoYNC3WSQGISHvgX0A6UAx81RjfPrzmecXASZwish5jTFZdxo0qvi9G/9nb\nWVsCm72pUVEAS3c48f9XtSwloew0Z2VbNhVqSxoC3LsWM35kDuMBM/0WTnksDDatpAKpbAZffALK\nRjtlHAo+hC3zqNmEu4nes0gRgwJLiSx13QE8BnxsjHlCRB7zPf5JLedeZ4w5VMtzjYu0bAYM6Qgl\nm9iy72TDju2LzFhW5HQkG3FyAbh9Rm3L3fS26uGShpJTajooh9yB2JW0kqAVvqciUOGyuoyD2+nD\nCk7YaGiGqqLEGXVVAJOAMb6/XwUWUrsCaFIM6NoagM37TzTcoEGRGUvL/wx0ZqS1IfB8n7FNc8UW\nGppZow57OZw64Gu156svgzi/d3zqFMbzl4iwgdRhTvGuWKhGqShRpq4KoLMxZp/v7/1A51rOM8AC\nEfECzxtjXqjjuFGnv08BbN13HM9nT+LObABB4hN8e+y27DSdacVpLpPiwPPbFwTi35saoeYIX2QQ\nGOe6b/6jE/+fnALLn4XSrVS32oOze6yqeUNRzl8OWkQWiMiGMD+Tgs8zxr/MCsvVxpghwM3A/SJy\nzTnGmyoi+SKSX1paejHX0qC0SU4gtaVFhReKP3mlYcov+yIzlngvB2C4tRm3BNm0/UlpTZ20bCdD\n1x8dZHsc4e9vNl66NeQFApljmp6DXFHqyHkVgDFmrDHm8jA/7wIHRKQrgO/3wVreY4/v90HgbaDW\nb6Ex5gVjTJYxJqtjx46Xck0NxoDEAwBssruHr6UeaXwO6KUtnKSvUdYm57i46lbTvDESXKc++Lo3\nv3v2ue6kxtMgR1EakLo2hJkNTPb9PRk469snIi1EpJX/b+AmYEPoeY2OkjwGnFwKwGa7h3MsOaXe\nhzXdr2KJpx8Ao1ybwZ0ME/5SeyOKpoa/sQaEb8AxYFLN8/tPjI/7oiiXQF19AE8Ab4rIvcBO4KsA\nItINeMkYMx7HL/C2OFk4bmCGMWZ+HceNPsW5DJQiADaa9AaLKS9cv5T9p21SEqro16sXDHw4fhKZ\nwpUn8Jt9/PjvRXCPW0VRwlInBWCMOQzcEOb4XmC87+9CiJn2uZEjPYcr3C9DJay3MzHGIPUdU16S\nx5L/PA18g5HefKRoIexaGj+JTBdaniBrigp+RbkAtCfwpZKWTerNj9KeExylFbtNh/qPwy/OZWmV\nz/xjbazZoSoeiEQPV0VRqtFSEHVAyg9zhVXBZ/Zg1ptepA2t35W43fNqltqO2WmUtfHs0MamjpYn\nUJSIogqgLqTnMMj1Mp/Zg1lHH8b7ywnXE1sOezlKK7olnKbngKugVaem1azkQtD4fUWJGGoCqgtp\n2Qy68esArG831lmZ1lcuQEkeS95+DoBRdj6ybS6smVk/YymKEheoAqgjg4ZcBcC6gx7Mx7+rv4Sw\n4lyWePoCcWr/VxQl4qgCqCOdWyfRKbGCkzRnp92h3oRyZbMUltv9ARhpbYo/+7+iKBFHFUBdKclj\nkNfJyF1rMp1jkU4IK8lj5dz/5TTJ9JXddL36brj+55rgpChKnVAFUFeKcxkkOwBYbfcJJIRF0gy0\ndiYLqwYCcK21BipOOAlQKvwVRakDqgDqSnoOWe5CAPLtftRoMhIxDJ/ZTi7dGGsttdfcUxRFuXBU\nAdSVtGyG3P1HXNhsMj05RYuI2+b39/oaW0xPkiknK6EI6jncVFGU+EAVQARonjmcy7u3w8ZizeD/\njrhtftGp7gCM6uSh2T3vqulHUZSIoAogQlzZsz0AK1reEBkB7a96WZLHwm1Ole0xI4er8FcUJWJo\nJnCEuCq9HS9/XsTKnUfr/mZBVS89VhK5VS8BcG3fTnV/b0VRFB+6A4gQV6a3A2DVrqN4di6vXr1f\nEkFVL1dX9eBkpSGzQwt6pDSP4IwVRYl3dAcQITq1SqJnSnN2Hj7DppfvZ5C1A1zNLs0f4K966a1k\ngckC4Np+sd0dTVGUxofuACLI8E5Of97Pvf19pRoqLi0c1Ff10gy7mw9c1wLwBe9n9d9zWFGUuEIV\nQATJSXISwhbZg5wDYtUpHHTbqlyKy5uTwnGuWvPzhmk8ryhK3KAKIIJcPWwQgk2+3Y8zJMPIBy49\naqc4l/lVTvLXWNcqXHi1+JuiKBFFFUAEaZdkMaj5Eapws8weAMufv/QVe3oOH9iO/f8L1gpA6r/j\nmKIocYUqgEjhC928tmIRAIu8l9dpxb7ruJdNdk9acsYp/4xBS0AoihJJ6qQAROQ2EdkoIraIZJ3j\nvHEislVECkTksbqMGbP4Qjevca0DfH6ACy0JEZT05eedvK0AjLVWkSRVzkHboyYgRVEiRl3DQDcA\nXwaer+0EEXEBfwduBHYDK0RktjFmUx3Hji18oZtDTBGtOEOh6cbOL71Nz/P5AIKSvnAlOtE/3a/i\n7QNdAC+3uhYHzq2jU1lRFCWYOu0AjDGbjTFbz3NaNlBgjCk0xlQCs4BJdRk3JvGFbrpvmMZ1fZyk\nsHmHLyBzNyjpy28yWl1yjKLjXjomeRnt3kK1/X/8k1oKQlGUiNEQiWCpQEnQ493A8AYYt+HxNSwf\n324/s7evZN76fXz32l7nfk1Q0pffZPT2yj0ATMrqjXvwHEdJpOeo8FcUJaKcVwGIyAKgS5infmaM\neTfSExKRqcBUgB49ekT67RuEMf060jzRxdrdxynZuJS0I0tqF+Bp2TDuCdj8LgyYRLnH8N7KQsDF\nrcNSoVsbFfyKotQL51UAxpixdRxjD5AW9Li771ht470AvACQlZXVKMNekhJcXN+/E++v28c7//pf\nvu96y1ndj3sCyg7XVAYleU4HMW8lFC/m/aqrOVb1LS6zdjLQ04GmullSFCX6NEQY6Aqgj4hkiEgi\ncDswuwHGjSq3ZTk6719VV2PbNngqYO4j8PFv4ZWbIX+6c2KQD8B4qni16noAJrs+QHYuruXdFUVR\n6k5dw0BvFZHdwEhgjoh84DveTUTmAhhjPMADwAfAZuBNY8zGuk079sk5PodU93F2m44sMVeAZYHt\nBWwnnHPOw/D+D50G8q5EEBcrGMB6k0k7TnJLwgo4vltLPyiKUm+IMbFrZcnKyjL5+fnRnsbFkz8d\n3n+Qpz238hfPbVzXeg+v3JTo7ABsT81zLbdTMqLiOHcu6cYS+zJ+4HqLh5u96xSU84WGqh9AUZQL\nQRORXdsAAAZjSURBVERWGmNqzcsKRjOB64PNjm/8664FJFPOpydS2WD1ccI4LTcggXNtDyx9hmX7\nhSX2ZbTiNPe65zrHjdb/URSl/lAFUB8McNIcUuQkX3d9DMD//N+nznPjn4QOfWucXuUVfrXLKfx2\nr3sebeQMiIC4It5gXlEUxY8qgPogawr0nwjAVPccWnKGj+2hzJ89A+b9GA5tq3H6i97xbKnsSJoc\n5Duu952D4oIrJ6v5R1GUekMVQH0x+kFwJdJJjvET9ywAfl45hZ1VbQku6rbIewV/9twGwG/c00mW\nSucJY0Ob7ir8FUWpN1QB1Bdp2TBlDvSfwNfdC8mx1nOINny96qestvvgMS5mecbwrapHsLH4vutt\nrnOt8b1Y1PSjKEq9oz2B65O0bLh9BlZJHv8o+Jy71rhZc6Ajt1b+Ggsb26d/73R/yg8T3wPjAssF\nfW6ClhdQR0hRFKUOqAJoIFq6bWZOastTC7byZmEiR2lFL9nD/d228+UvfRX4qhPtk5wSyAxeM1N9\nAIqi1BuqAOqboHLPya5Epo17gscOPESV10ui2wVfChLwadlOX4CQ6qCqABRFqQ9UAdQ3oeWeyw4j\nU2aTWFuFzzDVQRVFUeoDVQD1TTiB7isbHRZfXwEtAa0oSn2jCqC+uRSBfi4FoSiKEiFUATQEfmG+\ndgasnQmD71ABryhK1FEF0BCU5MH0CY4ZCGD1GzDlfUcJlOSpuUdRlKigCqAhKM4Fb1XgcXCBt5CG\n8KoEFEVpKDQTuCFIzwFXQuCx3xkcpiG8oihKQ6E7gIbAXxZi7QxAavoANORTUZQooQqgoQgX2aMh\nn4qiRBFVANFGQz4VRYkS6gNQFEWJU1QBKIqixCl1UgAicpuIbBQRW0RqbUIsIsUisl5E1ohII+zy\nriiK0vSoqw9gA/Bl4PkLOPc6Y8yhOo6nKIqiRIg6KQBjzGYAEYnMbBRFUZQGo6F8AAZYICIrRWRq\nA42pKIqinIPz7gBEZAHQJcxTPzPGvHuB41xtjNkjIp2Aj0RkizFmUS3jTQX8SuKUiGy9wDFC6QDE\nm8lJr7npE2/XC3rNF0vPCz1RjDGXOEbQm4gsBH5kjDmvg1dEfgWcMsb8uc4Dn3ucfGNMrY7ppohe\nc9Mn3q4X9Jrrk3o3AYlICxFp5f8buAnHeawoiqJEkbqGgd4qIruBkcAcEfnAd7ybiMz1ndYZWCwi\na4E8YI4xZn5dxlUURVHqTl2jgN4G3g5zfC8w3vd3ITC4LuNcIi9EYcxoo9fc9Im36wW95nojIj4A\nRVEUpfGhpSAURVHilCanAERknIhsFZECEXks2vNpCETkZRE5KCJx4VwXkTQR+VRENvlKkTwY7TnV\nNyKSJCJ5IrLWd82/jvacGgoRcYnIahF5P9pzaQgasnROkzIBiYgL2AbcCOwGVgB3GGM2RXVi9YyI\nXAOcAl4zxlwe7fnUNyLSFehqjFnlizBbCXypKf+fxUm3b2GMOSUiCcBi4EFjzLIoT63eEZGHgSyg\ntTFmYrTnU9+ISDGQ1RClc5raDiAbKDDGFBpjKoFZwKQoz6ne8SXVHYn2PBoKY8w+Y8wq398ngc1A\nanRnVb8Yh1O+hwm+n6azeqsFEekOTABeivZcmiJNTQGkAiVBj3fTxAVDvCMi6cBQYHl0Z1L/+Ewh\na4CDwEfGmCZ/zcBTwI8BO9oTaUAarHROU1MAShwhIi2B/wAPGWNORHs+9Y0xxmuMGQJ0B7JFpEmb\n+0RkInDw/7d3hyoRBHEcx7//ItgsBuEEDeJDnMmkBrNBk1FfwJfwDWyKImiwieBFgwgiiD6AV0w+\ngPAzzASjwdmBmd8Hjr3b9INl+c3Ozu1KeqqdZWBr+ThvAgd5ireI1gpgCiz++j3K+6wxeR78CjiT\ndF07z5AkfQETYKN2lsLGwHaeE78A1iPitG6k8iRN8/aT9D+rYu+Mba0AHoGViFiOiBlgB7ipnMn+\nWb4hegK8STqunWcIETEfEXP5+yxpocN73VRlSTqSNJK0RDqX7yXtVo5V1NCPzmmqACR9A4fALenG\n4KWk17qpyouIc+ABWI2Ij4jYr52psDGwRxoRPufPVu1QhS0Ak4h4IQ107iR1sSyyM4M+OqepZaBm\nZvZ3TV0BmJnZ37kAzMw65QIwM+uUC8DMrFMuADOzTrkAzMw65QIwM+uUC8DMrFM/Yvq4uviPnB4A\nAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x7f7daf6a2090>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig = figure()\n", | |
| "ax = fig.add_subplot(111)\n", | |
| "\n", | |
| "ax.plot(T, Z_nominal, linewidth=2, zorder=3)\n", | |
| "ax.plot(T, Z_data, '.')\n", | |
| "\n", | |
| "show()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In PyMC3, model initialization is quite different, and mostly for the better. Basically initialize a `Model()` context and populate it." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "model = pymc3.Model()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "with model:\n", | |
| " C = pymc3.Uniform(\"C\", lower=CMIN, upper=CMAX)\n", | |
| " \n", | |
| " Z = pymc3.Normal(\"Z\", mu=response_model_op(C), sd=SIGMA, observed=Z_data) # NOTE mu here" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "You can use `find_MAP()` to compute the maximum a priori estimate of the parameters, i.e. the optimum starting point. This is optional, you can also specify the starting point by setting `value` like in PyMC2.\n", | |
| "\n", | |
| "**NOTE** that Theano does a lot of magic automatic differentiation stuff. However, because we're using a non-analytic node (as generated by `as_op`) it can't do automatic differentiation. As a consequence, if you use `find_map` you must specify a gradient-free method. Powell's method is a decent choice, I think it does this by default." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.\n" | |
| ] | |
| }, | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Optimization terminated successfully.\n", | |
| " Current function value: -434.357388\n", | |
| " Iterations: 2\n", | |
| " Function evaluations: 47\n", | |
| "{'C_interval_': array(0.1994328510713635)}\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "MAP_estimate = pymc3.find_MAP(model=model, fmin=opt.fmin_powell)\n", | |
| "print(MAP_estimate)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The MAP estimate is actually worthless in this case, so I'm not going to use it. You can pass `start=MAP_estimate` to the `.sample` method. If you omit this it just randomly samples from the prior (I think). \n", | |
| "\n", | |
| "The `njobs` parameter sets the number of parallel chains." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "100%|██████████| 100000/100000 [00:13<00:00, 7596.55it/s]\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "with model:\n", | |
| " step_method = pymc3.Metropolis(vars=[C])\n", | |
| " trace = pymc3.sample(int(1e5), step=step_method)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "\n", | |
| "C:\n", | |
| "\n", | |
| " Mean SD MC Error 95% HPD interval\n", | |
| " -------------------------------------------------------------------\n", | |
| " \n", | |
| " 1.511 0.018 0.000 [1.476, 1.548]\n", | |
| "\n", | |
| " Posterior quantiles:\n", | |
| " 2.5 25 50 75 97.5\n", | |
| " |--------------|==============|==============|--------------|\n", | |
| " \n", | |
| " 1.476 1.498 1.511 1.523 1.548\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "pymc3.summary(trace)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7f7d910c4090>,\n", | |
| " <matplotlib.axes._subplots.AxesSubplot object at 0x7f7d8c893250>]], dtype=object)" | |
| ] | |
| }, | |
| "execution_count": 12, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAACICAYAAAALQhGXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8ZHW5+PHPk97bJtnd7G52s53tfWHpRSmCC0oRVJAi\nF7tX/SnqverVay8XFRVpAoJ0EFSKVEFYdmF7Y3tLtqRsek/m+f0xM9lJMjWZyUl53q/XvDI5M3Pm\nmTMzZ85zvt/v8xVVxRhjjDHGGGNM/8Q5HYAxxhhjjDHGDAeWXBljjDHGGGNMFFhyZYwxxhhjjDFR\nYMmVMcYYY4wxxkSBJVfGGGOMMcYYEwWWXBljjDHGGGNMFFhyZYwxxhhjjDFRYMmVMYOAiFwjIu+J\nSIOIHBGR50XkNKfjMsYYY/yx3y1j/LPkyhiHichXgNuAHwGjgWLgd8CHnYzLGGOM8cd+t4wJTFTV\n6RiMGbFEJBsoA65X1cedjscYY4wJxn63jAnOWq6McdYpQArwtNOBGGOMMWGw3y1jgrDkyhhnjQIq\nVbXD6UCMMcaYMNjvljFBWHJljLOqgHwRSXA6EGOMMSYM9rtlTBCWXBnjrFVAK3Cp04EYY4wxYbDf\nLWOCsOTKGAepai3wHeB3InKpiKSJSKKIXCgiP3M6PmOMMcaX/W4ZE5xVCzRmEBCRjwP/CZwE1ANr\ngR+q6tuOBmaMMcb4Yb9bxvhnyZUxxhhjjDHGRIF1CzTGGGOMMcaYKLDkyhhjzLAnIveKSLmIbAly\nn7NEZIOIbBWRf/ksv0BEdojIbhG5dWAiNsYYMxRZt0BjjDHDnoicATQAD6jqHD+35wBvAxeo6kER\nKVTVchGJB3YCHwBKgXeBq1V12wCGb4wxZoiwlitjjDHDnqq+ARwPcpdrgKdU9aDn/uWe5cuA3aq6\nV1XbgEeAlTEN1hhjzJA1JCaAy8/P10mTJjkdhjHGmChYu3ZtpaoWOB1HD9OBRBF5HcgEfq2qDwDj\ngEM+9ysFlvtbgYjcDNwMkJ6evnjmzJkxDdgYY8zAiOR3a0gkV5MmTeK9995zOgxjjDFRICIHnI7B\njwRgMXAukAqsEpF3IlmBqt4J3AmwZMkStd8tY4wZHiL53RoSyZUxxhgTY6VAlao2Ao0i8gYw37N8\ngs/9xgNlDsRnjDFmCLAxV8YMMi3tnbS0dzodhjEjzTPAaSKSICJpuLv+bcddwGKaiJSISBLwMeBZ\nB+McdJraOujodDkdhjHGDArWcmXMINDU1sGTa0v528YjrD9UTYdLKcpO5eplE7jp9MmkJMY7HaIx\nQ5qIPAycBeSLSCnwXSARQFXvUNXtIvICsAlwAXer6hbPYz8PvAjEA/eq6lYHXsKg9dK2Y+SlJ3H6\ntME2jM4YYwaeJVfGOKisppk/rzrAw2sOUtvczswxmdxwaglpSQmsO1jNL/65k0ffO8Q91y1l+uhM\np8M1ZshS1avDuM/PgZ/7Wf4c8Fws4houjje2OR2CGUIOVjURFwfjc9Mcef49FQ3UNrezqDjXkec3\nw5slV8YMIJdL2VBaw6vby3n1/XK2HakjTuCCOWO48bQSFhXnIiJd9397TyVffmQDl//hbe6+binL\nSvIcjN4YY0a2prYOWttd5KYnOR3KkLb+UDXgXHK1pawWwJIrExMxS65EZALwADAaUOBOVf21iOQB\njwKTgP3AlapaHas4jHGay6WsP1TDPzYd4bnNRzha10J8nLC4OJdvXDCTi+eNZUKe/x+YFVPyefIz\nK7juT2u47t41PHDjMpZOsgTLGGOc8NK2YwCsXDDO4UiMiYyqdjt5O5y1tHey/mANSyblkhg/8OUl\nYtly1QF8VVXXiUgmsFZEXgI+Bbyiqj8RkVuBW4FvxDAOYwaMqlLR0Mqe8kZ2VzSwubSGf+2s4Fhd\nK0nxcZw5o4BvzpvJWdMLyU5LDGudE/LSePTmU7jqzlVc/6d3efCm5SyYkBPjV2KMMcaY4aCxtYOX\ntx9jUXFuwJO5w8mOo/WU17dQVt3MpPz0AX/+mCVXqnoEOOK5Xi8i23FPxrgS96BigPuB17Hkygxh\nLpfyr10V/G3DYf69u5Ly+tau27JSEjh9WgHnnlTIebNGk5USXkLVU0FmMn+56WSu/OMqrr1nNX/5\n9MnMGZcdrZdgjDHGjDjPbjzM2TMKyAzx29zW4SIhToiLG5otP/UtHQAcrmkeEcmV0wZkzJWITAIW\nAquB0Z7EC+Ao7m6D/h7TNdN9cXFx7IM0JkIul/KPzUf49Su72F3eQG5aIqdNK2DhhBymFmYwtTCD\nsdkpUWuGH5Odwl8+vZwr71jFJ+9ZzYM3LWd2kSVYxhhjTF+oKqXVzZw0tnty1d7poqm1s6uHyfNb\njlCUkzosuuU3t3WyZv9xTp6cR3KCVSKOhZgnVyKSATwJfFlV63wPNFVVRUT9Pa7nTPexjtOYSPx7\nVyU/feF9NpfVMmN0JrddtYCL5o4lKSG2fXvH56bxl0+fzDV3vcPH7nyHP31qKUuGwc7eGGMGu7UH\nbHh4LFQ3trG3spHS6iZWTMmnIDM5ps+38VBNyPus2XecyoZWLplX1NVadbimOaZxxcreigY2ewp4\nAOwub6CmqY3S6mamFGQ4GFl0tXZ0UtvUTmFWitOhxHYSYRFJxJ1YPaSqT3kWHxORsZ7bxwLlsYzB\nmGjadriOT96zmk/cs5rjjW388or5PPel07l04biYJ1Zek/LTefwzK8jPSOaau1fzl9UHUbXzD8YM\nZ1UNrTyzoYzWjv5PMF7f0s6m0hpUlbqWdupb2qMQYXedLuWdvVUxWbdTSqubnA4h5upa2ikb4CTi\njV0VXdt2f1VjTJ+rpb2z13P4+/kcTlMLbDtSF9X1lde10OmK7jHH8cY2apr6t83f2XucVXuroh5b\nX8TsaFDcTVT3ANtV9Vc+Nz0LXOe5fh3wTKxiMCZaDtc089XHNvKh377JptJa/utDJ/HKV8/ko4vH\nE+9AH+xxOak8+ZkVLC/J41tPb+arj22kqa1jwOMwxgkikioiM5yOYyD9e3clAJUN/T/oW7PvOPsq\nG2lo7eC1993TQkRbVWMrx+paup0x74udx+pZe6CayobW0HceJFR1yCaVr71fznv7j8dk3XUt7dQ2\nD77tsqu8nn/trHA6jAHTn5EKtU3trNpbxdbD/r/XbR0uNpfW4oogwWnrcPHmrop+vwcNnnFlrkFw\nsjmWp9pPBT4JnCMiGzyXi4CfAB8QkV3AeZ7/jRmUVJW/rD7Iub/8F3/bdJibT5/MG//vbG46fTIp\nic72Vc5LT+K+65fxn+dN5+kNZVz6u7fYXV7vaEzGxJqIXAJsAF7w/L9ARJ51NqrBrb6lndfeL49K\nqxdAR6eLxtYwT+b08zhn+5E6SqubeGt3JS3t0Yk/UgeqGqnrkSwFOzu+q7yBV98v7/WYke6198t5\nfUf0EnmXS4Mmay6X0t7pCmtdPVtNBvLwfOOhGp7ZUBaz9Quhs6mqhlYOVoVumW3tdH8HvYlMTzuO\n1rO3soGDx8Nv5d0QpJtmp0upbRp636OYJVeq+m9VFVWdp6oLPJfnVLVKVc9V1Wmqep6qxuYUiTH9\n1NLeyRceXs+3nt7M4om5vPKVM/nmRSeFXUJ9IMTHCV86bxoP3LCMyoY2LvntWzyxttTpsIyJpe8B\ny4AaAFXdAJQ4GdBgt6u8gbqWdsrr+tf609TWwaHjTby9p4qXtx+LUnTh6+iR0KzZd5wXthwJcO/+\n83a33nCohtd6tO5tKatFVVl7oJqGHolmlad1sbmts+u+Pe8TrtqmdjYcqulX1+/a5vagien+ykY2\nl3Zviahtbu+aaDdc7Z0uqmPUnW7H0fper2HbkTpe31HeLdF/9f1j7KtspKPTxTv7qnhuc/ifjzd2\nVrB6b1XQ+9S1tPPCliNBt2dbR/eE7o2dFfxz69GA9/d2U4xm75OeMfjj+5H69+7Kromd+8PbahTJ\npzVYrBsO1fD6zvI+n1hpaO1gc2ntgA+dGPiZtYwZAmqb2rn2njX8Y/MRvn7BDB64YdmgLl96+rQC\nnv/S6cyfkM3XHt/ID/6+bVD0OzYmBtpVtedR34j5sPseJKhqRAcdqu4D/r52C3pzVyXrDlZTHcHY\niIog3fn2VjQE7F4UjiO1zbSGcRDZFw2tHTy78XDA8UeNrR0cb2yjtLqJ9QcDH5TWtbSzp6KBd/f1\n7Tzyqr1VHKhq7NfrfH1HOa9s791ipKq4XMrG0hr2Vjaw7fCJsTlv765kT0UDxxvbOBRmK8Safcd5\nY1dFyC5hb3u6uHpjCKW2qZ33j9b1KijiHRfV2NbBGzsrPOMHO9hUWsM/Nh+hor73Z6+0OvB4suqm\nNo7WtQSNZV+F+704Wuv/frvL63l+y5GuxNq73uYwvqcvbTtGY2sHb++pDNjKHM73vamtg+e3HGF3\neUOv2/pbvdi3FexYXQvPbCgLmci1dnSGlez5421R7HlixUtV2VxaG3C/tqeigb2VDTS2DWyrtyVX\nxvRQ29zOJ+5ZzYZDNfz26oV89qypQ2Jui9FZKTx443I+tWIS9/x7H599aG3UugEZM4hsFZFrgHgR\nmSYivwXedjqoWAp0sLqlrI4Xtx7ltTC7WlU0tPDPbUe75rwJx+7yBp7ZUMb6g9X96pZX2dDKy9uO\ndTvps7ms1u8BYKdLeXt3JVU9ErNYtYr44+1udiRActXh0q4D1eONbd26p4XbFS0aGlo7wmrx6HD1\njmnNvuP8bdPhrv93+elW/uauCtYFSR47Ol1d76k36Q6VLvkm3HsrQxewUM8aAx1gr9pTRXVTW6/W\nRX/6k8z7OlzbzItbj9Lp0m6Jw1ZPgtrXVqgtZbVU1LfywpbeLV0Hq5p4cevRkIUfmjyJxLEQiWIg\n4SS8Cuw65v7u+htfuOPoiUT9hS1HeTFIy53249xYZUMbeysbep3g6Pl5t5YrYxzU1NbBp/60hveP\n1vHHTy7m4nlFTocUkYT4OL734dl85+JZvLj1GJ95cJ1j4xSMiZEvALOBVuBhoA74sqMRxYiq0tzW\n2e0A2NeRWveBf12IIgHeU0M1EY5dUNWug9FIxlD489buShrbOsIq9LCrvJ6KhtauIh5e6w66u+CF\ns0/beaye5zcfYV+Ig/f2TlefDrxqemxz37FEvi174ZyxP1rbQnl9ZAfC7Z0untlQxivbj/HStsBd\nNHu+tqqG1q6xYKFaacLxj81Hgj6/N4adx/yPB46k62FNU1vA9Qy0ivpWWto7WbWniuc9XVN9W6s2\nHKrhzV3dCzQ8u/Ewz2wo45kNZWF9hnuOa6xsdCelPU+OuFzqt6dKOB9rVeV9n0TocE0zz248HHC8\nYLgNXz1bWr3dBTeV1vRq7epZmbGqoZXtEVY4rOqxjmAtlANhQCYRNmYo6HQpX3x4AxsP1fD7jy/m\n7JmFTofUZzecVkJyYhzffnoLN/95LXd+crHjBTiMiQZVbQK+7bkMa1sP17GnonvLjveAaXd5vd+u\nRmv2HaeqoZUL547tWhZOYtTW4eo1ncTG0r6f5fftPhQqwWjvdBEv0tVDoKMz8FHhvopG9lY2cNrU\n/K5l+ysb2VPRwLknjebQ8Sa2Hq7rarXfVFpDTmoiuelJvdbV0eniuc1HmFqYEfGE7KoaVgXDtzwJ\nor+z842tHaQnJ7B6n3ucz8oF42jvdLG7vIEZozOJi5Nulc8OVjVxrL6FpZPyulonfFU2tNLU2knx\nqBNd2Nf3KBbgTVgXTsgNGndbBK1vPXtIHK5pZkJeGq9sPxbxWLNQXQq3H6lj+ujMiNYZqdd3lAdM\nuHsW0KhqPPEZ8H2vGlo7oJVurcq+66xrdndhfHtP9xMIvs/68vZjrFwwrlcMPUN7bUc5Da0dFOel\nkZuWREZK90P7YC1Dx+pa2XH0RML6rqdSZG1TO1kpgceX+8YQ6C17duNh8jNOfO98i3Z4E7pmP59j\n72f0pLFZARNEl0u77f9cqsSHUbhjoFhyZYzHj5/bzsvbj/H9lbO5YM4Yp8Ppt48vn0hiXBzfeGoT\nN97/Lndfu5TUJEuwzNAmIq/hp+eRqp7jQDgx5W/MyLqD1RRkJrPtSPcz+O2e7lne1ix/gh3otnR0\ndiVXdS3tZCYncKCPcw6parcDuud9ik60tLt6nbV/bvMRctOSOHVqfsjWFO/Ze9+5ezaWuhOIQBXX\n/I1X2lvR0JWglFY3M7som+a2Tupa2tl6uJZphaEP4P2dXfc9gPY92G7rOPG6UxLjOXS8iXUHq1le\nMqrb43ccrWdPRQPpSQkUj0rr6mK481h9yFY4byLX1uni0PEmFk/KDTheqq/FC7xxnzOzkEyfg2/f\nrmrrDlZTXt/SpyIeR+taqGlyd7OcOCo96H37OXwooGAVCIONN/R3EiFQq7IC+yp6d4kN1JWvw/Oe\ngjsxKR6Vhqqyv6qpazsfPN7EweNNXdut3eWio9PlOdHh/iweqT0xcXBpdRMzxvj/nK87WE2nS5mU\nn05ts3t/oJwY89Xc3tnV9XHnsfquiZ9dPcaE+tuHeYX6PPuz9sBxJo5K53BNM/sqG1kyKa/XfQZB\nJfbwkisRmauqm2MdjDFOefCdA9z97318asUkrj1lktPhRM2VSyeQEC987fGN3Hj/u9xznSVYZsj7\nms/1FOCjwLCc5C1Q1xx/4xdee7/cb0tWsMlQfY9B1HNc6C3bnpvWu6UnlPK6FjJSEoJ2E/O20vRU\n3dTG5rJaDlQ1kuenlamnSCZ5Xb2vivNnj6GuuZ289CT2VzX1GnvT2tHJP7ed2K7eFsNIJtTdcbSe\nhtYT79kanyIWrR2urvdt5YJxXQfwPbeHt3tXz7l6Qh2I+ias3tfWcwzS2z26WUbCW5TAO/5qb0Uj\n8yfkdN3+r50V3eZ87E+3LO98R/UtHRRmJpMY378RLP/eVcnMsZkRzb3kK5yurO2dLt7YFdk8TeF0\ny2xo7SAlIY4tPsVGmts7eWHLEdo61W8Lm/ekSF1zO//wUy3R+1kJNefYxtIaXKpsLqulJD+duuaO\nrpY63zFllQ3uLpIpifG0B2l19hWoi3Kwk0OdLqW0upkjtS2kenrirDsQ3kmCgc63wm25+r2IJAP3\nAQ/5qdRkzJD1xs4KvvvsVs6ZWch/XzzL6XCi7iOLxiMCX3nMEiwz9Knq2h6L3hKRNY4EM4j0TKze\n23+cJZPyWLUneGlpr07PQVpLuzvLiqQioNeqEGWsQ/EesAXrFhhqfFkgB6qaeP9oHZkpCX4LevRM\nRPoy2a3v2JW+2FfZ2GssV7jCmT8qWOXGUF7dXt6tq+D+qu7JVSzsqWhgT0UDhZkp/VpPVWMrb+3u\n22uvbGwN6zsUScl3CP/79cr2Y6QnJdDYo0BGf6pH+nZnLQuRBHs/j+Ek93sqGoImR74CvX7fExLP\nbCgjPal3mtLp0q7WusEwYbA/YZ0OUNXTgY8DE4C1IvIXEflATCMzZgAcrGriCw+vZ1phBr+5emG3\nM2/DyWULx/OrK+ezam8VNz3wrt9+zsYMBSKS53PJF5HzgZADZkTkXhEpF5EtAW4/S0RqfSa9/47P\nbftFZLNn+XtRfDkxU1bTTF1Lu98qcV6VPl12Vu2pchd3iOAc7zMbyobMlA/exMdfYtXS3hmzku7+\nuFzqt+vSptKaru51wbrUldU09+paGev4wxmDFavPgr9iHy6XRtR62VfhnpyIlO84p1B6JlbRFEnL\nbCj+Kn/2l/e1q4aeEHogq3SGEvaYK1XdJSL/BbwH/AZYKO46pN9S1adiFaAxsdLc1sktD65FVfnj\nJxeTkTy8hyBetnA8qvDVxzdy0wM2BssMWWtx9/IQ3N0B9wE3hvG4+4DbgQeC3OdNVb04wG1nq2rf\n+1U5IFRp6s0+ldo6XO7iDnERDmSJtCtUMMHGZwwngao/+tpT0cCccf7PGby3v29zZkVbX0t999cz\nG8qYPz62rWbDWSSNPeHOceavwEo0vRpGmX0vJyY47yncMVfzgOuBDwEvAZeo6joRKQJWAZZcmSFF\nVfn2Xzez/Wgd9163NOTA2eHiI4vGA+4E69MPvMfd1y2xKoJmSFHVkj4+7g0RmRTdaIafSLvZ9LWb\nXtB1hjHGZSQI98DWKe/0sxtof2yJ0nxVI1EklSDD9e4gSfjBfwuWDnCjVrin6n8L3I27laqrDVFV\nD3tas4wZUh585wBPrSvjy+dNG9Il1/viI4vcLVhfe2IjX35kA7/7+KJh2x3SDB8i8pFgt0epB8UK\nEdkElAFfU9Wt3tUDL4tIJ/BHVb0zQIw3AzcDFBcXRyEcM5IFm7x3pBsq3VHN4NDa2QkELi0fbeEm\nVx8CmlW1E0BE4oAUVW1S1T/HLDpjYmDDoRq+//dtnD2jgC+eM83pcBzx0cXjqWlu5wd/38b3/7aV\n/1k5x+mQjAnlkiC3Kf3vQbEOKFbVBhG5CPgr4N1BnKaqZSJSCLwkIu+r6hu9gnAnXXcCLFmyxI7+\njDFmEKisb+t3YZRIhJtcvQycB3hHq6UB/wRWxCIoY2Kltqmdzz20jsLMFG67amHXpJUj0Y2nlXC4\nppl7/r2PWUVZXLXUzrSbwUtVr4/x+ut8rj8nIr8XkXxVrVTVMs/ychF5GlgG9EqujDHGDD67yuuZ\nVZQ1YM8XbnKVoqpdZUA8Z/bSgj3AmMFGVfnq4xspr2/h8VtWkJ02cE3Eg9U3L5zJjqP1/PczW5k5\nJivmpXWNiQYR+RAwG/c8VwCo6vf7uc4xwDFVVRFZhruabpWIpANxqlrvuf5BoF/PZYwxZvgKd2a2\nRhFZ5P1HRBYD0avfaMwAuPvNfby8/RjfvPAkFlgSAUBCfBy/uXohBRnJfPGR9TQGKf9rzGAgIncA\nVwFfwF0x8ApgYhiPexh3AaYZIlIqIjeKyC0icovnLpcDW0RkI+6KuB9T9wydo4F/e5avAf6hqi9E\n/YX1YNMlGGPM0BRuy9WXgcdF5DDuH7MxuH/cjBkS1h44zk9feJ8LZo/h+lMnOR3OoJKXnsT/XbWA\nq+5cxQ+f286PLpvrdEjGBLNCVeeJyCZV/R8R+SXwfKgHqerVIW6/HXep9p7L9wLz+xxtH1U1joyy\n5MYYM9yElVyp6rsiMhOY4Vm0Q1WtVqoZEo43tvH5v6ynKCeVn10xD4lwHpeRYFlJHjefPpk/vrGX\nD5w0esRVUDRDirfXRJNnOpAqYKyD8cREkLl/jTHGDGLhdgsEWArMAxYBV4vItbEJyZjocbmUrzy2\ngaqGNn7/8UVkpdg4q0C+8sHpzByTydef3MTxxjanwzEmkL+LSA7wc9wV/vYDf3E0ohjYdqQu9J2M\nMcYMOmElVyLyZ+AXwGm4k6ylwJIYxmVMVPz5nQO8vqOC/7r4pICz3Ru35IR4fnXlAmqa2vivv25G\nI5xM1JiBoKo/UNUaVX0S91irmar6HafjirbWDhtzZYwxQ1G4Y66WALPUjrbMELKvspEfP7+ds2YU\n8MmTQ453N8Csoiz+8wPT+dkLO3h+y1EumjvseluZIc4zye8jwKOqugewwUnGGGOCcrl0wKbfCbdb\n4BbcRSyMGRI6XcpXH9tAckI8P/2ojbOKxM2nT2bOuCy+88xWaptsaKUZdC4BOoDHRORdEfmaiNgk\nbcYYYwJqbBu4asjhJlf5wDYReVFEnvVeYhmYMf1x5xt7WXewhu+vnM3orIGblXs4SIiP46cfnUd1\nUxs/fG6b0+EY042qHlDVn6nqYuAa3GOB9zkcljHGmEFsIE+yh9st8HuxDMKYaHr/aB3/99JOLpwz\nhg/PL3I6nCFpdlE2N58xmT+8voeVC8Zx6tR8p0MypouITMQ9HchVQCfwdWcjMsYYY9zCarlS1X/h\nrsiU6Ln+Lu4qTcYMKm0dLr762EayUhP430vnWHfAfvjSudMoyU/nm09ttglNzaAhIquBp3H/fl2h\nqstU9ZcOh2WMMcYA4VcL/DTwBPBHz6JxwF9jFZQxfXX7q7vYeriOH142l1EZyU6HM6SlJMbz44/M\n5eDxJm57eafT4Rjjda2qLlLVn3gm+DXGGGOCGshT7eGOufoccCpQB6CquwCbZdQMKhsP1fC71/fw\nkUXjOH+21V+JhpMnj+LqZRO46829bC6tdTocY1DVHU7HYIwxxgQSbnLVqqpds4qKSAJgZdnNoNHU\n1sF/PrqB0ZnJfPeS2U6HM6zceuFJ5Gck8/UnN9He6XI6HGOMMcaYQSvc5OpfIvItIFVEPgA8Dvwt\ndmEZE5kfP/c+eysb+cWV88lOTXQ6nGElOzWR76+cw/Yjddz1pvXCMsYYY8zQMpBD8MNNrm4FKoDN\nwH8AzwH/FaugjInEazvK+fM7B7jptBJWTLGqdrFwwZwxXDhnDLe9vIu9FQ1Oh2NGMBFJE5H/FpG7\nPP9PE5GLnY7LGGOMgfCrBbpU9S5VvUJVL/dct26BxnHHG9v4+hObmDE6k6+dP8PpcIa1//nwbFIS\n4vjGk5vodNnX3zjmT0ArcIrn/zLgf50LxxhjjDkh3GqB+0Rkb89LrIMzJhhV5dYnN1Hb1M7/XbWA\nlMR4p0Ma1gqzUvjeh2fz7v5qbn91t9PhmJFriqr+DGgHUNUmBrYQlDHGmCFGBvBnItxJhJf4XE8B\nrgDygj1ARO4FLgbKVXWOZ1ke8CgwCfe8WVeqanVkIRvj9ud3DvDPbcf4rw+dxKyiLKfDGRE+smg8\nb+6q5Nev7OTkyXksnzzK6ZDMyNMmIql4iiqJyBTcLVnGGGOM48LtFljlcylT1duAD4V42H3ABT2W\n3Qq8oqrTgFc8/xsTsa2Ha/nfv2/nnJmF3HhaidPhjCg/uHQOE0el89mH1nHoeJPT4ZiR57vAC8AE\nEXkI92/J10M9SETuFZFyEdkS4PazRKRWRDZ4Lt/xue0CEdkhIrtFxH63jDHGBBRut8BFPpclInIL\nIVq9VPUN4HiPxSuB+z3X7wcujTRgYxpbO/jCw+vJTU/k55fPQwayBIwhIzmBu69bQluni5vuf4+6\nlnanQzIjiKq+BHwE+BTwMLBEVV8P46H30fuEX09vquoCz+X7ACISD/wOuBCYBVwtIrP6Fr0xxhgn\nDMZqgb9HEUNjAAAgAElEQVT0ufwYWAxc2YfnG62qRzzXjwKjA91RRG4WkfdE5L2Kioo+PJUZrr77\n7Fb2VTZy21ULGZWR7HQ4I9KUggz+8PHF7Klo4Po/vUtDa4fTIZlhzvckHzAROAIcBoo9y4IKcMIv\nHMuA3aq61zPf4yO4TxQaY4wxvYQ15kpVz472E6uqikjAkmOqeidwJ8CSJUusNJkB4On1pTyxtpQv\nnjuNU6bYeB8nnTYtn99evZDPP7yeG+57l/uuX0paUrjDOI2J2C+D3KbAOVF4jhUisgl3BcKvqepW\nYBxwyOc+pcByfw8WkZuBmwGKi4ujEI4xxpihJqwjIRH5SrDbVfVXYT7fMREZq6pHRGQsUB7m44xh\nb0UD3356C8sm5fHFc6Y6HY4BLpw7lttcypceWc9N97/HPdctJTXJqjaa6IvFSb4e1gHFqtogIhcB\nfwWmRbKCaJ4ULMhIpqLB6nQYY0w0DGRF6XC7BS4BPoP7DN444BZgEZDpuYTrWeA6z/XrgGcieKwZ\nwVo7OvnCw+tJSojj11cvICE+3I+uibVL5hfxyyvns2pvFZ9+4D2a2zqdDskMYyKSIiJfEZGnRORJ\nEfmyiKT0d72qWqeqDZ7rzwGJIpKPuxVrgs9dx3uWxZTt44wxZmgKtw/PeGCRqtYDiMj3gH+o6icC\nPUBEHgbOAvJFpBR3haefAI+JyI3AAfo2bsuMQD95/n22Hq7j7muXMDY71elwTA+XLRxPpwu+/sRG\nrvvTGu791FIykq2LoImJB4B64Lee/68B/ox7ipA+E5ExwDFPl/VluE8+VgE1wDQRKcGdVH3M85zG\nGGNML+Ee/YwG2nz+byNIMQoAVb06wE3nhvmcxgDw0rZj/Omt/Vx/6iTOmxX0Y2ccdPni8STGC195\nbCPX3rOa+25YRlZKotNhmeFnjqr6Vut7TUS2hXpQgBN+iQCqegdwOfAZEekAmoGPqaoCHSLyeeBF\nIB641zMWyxhjjOkl3OTqAWCNiDzt+f9STpRUNyZmDtc08/+e2MjsoixuvXCm0+GYEFYuGEdyQhxf\neHg9n7h7NQ/csIyctCSnwzLDyzoROVlV3wEQkeXAe6EeFOSEn/f224HbA9z2HPBcH2LtM3deZ4wx\nZqgJdxLhHwLXA9Wey/Wq+qNYBmZMR6eLLz+ygfYOF7dfs4jkBCuUMBRcMGcsd3xiMe8fqefqu1ZT\nZYPyTXQtBt4Wkf0ish9YBSwVkc2eSn/GGGOMYyIZMZsG1Knqr4FST/9zY2Lmh89tZ83+4/zvZXMo\nyU93OhwTgXNPGs3d1y1hb0UDV9/1DuX1LU6HZIaPC4AS4EzPpcSz7GLgEgfjiqpOa7kyxpghKazk\nSkS+C3wD+KZnUSLwYKyCMubRdw/yp7f2c8OpJVy2cLzT4Zg+OGN6Afddv4zS6mY+9sd3OFLb7HRI\nZhhQ1QNAHZANjPJeVPWA57ZhwcYrGmPM0BRuy9VlwIeBRgBVPUxkJdiNCdsr24/xrae3cPq0fL51\nkY2zGspOmTKKB25YRnl9K1fcsYr9lY1Oh2SGOBH5AbAJ+A3uiYV/CfzC0aBiYHZRltMhGGOM6YNw\nk6s2T9UkBRAR66NlYmLVnio++9A6Zhdl8YdPLLa5XoaBJZPyePjTJ9PY2sHld6zi/aN1TodkhrYr\ngSmqepaqnu25nON0UNEmIk6HYIwxpg/CPXJ9TET+COSIyKeBl4G7YheWGYne2FnB9fetYUJeGvdd\nv8zmSRpG5o7P5vFbTiEhTrjqj++w7mC10yGZoWsLkON0EMYYY4w/4VYL/AXwBPAkMAP4jqr+Nvij\njAnfMxvKuOn+9yjJz+CRm08mL93Kdw83UwszefyWU8hNS+Tae9aw4VCN0yGZoenHwHoReVFEnvVe\nnA7KGGOMgTDmuRKReOBlVT0beCn2IZmRxOVSfvXSTm5/bTfLJuVx57WLbV6kYWxCXhqP/scpXHHH\nKq67dw2P3HwyJ421sSUmIvcDPwU2Ay6HYzHGGGO6CdlypaqdgEtEsgcgHjOCNLZ28JmH1nL7a7u5\naskEHrxpuSVWI8DorBQeumk5aUnxfPKe1eytaHA6JDO0NKnqb1T1NVX9l/fidFAmOsZkpVCYmdKv\nxztp+ujo1fpaOCE3ausajEalJzsdgjExEe6YqwZgs4jcIyK/8V5iGZgZ3spqmrn8jlW8tO0Y/33x\nLH7y0bkkJVjxipFiQl4aD960HICP372aQ8ebHI7IDCFvisiPReQUEVnkvTgdVCycP3tM2PeNC6MA\nxsRRg78WVUFmMkkJvV9LUpjFjZZPHhX09tERJF+nhFhXT1MKMpg5JnRyVZSTGtb6ikelRfT8XtMK\nM7lo7lgWTMjpdxf7aCaLE/K6v57s1MR+JdImek6dmu90CMNKuEezTwH/DbwBrPW5GBOxtQeqWXn7\nW5Qeb+LeTy3lxtNKrDLWCDSlIIMHblhOY2sHn7hnNcfqbKJhE5aFwMnAjxjGpdgBUhLjI7p/fkYy\nUwoyAt4+f3zkHVAWFUe/9aQgM3CLRaDfgnNOKuz2v28vh7HZwZMV39u9a88No5dEYYStYKlJ8WH9\nli2dlMfKBeMiWnckEuOFxPg4Jo5KJznhxGdo5pjQXbB7xjW5IHRCnhzkxKjvZ3h8bu/3aUz2iW2c\nOEDVgbNTw5tD7pyZhSGTv4S46MZ8/uwxnD2zMPQdfSwryQvr8xzIygXjyM+IXSvipAhO6kQjmR8M\nxdCCfipEpBhAVe/3dxmYEM1w8tS6Uq6+8x3Sk+N5+nMrOGtGZDsRM7zMKsri/huWUVnfyifuXk1V\nQ6vTIZlBzqf8+tnDuRR7KP4Ohk6dms+ccYETqEAH/stLRlEQ4OBqQl4aKxeMi2oy4D0gXTChd9HH\nURlJTM4PnCB6jc9N5eyZhXxw1hiWleR1u+3COWO7/V+QmUx6UvcDrtSk8BLXSH6jvFt3WUkeH5g1\nOuzH9Yw/mnLSTiQSCfGxOYk5KsiB+VkzCrquh0pUZoTR6hdtp03N7/oc9kz6M8OYyHv+hOiOmElJ\njO82gXhifFzIkyy+Jw+WleR1JY/LS3q3vPb8rKWGWHeohDctKXQiMzY7vJMUKxeM69cYbG8svp95\np4RKuf/qvSIiT8Y4FjOMqSq/e203X3lsI4sn5vLXz57K1EKbh9rAwuJc7r5uKQePN3HtvWuobW53\nOiQzyInIh0Tk6yLyHe/F6ZgG0gVzxrBiSmRd1gI5Y1oBY7JTWBFGt6BALUSRnm2OC3KMn5WSSK5P\nV7ZgrSJZKYldSVJyQhzzx7sPkoN1MfceqI7OTOGMae4Df99xWj1fY3ZqIhfOGctFc7snbP54k9ex\n2amkJSVwSo/3qDgvjRVTem/nsdmpnNaP7d/zsb459ITcvnUthPDHRGUFSUJ8W8566pnM9CfWnsLt\nRpqXnsTEUemsXDCu2+fAm4yHSkh9W8FCnRSItJup1zkzC5kzLpvivMDbx9tylZmSyNJJeUwalc7o\nrMDvX7itd9MK3a8p0HvsngI3fOfPHsOsKBWx6nnCxPuaMpIHf3Ll+6maHMtAzPDlcik/+Pt2fv7i\nDi5dUMT9Nyzr9uNpzClTRnHHJxez81g919z1Dkdqm50OyQxSInIHcBXwBdy/UVcAEx0NKoZ6HoDm\npiWRnNC7+1m4Pat7tlBFsi9ePDG31wFiQUZyRGebT548inif7OqsGYXM9dPaNj43lYLM5K7WAyH4\nC7xgzlgm5YfufjS5IJ1FxblMyEslNz2JlQvGdRun5e9gMSkhrk9d1nxbamYXZbOwODdgl8i89CSm\nFWYGTWgWT+zdRXNsdmrQlqP+OG1aeONwQrV++HPxvCLGZKd0294i/hO6SAtd5Wckc/7sMQFbW33f\n4UCtud7Fc4q6fza9r3VqYQaT8zO6tW7NDdHtNlQ300DdARPj44J29wWYXZTFmdMLyEhOID05gfkT\ncvy+tp4f7/4OyUj0cyLjvJNGMzk/o1f3vPg4ISUxnml96Pq3eGJur+6F8T3O0niTq8EwlU+ovYUG\nuG5MWNo7XXz18Y3c+9Y+rj91Er+6coEVrjB+nT2jkDs/uYT9lY2svP0t1uw77nRIZnBaoarXAtWq\n+j/AKcB0h2OKmdljux+wnTG9IMA9/UtLSuDM6QVdXZ/SfQ54Ij1ojY+T3t3pPMc3Z80o5PzZY5jq\nOdNdlJPKND+9E3oWlMhOTWSynwPHxRPzWDEln6WT8jhlyqio/W7EiTAhL63XQeXiibkx7abu3S6+\nLpo7tqsLo4gwqyiL1KTur3PJJHc3rqzUxF4HkyX56Szxk3AFGgeUGOXxQb4iKb4CvQ+MwZ1ELJ/c\nvduaN2n4kE/L4enTCro+W/FxQkmPpHp8bipxwZpHPYJ1x/R+PHzjPHN6QbcD+FDJlC/v+xhIfkZy\n0BbAUOLixO/3uWehG29LnPd1+PtceuMByPZ0sctJS+TUqfmc2WP/Mz4ntVt3zkmj0klPTmDu+Oyu\npD8+Tjh7ZiHnnXSiq+yFc8Zy8bwiTp2a32u943wKvkzOzyA7NZHxuWnM99ONGGDFlHxOn1bA9NEZ\nnDo1v9sJjHC6LcZCqGedLyJ1uHefqZ7reP5XVbUJakxAzW2dfPahtby2o4L/d/4MPnvWFCtcYYI6\ne2YhT332VG7+83tcdecqPnPmFL583nRLyI0vb7Nmk4gUAVVA6D5bQ1R2WiJxIrgi7H7jdfq0fFIS\n47sOvHLSEt1bDPrVtTAjOYHRWSldZ9R9D9b2VTYypSCDo7X9L1KTlBDX1QKUFB9HW2d4U5ulJMbT\n0t5JZkoCRTkp5KQmsvNYPWkBxlqN93RJG5OdwtEBKq7jrzVs2uhMSqtPtNwHGpwfqFUmOSGOiT5V\nBtXnvHheRuzO6EdafMXLm7x4D/J7bpMCz3uf4LM8Lz2JvPQkOlwuJuSmUdXoHqubn5HM/Ak53bZZ\namI8ze2dZKcmUtvczswxWRz29IxISwx8COxtKU306RYY7GSEv2Idvvx1bx2dlcLkgnRW7anq1fKc\nnBBPa0dn0HXOH5/jN0n1dcGcMRypaWH9oWrA3Zq6qDiXopxUFgYpVuNt4U1NiuesGYVkJif0Slgv\nnDOWpIQ4Wto72XG0vtc65o7LJj8jyW/Lqvc33d/Y0SWT8ijbUOZeR5AE1vvZTk6M60pMe67PqSPO\noMmVqvbt22JGvJqmNm68/z3WH6zmR5fN5ZrlxU6HZIaIGWMy+ccXT+cHf9vG71/fwxu7Kvi/Kxf0\nqSuBGZb+LiI5wM+Bdbh7VdzlbEixtawkj3f2VoV9/7OmF/L6znKg90HvxFHpbDhUA/SvOpsIfotn\nJCfEc/G8IoCu5CozJYE5Rdn9rgZ3zkmFrD1QTUV96MI3Z04voLmts6vbY3JCfMgy7dB9+zghKyWR\nS+YV0eHqWzI9KT+920lMb06emhgfdhW1hRNySfRTDr8wM4Xy+pau68mJcX6n0fC+/+GYkJtGS7uL\nKT5VCU+ZMoq65vaQ47LnecbYVXoKIeWkJfZ6jbnpSTTXNJOblkRtc3v3Low9XmJxXhoV9a0crWvp\nSvbCOSF88byirnGEk/Mz2FsZeu7GjOQEZo7JpLXD5Qml+/OcPbOAF7YcDbqO3LSkrpalQBLj4yge\nldaVXEHvkvheM8ZkdiVJwolWan/jszJTEkKe9IyPk66TFrE0GE/ZO1+v0Aw7x+pauPaeNeyrbOR3\n1yziwjAGAhvjKyM5gZ9ePo+zZxbyzac2cdFv3uS6UybxhXOnhT0Q1wxPqvoDz9UnReTvQIqq1oZ6\nnIjcC1wMlKvqnCD3WwqsAj6mqk94lu0H6oFOoENVl/TvVURmdFYKs4uyaGkP3GrjOx4h1AFXIKdO\nzeet3ZV9emww43PTuo03mVWU1bU8EskJ8WQmJ1JR3xrygColMb7PrSmRWDwxl2N1rZRW+5+rLy89\nieONbRGtMy5OSAqjW5uv9KQEGts6ei33bgNv160547JJiBOK89J4duNhwP358p0Kw3d+LW9ylpwQ\nxylTRnGgqpENh2pITYrvNT4tJy2Jmqa2oK0po9KTu1qavK+1Z5XAwkz/E0kX5aRyuKb3eFxvAu1v\nvNaCCTkU5aRSlJ1CfmYy43JSSUmMY/2hml6JWEJ8XMAEPNg8cr6vtyAzOWhyVZyXxsHjTZzr6SJX\n7fls9PxdS4p3t9hOKYzO3HQfnDWGzhCt3zPHZFHb1B6y5fbieUUBv3/hVgYM5czpBd1aK71OGpvF\n+0frOXlyHptLA+/2T52aT1pSPGsPVNPYRrcuiQPBkisTVRsP1fCZB9dS29zOfdcvDasClTGBXDBn\nDIsn5vKLF3dwz1v7eHp9GV/94AyuXDLe747XDF+epOeQqh71/H8t8FHggIh8T1VDDdK7D7gdeCDI\nc8QDPwX+6efms1U1+plHmIKdxb9kXlHYBS2Cyc9IZnxuKrlpSWwu83/gkp6UQEFGMjP7UfErOSE+\naJekoSA3LYnlk/NIToinqiFw8rRiSj6dfWyFAveJpuSEOGaH2N7j81LZcbS+VwtIfJx060LorzBC\nz+Sqr84IowBGuEUy/FkyMRf187nJz0ju6qLWU2J8XNcYHu/fwqyUiMaIfWju2K4WrMkFGRytawk4\nr9SY7BQ+OGsMqUnxVNS38vYe9y7D+64sLM7t9tnPTU/itKn5vYowiEivapO9Wn4j+M6HO/VAOHom\nz97EMy89KeK54QIJ1AVz+ujMruqkcZ6RSv5aF73dA5dOyuNYXUu3saYDwY5OTFS4XMo9/97HFXes\nQkR49D9OscTKREVBZjI/vXwef/v8aUwpyOBbT2/mwl+/yUvbjkVcBtYMaX8E2gBE5AzgJ7gTpVrg\nzlAPVtU3gFAJ2BeAJ4HyfkU6QDKSE5g3Poe4OOl1gDF9dKbfeW5CWTwxz2+BCa+4OGGFn4NBfyJJ\n+M6eWcjSEIP+B5OEeAlaZtwrPk76NWY0Pk64YM7YroPWnLQkMlN6HyiOz00jToRxIcb++BNJfONz\n0yjOS+Oksb2TfZHen8NoEpGAhSpiOS43IT6uK6EoyExm5YJxQVtFvYlMQWZyWN+TURnJYW23k8Zm\nMW98Tr8KX8RCUkIcy0tGxXS+Nn+WleQxY0xm0O6uKYnxTIxgEuNosZYr029bD9fyvWe38u7+as6d\nWcgvrphvpdZN1M0Zl82j/3EyL249ys9e3MGnH3iPJRNzufXCmSErMZlhId6ndeoq4E5VfRJ398AN\n/V25iIwDLgPOBpb2uFmBl0WkE/ijqvpN5kTkZuBmgOLi2I8zPTdIV5dQ5dHDGf8UTvIQLVkpiWEd\nNCYnuuPu7/itoapntTavjOQELpkf/lgncHcl3V/VyLicVFzFud2qtPXmPviPj5OulhdvkpeSFPq9\nCFTBcLgbk5XC8ca2qFSt81ZGTIwX1h6oDlicxQljotQdMBLpyQnMHDM46+pZcmX67EBVI7/8506e\n3XiYnLREfnnFfD6yaJxVBDQxI+I+i3veSaN57L1Sbnt5J5ffsYoPzBrN18+fYUUvhrd4EUlQ1Q7g\nXDxJjEc0fstuA76hqi4/+7DTVLVMRAqBl0TkfU9LWDeepOtOgCVLlgzaZtUL54wN2ark2xVqMJla\nkEFKQnzI6myxlukzUenorBT2VzWSF2Fpe6fNn5DTVd46UJED72S8/lqqphRkkJvmvxqcrw/OGhOy\nqt1wNbUwg+JRaVE9UTE+N21ACkWYvrPkykSsvL6F21/dzV9WHyQhXvjc2VO4+YwpVmjADJiE+Diu\nWV7MpQuL+NNb+7nj9T2cf9sbXLF4Al88b1qIM7BmiHoY+JeIVOIux/4mgIhMxd01sL+WAI94Eop8\n4CIR6VDVv6pqGYCqlovI08AyoFdyNVSE04UqGmMaE+JOtHZES1ycdCu4MNCyUhOZXZRFvk/xhDHZ\nKVwyryisuZWGmrgeY7Z8iUhYExhHc7zPUCMSXvdRM7xYcmXCVtfSzl1v7OXuN/fR3uniY8sm8MVz\npkVtAKMxkUpLSuBzZ0/l6mXF3P7qbh585wBPry/j6mUT+NzZU+2zOYyo6g9F5BXcc1r9U08MuIvD\nPVaqv+sv8V4XkfuAv6vqX0UkHYhT1XrP9Q8C3+/v840EUwoyEBFKHBjzEAveKmn+kqjhmFgZY/rG\nkisTUkt7Jw++c4Dfvbab6qZ2Lp43lq99cAaT8ofHD6YZ+vLSk/jOJbO46fQSfvvqbh5afZBH3j3E\n5YvH88lTJg7aftkmMqr6jp9lO8N5rIg8DJwF5ItIKfBdINGzjjuCPHQ08LSnRSsB+IuqvhBZ5CNT\nXJx0zRc0HIzUrm1m5MhOS+RoXQspiSNzjFy0WHJlAmps7eDhNQe56829HKtr5fRp+Xzjgpl+J440\nZjAoyknlxx+Zyy1nTub3r+3hibWlPLT6IDPHZPKBWaNZVJzLggk5VnBlBFLVqyO476d8ru8F5sci\nJmOMGUxmjM5kTFZKwFLoJjyWXJleymqaefTdQzywaj81Te2cMnkU/3flAiutboaMiaPS+enl87j1\nwpn8dUMZ/9h0hN+9thvvdDOTRqUxe1w2s4uymFOUzfzxOX2eeNUYY4wZDkTEEqsosOTK0NrRydbD\ndazaU8XrO8p5d3814J7R+rNnT2HREJ/s0YxcuelJXH9qCdefWkJjaweby2pZf7CGDYeq2Xiohn9s\nOgK45+OZMTqTkye75+pYOimPgszQA7WNMcYYY3xZcjXCqCrH6lrZcKiadQdrWHugms1ltbR1uACY\nMy6Lr3xgOpcuGOdoRSZjoi09OYGTJ4/i5MknJlatbWpn6+Fa1h6oZs3+4zz67iHue3s/AJML0lle\nMorlJXksK8mjyCoQGmOMMSYES66GMZdLOVLXwvbDdWwqq2VLWS2bSmupbGgF3PNXzB2fzXWnTGTx\nxFwWT7Sz9WZkyU5LZMXU/K4ur+2dLraU1bJ633HW7DvO3zcd5uE1BwGYkJfKwgm5TCvMYKrnMnFU\nelhlrY0xQ8cZ0wqIj7fiFcaYvnEkuRKRC4BfA/HA3ar6EyfiGIo6Ol3UNLdT09RGdVM71Y1t1DS1\nU93U1rW8sqGNA1WNHKhqotXTIhUn7snszpiez9xx2cwbn8OccVk2/4IxPhLj41hYnMvC4lxuOXMK\nnS5l+5E61uw7zup9Vaw9UM2zGw933T8+TijOS2N8bqrnksa4nFSKclLJSUskO9V9SUm075kxQ4UV\nvDHG9MeAJ1ciEg/8DvgAUAq8KyLPquq2gY4F3N3kvLOlqM8y7/8nbvMsU9/Hnrit06V0dCodLqXD\n5TpxvdPl+etZ7lLaO13d7t/a0UlTWydNrR00tXfS1NpJQ2sH1Z4EqqbpRAJV39IR8LUkxLkHIual\nJ1Kcl86Z0wsoyc9g+ugMZhVlkZZkDZXGRCI+TpgzLps547K54TT3NEiNrR3srWhkd0U9u8sb2FfZ\nSFl1M/88XEdVY5vf9SQnxJGVmkhqYjwpiXGkJMaTkhBPsud6coJnWWIcKQnxJ64nxpOcGE9KQs/7\ndV9PSmIcSQlxiAhxwom/CCLuMWVxEvhMvO9+LZBu6w2yLmOMMWYkc+Joexmw21PeFhF5BFgJxCy5\nOueXr7OvshEI7yDCaSKQkZRATnoiuWlJ5KQlUZKf7rme2O2v93pOWiIZyQl20GNMjKUnJzB3fDZz\nx/eekqC5rZOymiaO1LZQ29xOTVM7tc3uS11zO83tnbS2u2jp6KSlvZP6lg4q6ltp7XDR0u5e5r3u\nGqT7qmuWF/Ojy+Y6HYYxxhgzKDmRXI0DDvn8Xwos73knEbkZuNnzb4OI7BiA2PzJByodeu5IWayx\nYbHGhsUafTGP88eeSz9N7P8qBre1a9dWisiBfq5mqHzuBoptj+5se3Rn26M72x4nRGNbhP27NWj7\nianqncCdTschIu+p6hKn4wiHxRobFmtsWKzRN1TiHAlUtaC/67D3szvbHt3Z9ujOtkd3tj1OGOht\n4USZqzJggs//4z3LjDHGGGOMMWbIciK5eheYJiIlIpIEfAx41oE4jDHGGGOMMSZqBrxboKp2iMjn\ngRdxl2K/V1W3DnQcEXC8a2IELNbYsFhjw2KNvqESpwmPvZ/d2fbozrZHd7Y9urPtccKAbgvRoVA+\nzxhjjDHGGGMGOSe6BRpjjDHGGGPMsGPJlTHGGGOMMcZEwYhNrkTkXhEpF5EtIe63VEQ6ROTyHsvj\nRWS9iPw9tpH2L1YRyRGRJ0TkfRHZLiKnDOJY/1NEtorIFhF5WERSnIxVRM4SkVoR2eC5fMfntgtE\nZIeI7BaRW2MZZ39iFZEJIvKaiGzzbNsvDdZYfW4fNN+tEJ+BAftu9TPOAf1emegY6H3MQAm0TxKR\nPBF5SUR2ef7m+jzmm57tsENEzvdZvlhENntu+42IiGd5sog86lm+WkQmDfTrjETPfd5I3hbgf986\nUreJv/33SNoW/n77Bur1i8h1nufYJSLXRRS4qo7IC3AGsAjYEuQ+8cCrwHPA5T1u+wrwF+DvgzlW\n4H7gJs/1JCBnMMaKe3LpfUCq5//HgE85GStwlr/31xP/HmCyZ5tuBGYN0ljHAos81zOBnYM1Vp/b\nB813K1isA/nd6sf7P+DfK7tE5f0e8H3MAL42v/sk4GfArZ7ltwI/9Vyf5Xn9yUCJZ7vEe25bA5wM\nCPA8cKFn+WeBOzzXPwY86vTrDrFNuu3zRvK28MTZa986ErdJoP33SNoW+PntG4jXD+QBez1/cz3X\nc8ONe8S2XKnqG8DxEHf7AvAkUO67UETGAx8C7o5NdN31NVYRycb9wbzHs542Va2JVZye5+jzdsVd\nvTJVRBKANOBw9CM8IcxY/VkG7FbVvaraBjwCrIxqcD30NVZVPaKq6zzX64HtuHfYMdOP7TpYv1u9\nDPR3qz/blAH+XpmoGPB9zEAJsk9aifugGs/fSz3XVwKPqGqrqu4DdgPLRGQskKWq76j7aOiBHo/x\nrjS5djgAAARHSURBVOsJ4FzvmerBJsA+b0RuCwi6bx2p28Tf/nvEbIsAv30D8frPB15S1eOqWg28\nBFwQbtwjNrkKRUTGAZcBf/Bz823A1wHXgAYVQJBYS4AK4E+eLgd3i0j6gAfoI1CsqloG/AI4CBwB\nalX1nwMfYS8rRGSTiDwvIrM9y8YBh3zuU0qME5Yw+Yu1i6e5eyGweqAD8yNQrIPqu+XhL9ZB993C\nT5yD+Htlghus+5io6rFPGq2qRzw3HQVGe64H2hbjPNd7Lu/2GFXtAGqBUVF/AdHhb583UrcFBN63\njrhtEmT/PeK2RQ8D8fr7tQ+25Cqw24BvqGq3gzwRuRgoV9W1zoTll99YcZ/xWAT8QVUXAo24m1Cd\nFGi75uI+g1ACFAHpIvIJB+LztQ4oVtV5wG+BvzocTzBBYxWRDNythV9W1ToH4vPlN9ZB+t0KtF0H\n23cr0DYdjN8rY4Lukzxnl4f9PDHh7PNGyrbwEXLfOlK2STj775GyLQIZrK/fkqvAlgCPiMh+4HLg\n9yJyKXAq8GHP8keAc0TkQceidAsUaylQqqrelooncO+0nBQo1vOAfapaoartwFPACufCBFWtU9UG\nz/XngEQRyQfKgAk+dx3vWeaYILEiIom4D2IeUtWnHAwTCBrroPtuBYl1UH23gsQ56L5XJiyDbh8T\nTQH2Scc83Xfw/PV2Gw+0Lco813su7/YYT3eqbKAq+q+k3wLt80bitvAKtG8didsk0P57JG4LXwPx\n+vu1D7bkKgBVLVHVSao6CfeX+7Oq+ldV/aaqjvcs/xjwqqo6eiY4SKxHgUMiMsNz13OBbU7FCYFj\nxd3sfbKIpHn6u56Luy++Y0RkjE9FmWW4vy9VwLvANBEpEZEk3J+DZ52LNHCsnmX3ANtV9VdOxugV\nKNbB+N0KEuug+m4F+awOuu+VCcug28dES5B90rOAtyLXdcAzPss/Ju6qXiXANGCNp1tQnYic7Fnn\ntT0e413X5bj3JYPu7HaQfd6I2xZeQfatI3GbBNp/j8Rt4WsgXv+LwAdFJNfTgvhBz7Lw6CCoBuLE\nBXgYdx/WdtxnSm4EbgFu8XPf++hRLdCz/CwGpqJZn2MFFgDvAZtwdxUKu9qJA7H+D/A+sAX4M5Ds\nZKzA54GtuKvPvAOs8HnsRbirXO0Bvu30ZyBQrMBpuJvMNwEbPJeLBmOsPdYxKL5bIT4DA/bd6mec\nA/q9skvU3vMB3ccM4Ovyu0/CPc7hFWAX8DKQ5/OYb3u2ww48Vb48y5d4Ptd7gNsB8SxPAR7HPaB9\nDTDZ6dcdxnbp2ufZtui9bx2p28Tf/nskbQv8//YNyOsHbvAs3w1cH0nc3pUbY4wxxhhjjOkH6xZo\njDHGGGOMMVFgyZUxxhhjjDHGRIElV8YYY4wxxhgTBZZcGWOMMcYYY0wUWHJljDHGGGOMMVFgyZUx\nxhhjjDHGRIElV8YYY4wxxhgTBf8f+2BlGcKgsqcAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x7f7d98d58210>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "pymc3.traceplot(trace)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 2", | |
| "language": "python2", | |
| "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 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment