Skip to content

Instantly share code, notes, and snippets.

@tillahoffmann
Last active May 4, 2016 07:55
Show Gist options
  • Save tillahoffmann/8e7c69e9cdde3a438f0bcfaeda2ca74c to your computer and use it in GitHub Desktop.
Save tillahoffmann/8e7c69e9cdde3a438f0bcfaeda2ca74c to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is an implementation of [Bayesian linear regression](https://en.wikipedia.org/wiki/Bayesian_linear_regression#Posterior_distribution) to illustrate the overestimation of precision if the number of features $p$ is comparable with the number of observations $n$. We first consider a simulation with $n\\gg p$ to test the inference algorithm, and then consider a simulation with $n \\gtrsim p$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"from scipy import special, stats\n",
"np.seterr('raise')\n",
"\n",
"def simulate(n, p, precision):\n",
" \"\"\"\n",
" Simulate data for linear regression.\n",
" \"\"\"\n",
" # Generate features\n",
" X = np.random.normal(0, 1, (n, p))\n",
" # Generate coefficients\n",
" theta = np.random.normal(0, 1, p)\n",
" # Generate observations\n",
" y = np.dot(X, theta) + np.random.normal(0, 1 / np.sqrt(precision), n)\n",
" \n",
" return X, theta, y\n",
"\n",
"n, p, precision = 1000, 10, 1.5\n",
"X, theta, y = simulate(n, p, precision)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def infer(X, y, **prior):\n",
" n, p = X.shape\n",
" default_prior = {\n",
" 'precision_theta': 1e-3 * np.eye(p),\n",
" 'mean_theta': 0 * np.ones(p),\n",
" 'shape_precision': 1e-3,\n",
" 'scale_precision': 1e-3,\n",
" }\n",
" \n",
" if prior:\n",
" default_prior.update(prior)\n",
" \n",
" # Compute posterior parameters for the Gaussian part of the posterior\n",
" precision_theta = np.dot(X.T, X) + default_prior['precision_theta']\n",
" cov_theta = np.linalg.inv(precision_theta)\n",
" estimate_theta = np.dot(cov_theta, np.dot(X.T, y) + \n",
" np.dot(default_prior['precision_theta'], default_prior['mean_theta']))\n",
"\n",
" # Compute parameters for the gamma part of the posterior\n",
" shape_precision = default_prior['shape_precision'] + 0.5 * n\n",
" scale_precision = default_prior['scale_precision'] + 0.5 * \\\n",
" (np.dot(y, y) - estimate_theta.dot(precision_theta).dot(estimate_theta))\n",
" \n",
" return estimate_theta, cov_theta, shape_precision, scale_precision\n",
"\n",
"estimate_theta, cov_theta, shape_precision, scale_precision = infer(X, y)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAGsCAYAAAAv5DqZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcVOWV//HPaUBwYRMEBQEFQQXcUIkbUopGjRqSqDFO\nlkkwGTMx20STmN84093ZzSSZmKhRE+NE4yROEvcl0URbMRp3UBBstkZFBJWtWaIs5/fHrZKy6equ\n7rq3nltV3/frVa+uqnvr3lMUXE4/dZ7zmLsjIiIiIiLb1YUOQEREREQkbZQki4iIiIi0oSRZRERE\nRKQNJckiIiIiIm0oSRYRERERaUNJsoiIiIhIG0qSRURERETaUJIsIiIiZWNmc8zs+FL3EUmaaTER\nERERMbMWYAiwBdgA/Am40N03hoxLJBSNJIuIiAiAA6e7ez9gEnAEcGnbncysR7kDEwlBSbKIiIjk\nGIC7LwfuBSYCmNkSM/uamc0G1ptZnZntZWZ/MLOVZrbIzL7wzkHM9jazP2a3vW5mP83btsTMTsze\n/7qZvWJm68xsnpmd0Haf7OMDzOxBM1ttZs+b2ZltjneRmc3Obv+tme2U8J+T1AAlySIiIvIuZjYC\neB/wTN7THwFOAwYQjTrfCTwL7AVMA75kZiebWR1wF7AEGAkMB37XzjnGARcCh2dHr08BWtrZr2f2\nXH8C9gC+CNxkZmPzdjsHeC+wL3AI8MnuvXOR7XqGDkBERERS4zYz2wKsJUp0v5e37XJ3fxXAzCYD\ng939O9ltLWb2S6JEupUocf6au2/Lbn+0nXNtBXYCJprZm+7+UoGYjgJ2dffLso8fNLO7gPOAb+bF\ntiIb253AoW0PYmb9gMuBQUTJ9BLgbeBj7v6Pgn8iUrOUJIuIiEjOdHd/sMC2V/LujwKGm9mq7GMj\n+nZ6JjACWJqXILfL3ReZ2ZeBBmC8mf0Z+Iq7v9Zm12HAy22eW0o0Qp2zIu/+RqIkva1JwKezr8u4\n+w0dxSeicgsRERHJsQ625bfDehlY7O67Z28D3b2/u5+R3TYyW3bRIXf/nbtPIUq6AS5rZ7dXiRLv\nfCOBZZ0dv825mtx9K3AW8GRXXiu1SUmyiIiIdNUTQGt2Ml8fM+thZhPM7IjstuXA981sFzPrbWbH\ntD2AmY0zsxOyk+zeBjYRlWC09TiwMXuunmaWAc4AftvN2N/r7vO6+VqpIUqSRUREBN49Utzhtmwp\nxRlEtb9LgJXAL4B+2W1nAmOBl4hGlj/czrF6A98HXicaLd4D+H9tz+fum7PHex/wBnAF8HF3X1BE\n3O9iZrsRlWOIdEqLiYiICGbWH/glUcuvbcAMd388bFQiIuFo4p6IiEA06/8edz8n23Jrl9ABiYiE\npJFkEZEal22N9ay7jwkdi4hIWqgmWURE9gXeMLPrzewZM7vWzHYOHZSISEhKkkVEpCdRD9kr3X0S\n0cSmS8KGJCISlmqSRUTkFeBld38q+/gPwNfb7mRmqs8TkYrj7h31/y5II8kiIjUuu5zvy2Y2LvvU\nNOCFAvum7lZfXx88hkqJq37q1OAxVMqfleKqjrhKoZFkEREB+CJwk5n1AhYDnwocjyThoYdCRyBS\nMZQki4gI7j4bODJ0HCIiaaFyCxERqWiZTCZ0CO1KY1yZ0AEUkMY/K1BcXZXWuLpLfZJFRKQoZub6\nP6PCmYE+Q6khZoZr4p6IiIiISDyUJIuIiNSK+vrQEYhUDJVbiIhIUVRuISKVRuUWIiIiIiIxUpIs\nIiIiItKGkmQRERERkTaUJIuIiIiItKEkWUREpFY0NISOQKRiqLuFiIgURd0tqoAWE5Eao+4WIiIi\nIiIxUpIsIiIiItKGkmSpCWY2zsyeNbO1Zvb57ONZ2cdfMLPnzez4Io4zp5j9REREpLKpJllSxcz+\nCfg34ABgHTAL+K67/63E4/4SWOvuF7X3uNzMbAlwvrs/EOL8It2hmuQqoJpkqTGqSZaqYGZfAX4M\nfBsYAowErgTOjOHwo4C5HTyuKGZ2jJnNMLOLzGxA6HhEpELU14eOQKRiKEmWVDCzfkAj8Dl3v93d\nN7n7Vne/x90vye5zgJk9aGars+URZ+a9fi8z+4OZrTSzRWb2hbxtfwVOAK4ws3V5j6/MPh5rZkvM\n7MTs/nub2R+zx3rdzH6ad6z8/QqeM2/fi8xstpmtMbPfmllvM7uB6BeAO7Pnvzi7/9fN7JXsc/PM\n7IQCf1ZjgE+6+6+Al4CzS/4ARKQ2qAWcSNGUJEtaHA30Bm5rb6OZ9QTuBP4E7AF8Ebgpm+Badtuz\nwF7ANOBLZnYygLtPA2YCn3f3fnmPL8w+XpB3njrgLmAJUSI7HPhdO/F0eM485wDvBfYBDgH+2d0/\nQZTcnpE9/w/NbBxwIXC4u/cDTgFaCvxZXQZck70/DhhcYD8RkaC2bYPLL4cf/xjeeit0NCJdoyRZ\n0mIQ8Ia7byuw/ShgV3e/zN23uPuDRMnsecCRwGB3/0529LkF+GV2W1dNJkp6v+bu/3D3t9390Xb2\nK/acl7v7CndfQ5RUH5q3Lb9GaiuwEzDRzHq6+0vuvqTtSc1sr2yMh5vZBcDHiBJuEZFUWb4c3vte\n+MMfoKkJJkyA29odBhFJp56hAxDJehMYbGZ1BRLlYcDLbZ5bSjTSOwoYbmarss8b0S+AD3cjjhHA\n0g6S9Zxiz7ki7/5GogR8B+6+yMy+DDQA483sz8BF7r68za4nAne5+7VmtgvwA+C+TmIVESmr5maY\nOhUuuAAuvRR69oT774fPfCaaOzh9eugIRTqnkWRJi8eAt4APFNj+KlECm28ksIwoeV7s7rtnbwPd\nvb+7d2fC38vAyGzZRWf7lXLOHaaXu/vv3H0KUQIO8P12XjcCmJe9Px24x93fKPKcIiJl0dAAX/hC\n9LNndjju5JPhqqvgkktgy5aQ0YkUR0mypIK7rwPqiSbTTTeznc2sp5mdZmbfBx4HNprZ17LPZ4Az\ngN8CTwCt2W19zKyHmU0wsyO6EcoTwHLg+2a2S3ai3TEF9ivlnCuA0bkH2b7NJ5jZTsDbwCagvdHs\n14lGpAE+Dny9yPOJiJRl4t78+fCXv0RJclunnQZ77gn/8z+JhyFSMiXJkhru/mPgK8ClwEqiWtvP\nAbe5+2aiVnDvA94ArgA+7u4LsqURZxDV+y7JvvYXQL/8w7c9XXuPs8c6ExibPf/LwIcL7NfVc+b7\nHvAfZrYq2/quN9HI8etEo+Z7AN9o53W/Bw42s08D33R31SOLSPEaGxM/xXe+A1/6EvTtu+M2M7js\nsihX37hxx+0iaRJ8MREzu44o2Vjh7ge3s30qcDuwOPvULe7+7TKGKCIiaDGRqpDwYiLNzXDssbBw\nIfTvX3i/D38YDjsMvtHeUIBIjEpZTCQNSfJxwHrghg6S5Ivc/f1lD05ERN6hJLkKJJwk//M/w5gx\n8J//2fF+uWT61VehV6/EwhEpKUkO3t3C3R8xs1Gd7NatNyciIiLlsXYt3HorvFREEdi4cTB6NDz0\nEJx0UvKxiXRHpdQkH21ms8zsbjMbHzoYERERebd774UpU2DAgOL2P/vsqIeySFpVQpL8NDDS3Q8l\nmqylVuQiIiLdUV+f2KHvuKNr/Y/POisaed66NbGQREoSvCYZIFtucWd7Ncnt7LuEaOneVe1sC/9m\nRES6qLv1cuWmmmQpZPNmGDIE5s6FYcOKf92kSdGS1ZlMYqFJjSulJjktI8lGgbpjMxuad38yUWK/\nQ4Kc4+4Vf6uvrw8eg95Hdb4XvY/03USqwcMPw9ixXUuQISq5+OMfk4lJpFTBk2Qz+1/gUWCcmb1k\nZp8yswvM7F+yu5xtZnPM7FngJ8C5wYIVERGRHdx+e/eWms4lydvaWzpJJLA0dLf4p062XwlcWaZw\nREREpAvco3rkO+/s+mvHjYPBg+Gxx6KWcCJpEnwkWXaUqZLirGp5H1A970XvQ0Ti9vzzUFcHEyd2\n7/W5CXwiaZOKiXtx0aQSEak0pUwqKTddY6tAQ0N0i9G3vgVvvgk/+Un3Xv/II/DlL8NTT8UalghQ\n4SvuxUkXcBGpNEqSpawSWHHvhBPga1+D007r3uvfegsGDYLly6Fv31hDE6mK7hYiIiJSYTZvjkaA\njz66+8fo3RsOPxwefTS+uETioCRZREREuuW552DUqOJX2Svk+ONh5sx4YhKJi5JkERER6ZZHH4Vj\njin9OFOmKEmW9FGSLCIiIt0SV5J89NHw9NNRfbJIWihJFhERqRX19bEeLq4kuW9fOPBAePLJ0o8l\nEhclySIiIrUixvZvy5bBhg3RctRxmDIlWt5aJC2UJIuIiEiXPfZYNIpsMTUw1OQ9SRslySIiItJl\ncZVa5Bx3XHTMrVvjO6ZIKZQki4iISJfFnSQPHgzDh8Ps2fEdU6QUSpJFRESkSzZtguefhyOOiPe4\nRx2lyXuSHkqSRUREakVME/eefhomTIBddonlcO+YNAmeeSbeY4p0l5JkERGRWtHYGMthnn4ajjwy\nlkO9y2GHwbPPxn9cke5QkiwiIiJdMns2HHJI/Mc95BCYMwc2b47/2CJdpSRZREREuiSpJHm33WDU\nKJg3L/5ji3SVkmQREREp2pYtMH8+TJyYzPFVlyxpoSRZREREivbii7D33rDrrskcX3XJkhZKkkVE\nRGpFfX3Jh5g9Gw4+OIZYCtBIsqRFz9ABiIhIOphZC7AW2AZsdvfJYSOS2MXQAi6peuScww6DWbNg\n2zao01CeBKS/fiIikrMNyLj7YUqQpZCkk+SBA2GPPWDBguTOIVIMJckiIpJj6P8F6UTSSTKoLlnS\nQRdDERHJceB+M3vSzD4TOhhJn5Ur4R//gBEjkj2P6pIlDZQki4hIzrHuPgl4H3ChmR0XOiBJl9wo\nslmy51GSLGmgiXsiIgKAuy/P/nzdzG4FJgOP5O/TkDfxK5PJkMlkyhihlKyhoaTJe+UotYDtSbJ7\n8gm5VJempiaamppiOZa5eywHSgMz82p6PyJS/cwMdw+eBpjZLkCdu683s12B+4BGd78vbx9dYyud\nWZR5dtPHPw4nnAAzZsQYUwFDhkRdLoYNS/5cUr1Kucaq3EJERACGAo+Y2bPA34E78xNkESjfSDLA\n+PEwd255ziXSHpVbiIgI7r4EODR0HJJeb78dtWWbMKE855swIUqSTz65POcTaUsjySIiItKphQuj\nrhZ9+pTnfLkkWSQUJckiIiLSqfnz4cADy3c+JckSmpJkERGRWlFf3+2Xzp8PBxwQYyydmDABXnih\npHmGIiVRkiwiIlIrSmj/Nm9eeUeSBw+G3r3h1VfLd06RfEqSRUREpFPlHkkGlVxIWEqSRUREpEPu\nSpKl9ihJFhERkQ4tWwa77QYDBpT3vOqVLCEpSRYREZEOzZtX/lFk0EiyhKUkWUREpFZ0c+Jeudu/\n5ajDhYSkJFlERKRWNDZ262Uh6pEBBg2CnXeGV14p/7lFlCSLiIhIh0KVW4BKLiScnqEDEBHpjqam\n6Ja7n8lE9zOZ7fdFJB6hyi1ge8nFqaeGOb/ULvMqKvQxM6+m9yMinWtthX79YN066Ns3dDRdZ2a4\nu4WOoxi6xlYBsy4X+K5dC8OHR//WLMDf1GuugSeegOuuK/+5pfKVco0NXm5hZteZ2Qoze66DfX5q\nZgvMbJaZHVrO+EQkvVpbYcqU6P6UKdFjEYlXrh45RIIM0bnnzw9zbqltwZNk4HrglEIbzew0YIy7\njwUuAK4uV2Aikm4L72nmuOd/DkRfx6puUaQT9fVdfknIemSA/feHF18Md36pXcGTZHd/BFjdwS7T\ngRuy+z4O9DezoeWITURSrLmZQy6axqBhvYFo0YEJEwLHJJJ23WgBF7IeGWDoUNi8Gd58M1wMUpuC\nJ8lFGA68nPd4WfY5EalV27bBOedQ981GLn5hBgAzZ1ZmTbJI2s2fH43mhmIG48ZBc3O4GKQ2qbuF\niFSeujp46CEYMIBcXqwEWSQZzc1hk2SIzt/cDEcfHTYOqS2VkCQvA0bkPd47+1y7GvK+SspkMmTU\nC0qkKjXNGvBOC7ipU7d/i5z2FnBNTU005QIXSbmtW2HJEhgzJmwc48apLlnKLxUt4MxsH+BOdz+o\nnW3vAy5099PN7CjgJ+5+VIHjqD2RiFQUtYCTNHvpJTjmmPAr3t18M/z+9/CHP4SNQypPpbeA+1/g\nUWCcmb1kZp8yswvM7F8A3P0eYImZLQSuAT4XMFwRCUG93UTi0cWJewsWwH77JRNKV6jDhYSQipHk\nuGiUQ6QKNTfDSSfB/feHL4xMgEaSpay6uJjINdfAk0/CL3+ZYExF2LABBg+G9euhR4+wsUhlqeiR\nZBGRgpqbYdq0aPSrChNkkbRbuBDGjg0dBey6a5Qkv/xy5/uKxEVJsoikUy5BbmyEGTNCRyNSkxYu\nTEe5BWjynpSfkmQRSR8lyCKpkJaaZNjeBk6kXCqhBZyI1Jq334bvfAc+8YnQkYjUrG3bYPHidCXJ\nGkmWclKSLCLpM3FidBOReNXXF73rq6/CgAFRPXAajBsHd98dOgqpJSq3EBERqRVdaAGXplIL0Eiy\nlJ+SZBEREdlBWjpb5IwaBStXwsaNoSORWqEkWUTCam4O34RVRHaQps4WEPVHHj06ikukHJQki0g4\nuS4WdboUiaRN2sotQG3gpLz0P5OIhKE2byKplrZyC4iSZI0kS7koSRaR8lOCLBJGkRP33GHRIhgz\nJtlwumq//ZQkS/koSRaR8nKHj35UCbJICI2NRe326qvQt290SxMlyVJO6pMsIuVlBg88kL7/fUXk\nHWkstYBoZFtJspSLRpJFpPyUIIuk2sKF6Su1ANh7b1i1Sm3gpDyUJIuIiMi7LFqUvs4WEDXC2Xff\naLlskaQpSRaRZK1ZEzoCEemiJUuinsRppLpkKRclySKSnOZmOOggmDMndCQiAlBfX9RuixdHI7Zp\npCRZykVJsogkI7/N28SJoaMRESi6BdzixRpJFlGSLCLxUx9kkYrV2hpNjBsyJHQk7VOSLOWiJFlE\n4qUEWaSi5eqRzUJH0j4lyVIuSpJFJF5m8P3vK0EWqVBpLrUAGDkSXnsN3nordCRS7bSYiIjEa+zY\ndK5CICJFSfOkPYCePaNEeckSOOCA0NFINdNIsoiISK0oYuJe2keSQSUXUh5KkkVERGpFY2OnuyhJ\nFokoSRaR7mtuhssvDx2FiMQozQuJ5ChJlnJQkiwi3ZPrYtG3b+hIRCQm27ZBSwvss0/oSDqmJFnK\nQUmyiHSd2ryJVKXly2HAANhll9CRdGy//WDRotBRSLVTkiwiXaMEWaRqVUI9MkQj3S+/DJs3h45E\nqpmSZBEpnjucf74SZJFKVV/f4eZKSZJ32gn22gteeil0JFLN1CdZRIpnBvfdBzvvHDoSEemOTlrA\nVcKkvZzRo6OkfsyY0JFItdJIsoh0jRLkqmVmdWb2jJndEToWCSPtC4nkyyXJIklRkiwiIjlfAl4I\nHYSEUynlFhCNIGvyniRJSbKIFPbGG1EdslQ9M9sbeB/wy9CxSDiVlCRrJFmSpiRZRNrX3AyHHQbP\nPBM6EimP/wa+Cui3ohq1aROsXg3DhoWOpDgaSZakKUkWkR3lt3k7/PDQ0UjCzOx0YIW7zwIse5Nq\n1MHEvZYWGDkS6iokM8iNJOvLLkmKuluIyLupD3ItOhZ4v5m9D9gZ6GtmN7j7J9ru2JCXZGUyGTKZ\nTLlilDg0NhZMlFtaKmfSHsDAgVFCv2oVDBoUOhpJi6amJpqammI5lnkV/QpmZl5N70ek7BYsgBNP\nVIJcRmaGu6dm5NbMpgIXufv729mma2ylMys49HrVVfDcc3D11WWOqQRHHBHFPXly6EgkrUq5xlbI\nlyoiUhZ9+sB//ZcSZJEa1NISrWRXSTR5T5KkJFlEthsxAj7ykdBRSEDu/lB7o8hS/ZYsqaxyC9Dk\nPUmWkmQRERHRSLJIG0qSRUREakV9fcFNlTZxD6KRZCXJkpTgSbKZnWpm882s2cy+3s72qWa2JrtU\n6jNmdmmIOEWqzoIF8N3vho5CRMqpQGeL9ethwwbYY4/yhlOq0aNVbiHJCZokm1kdcAVwCjABOM/M\nDmhn14fdfVL29u2yBilSjXJdLIYODR2JiKRArtTCUtNnpTh77w0rV8Jbb4WORKpR6JHkycACd1/q\n7puB3wHT29mvwv7ZiqRYLkFuaIDzzw8djYikQCVO2gPo2TOab9zSEjoSqUahk+ThwMt5j1/JPtfW\n0WY2y8zuNrPx5QlNpAopQRaRdlTipL0cTd6TpFTCintPAyPdfaOZnQbcBowLHJNIZfrc55Qgi8gO\nKjlJ1uQ9SUroJHkZMDLv8d7Z597h7uvz7t9rZleZ2e7uvqq9A2rJVJEO3HUX9O4dOoqaFueSqSJd\n1tDQ7uS9JUvg6KPLHk0sNHlPkhJ0WWoz6wG8CEwDlgNPAOe5+7y8fYa6+4rs/cnA/7n7PgWOpyVT\nRaSipG1Z6o7oGlsFCixLPWkSXHtttMxzpbnlFrjhBrjtttCRSBqVco0NOpLs7lvN7PPAfUT10de5\n+zwzuyDa7NcCZ5vZvwKbgU3AueEiFhERqT5LllR2uYVGkiUJQUeS46ZRDpE8r70WtXirtJ5ONUYj\nyVJW7Ywkr1kTdYhYt64yLxfr1sGwYdDaWpnxS7JKucaG7m4hIkloboYjj4S//z10JCKSckuXVmaP\n5Jx+/aBPH3j99dCRSLVRkixSbZqbYdo0aGys3Jk4IlI2lVxqkaM2cJIEJcki1SQ/QZ4xI3Q0IpI2\n9fU7PNXSUpkLieTbd98o2ReJk5JkkWqxYIESZBHpWDvt3yq5R3KORpIlCUqSRapF//7w4x8rQRaR\nLqnUJanzKUmWJChJFqkWQ4bAOeeEjkJEKkxLC4waFTqK0ihJliQoSRYREalR7tVTk6wkWeKmJFlE\nRKRGrVkT/RwwIGwcpRoxImoN//bboSORaqIkWaQSNTfDpZeGjkJEKk2biXu5UotK7ZGc06sXDB8O\nL70UOhKpJkqSRSpNrs3b6NGhIxGRStPY+K6HuYVEqoHqkiVunSbJZvYDM+tnZr3M7K9m9rqZfawc\nwYlIG+qDLCIxqob2bzmqS5a4FTOS/F53XwecAbQA+wFfTTIoEWmHEmQRiVk1JcmjR2tBEYlXMUly\nz+zP04Hfu/vaBOMRkUIuukgJsojEqhrav+Wo3ELi1rPzXbjLzOYDm4B/NbM9gH8kG5aI7OCWW6LZ\nKSIiMVFNskhhnY4ku/slwDHAEe6+GdgITE86MBFpQwmyiJSqvv5dD6up3EI1yRI3c/eOdzDbBfgK\nMNLd/8XMxgL7u/td5QiwK8zMO3s/IiJpYma4e0U04NI1trqsXQt77w3r1lV+CziIFkbp3z8aHR84\nMHQ0khalXGOLqUm+HnibaDQZYBnw7e6cTESKtGwZbN0aOgoRqWJLl1ZHj+QcM03ek3gVkySPcfcf\nAJsB3H0jUCX/pERSqLkZjjoKHnkkdCQiUsWqqdQiRyUXEqdikuS3zWxnwAHMbAzwVqJRidSq/DZv\nU6eGjkZEqlg1JsmjR8OiRaGjkGpRTHeLeuBPwAgzuwk4FvhkkkGJ1CT1QRaRMqrWJPn550NHIdWi\nmO4W9wMfIkqMf0vU5aIp2bBEasyiRUqQRSR5DQ3v3M3VJFcT1SRLnIrpbnF8e8+7+8OJRFQCzbyW\nirVmDTz4IE0DP0hTU/RUUxNkMtH9TGb7faku6m4hZWUWtYEADj8crr4ajjwycEwxmj8fzjwTFiwI\nHYmkRSnX2GKS5DvzHvYBJgNPu/uJ3TlhknQBl2rR2gr9+kWtmfr2DR2NJElJspRVXpI8aFCUVO6x\nR+CYYvSPf0Rt4DZuhB49QkcjaZBoCzh3PzPvdjIwEVjdnZOJSOdaW2HKlOj+lCnRYxGROLW2Rgnl\n4MGhI4lXnz5R0v/KK6EjkWpQTHeLtl4BDow7EBGJzJkDc+dG9194Yft9EZG4VFuP5HyqS5a4dJok\nm9nPzOyn2dsVwEzgmeRDE6lSzc3wb//2zleebU2cCBMmRPfHj99+X0QkLtXY2SJn9Gj1SpZ4FNMC\n7qm8+1uA37r73xKKR6QqNTVFt0FvNnPuL6bx2CmNPNto7U7I69sXZs6MapJnzlRNshTPzOqAs939\n/0LHIilVXw9Ud5KsBUUkLp1O3KskmlQiqZbtgzzjlUZ+5e23ecsl07n76m5R/eKeuGdmT7n7EXEd\nr82xdY2tEhdfHNXufv3roSOJ3403wr33wv/+b+hIJA1KucYWHEk2s+fJrrLXdhPg7n5wd04oUpOy\nCfKmbzRy/YUzuLy1/RFiJcMSg7+Y2cXAzcCG3JPuvipcSJI2LS3V1fotn2qSJS4dlVucUbYoRKrd\npZey6RuNHH1tNII8ZYpKKSQx52Z/Xpj3nAOjA8QiKbV0afWWW6gmWeKicguRBOXKJ2zbVm67swfP\nPQfbtkHPnlGSfNRRoSOU0NQnWULYY4+ok87QoaEjiZ877LorrFwJu+0WOhoJLenFRI4CfkbU9m0n\noAewwd37deeESdIFXNLMDA45BGbPjn5qJFkgkZrkXsC/ArnVUpuAa9x9cwzH1jW2CqxfD0OGwIYN\n1dkCDqLOQDffDAcdFDoSCS3RxUSAK4DzgAXAzsCngSu7czKRWpVbEOSee6KfSpAlQT8HDgeuyt4O\nzz4nAg0NVd0jOUd1yRKHYlrA4e4LzayHu28FrjezZ4FvJBuaSIV66SXYay/o1QuIEuTDDos2TZoE\nxx4LP/pR9FgT9SQBR7r7IXmPHzCz2cGikXRpbKTlyIaqrUfOUV2yxKGYJHmjme0EzDKzHwDL6d5K\nfSLVL9vFguuvh5NOAqK6v6VLo82rVsEPf6haZEnUVjMb4+6LAMxsNLA1cEySItXcIzlHvZIlDsUk\nux/P7vd5onZCI4CzkgxKpCLlEuTGxncSZNAKelJ2XwUeNLMmM3sIeAC4KHBMkiItLVG5RTXTSLLE\noZiR5MOBu919HdCYcDwilSk/QZ7x7oVCtIKelJO7/9XMxgL7Z5960d3fChmTpMvSpXD44aGjSJZq\nkiUOxSTb6D7YAAAgAElEQVTJZwL/bWYPEzWn/5O7b0k2LJEKsnhxwQQ5fwW9qVNViyzJMbMT3f0B\nM/tQm037ZWd33xIkMEmdWim3WLIkagdXzRMUJVlF9UnOthQ6jahJ/XHA/e7+6YRj6zK1J5Ig1q+H\nBx+EM88MHYlUoLhawJlZo7vXm9n17Wx29wJroXftHLrGVrqGBoZc1cBzz8Gee4YOJllDhsCsWTBs\nWOhIJKRE+yTnnaQXcCrwKeB4dx/cnRMmSRdwEak0cfZJNrM64Gx3/784jtfO8XWNrXAbNsDgwbBx\nY/WPsB51VPTt3bHHho5EQkq0T7KZnWZm/0PUJ/ks4JdAlf/+KSJSedx9G/C10HFIetVCj+ScMWNg\n0aLQUUglK6Ym+RNEtcgXaPKH1Lr8GuOmpu11xaoxlhT5i5ldTHTd3pB70t1XhQtJ0qIW6pFz1OFC\nStVpkuzu5yUZgJmdCvyEaFT7One/rJ19fkpUE70B+KS7z0oyJpFCMsOaySz/EVx9NdZo7yTMIily\nbvbnhXnPOTC60AvMrDfwMLBT9na7u/+/xCKUYGqh/VvO6NHRdBGR7gq6KEi2fu4K4BRgAnCemR3Q\nZp/TgDHuPha4ALi67IGKwPY2b+95D63ro+8qc8tNi6SFu+/bzq1ggpx9zVvACe5+GHAwcKKZqZKz\nCi1dWjsjySq3kFKFXjlvMrDA3Ze6+2bgd8D0NvtMB24AcPfHgf5mNrS8YUqtaWqChobolsnAz77Q\nzNojpzH/vEZaz5nBlCnRflOmKFGWdDGzXczsUjO7Nvt4rJmd0dnr3H1j9m5vov8bVicYpgTScvec\nmkmSVW4hpQqdJA8HXs57/Er2uY72WdbOPiKxymSiBPmii+DVh5q58NZp9P/vRg74wQzmzIG5c6P9\nXnhh+32RlLgeeBs4Jvt4GfDtzl5kZnVm9izwGtDk7i8kF6KE0jJ3Q80kycOGwZo1UUcPke4oWJNs\nZs8T1bG1y90PTiSiEjU0NLxzP5PJkNFsKumm1tZopPjf+Q6N1sjF58ygL9uXmZ49W8tMS9c1NTXR\nlGwx+xh3P9fMzoNohNis814G2c4Yh5lZP+A+M5vq7g+13U/X2MrWwj41kyTX1UWlJUuWRNdtqQ1x\nXmML9kk2s1xpf27yx43Znx8FcPdLSj652VFAg7ufmn18SXTo7ZP3zOxq4EF3vzn7eD4w1d1XtHM8\n9fCU2Dz2GBx/PGzdso2evep4+OGo7yZECXS/frBunZaZltLE2Sc5e7xHgWnA39x9kpmNAX7r7pO7\ncIz/ADa6+4/aPK9rbAXbuBEG7bqJDVt3pi7098hlcsYZ8JnPwPS2hZxSM0q5xhYcSXb3pdmDn5yd\nzJFziZk9A5ScJANPEi2ZOgpYDnwEaNtN4w6iRP3mbFK9pr0EWSRu20eM694ZMdYy01IBGoA/ASPM\n7CbgWKJFoAoys8HAZndfa2Y7AycDjUkHKuW1dCmM5CXq6vYPHUrZqC5ZSlFMn2Qzs2Pd/W/ZB8cQ\nUy2zu281s88D97G9Bdw8M7sg2uzXuvs9ZvY+M1tI1AKuw4u9SFz69oWZM6MR45kzo8dKhiXt3P0+\nM3saOAow4Evu/kYnL9sL+HW2LKMOuNHd/5pwqFJmLS0wiqVAbSXJ6nAh3VVMknw+8Csz6599vAaY\nEVcA7v4n2vyLdfdr2jz+fFznE+lQSwvsuSf06QNsL6VQSYVUCjP7q7tPA+5u57l2ufvzwKRyxCfh\nLFkC+07aPXQYZTVmDNx3X+gopFIVs5jI08AhuSTZ3dcmHpVICLk+yNdeS9POp72rrCI3V0kjyZJW\nZtYH2AUYbGYDiUaRAfqhjkBCNAaw7zlHhA6jrFRuIaXoNEnO9iT+LjDM3U8zs/HA0e5+XeLRiZRL\nLkFubITTTiODkmGpOBcAXwaGAU+zPUleR7Rok9S4JUtgUo19X7DvvtEvB1u3Qo8eoaORSlNMbfH/\nAH8muvACNBNdiEWqQ36CPCO2SiKRsnL3y919X+Bidx+dt9reIe6uJFmicot9Q0dRXrvsArvvDq++\nGjoSqUTFJMmD3f3/gG0A7r4F2JpoVCLl0tKiBFmqzWtm1hcgu/LeLWZWY+OH0p6WltpLkkGT96T7\nikmSN5jZILILi2TbsKkuWarDnnvCL36hBFmqyX+4e6uZHQecBFwH/DxwTBJYayts2gR77BE6kvIb\nM0Z1ydI9xSTJXyHqVTzGzP4G3AB8IdGoRMqlTx849dTQUYjEKfdN3+nAte5+N7BTwHgkBZYsiVaf\ns8aG0KGUnSbvSXd1mCSbWR3QB5gKHEM0MWSCuz9XhthERKTrlpnZNcC5wD1m1puYettL5WppiZJk\nGmtvjZgxY1RuId3T4YXT3bcBV7r7Fnef6+5z3H1zmWITEZGu+zDRZOtT3H0NsDvw1bAhSWi1OGkv\nRyPJ0l3FjC781czOyq7EJFK5mpvh4x+HbdtCRyKSGHffCCwCTsmuaDrE3bWcQo2r1Ul7EI0kL1wY\nOgqpRMUkyRcAvwfeMrN1ZtZqZusSjkskXrk2byecAHX65lmql5l9CbgJGJK9/cbMNI+kxuVqkmvR\nkCHw9tuwenXoSKTSdLiYSHb0eIK7v1SmeETipz7IUlvOB97j7hsAzOwy4DHgZ0GjkqBqudzCbHtd\n8hG1teCglKizmmQH7i5TLCLxU4Istcd4dy/7rWxffU9qkHteuUV9fehwgthvP03ek67rdFlq4Bkz\nO9Ldn0w8GpG4/fd/K0GWWnM98LiZ3Zp9/AGiXslSo3JlBgMGAA0NIUMJZr/9VJcsXVdMkvwe4KNm\nthTYQDQi4e5+cKKRicThqqui79pEaoS7/9jMmoDjsk99yt2fDRiSBJYrtajlS+GYMfDoo6GjkEpT\nTJJ8SuJRiCSllv9XkJpiZn2AzwL7Ac8DV7n7lrBRSRq80yO5hu23H9x4Y+gopNJ0Os3f3ZcCI4AT\ns/c3FvM6EREpq18DRxAlyKcBPwwbjqRFLU/ay1G5hXRHpyPJZlZPdOHdn6jWrRfwG+DYZEMT6aLF\ni2GPPaBv39CRiIQw3t0PAjCz64AnAscjKdHSAuPGhY4irOHDo9rsDRtg111DRyOVopgR4Q8C7yeq\nR8bdXwWUhUi6NDfD1Knw4IOhIxEJ5Z3VUFVmIfne1SO5Rifu1dVFo+laeU+6opgk+e1sKzgHMDP9\nDibpkt/m7f3vDx2NSCiHZBd8WmdmrcDBWgBKoE25RWNj0FhCUsmFdFUxE/f+z8yuAQaY2WeAGcAv\nkg1LpEjqgywCgLv3CB2DpI87LF2qiXug5aml6womyWbW293fcvcfmtnJwDqiuuT/dPf7yxahSCEv\nv6wEWUSkAytWRDW4mqoRjSQ//3zoKKSSdDSS/BgwycxudPePA0qMJV323BOuvx5OOil0JCIiqbR4\nMYweHTqKdNhvP7j11s73E8npKEneycz+CTjGzD7UdqO735JcWCJF6NVLCbKISAeUJG+ncgvpqo6S\n5M8CHwUGAGe22eaAkmQREZEUW7SoTZJcXx8sltBGjYLly+Gtt6B379DRSCUomCS7+yPAI2b2lLtf\nV8aYREREJAaLF0fdMd9Roy3gIPryccSIqG/0/vuHjkYqQTEr7l1nZseY2T+Z2Sdyt3IEJ/KO5mY4\n+2zYujV0JCIiFUPlFu+mkgvpimJW3LsRGAPMAnIZigM3JBiXyHb5bd56qMuViEixlCS/m3olS1cU\n0yf5CKLlTj3pYER2oD7IIiLdsmkTvPFGtCSzRJQkS1cUs+LeHGDPpAMR2YESZBGRbmtpiSar6Qu4\n7caOhQULQkchlaKYJHkw8IKZ/dnM7sjdkg5MhGuuUYIsItJN7ZZa1PDEPYBx46LxF5FiWGdVFGY2\ntb3n3f2hRCIqgZmpKqSauINZ6ChEEmVmuHtF/EXXNbay/OxnMG8eXHVV3pNm0bW1Rm3eHK0+uHat\n2sDVilKusZ3WJKcxGZYaoQRZRKTbNGlvR7k2cIsXw4EHho5G0q5guYWZtZrZunZurWa2rpxBioiI\nSNcoSW7fuHGqS5bidLSYSN9yBiI1buFC2H336CYiIiVTkty+sWNVlyzFKWbinkiympvhhBOgqSl0\nJCIiVcE9SpL33Td0JOmjyXtSLCXJElZ+m7cPfSh0NCIiVWHlSth5Z+jfv82G+vog8aSJyi2kWEqS\nJRz1QRYRSUTBUosabwEHKreQ4ilJljBefVUJsohIQhYvhjFjQkeRTiNGwKpVsH596Egk7ZQkSxh7\n7gk33aQEWUQkAZq0V1hdnZanluIoSZYw6urg+ONDRyEiUpUWLdKkvY6o5EKKoSRZRESkyixaFI2W\nSvs0eU+KESxJNrOBZnafmb1oZn82s7ZzcHP7tZjZbDN71syeKHecIiIilWbBgmi0dAeauAeoDZwU\nJ+RI8iXAX9x9f+AB4BsF9tsGZNz9MHefXLboJD7NzXD66bB5c+hIRESqXmsrrFsHe+3VzsbGxrLH\nk0Zjx2okWToXMkmeDvw6e//XwAcK7GeoLKRy5dq8nXUW9OoVOhoRkaq3aFHU2aJO/3MWpJFkKUbI\nf0JD3H0FgLu/BgwpsJ8D95vZk2b2mbJFJ6VTH2QRkbJbuFD1yJ0ZMgTefjtqBSdSSM8kD25m9wND\n858iSnovbWd3L3CYY919uZntQZQsz3P3R2IOVeKmBFlEJAglyZ0z2z6afNRRoaORtEo0SXb3kwtt\nM7MVZjbU3VeY2Z7AygLHWJ79+bqZ3QpMBgomyQ15kxIymQyZTKZ7wUtpbrpJCbJIO5qammhqagod\nhlSxhQthsmbwdGr//ZUkS8fMvdAAbsInNrsMWOXul5nZ14GB7n5Jm312Aercfb2Z7QrcBzS6+30F\njumh3o+ISHeYGe5uoeMohq6xlWHqVKivhxNPbGdjQ4M6XGR961uwcSN873uhI5EklXKNDVmTfBlw\nspm9CEwDvg9gZnuZ2V3ZfYYCj5jZs8DfgTsLJcgiIiLSSbmFEuR3HHggzJ8fOgpJs2AjyUnQKIeI\nVBqNJEucNmyAwYOjn+pu0bE5c+Ccc2DevNCRSJIqdSRZqsXChfDaa6GjEBGpeYsXR8tRK0Hu3H77\nwZIlauEvhemfkZRmwQI44QR4+OHQkYiI1LyFCwustCc76NMH9t47+sVCpD1KkqX7FiyIZoY0NMCH\nPxw6GhHpJjPb28weMLO5Zva8mX0xdEzSPWr/1jUHHKC6ZClMSbJ0T36CfP75oaMRkdJsAb7i7hOA\no4ELzeyAwDFJNyxY0EmSrIl773LAAapJlsKUJEvXrVihBFmkirj7a+4+K3t/PTAPGB42KumOTkeS\nGxvLFksl0EiydERJsnTdkCHw+98rQRapQma2D3Ao8HjYSKQ7VG7RNWoDJx1RCzgRkYDS1ALOzHYD\nmoBvufvt7WzXNTbFNm2CgQOj9m89ehTYyQz0Gb7jzTdhzBhYvTr6o5HqU8o1NtFlqUVEpDKYWU/g\nD8CN7SXIOQ15Na2ZTIZMJpN4bFKcJUtgn306SJBlB4MGQa9eURXhnnuGjkbi0NTURFNTUyzH0kiy\niEhAaRlJNrMbgDfc/Ssd7KNrbIrdfjtcey3cfXcHO2kkeQdTpsA3vxl1M5Xqo8VEJDnNzTBtGrz1\nVuhIRCQhZnYs8FHgRDN71syeMbNTQ8clXdPcDPvv38lO9fVliaWSaPKeFKJyCykslyA3NkLv3qGj\nEZGEuPvfAH1JX+Hmz4fJkzvZSS3gdqAkWQrRSLK0Lz9BnjEjdDQiItKJF18sYiRZdqAOF1KIkmTZ\nkRJkEZGKoyS5ezSSLIVo4p7s6Pvfj3ohK0EWSVxaJu4VQ9fY9Fq1KupssXatWpl11dat0Ldv1OGi\nb9/Q0Ujc1AJOuq2pKbrl7kfdnC4hMxoygWISEZGuyY0iK0Huuh49otHkF16A97wndDSSJkqSa1wm\nE91aW6Pqijvv1G/SIiKVpuhSi4YGTd5rx8SJMHeukmR5N9UkC62tUZ9IiH62toaNR0REuqboJLmx\nMfFYKtGECTBnTugoJG2UJNe65maa71/K3LnRwxde4J37IiJSGTRprzS5kWSRfCq3qEG5OuRBbzZz\n7i+mseCUHzBo0ChWrIDx46PfqEVEpHIoSS6NRpKlPepuUauam9l24jQ+vayRy9dFXSz69YN161ST\nLFJO6m4hpdq6FXbbDd58E3bZpZOdtSx1u7Ztg/794aWXYODA0NFInNTdQrommyA3WiPXM4N7xsKn\nPgVTp8KPfhTtkpvQJyIi6dbSEnXt7DRBloLq6qJvUufOheOOCx2NpIWS5Frz5pswbRqL/7mR7/4g\nGkFetQqmT4fvfS9wbCIi0mVdKrWor080lkqWq0tWkiw5SpJrzaBBcMcdDN3vMCbcDbNnqw5ZRKSS\nvfhi1Oe3KGr/VpDqkqUtdbeoRYcdRt++MHNm9HDmTNUhi4hUqvnzNWkvDupwIW0pSa5hucRYCbKI\nSOVSZ4t4aCRZ2lJ3i2rnvsM6pe0vRa3JeiIhqLuFlGqvveCJJ2DEiNCRVDZ32H336JeOIUNCRyNx\nUXcLaV9zc9S24v773zXtWcmwiEh1WLs2WiV1+PDQkVQ+s2g0ee5cJckSUblFtWpuhmnT4Pzz1RdI\nRKRKvfACHHhg1MKsKJq41yHVJUs+JcnVKJcgNzbCjBmhoxERkYTMndvF7kSNjYnFUg1Ulyz5lCRX\nGyXIIiI1o8tJsnTo4IPhuedCRyFpoSS52tx3nxJkEZEaMXdu1Ote4nHIIfD889Ey1SLqbiEiEpC6\nW0gphg+Hv/0N9tmnyBeYRW0cpKB9943Gm8aODR2JxKGUa6xGkkVERCrQmjWwbh2MHBk6kupy6KEw\na1boKCQNlCSLiIhUoLlzu9jZAqC+PrF4qsUhhyhJloiS5ErW3Bx1PRcRkZrTrUl7agHXKY0kS46S\n5EqV62Lx5JOhIxERkQDU2SIZSpIlR0lyJcolyA0N8LGPhY5GREQCUJKcjFGjYMMGeP310JFIaEqS\nK01+gnz++aGjERGRQJQkJ8MsGk2ePTt0JBKakuRKsnatEmQREWHVKli/HkaMCB1JdVLJhYCS5MrS\nvz/ce68SZBGRGpdbRMS62v1VE/eKog4XAkqSK8/EiaEjEBGRwF54oZulFo2NscdSjTSSLBAwSTaz\ns81sjpltNbNJHex3qpnNN7NmM/t6OWMUERFJI9UjJ2v8eFi0CP7xj9CRSEghR5KfBz4IPFRoBzOr\nA64ATgEmAOeZ2QHlCS8FtHSoiIi0Y84cJclJ6t0bxo2LfhmR2hUsSXb3F919AdBRRdVkYIG7L3X3\nzcDvgOllCTC05mY48khobQ0diYiIpIh71Hnh0ENDR1LdDj0UnnkmdBQSUs/QAXRiOPBy3uNXiBLn\nitfUFN1y9zOZ6H4mA5lh2TZvjY3Qt2+Q+EREJJ1eeQV69YI99wwdSXU74ohova7PfCZ0JBJKokmy\nmd0PDM1/CnDg3939ziTPnXaZTHRrbY1y4TvvzObDzXkJ8owZgaMUEZG0mTWrhFHk+vpYY6lmkyfD\nr34VOgoJKdEk2d1PLvEQy4CReY/3zj5XUENee5tMJkMmN0SbQq2tMGVKdH/KFHjkV83sNl0Jskg1\na2pqoin3NZJIN5SUJKsFXNEOOQRefBE2bYKddw4djYRgHnhymJk9CFzs7k+3s60H8CIwDVgOPAGc\n5+7zChzLQ7+frnjsMTj+eNiyJfrqbN5Xf8WYMShBFqkhZoa7d7XbbRCVdo2tVmedBWefDeedFzqS\n6nf44XDFFXD00aEjke4q5RobsgXcB8zsZeAo4C4zuzf7/F5mdheAu28FPg/cB8wFflcoQa5EEydu\nn508fjwMuWSGEmQREemQJu2Vz5FHRnXJUpuCjyTHqRJHOVpboV8/WLdOc/REapFGkqUr1q2DvfaK\nfvboETqa6nfdddHk+htvDB2JdFcp19i0d7eoWvndLaZOhR/9KLqfm9AnIiLS1nPPRd9CKkEujyOP\nhB/8IHQUEoqS5EAyw5rJfHBTNDNARESkCCVN2oNo4p4m7xVt/HhYtgzWrIEBA0JHI+UWcsW92pVr\n8zZ7duhIRESkgpScJDc2xhZLLejZEw47DJ7eobWA1AIlyeWW3wf5E58IHY2IiFSQkpNk6TJN3qtd\nSpLLSQuFiIhIN23ZAvPmwUEHhY6ktihJrl1Kkstl/Xo46SQlyCIi0i0vvgh77w277RY6ktpy5JHw\nxBOho5AQNHGvXHbbDf76Vxg7NnQkIiJSgVRqEcaYMbBxI7z6KgwbFjoaKSeNJJeTEmQREemmp56C\nSZNKPEh9fSyx1BIzOOYYePTR0JFIuSlJFhERqQCPPw7veU+JB1H7t2457jh45JHQUUi5KUlOyrZt\noSMQEZEqsXlz1DX08MNDR1Kbjj1WSXItUpKchObmqHBs9erQkYiISBV4/nnYd1/o2zd0JLXpiCOi\nziLr14eORMpJSXLccm3evvxlGDgwdDQiIlIFYim1kG7r0ydaVOTxx0NHIuWkJDlO6oMsIiIJeOIJ\nmDw5dBS1TXXJtUdJclyUIItIBTOz68xshZk9FzoW2VFsSbIm7nWb6pJrj7l76BhiY2Ye7P3cfDNs\n2KAEWUS6xMxwd0tBHMcB64Eb3P3gAvuEu8bWsHXrov68q1dDr14lHswM9Bl2y5tvRnXhq1ZBT60y\nUTFKucbqY47LueeGjkBEpNvc/REzGxU6DtnRU09Fc8FLTpClJIMGwYgR8NxzMfSrloqgcgsREZEU\nUz1yeqjkorYoSRYREUmxxx9XkpwWxx0HM2eGjkLKReUW3dHcHBUlHXVU6EhERMqqIW/iVyaTIZPJ\nBIulVjzxBPz4x6GjEICpU+Hii6P1wuo0zJhKTU1NNDU1xXIsTdzrqlwXi299Cz75yWTPJSJVLy0T\n9wDMbB/gTnc/qMB2Tdwrs1deifrzrlwZzbkrWUODOlyUaOxY+MMf4JBDQkcixSjlGqvfg7oiv82b\nEmQRqSJm9r/Ao8A4M3vJzD4VOiaJvto/7riYEmRQghyDk0+Gv/wldBRSDkqSi6U+yCJSxdz9n9x9\nmLv3dveR7n596JgEHnoo+opf0uOkk5Qk1wqVWxRj0yYYPx7+4z+UIItIrNJUbtEZlVuU3/jx8Jvf\nqOVYmqxeDaNGweuvQ+/eoaORzqhPctJ23hkefBD22Sd0JCIiUiNWroRXX1Xta9oMHAgHHAB//7tG\n+audyi2KpQRZRETKaObMqC9vjx6hI5G2VHJRG5Qki4iIpNBDD8Hxx8d8UE3ci4WS5NqgmuT2bNmi\nhdlFpCxUkyyFHHooXH11zC35zUCfYcn+8Q/YY4+oRV///qGjkY6oBVycmpvhoIOiinwREZEAVq+G\nRYvg8MNDRyLt6dMHjj4aYlqzQlJKSXK+XJu3r341+hVRREQkgEceiUaQe/UKHYkUcsopcM89oaOQ\nJNV0TUFT0/bfAhfe08yVL06jKdNI/9EzyASMS0REatvDDydQjyyxev/7o+4WP/+5lqiuVqpJBtY/\n08zqw6cx+MpGdv6c+iCLSPmoJlnac+SR8MMfJtBiTDXJsTrwQLjhhujzknRSTXIJWlvh38+aTz2N\nHH3tDFpbQ0ckIiK17PXXo+q/o49O4OD19QkctHZNnw533BE6CklKzY8kP/ZY9JXWli1R7dfDD8c8\nk1hEpAMaSZa2broJfv97uO220JFIZx59FD77WXjuudCRSCEaSS7BxIkwYUJ0f/z47fdFRERCuPde\nOO200FFIMd7zHnjtNViyJHQkkoSaH0mGqOSiXz9Ytw769k0gMBGRAjSSLPm2bYOhQ+Hpp2HkyNDR\nSDFmzIh6Wn/xi6EjkfaUco1Vd4um6P7UqfCjH0X3M5noJiIiUk5PPQVDhihBriTTp8NPf6okuRpp\nJFlEJCCNJEu+xsbo280f/jB0JFKsDRtgr71g6VIYODB0NNKWapJFRESqQOL1yA0NCR68Nu26K5x0\nEtx6a+hIJG4aSRYRCUgjyZLz5pswejSsXAm9eyd0EvVJTsQf/whXXgkPPBA6EmlLI8kiIiIV7r77\novkwiSXIkpjTT4dZs2DZstCRSJyUJIuIiKTA7bfDGWeEjkK6o08f+NCH4Le/DR2JxEnlFiIiAanc\nQiCa/DV8OCxcCIMHJ3gilVsk5sEH4StfgWefDR2J5KvIcgszO9vM5pjZVjOb1MF+LWY228yeNbMn\nyhmjiIhIOdxzT7QwRaIJsiRq6lR44w2YOzd0JBKXkOUWzwMfBB7qZL9tQMbdD3P3ycmHFV5Trnlz\nhauW9wHV8170PkTS6eab4dxzy3Ci+voynKQ21dXBeedFy4pLdQiWJLv7i+6+AOhsCNyosdrpakkA\nquV9QPW8F70PkfRpbYX774cPfKAMJ1MLuER97GPwm9/A1q2hI5E4VELy6cD9ZvakmX0mdDAiIiJx\nuusuOO442H330JFIqQ4+OKotv/vu0JFIHBJdltrM7geG5j9FlPT+u7vfWeRhjnX35Wa2B1GyPM/d\nH4k7VhERkRBuvhk+/OHQUUhcLrwQrrgC3v/+0JFIqYJ3tzCzB4GL3P2ZIvatB1rd/ccFtmvKrohU\nHHW3qF3r1sGIEdGSxgMGhI5G4vDWWzBqFDz0EOy/f+hopJTuFomOJHdBu8Gb2S5AnbuvN7NdgfcC\njYUOUin/0YiIiEA0inziiUqQq0nv3nD++XDVVXD55aGjkVKEbAH3ATN7GTgKuMvM7s0+v5eZ3ZXd\nbSjwiJk9C/wduNPd7wsTsYiISHzc4ec/h89+town1cS9srjggmgC3/r1oSORUgQvtxARkcqgcot4\nPfEEfOQj0QIideUastJiImXzwQ/CKaeU+Zcg2UFFLiZSqmpajKQL7+VUM5tvZs1m9vVyxlgMMxto\nZveZ2Ytm9mcz619gv1R+JsX8+ZrZT81sgZnNMrNDyx1jsTp7L2Y21czWmNkz2dulIeLsjJldZ2Yr\nzHYxiXgAABGKSURBVOy5DvZJ/WfS2fuolM9D4nX11dGIY9kSZCmriy6C//ov2Lw5dCTSXZX8T7Oa\nFiPp9L2YWR1wBXAKMAE4z8wOKE94RbsE+Iu77w88AHyjwH6p+0yK+fM1s9OAMe4+FrgAuLrsgRah\nC39XHnb3Sdnbt8saZPGuJ3of7aqUz4RO3kdWJXweEpPVq+GWW+BTnwodiSTluONgn32isgupTBWb\nJFfTYiRFvpfJwAJ3X+rum4HfAdPLEmDxpgO/zt7/NVCoNX4aP5Ni/nynAzcAuPvjQH8zG0r6FPt3\nJfUTXbPtHld3sEtFfCZFvA+ogM9D4nPjjXDaaTBkSOhIJEn19fCd78CWLaEjke5IW6KShGpZjGQ4\n8HLe41eyz6XJEHdfAeDurwGFLv9p/EyK+fNtu8+ydvZJg2L/rhydLVG428zGlye02FXKZ1KMavg8\npAjbtgWYsCdBHH981OJPo8mVKS0t4NpVTYuRxPReguvgfbRXQ1lodkgqPpMa9zQw0t03ZksWbgPG\nBY6plunzqCG33AK77holUGVXXx/gpLWtvh4+/eloyeqeqc66pK1Uf1zufnIMx1ie/fm6md1K9FV0\n2ROyGN7LMmBk3uO9s8+VVUfvIzsxaai7rzCzPYGVBY6Ris+kjWL+fJcBIzrZJw06fS/uvj7v/r1m\ndpWZ7e7uq8oUY1wq5TPpUBV9HtKJbdvgm9+E7343ajRRdmoBV3aZDIwcCb/4Bfzrv4aORrqiWsot\nCi5GYma7Ze/nFiOZU87AuqHQZfNJYD8zG2VmOwEfAe4oX1hFuQP4ZPb+PwO3t90hxZ9JMX++dwCf\nADCzo4A1ufKSlOn0veTX7ZrZZKJ2kGlNyIzC/y4q5TOBDt5HhX0eUoLbboOddoLTTw8diZTTT34S\njSi/8UboSKQrKjZJrqbFSIp5L+7+/9u78yipyjuN499HBVd2EUfWCVECBJBOcEWj4kzco0EHcUeT\n6CiBEzJx3CbGI57IuEQdjRGPRkLiuOEQJio6aohiooC00CoaBXcDniAobqy/+eO93RYl3V0gXbea\nfj7n3NO37r116+mqXn711nvfdy0wGngEeAG4KyIW5JW5HhOAf5L0MjAMuBKax2tS3/Mr6WxJP8iO\neRB4TdKrwC3AubkFbkAp3wtQO+xgNXAdMCKnuA2SdCfwZ2APSW9KGtUcX5PGvg+ayethX05tK/JP\nf5pTK7LlZuBAGDkSLroo7yS2MTyZiJmZlUSeTORLmTo1FcnPPusiuSVavhz69oVp02DIkLzTtBxq\niZOJmJmZNRcrV8IFF6Qi2QVyy9S+PUyYAOee6yHhmgsXyWZmZk3sqqugTx846qicg/jCvVydeip0\n6ACXX553EiuFu1uYmVlJ3N1i0yxcCHvvnbpZ9OyZcxgJ/BrmavFiGDwY7r47p2EAWxh3tzAzM6tA\nETB6NJx/fgUUyFYRdt0VbrstjZv8vsewqWguks3MzJrI5Mnw1lvwox/lncQqyRFHwPHHp+4X7p9c\nudzdwszMSuLuFhunpgYOOQQeeywNAVYR3N2iYqxeDUcfnSYaueUWX9DZVNzdwloMSR0lVUuaK+lv\nkt4uuJ3bDJKShmWzB5qZ8cEHMHw4XHttBRXIVlFatYJ770191cePzzuNbUhFT0ttViybhWwwgKSf\nAh9FxLXFxymfJi83z5gZ69bBqFFw6KHp4/SKcumleSewAm3awAMPwH77Qbt2MGZM3omskFuSrTmr\n+/hEUm9JL0j6raTnge6SlhXsHyHp1mx9F0lTJM2S9HQ2DfD6J5ZmS9q94PaTkgZK2lvSnyU9m23r\nvYH7Xi5pTMHtBZJ2y9ZPk/RM1vJ9Y7Zta0m/kTRP0nxJozfT82NmZbZuHZx9NixdCr/4Rd5pNsBD\nwFWcXXeFxx+Hm26CSy5xb5hK4pZk25L0AU6JiGpJW/PFlt3a2zcAEyJilqSewB+AAUXH3kWaGni8\npK5Ah4iYL6kNMDQi1kn6NnAFcGIjuQJAUn/gOGDf7P63SDoRWATsHBGDsuPabsL3bmY5i4DzzoMF\nC+Chh2DbbfNOZM1Fr14wc2a6oG/JklQwt26ddypzS7JtSRZGRHUJxx0K/EpSNTAVaCep+N/ZvcDx\n2fqI7DZAB+B+STXA1UC/jch3KPBNYE722AcCvYFXgT0kXSfpnyPiw404p5lVgFWrUgvyvHnw4IPp\nY3SzjdG5c2pRXrwYDjgAFi3KO5G5SLYtyccF6+tY/+d7u6Jjh0TE4GzpERErC3dGxJvAR5L6kork\nu7NdVwDTI2IAcOwGzguwpuixt8++Crg9Iqqyx+0bEVdk/awHAk8C50q6peTv2Mxy9+67cPDBqbiZ\nPh3a+rMg20Rt2sC0aXDSSWkCmsmT3f0iTy6SbUtS10c5u2jv/ayv8lakbg61HgV+WHcnaVA957sb\nuBBoHREvZdvaAu9k66Pqud/rwDeyc+8FdC943H+R1Cnb11FSd0k7A1tFxBTgUrILE82s8k2fDkOG\npI/Jp051gWxfngRjx8LDD6fRUQ48EObOzTtVy+Qi2bYkxe+3LwAeAWYCbxVsHw3sn10o9zzwvXrO\ndx8wks9bkQH+E7ha0pwNPF6te4FdJc3Pzr0QICKeBy4DHpU0D3gY2IVURD+RdcG4nVSYm1kFe/11\nOO64NJvepElw8cWwVXP4j+oL95qNqiqYMwdOPz29CTv5ZHjuubxTtSyeTMTMzJB0GHAdqfHktoiY\nsIFjWvxkIn/9axq14p57YNw4+PGPYbsNdbqqVJ5MpFn68MM04cj110PfvnDWWXDMMbDDDnknq3ye\nTMTMzDZZ1iXpRuDbQH9gpKSv5ZuqdDNmzGjS8y9fnvqGHnlkuqBql13gxRdT63FDBXJT59oUM/IO\nUI9KfK6gcnK1bQs/+Um6mO+MM+Caa2bQtSuccgr87ndpRIxKUCnP1+biItnMzPYCXomINyJiNWkI\nxO/knKlkm/sf8/LlaZSByy6DYcOgZ0+YMgVGjkzdLC67DLp0KX+uzWFG3gHqUYnPFVRertatU7eL\nI4+cwUsvwdCh6WezTx/o3x/OPBNuvhmeeiqN1V1ulfZ8fVkeJ9nMzLqyfr/9t0mF8xZlzRpYsSJN\nGb1sGbz3XlrefhveeANeey21EC9fDgMGpFbjcePShVMe0s0qTZcucM45aVm9GmpqYNYsmD0b7rgj\njde97bZpDOaePaF793SfLl2gUyfo0CEtbdrATjulpXXr1CPHEhfJZma2WY0ZkwrO+tTXJbZ2e+H+\niPW31y7r1n3+ddEieOSRtL56dSqGV61Ky8qV8Omn8PHHaV+bNmn63/btU7HQuTN06wYDB6Y+nv36\nQY8ezeQiPLNMq1bpQr+qqlQ0Q/r9WLw4vQF88830ZnDJklQ8L12a3gwuW5beOK5YkX5H1q6F7bdP\n3Yhat05Lq1afL1tv/cVFSr8vUvq9nzkzrRcusP567W1Ird/Dh5f3+SqVL9wzM2vhJO0D/CwiDstu\nX0AaSXFC0XH+h2Fmzc6mXrjnItnMrIXLpnF/GRgG/A2YBYyMiAW5BjMzy5G7W5iZtXARsVbSaNK4\n4rVDwLlANrMWzS3JZmZmZmZFfGmCmZnVkXSbpCXZjJEb2n9SNlvlPEkzJQ2okFzHZJmqJc2RdEgl\n5Co4boik1ZK+Wwm5JH1L0nJJc7PlkrwzZccclL2Gz0v6Y1NnKiWXpH/LMs2VVCNpjaT2FZCrk6SH\nJD2X5TqjqTOVmKu9pPuz38enJfUrQ6Zukh6X9EL2XIyp57gbJL2SPWd7NnZeF8lmZlbo16RJReqz\nCDgwIgYB44Fby5Kq8VyPRsSgiBgMjAImlidWo7lqJ2u5kjQVfbk0mgt4IiKqsmV83pkktQNuAo6K\niK8DJ5QhU6O5IuLqiBgcEVXAhcCMiFiedy5gNPBcROwJHAxcI6kc3Wgby3URUJ39jTgduKEMmdYA\n4yKiP7AvcF7xhEiSDgd6R8TuwNnArxo7qYtkMzOrExEzgWUN7H86Ij7Ibj5NGmO5EnJ9UnBzJ+Dv\nTR6KxnNlfgjcB7zX9ImSEnOVdUTcEjKdBEyJiHey4yvpNaw1EvjvJoxTp4Rci4HaEbzbAEsjYk0F\n5OoHPJ4d+zLQS1LnJs60OCKey9Y/Ahbwxb9N3wF+kx3zDNBOUoPTArlINjOzTfU94KG8Q9SSdKyk\nBcCDwAY/bi03SbsBx0bEzZS5KC3BvtnHzg+U4yPxEuwBdJT0R0mzJZ2ad6BCkrYHDgOm5J0lcyvQ\nX9K7wDxgbM55as0DvgsgaS+gB9CtXA8uqRewJ/BM0a7iSZPeoZE3+R7dwszMNpqkg0ndGobmnaVW\nREwFpkoaCkwG+uQcCeA64N8LbldKofws0CMiPsk+hp5KKlLztA1QBRwC7Aj8RdJfIuLVfGPVORqY\nWaauFqW4EJgXEQdL6g38n6SBWUtqnq4Erpc0F6gBqoG15XhgSTuRPrUZuzmeBxfJZma2USQNJPX5\nPSwiSv2YumwiYqakbSR1ioilOcf5JnCXJAE7A4dLWh0R0/IMVVhARMRDkn4pqWNEvJ9jrLeBv0fE\nZ8Bnkp4ABgGVUiSfSJm6WpRof+AKgIhYKOk14GvAnDxDRcQK4Mza21muRU39uFl/7PuAyRHx+w0c\n8g7QveB2t2xbvdzdwszMiol6Wjwl9SB93HxqRCwsa6qGc/UuWK8CKGOBXG+uiPhKtvwj6R/4uWUs\nkBt6vroUrO9FGhK2HAVyvZmA3wNDJW0taQdgb1Lf0nJoKFftRYXfImUsp4ZyLQAOhbrXcw/KUIw2\nlktSO0mtsvXvA38qU+v27cCLEXF9PfunAadlufYBlkfEkoZO6JZkMzOrI+lO4CCgk6Q3gUuB1qRp\nqicC/wF0BH6ZtY6ujoi9KiDXcEmnAauAj4ERTZ2pxFyFyjYxQQm5jpf0r8Bq4FPK8Hw1likiXpL0\nMDCf9PH8xIh4Me9c2WHHAg9HxKdNnWcjcv0c+LWkeaSC9fxyvNEpIVdfYJKkdcALwFllyLQ/cDJQ\nI6ma9Lt2EdCTz3++HpR0hKRXSX8jRjV6Xk8mYmZmZma2Pne3MDMzMzMr4iLZzMzMzKyIi2QzMzMz\nsyIuks3MzMzMirhINjMzMzMr4iLZzMzMzKyIi2QzMzMzsyIuks3MzMzMirhINjMzs3pJWitprqQa\nSXdL2m4znHPml9lvVg6ecc/MzMzqJenDiGibrf8WmBMR1xUdo3BBYVsYtySbmZlZqZ4Eviqpp6SX\nJE2SVAN0k3SypGeyVuebJQlA0mmS5kmqljQp27Yi+7qDpD9k++ZLOqFwf7Y+LmvFni9pbLatp6QX\nJU2U9Lyk6ZK2LfeTYVs2F8lmZmbWkNpidxvgcKAm2747cGNEDAB2BEYA+0VEFbAOOFlSP+Bi4KCI\nGAyMze5b2+p8GPBORAyOiIHA9ML9kr4BnA4MAfYFvi9pUHbMV4H/ioivAx8Aw78QXOqVdRGZI2ma\npCkupq1U2+QdwMzMzCra9pLmZutPArcBXYHXI2J2tn0YUAXMzlqQtwOWAO2BeyJiGUBELC86dw1w\ntaSfAw9ERHFf5P2B/4mIzwAk3Q8cAPwv8FpE1BbszwK9NpC9a0SMkPSDiJi4Cd+7tWAuks3MzKwh\nn2Stw3WynhQfF24CJkXExUXHjW7oxBHxiqQq4AhgvKRHI2J8iblWFqyvJRXmxed/KlvdrcRzmtVx\ndwszMzNriErY/hhwvKTOAJI6SOoBPA6cIKlj7fbC+0r6B+DTiLgTuIrUGl147ieBYyVtJ2lH4Lhs\nW0O51g8pfYX1C3qzkrgl2czMzBpS36gVddsjYoGkS4BHJG0FrALOi4hZkq4A/iRpDVANnFlw3wHA\nVZLWZfc5p/DcEVEt6Q5gdrZtYkTMk9SzgVzF9snub7ZRPAScmZmZmVkRd7cwMzMzMyviItnMzMzM\nrIiLZDMzMzOzIi6SzczMzMyKuEg2MzMzMyviItnMzMzMrIiLZDMzMzOzIi6SzczMzMyK/D93TkKm\nwpmftAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x110236450>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def plot(estimate_theta, cov_theta, shape_precision, scale_precision):\n",
" # Show the results\n",
" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 6))\n",
"\n",
" ax1.set_aspect(1)\n",
" ax1.set_title(r'Coefficients $\\theta$')\n",
" ax1.set_xlabel('True values')\n",
" ax1.set_ylabel('Inferred values')\n",
" ax1.errorbar(theta, estimate_theta, np.sqrt(np.diag(cov_theta)),\n",
" ls='none', marker='.')\n",
" tmin, tmax = np.min(theta), np.max(theta)\n",
" ax1.plot((tmin, tmax), (tmin, tmax), ls='--', c='r')\n",
"\n",
" ax2.set_title(r'Precision $\\tau$')\n",
" ax2.axvline(precision, ls='--', c='r')\n",
" ax2.set_ylabel('Posterior')\n",
" ax2.set_xlabel(r'Precision $\\tau$')\n",
" mean = shape_precision / scale_precision\n",
" std = np.sqrt(mean / scale_precision)\n",
" sigma = 5\n",
" linx = np.linspace(max(1e-5, mean - sigma * std), mean + sigma * std, 100)\n",
" pdf = np.exp(shape_precision * np.log(scale_precision) \n",
" - special.gammaln(shape_precision) + (shape_precision - 1) * \n",
" np.log(linx) - scale_precision * linx)\n",
" ax2.plot(linx, pdf)\n",
"\n",
" fig.tight_layout()\n",
" \n",
"plot(estimate_theta, cov_theta, shape_precision, scale_precision)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"One-sided p-value: 4.25623504633e-05\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAGsCAYAAAAv5DqZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xuc1XW1//HXwpFBERQBb4yAXBVQBLygSGw1L3jvl1qd\n8qRZekotj3axsgY6ZcbE0ZPZZSbTo0cr09LUREtnK6FpgRAXRe4iCnIRZgAZbuv3x2dvZzPMMJuZ\nvfd3X97Px2M/5rv39zP7u5jyO2s+e33Wx9wdERERERFp1CHqAERERERE8o2SZBERERGRJpQki4iI\niIg0oSRZRERERKQJJckiIiIiIk0oSRYRERERaUJJsoiIiIhIE0qSRUREJCfMbI6ZfaS9Y0RywbSZ\niIiISGkzs6XAIcB2YBMwBbjO3TdHGZdIlDSTLCIiIg6c7+5dgZHACcCtTQeZ2T65DkwkKkqSRURE\nBMAA3P1d4GlgGICZLTGzr5vZLGCjmXUws8PN7BEze8/MFpnZDR++iVmFmT2aOLfazH6Scm6JmZ2R\nOP6Gmb1tZnVm9rqZnd50TOL50WZWa2bvm9lsM7uwyfvdbGazEud/Y2Yds/xzkhKhJFlEREQ+ZGZH\nAucBM1Je/iQwHjiIMOv8BPAacDhwJvAVMzvLzDoATwJLgN5AL+C3zVxjEHAdMCoxe30OsLSZcWWJ\na00BegJfBh40s4Epwy4DzgaOAoYDV7btXy6yq7KoAxAREZG88JiZbQc2EBLdH6ac+x93fwfAzE4C\nerj7DxLnlprZrwiJdD0hcf66u+9MnH+pmWvtADoCw8xsrbu/1UJMo4HO7v6jxPNaM3sS+BTwvZTY\nViViewI4vumbmFlX4H+A7oRkegmwFfiMu29p8SciJU1JsoiIiABc7O61LZx7O+W4D9DLzNYlnhvh\nk+mpwJHAspQEuVnuvsjMbgQmAEPM7BngJndf2WToEcDyJq8tI8xQJ61KOd5MSNKbGgl8PvF9MXe/\nf0/xiYDKLURERCSwPZxLbYW1HFjs7gcnHt3c/UB3vyBxrnei7GKP3P237j6WkHQD/KiZYe8QEu9U\nvYEVrb1/k2vF3X0H8HHgH3vzvVK6lCSLiIjI3ngVqE8s5utkZvuY2VAzOyFx7l3gdjPb38zKzezU\npm9gZoPM7PTEIrutwAeEEoymXgE2J65VZmYx4ALgN22M/Wx3f72N3yslRkmyiIiI7GnThF3OJUop\nLiDU/i4B3gNqgK6JcxcCA4G3CDPLlzfzXuXA7cBqwmxxT+BbTa/n7tsS73cesAb4KXCFuy9II+5d\nmNkBhHIMkbRoMxERkSJjZuXAi4SFUR2Bx939W03GjAMeBxYnXvqDu38/p4GKiOQxLdwTESky7t5g\nZqe7++bE5g/TzGyMu09rMvRFd78oihhFRPKdyi1ERIpQynbC5YR7/fvNDNvTQi0RkZKmJFlEpAgl\ndkV7DVgJxN19XjPDTjGzmWb2lJkNyXGIIiJ5TTXJIiJFLLGJwrPAN9z9hZTXDwB2JkoyxhM2ZBgU\nVZwiIvlGSbKISJEzs+8Am9198h7GLCFsEbyumXP6RSEiBcvd21RapnILEZEiY2Y9zOzAxPF+wFnA\nzCZjDk05PokwabJbgpzk7gXxqKyszM57Z+l9sxZvicdaaPEWUqyFFm97KEkWESk+hwO1iZrkvwN/\ncvfnzOxaM7smMeZSM5uTGHMn8Imogi0IEydGHYGI5JhawImIFBl3nw2MbOb1X6Yc3w3cncu4REQK\niWaSRUSkaMRisahD2CuFFG8hxQqFFW8hxQqFF29baeGeiIjskZl5yf+uMINS/xmIFCAzw7VwT0RE\nREQkM5Qki4iItKayMuoIRCTHVG4hIiJ7pHILESlUKrcQEREREckgJckiIiIiIk0oSRYRERERaUJJ\nsoiIiIhIE0qSRUREWjNhQtQRiEiOqbuFiIjskbpboM1ERAqUuluIiIiIiGSQkmQpeGY2yMxeM7MN\nZnZ94vnMxPMbzGy2mX0kjfeZk844ERERKX4qt5CcMbN/A/4TOBqoA2YCt7n7tHa+76+ADe5+c3PP\nc83MlgBXu/vzUVxfJNNUboHKLUQKlMotJO+Z2U3AfwPfBw4BegN3Axdm4O37AHP38LygmNmpZvY5\nM7vZzA6KOh4REZFSpCRZss7MugITgS+5++Pu/oG773D3P7v7LYkxR5tZrZm9nyiPuDDl+w83s0fM\n7D0zW2RmN6Scew44HfipmdWlPL878XygmS0xszMS4yvM7NHEe602s5+kvFfquBavmTL2ZjObZWbr\nzew3ZlZuZvcT/gB4InH9rybGf8PM3k689rqZnd7Cz6o/cKW7/xp4C7i03f8DiEj7VVZGHYGI5JiS\nZMmFU4By4LHmTppZGfAEMAXoCXwZeDCR4Fri3GvA4cCZwFfM7CwAdz8TmApc7+5dU55fl3i+IOU6\nHYAngSWERLYX8Ntm4tnjNVNcBpwN9AWGA591938nJLcXJK7/YzMbBFwHjHL3rsA5wNIWflY/An6Z\nOB4E9GhhnIjkklrAiZQcJcmSC92BNe6+s4Xzo4HO7v4jd9/u7rWEZPZTwIlAD3f/QWL2eSnwq8S5\nvXUSIen9urtvcfet7v5SM+PSveb/uPsqd19PSKqPTzmXWv+0A+gIDDOzMnd/y92XNL2omR2eiHGU\nmV0LfIaQcIuIiEiOlUUdgJSEtUAPM+vQQqJ8BLC8yWvLCDO9fYBeZrYu8boR/rh7sQ1xHAks20Oy\nnpTuNVelHG8mJOC7cfdFZnYjMAEYYmbPADe7+7tNhp4BPOnu1Wa2PzAJeLaVWEVERCQLNJMsufAy\n0ABc0sL5dwgJbKrewApC8rzY3Q9OPLq5+4Hu3pYFf8uB3omyi9bGteeauy2Bd/ffuvtYQgIOcHsz\n33ck8Hri+GLgz+6+Js1rioiISAYpSZasc/c6oJKwmO5iM9vPzMrMbLyZ3Q68Amw2s68nXo8BFwC/\nAV4F6hPnOpnZPmY21MxOaEMorwLvAreb2f6JhXantjCuPddcBfRLPkn0bT7dzDoCW4EPgOZms1cT\nZqQBrgC+keb1REREJMOUJEtOuPt/AzcBtwLvEWptvwQ85u7bCK3gzgPWAD8FrnD3BYnSiAsI9b5L\nEt9bA3RNffuml2vueeK9LgQGJq6/HLi8hXF7e81UPwS+Y2brEq3vygkzx6sJs+Y9gW82832/B44z\ns88D33N31SOL5Ast3BMpOZFuJmJm5YQ6z46Jx+Pu/q3IAhIRkd1oMxG0mUhE3nkH7roLBg2Cq66K\nOhopRAW7mYi7NwCnu/sI4DjgDDMbE2VMIiIiEq26OvjSl2DYMNi4Eaqq4JproKEh6siklERebuHu\nyRrMckI870cYjoiIiETsO9+BVavgjTfCTPIrr8CaNRCLQX191NFJqYg8STazDmb2GrASiLv7vKhj\nEhERkWi88QY89BD84hdwyCHhtS5d4JFHoFcvuPvuaOOT0hFpTXKqxNbFzwLfcPcXoo5HREQC1SSj\nmuQcuuCCMGP81a/ufm7ePDj9dFi8GDp3znloUoDaU5OcN5uJuHudmT0FnADskiSbme5MIlLQ2nqT\nljxRWRl1BCXhL38JM8mPPtr8+SFDYNy4MMt88825jU1KT6TlFmbWw8wOTBzvB5wFzGxurLsXxaOy\nsjLyGPRvKc5/h/4t+fuQIqAWcFnnHhLfqiooL2953K23wo9/DB98kLvYpDRFXZN8OFCbqEn+O/An\nd38u4phEREQkx/75z9C94pKW9mZNOO44GD0aampyE5eUrkjLLdx9NjAyyhhEREQkeg8/DJddFsq/\nW/Ptb8Oll8INN6Q3XqQtop5JLjmxWCzqEDKmWP4txfLvAP1bRKQwucPvfx+S5HSMGgUdO8Jrr2U3\nLiltedPdYk+0slpECll7VlfnA92DJdv+8Q/49Kdh/vz0Z4a//nXo1Am+973sxiaFrWB33BMRESkI\nWriXVclZ5L0pnbjkEnjssezFJKKZZBGRLNNMchFQn+SscYd+/eCPf4Tjj0//+3bsCJuLTJsG/ftn\nLz4pbJpJFhERkYI0fTqUlcHw4Xv3ffvsAxddBI8/np24RJQki4iISGTaUmqRpJILySYlySIiIhKZ\nZ58NW1G3xRlnwKxZ8N57mY1JBJQki4iISETWr4cFC+CEE9r2/Z06wdlnw5NPZjYuEVCSLCIi0rrK\nyqgjKEovvQQnnRR6HrfV+efDM89kLiaRJCXJIiIirVELuKyYOhXGjm3fe4wdC3/7m5qPSOYpSRYR\nEZFIvPhi+5Pkfv1g505YtiwzMYkkKUkWERGRnPvgA5g5E0aPbt/7mMFpp4XZZJFMUpIsIiIiOffq\nqzBsGBxwQPvfS0myZIOSZBGRImNm5Wb2ipm9ZmZzzey2Fsb9xMwWmNlMM9uLvc5E2i8T9chJSpIl\nG5Qki4gUGXdvAE539xHAccAZZjYmdYyZjQf6u/tA4FrgF7mPtIBo4V7GZaIeOWn4cHjrLVi3LjPv\nJwJKkkVEipK7b04clhPu9e83GXIxcH9i7CvAgWZ2aO4iLDATJ0YdQVHZvh3+/vcwA5wJZWVw8smh\npZxIpihJFhEpQmbWwcxeA1YCcXef12RIL2B5yvMViddEsm7mTOjdG7p3z9x7JlvBiWRKWdQBiIhI\n5rn7TmCEmXUFnjWzce7+Qlvfb0JKuUEsFiMWi7U7RildL70EY8a0Pm5vnHYafPe7mX1PKTzxeJx4\nPJ6R9zIvgO7bZuaFEKeISHPMDHe3CK//HWCzu09Oee0XQK27/y7x/A1gnLuvaub7dQ82024VGXTV\nVaH127XXZu49N22CQw+FNWvCdtUi0L77r8otRESKjJn1MLMDE8f7AWcBM5sM+xPw74kxo4H1zSXI\nItnw2mswYkRm37NzZzjmGPjnPzP7vlK6lCSLiBSfw4HaRE3y34E/uftzZnatmV0D4O5/BpaY2ULg\nl8CXogu3AFRWRh1B0WhogPnz4dhjM//eJ54I06dn/n2lNKncQkQky6Iut2gv3YMlk2bMgH//d5gz\nJ/PvXVMD06bBffdl/r2lMKncQkRERApCNkotkkaODEm4SCYoSRYREZGcee01OD5L+zsOGwYLF8KW\nLdl5fyktSpJFREQkZ7I5k1xeDgMHZqeUQ0qPkmQRERHJiR074F//yt5MMqjkQjJHSbKIiEhrUjZT\nkbZbuDDssnfwwdm7xogRYbZapL2UJIuIiLRm4sSoIygKM2dmr9QiacQIzSRLZihJFhERkZzIZj1y\n0vHHh5rk7duzex0pfkqSRUREJCdykSR36QK9eoUNS0TaQ0myiIiIZJ17bpJkUMmFZIaSZBEREcm6\nd9+FnTvDLG+2jRypxXvSfkqSRUREWlNZGXUEBW/ePBg6FCwHG7Srw4VkgpJkERGR1qgFXLu9/joM\nGZKbayWTZPfcXE+Kk5JkERERybp583KXJPfsCZ06wdtv5+Z6UpyUJIuIiEjWzZsHxxyTu+sNHQpz\n5+buelJ8lCSLiIhI1uWy3AKUJEv7KUkWERGRrFq9GrZuhcMPz901lSRLeylJFhERaY0W7rXL66+H\nUotcdLZIUpIs7aUkWUREpDUTJ0YdQUHLdakFhCR53jx1uJC2U5IsIiIiWZXLzhZJ3bqFLarfeiu3\n15XioSRZREREsipZbpFrKrmQ9og0STazCjN73szmmtlsM/tylPGIiIhI5kUxkwxKkqV9op5J3g7c\n5O5DgVOA68zs6IhjEhERkQzZsAHWr4fevXN/bSXJ0h6RJsnuvtLdZyaONwKvA72ijElERGQ3lZVR\nR1Cw3ngDBg+GDhFkHEqSpT2inkn+kJn1BY4HXok2EhERkSbUAq7Noiq1gHDd11+HnTujub4UtrKo\nAwAwswOAR4CvJGaURaRExePh0dAAU6bAuedCeTnEYuEhIoUlivZvSQcdFB7LlsFRR0UTgxSuyJNk\nMysjJMgPuPvjLY2bkPJXfCwWI6bfliJFKRaDUaNg7FiYNSv0OJ06NbRyKhTxeJx4PB51GCJ5Yd48\n+Pzno7t+suRCSbLsLfOIu2yb2f3AGne/aQ9jPOo4RSR3Xn4ZPvIR2L4d9t0XXnwRRo+OOqq2MzPc\nPYd7jWWW7sHSHgMGwJNPwtERLcu/6SY49FD4xjeiub5Eqz3330hnks1sDPBpYLaZvQY48C13nxJl\nXCKSXcmSiuRxLAZjpk2i4qQjOOYHn2HYsDD7M2tW+Jh26NDoYhWRttu2DZYvh379ooth6FB44YXo\nri+FK/KZ5HRoFkOkONXXQ9eusOV7kyi/vyZkzL167XKurq6wSi2ao5nkIjBhghbvtcGbb8L48bBo\nUXQxvPQS3HgjvPpqdDFIdNpz/1WSLCKRqK8Pdcdnz5rE9R1rOPhfcQ4Y3KvZWWYo7IV7SpKLgFko\nkJe98tRTcNddYRFuVNauDTPZ69eH/xmltBRsuYWIlK45c+Dc2ZO4mho+sjPOb9/vxWgKOxkWkV0t\nWAADB0YbQ/fuYW3DqlVw2GHRxiKFRUmyiERi2OFr2a/T48Q2x+k5tJfqjkWK0IIFYSORqA0eDPPn\nK0mWvZM3m4mISOmIx2Hyfd353Q1/49ARvRg/HiZPbiyzEJHi8Oab0c8kQ2OSLLI3NJMsIjnXWFJh\n/PD2aGMRkezJh3ILUJIsbaOZZBERkdZUVkYdQcHZsgVWroS+faOOREmytI1mkkUk45rrUHH4uzMY\n/KmRWpQnhUnt3/ba4sXQpw+U5UGmoSRZ2kIt4EQka5K9jj/4ryo63V8DM2fC/vtHHVbOqQWclKLH\nH4eamrDbXtS2bm3su96xY9TRSC615/6rcgsRyYpkH+SvUsV7/1XNxj89X5IJskipypd6ZAiJ8ZFH\nRrupiRQeJckikhVz5sA5s6u4hmrG7axlzvqKqEMSkRzKpyQZVHIhe09JsohkxYiX7ua6smpOp5YD\nh1aoD7JIicm3JHnQICXJsnfyoJxeRApZc4v0li6FwzZ8lAM/fzGHvFzBOeeEPsjaTU8K1oQJWry3\nl/ItSR48GF55JeoopJBo4Z6IZERykV5dHXTpEnU0+UUL94qAGZT6z2AvbN4ctoPeuBH22SfqaIJ4\nHL79bZg2LepIJJe0cE9EIpVcpAfha319tPGISLQWLoSjjsqfBBlUkyx7T0myiLTbnDkwd244njev\n8ViiY2YVZva8mc01s9lm9uVmxowzs/VmNiPxuDWKWKX4LFgQaoDzyWGHhVZwa9dGHYkUCiXJItJu\nI56r4ruH/RKAIUPQIr38sB24yd2HAqcA15nZ0c2Me9HdRyYe389tiFKs8q0eGULFjGaTZW8oSRaR\ndln0H1VsuqOasovOZ8QIGD8+LNJLLuaTaLj7SnefmTjeCLwO9GpmaMHWSkv+WrQI+vePOordKUmW\nvaHuFiLSdlVV9H+uGmbV8s2KCr4ZdTzSLDPrCxwPNLe2/xQzmwmsAL7m7vNyGFrhqKyMOoKCsngx\nXH551FHsbuDAUC8tkg7NJItI21RVQXU11NZChTYKyVdmdgDwCPCVxIxyqulAb3c/Hvgp8Fiu4ysY\nav+2VxYvhn79oo5idwMHhlIQkXRoJllE9t6GDfD000qQ85yZlRES5Afc/fGm51OTZnd/2sx+ZmYH\nu/u6pmMnpCSJsViMmBpeSwu2bYN33oHevaOOZHeaSS5+8XiceIbq/dQnWUQky6Lqk2xm9wNr3P2m\nFs4f6u6rEscnAQ+7e99mxukeLGlbuBDOPjvMJueb998PyXtdXVjIJ8WvPfdfzSSLiBQhMxsDfBqY\nbWavAQ58C+gDuLtXA5ea2ReBbcAHwCeiileKR76WWgB06wbl5fDee3DooVFHI/lOSbKItCi55XRD\nA0yZAueeG37BaHvp/Ofu04A9buXg7ncDd+cmIikV+ZwkAwwYEOqSlSRLa5Qki0iLYjEYNQquH/kS\nMxeeijtMnaptp6UETZigxXtpytf2b0nJxXunnRZ1JJLv1N1CRPZo3S2TuHXhZ+lCnXbTk9I1cWLU\nERSMQphJ1uI9SYeSZBFp2aRJHDmlhuuOqaWertpNT0Rale9JstrASbpUbiEizZs0CWpq6PBCLY8e\nWEHXriq1EJE9cw9Jcr6XW2gmWdKhFnAisrtf/AImT+bl22p5Zm7ogxyPNy7W08K9vRNVC7hM0T2Y\n0C+s1H8GaVizBgYNgnW7ddrOH2oDV1rUAk5E2qSl7hXnHn0uo2sv4JSKCk65LOooRaRQ5HupBagN\nnKRPSbJIiWtogHvvhVWrQp3esGHw4IN9ufRStXsT+VBlZdQRFIRCSJJBbeAkPUqSRUpYLBYS4R//\nODzfuhXWrg31ek8/rRpkkQ+p/Vta8r39W1KyLllt4GRP1N1CpMT16QMdEneC7dsbF7TMng0PPBBd\nXCJSeAptJllkT5Qki5S4bT+YxPXb7gB2XZc0ZAhccUVEQYlIQSqUJFkdLiQdKrcQKWGLrp3Egb+v\n4bnucViz67k771SphYjsnUIpt9BMsqRDLeBESkyyo8XJL0xi6Es1PPiFOHVdenH77dCzJ6xeDZ07\nww03aOFepqgFnJSChgbo2hU2bYKyPJ+CUxu40qEWcCLSoubavH3k75M4emoNp+yIs+G+Xhx8cKhN\nrqgISfLDD8N550UduUgemTBBi/dasWwZHHlk/ifIoDZwkp4C+L+yiLRFanL81FOwcmVIgHds2MjN\nPWoZRZx36MW+W2H5crjllvBLo6wMXn01PDSLLJIwcaKS5FYsXgxHHRV1FOnr3z+UhyhJlpYoSRYp\nUrEYjBoFY8eGThVJb7x9ALNqnqbbzfDWLNhnH9i2rbHlm/IAEWmLpUsLL0leuBBOPTXqSCRfqbuF\nSBGbMwfmzt31td694a9/hTPOCM+3bAlf587dfayISLqWLCm8JHnRoqijkHymJFmkiA0bBkOHNh4D\n/OQn4etDD4WvyfrBPn0ax4qI7C0lyVJslCSLFLEuXeAXn4wDzgUXwLhxodb4vfdCfTLAzp0waBBc\ndhlMnhzqmEVE9tbSpdC3b9RRpG/AACXJsmeqSRYpUvE4bP1BFaNmVHP+qa9SXt7tw4V4o0bB9Okw\naxYce6y2nxZpVWVl1BHkPc0kS7GJvE+ymd0DXACscvfjWhijHp0iaUp2tTjpxSqGTavmwc/X0tCz\nYrdOFfX1oadpXZ0S5GxTn2Qpdhs3wiGHhB7JhdJ32B0OOADefTfcC6U4FXqf5HuBu4D7ow5EpNAk\nE+LkcTIJ/sTyKjpPq+bUrbX0mFaxy0xx6veMGxdKLEDt3kSk7ZKlFoWSIEOItV+/MJs8YkTU0Ug+\nijxJdve/mVmfqOMQKTSpye5f/wrTpsEpp8CV239F72erOWZHLSuo4L15oWvF6NFhrJJhEcm0Qiu1\nSEqWXChJluZo4Z5IgYrFQk/ja64JCTKEXse9rjmfHX+p5aBhFQAMGaKuFSKSXUuWFNaivSTVJcue\nKEkWKWD19aFkImnePJiz9nAOODqUWIAW5YlI9hXaRiJJ6nAhexJ5uUW6JqRsAxaLxYjp82IR5swJ\nMzhJffrAunWNu+ap5jga8XicuHrpFZcJE7Qd5R4sWQJjxkQdxd7r3x8eeSTqKCRfRd7dAsDM+gJP\nuPuxLZzXymqRZvz5z/C5q5xV7xmdOsEXvxhWaSshzi/qblEEzEI7BGnW8cfDPfeE9pKFZNGisPvo\nsmVRRyLZ0p77b+RJspk9BMSA7sAqoNLd720yRjdokSaSfZBpaOCWjbdyzjlQXq4EOR8pSS4CSpL3\n6KCDYPFiOPjgqCPZO9u2hTZwdXXh/inFp6BbwLn7v0Udg0ghiv2jCpZWQ20tZ1dEHY2IlKr33w9/\nP3TrFnUke2/ffaGiItRUDx4cdTSSb7RwT6RAxONw5ZVhlvgHB1Wx+JvVfOygWq68tUJbSYtIZJKd\nLQqpR3IqLd6TlkQ+kywiLbvzTnjsMdi+PWwj3bMnXLOhisvrqvkItWxfXsGvJ6q8QkSiU6idLZLU\nBk5aoplkkTx2443wxBOwYQNs2QK++QOGb36J0wkbhaxeHX5BiUiWVVZGHUHeKtQeyUn9+8PChVFH\nIflISbJInnvlldDqDeDttftx0fY/soJQhNy3L1xxRXSxiZQMtX9rUaHutpekmWRpiZJkkQJ2223a\nKEREoqVyCylWSpJF8tzJJ8OwYeG4c+ddz912W1jMp4V7IhKVQi+36NcvJPo7d0YdieQbJckieeyP\n//EMF12wkwMPhE6dQi/SssRy25494bzzGjteiIjkmntIMAs5Se7cOdxb33kn6kgk36i7hUi+qqri\nY89V87GXX4YePYjHw4xxQwM88wwfbh4iIhKVNWvCH/Bdu0YdSfv06xdKLirUc15SKEkWiVAy8U0e\nJ2eEr1hZRf/nwkYh9OgB7LqT3g9/mNMwRWTCBC3ea0ahzyIn9e8fdgwcNy7qSCSfRL4tdTq0JaoU\nu/r6MBNTVwddflEF1YkEWdMaRUHbUhcBbUvdrN//Hn77W3j00agjaZ+JE8MW1d//ftSRSKa15/7b\nak2ymU0ys65mtq+ZPWdmq83sM225mIjsrr4exo4Nx5OG/i87f6EEWUQKQ7HMJCfLLURSpbNw72x3\nrwMuAJYCA4CvZTMokVIyZw7MnRuOq9+9kNfuiCtBFpGCUCxJcrLcQiRVOjXJyTHnA7939w1WqBu0\ni0SsuRrkhgbo0yfMYhw+9GAGnR5dfCIie2PpUjj33KijaD/1SpbmpJMkP2lmbwAfAF80s57AluyG\nJVKckgvvnnkGXngBXn8dBg6ErVvD6+PHw/TpaukmIoVh2bLwR36hO+QQ2LIFNmyAAw+MOhrJF62W\nW7j7LcCpwAnuvg3YDFyc7cBEilE8HhLke+4Jz997z5kxA5YvhzFjQku31NlmEckTlZVRR5B3kj2S\niyFJNgt1ySq5kFTpLNzbH/gS8PPES0cAJ2QzKJFiFYvBRRfBunXwNSbxA77NBx+Ecxs3ws03hy5T\nmkmW9jKzCjN73szmmtlsM/tyC+N+YmYLzGymmR2f6zgLhtq/7WbtWujYsXhmXrV4T5pKZ+HevcBW\nwmwywAqvNk6iAAAgAElEQVRATVJE2iAeDy2TvrpzEl+ghru57sNz8+Y1LuATyYDtwE3uPhQ4BbjO\nzI5OHWBm44H+7j4QuBb4Re7DlEJVLIv2klSXLE2lkyT3d/dJwDYAd98MaOWeSBvEYvDF+klc7TXE\niPMOvT78JTNkCAwdGmV0UkzcfaW7z0wcbwReB3o1GXYxcH9izCvAgWZ2aE4DlYJVjEmyyi0kVTpJ\n8lYz2w9wADPrDzRkNSqRIrXo2kns/2BjggyhzKJbt5AgT56semTJPDPrCxwPvNLkVC9gecrzFeye\nSIs0q9iSZJVbSFPpdLeoBKYAR5rZg8AY4MpsBiVSDJIL8BYsgOefh/5HbuXGmTP5aZ8469/pBZuh\nZ0+4+mo45xzVIUt2mNkBwCPAVxIzyiIZsXQpDBoUdRSZo5lkaarVJNnd/2JmM4DRhDKLr7j7mqxH\nJlLgYjEYNSrsprdyJaxf35HLtj3E8M4hce7VK8xadOkSdaRSrMysjJAgP+DujzczZAVwZMrzisRr\nu5mQsnAtFosRK7W/6iZM0OK9JpYuhbPPjjqKzOnTB1asCNtT77tv1NFIW8XjceIZ+kjWvJW96M3s\nI8297u4vZiSCNJiZtxanSBSSs8UNDTBlSmiqX14OBx0E69eH1m6//vWu39OhA1x1FSxc2Dh7HItp\nJrmYmRnunvO1HGZ2P7DG3W9q4fx5wHXufr6ZjQbudPfRzYzTPdgs9DyTDw0bBg89BMcdF3UkmdO3\nL/z1rzBgQNSRSKa05/6bTpL8RMrTTsBJwHR3P6MtF2wL3aAln9XXh9niWbNg+HCYOrVxdvidd8KM\nMUCnTqFZfdMxUvyiSJLNbAzwIjCbsKbEgW8BfQB39+rEuJ8C5wKbgKvcfUYz76V7sJLkXbiHe9iK\nFcXTAg7gzDPhG98orhnyUtee+2865RYXNrnYkcCdbbmYSDGaM6exddvcufDJT4Zdm46c+QSPbBoP\nlLHPPmGW+bHH4LbblCBL9rn7NGCfNMZdn4NwpMgUW4/kJLWBk1TpdLdo6m3gmEwHIlKohg1rbN02\ndCjU1MDHF1fxvU030Y11QJh16d4dxo2DV18NpY3qYiEihWrZsuLqbJGkDheSqtWZZDO7i0T7N0JS\nfTyw28dxIqWqS5dQPtG1a/i67ptVXPBuNadTy2oOAUIP5Dvu0AyyiBSHYtmOuqn+/eGVpo0SpWSl\n0wLunynH24HfJD7GEyl5yYV7EGaJX/5/VYycXs2lh9WyYmXFh+PuvFMJsuw9M+sAXOruD0cdS8mr\nrIw6grxSbD2Sk9QGTlK1unAvH2jRiOSLlrpZxGIQW/FgqKOoreWRv1dw2WXhezp3hhtuSBkXiyp6\niUq7Vleb/dPdT8h0THsZg+7BsosbbggdIL7ylagjyawNG8Ji6/r6sFZTCl9WuluYWXJF9G6nCCuj\nc9b0RTdoySd//jN87nOwalXoWPHFL4ZSizNPrGOfDzZS/eQRLF0aZiO2bQszEwMGwJVXKkEuVe1M\nkm8H1gC/I3SgAMDd12UovHRi0D1YdnHhhWEjpEsuiTqSzOveHebNg0O1QXtRyFZ3iwvaGI9IUevW\nLazshtDS7YknYMYM6NKlK9CVUy+NNDwpPp9IfL0u5TUH+kUQiwgQyi2OOirqKLIj2eFCSbK0mCS7\n+7JcBiKSr+JxuO++sPnHokWwzz5QVgbbt4fzS5eG1m+jd9uCQaT93L1IUxEpVO7Fu3APGpPkU0+N\nOhKJWqst4MxstJn9w8w2mtlWM9thZnW5CE4kH8RicNddsHFj2F56xQqYNg2MnUBo+5ZsASeSaWa2\nr5l92cweSTyuNzNtmiuRWbcuTBQcdFDUkWSHFu9JUjp9kn8KfApYAOwHfB64O5tBieSb1A1DAB4/\nrYqHDgs7/Y4fD9OnRxSYlIKfA6OAnyUeoxKvSS5NmBB1BHmjWDtbJGlDEUlKazMRd18I7OPuO9z9\nXsIWpiIlY9iwxvq0r1LFZz6o5r7uX6VPH3j33V1bwYlk2Inu/ll3fz7xuAo4MeqgSs7EiVFHkDeU\nJEupSKdP8mYz6wjMNLNJwLu0bac+kYI1fTpcfjnsc0cV1xA2Crn87ApuuUgdKyTrdphZf3dfBGBm\n/YAdEcckJazYk2TtuidJ6STJVxCS4uuB/wSOBD6ezaBEonbnnfDYY2Fx3vTp0LMnXFtXxeWJBHkF\nFTz2mCaXJCe+BtSa2WJCC84+wFXRhiSlbNmykEgWqyOOCP2SN26EAw6IOhqJUjozwqMIfZHr3H2i\nu9+UKL8QKVo33hhau23YENq8le+znbE9XueMRIIM4RfFAw9EHKgUPXd/DhgIfBm4ARjs7rXRRiWl\nrNhnkjt0CO3ttHhP0kmSLwTeNLMHzOwCM0tn9lmk4L3ySliwB7BwaRkb7vg1Bw1r3Gp6yBC44oqI\ngpOiZ2ZnJL7+P+B8YEDicX7iNZFIFHuSDOpwIUGrSXJikcgA4PeELheLzOxX2Q5MJN8sWAAXJLbY\nOe64cDx5shbsSdaMS3y9sJmHNnvKtcrKqCPIC8keyaWQJKsuWdKaFXb3bWb2NGGXp/2ASwit4ESK\nQmp3ing8LMZraAgfuS1ZEr6uXg3l5TBuXONivVhMC/ckO9y90sw6AE+7+8NRx1Py1AIOgPffD+UI\nxdojOalfP3jjjaijkKils5nIeDO7j9An+ePAr4DDshyXSM4kE+SGBvjDH+CFF+Cg5x5l9dsNfCKx\nIfAnPhES5FgsjJ0wITyUIEs2uftO4OtRxyGSVAqzyKCZZAnSmUn+d+B3wLXu3pDleERyKpkg19XB\n3XfD1q2hD/L5f6vm+YYxLFx+GMOHhwRZJCJ/NbOvEu7Dm5Ivuvu66EKSUqUkWUqJuXu0AZidC9xJ\nmNW+x91/1MwYjzpOKV719TByJCxcGBLkZB/k//1rBWeeGXV0UgzMDHe3Nn7vkmZednfPWRMu3YMl\n6b//G956K7TJLGYNDdC1K2zaFLbglsLVnvtvpJuCJOrtfgqcAwwFPmVmR0cZk5SeOXNC3XFqgryC\nCm68MSTQIlFy96OaeRRxl1rJZ6Uyk1xeHnZZXb486kgkSlHvnHcSsMDdl7n7NuC3wMURxyQlJB6H\nP/0Jrih/eJcEGWDePPVBluiZ2f5mdquZVSeeDzQzdbfINS3cA0onSYZQcrFQu0KUtKiT5F5A6t9p\nbydeE8mZt96CP/mFfIQXWblPBfvvH14/6ij1QZa8cC+wFTg18XwF8P3owilR2l4TCElynz5RR5Eb\nqkuWFittzGw2oeVbs9z9uKxE1IIJKX/Fx2IxYmorIBkQi4WP1R5+eD+2sx8dHD7+8VCCcc45oQ+y\n2rzJ3orH48Qz10C7v7t/wsw+BeDum82sTfV1Iu3hHnYaLZWZ5AEDlCSXuj2Voyc/zrsu8TX5wfOn\nM3j9FUDvlOcVidd2M0EfdUkGtNQPuU+fcDM89tjQ5aJLl+hilMLX9A/5ie2bhdxqZvuRmLQws/6A\nOg1Jzq1fH74We4/kpAEDws6rUrpaTJLdfRmAmZ3l7iNSTt1iZjOAWzJw/X8AA8ysD/Au8EnCrn4i\nGRePw333hRqzBW/s4L21+7BlCxx9NPzkJ3D++TB1qhJkyTsTgCnAkWb2IDAGuCrSiKQkJeuRS+Vz\nDJVbSDqNTczMxrj7tMSTU8lQLbO77zCz64FnaWwB93om3lukqVgMRo2C6gGTOGDtIv6DXzJjBpx6\nKrz6athJb/LkxrEqsZB84O7Pmtl0YDRgwFfcfU3EYUkJKqV6ZGhMkt1L5w8D2VU6SfLVwK/N7MDE\n8/XA5zIVgLtPAQZn6v1EUjUtr7hu0yQueq+GGOHFbdtgyBD4vDZZlzxlZs+5+5nAU828JrlSWRl1\nBJFbsiQsaC4VXbtC587w7rtwxBFRRyNRaDVJdvfpwPBkkuzuG7IelUiGJGeE6+th88RJfKxvDWcN\njPPOgtBEpbwczjsv0hBFmmVmnYD9gR5m1o0wiwzQFXUByj2ti2HJkjC7WkqSi/eUJJemVssmzOxQ\nM7sH+K27bzCzIWZ2dQ5iE2mTeDz8PvvmN2HECLjpJrjziEl8gRrGeZwRFzTmF1/6ElRXN842i+SR\na4HpwNGJr8nH44RNmERyqpR6JCepV3JpS6fc4j5Cn85vJ56/CfwOuCdLMYm0W0MD3HMPrF4Nby3Z\nwX9tXEaMOCuX9+KDeGNd3YwZpXfTl8Lg7v8D/I+Z3eDud0Udj0iplVuA2sCVunSS5B7u/rCZfRPA\n3beb2Y4sxyXSZrEYbN8Ot98enq/bsA+T+93NO4th+LHwwgvqYCEFZaWZdXH3ejO7FRgJfN/dZ0Qd\nmJQO99KcSR4wAJ54IuooJCrpJMmbzKw7jT06RwOqS5a8kFyY19AAU6bAueeGOuOmSfCdd8JFF6nF\nmxSk77j7783sNOCjQBXwc+DkaMOSUrJmDXTsCAce2PrYYqJyi9KWTpJ8E/AnoL+ZTQN6ApdmNSqR\nNCXbuo0dC7NmhdmOqVPDufvuCzvn9ewJL72kFm9SsJKf3J0PVLv7U2ambalzbcKEkl68t3Rp6ZVa\nQJhJXrhQbeBKlbm3uPM0ZtaB0JvzVUKbNgPmu/u23IT3YRy+pziltL38MnzkI6HEYt994cnP/IZ/\nHHYRG70zt98Ot9wSZpeVGEtUzAx3b9OvWDN7krAT6VmEUosPgFfdfXgGQ2wtBt2DzUKmVKIefhh+\n9zt49NGoI8kt97DD4OLF0L171NFIW7Tn/rvH7hbuvhO42923u/tcd5+T6wRZpDXDhsHQoeF4Uo9J\njHnmu3Rq2EB5eZg9Li+PNj6RdroceAY4x93XAwcDX4s2JCk1pbhoD8LfRlq8V7rSKbd4zsw+DvxB\nUwmST1LrkXfuhK8xic9sqeFfv6zl5svU1FKKg7tvNrNFwDlmdg4w1d2fjTouKS1LlzZORpSaZF3y\nSSdFHYnkWjrbS18L/B5oMLM6M6s3s7osxyXSqlgs3LTuvRfOnT2Ja6yGmk/V8szcCvU9lqJhZl8B\nHgQOSTz+z8xuiDYqKTWlOpMMmkkuZXucSTYzA4a6+1s5ikekValbTT/+OIx57498gRo+2qGW315R\nwejRUUYnknFXAye7+yYAM/sR8DKg3smSM6W6cA9CkvzCC1FHIVForSbZgadyFItIWpIL8BoaYOZM\neHbf8zmNv7GlRwXr1kUdnUjGGY0dLkgca519rlVWRh1BZHbuhGXLGjdhKjVqA1e60qlJnmFmJ7r7\nP7IejUiaRo2CGxIfOG/c2pGNHMqKGXCESpGl+NwLvGJmf0w8vwTteJp7Jdz+bdUq6NoVOneOOpJo\nJNvASelJJ0k+Gfi0mS0DNhFmMNzdj8tqZCJ7MGcOzJu362vf/S5UVKjVmxQXd/9vM4sDpyVeusrd\nX4swJCkxS5aU3k57qQ4/HDZuhLq68MeClI50kuRzsh6FSDNSu1c88AB06BAe61dvY9iIfenUCTZv\nDmOHDYM77tBuelI8zKwT8B/AAGA28DN33x5tVFKKSnnRHoTfO8mSi5Ejo45GcqnV7hbuvgw4Ejgj\ncbw5ne8Taa9YDG6+GZ5+GlasgK1b4ceHVvGTzVfT0AAjRkC/fmHsBRfA9OmRhiuSaf8LnEBIkMcD\nP442HClVS5eW9kwywMCBsGBB1FFIrrU6k2xmlYQb9WBCbdy+wP8BY7Ibmkgoq5g7Nxx/dnUVw1dV\nczq1rPgnXHcdfPSjYba5vLyx44VKLaRIDHH3YwHM7B7CzqciObdkCZx4YtRRREtJcmlKp9ziY8AI\nYAaAu79jZvpQW3IiuZveWbOquNaq+Qi1rKACCKUVJbyWRorfh7ubuvv20JEzfYnE+gJgVXNrSMxs\nHPA4sDjx0h/c/fttD7fITZhQsjecpUvh8sujjiJaAwfCiy9GHYXkWjpJ8lZ3dzNzADMr0fWtkk2p\n9cePPAKHHgplZbB+Pdx2cBUDqeblH9TyzjcrwMNWoVddFXXUIlk1PGXjJgP2SzxPLp5ubQnRvYRe\nyvfvYcyL7n5R+0MtARMnlmySXOoL9yAkyfeop0zJSae2+GEz+yVwkJl9AfgrUJPdsKTUpNYfL1wY\nVhJffz3MmuUcsGkVV/Wp5ds/r+Cww8L4UaPgmmvgzjsjDVska9x9H3fvmnh0cfeylONW19i7+9+A\n91sZpn7Lskc7dsDbb5duj+QklVuUphZnks2s3N0b3P3HZnYWUEeoS/6uu/8lZxFKyfjOd2DWrHA8\naxZ85jMAxqfe/jE1NfBqoiIzHm+sOz7++NzHKVJETjGzmcAK4GvuPq+1b5DS8vbb0KNHWPdRyg47\nDLZsCZ9uHnRQ1NFIruyp3OJlYKSZPeDuVwBKjCWr/uu/4LnnwmI9CKUXACtXwh/+EJJmLcoTyZjp\nQG9332xm44HHgEERxyR5ZvHi0P6s1JmFTUUWLNAixlKypyS5o5n9G3Cqmf2/pifd/Q/ZC0uKVbL2\nOHkci4VFIUnLlzf/feeeqwRZJJPcfWPK8dNm9jMzO9jdm93cfUJKPW4sFiOm/yBLwpIlja02S12y\n5EJJcn6Lx+PEk4lGO5m7N3/C7DTg08DlwJ+anHZ3/1xGIkiDmXlLcUphqq8POxfV1YUuFfE4PPMM\nvH37AzzBhWyg8fOssjKYOhVGj44uXpH2MDPcPef1v2bWF3gi2UquyblD3X1V4vgk4GF379vC++ge\nXKLdLW69FTp2DDualrpvfzv8LCoro45E9kZ77r8tziQnFn38zcz+6e5a0ykZU18PY8eG4xEj4LLL\nQmlFp7uq+C7VPM8ZuyTJ3brB44+HejBNXomkx8weAmJAdzN7C6gEOhImOaqBS83si4RWcx8An4gq\n1oJQggkyhHKL886LOor8MHAg/PWvUUchudRqCzh3v8fMTgX6po539z21FRJp0f33w+zZ4Xjx4rCb\n3lGPVPGZ7WGjkHfotcv4W26Bm26KIFCRAubu/9bK+buBu3MUjhSoxYtVbpE0cCD8/OdRRyG51GK5\nxYcDzB4A+gMzgR2Jl93dv5zl2FJj0Ed9RSQ5k5zsZDHns1V0/N/ETnqJjUKSyspCE/uBA8MssmaS\npRBFVW6RKboHl65DDoF//YsP22+Wsvfeg6OPhnXNVu1LvspKuUWKEwjbo+oOKe2WXLh3xhkwfz6c\nseUpOj1QzUVda1lRV0FZWWjr9tZb4Yb01a/COecoORYRybX6eti0KWzuJNCzZ+gbvXYtdO8edTSS\nC+kkyXOAw4B3sxyLlIBksvvAA6EOeQrncsrOlzj34p7MeyDUHqf2Qy4vb+yGoURZRCR3liyBo44K\n7c8k/BySHS6UJJeGdJLkHsA8M3sVaEi+qK1MJR3xONx3X2jzlmz1ltze1Ax2+j6spidvvhl2dHr4\n4XA+FivZdTIiko9KsLuF6pF3l0yS1W2pNKSTJE/IdhBSvJJ1xMmWbwA/+xk8+mho65Y0aBDcfrtm\ni0UkT02cqCRZtD11iUmnu8ULuQhEilOy/3FNTXjekQauvLKcn/0slFUkd9d77TUYNSqyMEVEpInF\ni8MEhjQaNAiefDLqKCRXOrR0wszqzayumUe9mdXlMkgpXLEYnHlmWOjwNSbxf3yGtWthypTQIzlp\n3rxQpywiIvlBM8m7GzQoLDqX0rCnzUS65DIQKW5fYxJfoIYYcfr2hTvuCK//61+hFdyxx8IVV0Qa\nooiIpNCW1LsbPBjefBN27oQOLU4zSrFotU9yPlCPzsKRbPG2YAE8/zz06AEXz5/EZ7eFBHllh16c\neCIcdFDou3n44aEc45xzQicL9UKWYqQ+yUXADEroZ7BzJ3TuHD4F3H//qKPJL4cdBv/4Bxx5ZNSR\nSDqy3SdZJG2xWKgtHjsWVq6Er/kkLtxWQ5d/xnnnhF7UrYcuTT6j+OEPIwlVRCR9lZVRR5BT774L\n3bopQW7O4MGh5EJJcvFTkiwZN2cOzJ0L4GxctZEYcXac34sxY2Dy5DBGM8YiUlDU2UISkknyRz8a\ndSSSbUqSJWOSpRZ1dWE76e3bjUq+B8C+6+DHP1ZvSRGRQrB4cdhIRHaXTJKl+KnsXDImFoObbw61\nyFu2hNeOOSZ8HTIEhg6NLDQREdkLmklumZLk0qEkWTLq/vth9uzG5/vuGxY5DB0aSi2SW0yLiEj+\nUpLcMiXJpUPlFpI5v/41nz3vQmpqejJrFgwfHnbVa7pQT0RE8tuiRXDNNVFHkZ+OOiosTP/gA9hv\nv6ijkWyKbCbZzC41szlmtsPMRkYVh2TIpEnwwx9yQMetH243rQRZRIpGiS3cW7QIBgyIOor8VFYW\nEmVtT138IuuTbGaDgZ3AL4GvuvuMPYxVj848FY/Dth9MYuSMGj51WJytPXsBsH49XHJJGKNOFlLq\n1Ce5CJRQn+S6OjjiCKivD/9s2d0ll8CnPw2XXRZ1JNKaguyT7O7zAcz0n2AhSnayOPmFSQx9qYaf\nfjbOX2p68dRTcN55UUcnIiJttWgR9O+vBHlPVJdcGlSTLHslmRw3NMDa3/6FwctrOGVHnDX/G2aQ\nv/WtsJGIyixERArTwoUqtWjN4MFQWxt1FJJtWa1JNrO/mNm/Uh6zE18vzOZ1JTtS+yDfeSfULP0o\nJ+74O+/Qi61bw5jZs+GBB6KMUkRE2kNJcusGD4Y334w6Csm2rM4ku/tZmXqvCSmLJmKxGDEVuWZc\nMglOHid/xMma4uSW0yNHJvsgG+vovst7DBkCV1yRm3hF8lU8HieufodSoBYu1MZPrUmWW7irLKWY\nRbZw78MAzGoJC/em72GMFo3kUH09dO0aZoyblk28/HIop9ixIzzv3h3Wrm08f8UVobemFuuJNNLC\nvSIwYULJdLgYNy78U08/PepI8lv37jBvHhx6aNSRyJ4U5MI9M7sEuAvoATxpZjPdfXxU8UhQXx+S\nYIARIxpX7k6ZAheeuZlt++5P375hYQfAJz8J990HmzaFG8UnP6mFeyJShEokQQaVW6QrOZusJLl4\nRT6TnA7NYuRGPB7qie+7D3buhA4dYNgwWL4crn6/ijP2ncrXB/8JgDVrwm56/fpBr17w0EPNzzyL\niGaSpXBs2gQ9eoSvHbQn7x5dfTWcdBJce23UkcieFORMsuSfZM3xq6/CnDmhvnjyZPjLOVV8nmrO\n9lp+U9NYq5ZawzxuXBibfB+VWoiIFJ7kdtRKkFt3zDHw+utRRyHZpCRZPhSPwzPPwKpV4fmqVVD/\n3Sq+sLOaGLWsLaugd+/G8UqGRUSKi0ot0jdkCDz7bNRRSDbpb0X5UCwGF10E778fnn9ubRVnLKrm\nzA61rKCC7dvhrbciDVFERLJISXL6hgwJC/ekeClJll0MGwZDhwI43Q/awZ0X1dLQswIIK3kff7yx\nxEJEpGSUyMI9Jcnp690b1q8P63GkOClJll106QJTpwIYQx+4Be9VwVVXhU4XV10F5eVRRygiEoGJ\nE6OOICeUJKevQwc4+mjVJRcz1STLh+68Ex57DLZvh86d4bbboKwMLrkEZsyIOjoREck2Jcl755hj\nQsnFySdHHYlkg5Jk+dCNN4aWNmPHhvY/GzeGWWW1dRMRKX5btoQF20ceGXUkhUN1ycVN5RYS1NTA\n228zZw7MnRtemjev8VhERIrbkiXQp0/4BFHSM2SIyi2Kmf5TKGHJPscnvVjFsGnVPPj58dR1DTfJ\nRYvCf/xhEZ+IiBQ7lVrsPc0kFzclySUsFoPRf6vivWnVnLq1lh7TKpg6Fb71LejaVaUWIiIfqqyM\nOoKsW7gQ+vePOorCctRRsHJlKFHs3DnqaCTTlCSXkOTMcUMDTJkCP+hWxajp1YzdHvogr5oL//mf\nUFGhHfRERHZRAi3g5s+HY4+NOorCUlYWZt/nz4eRI6OORjJNSXIJSW47PXYsdJk1lWEdq/FptXT7\nfAVvzQqlFXfcodljEZFSNH8+fPzjUUdReJIlF0qSi48W7pWY5MK8v3EaI3f+k6XbKxJ9kVVeISJS\nyubPh8GDo46i8GjxXvHSTHIJSC2zeOop6NYNVq82yg85kMcfD6UXKq8QESlddXWwYUMot5O9M2QI\nPPhg1FFINihJLgGpZRazZ4etp1evhjfe0MyxiIjAggUwcGDYRU72jjpcFC/951Ai5r1Sz5w5ieN5\nYZvpyZPDWpR4PMrIREQKQJEv3FOpRdsNGABvvRU2Y5HiYu4edQytMjMvhDjzVlUV25/+Cyese5ZZ\ns2D4cNUfi+SSmeHuFnUcbaV7MGAGRfwzqKwM/7zvfS/qSArT0KGh5OL446OORJpqz/1XM8nFrqoK\nqqspu//XWqAnIiLNmj8fBg2KOorCdeyxoZxRiotqkotE0x7I554Lp/+zitNer2bWHbU886uwGkML\n9EREpKn58+Hmm6OOonAdd5yS5GKkJLkIJBPkujr4+c9DXdTHFlYxpms1+79SyykVFZxyWdRRikgu\nmdk9wAXAKnc/roUxPwHGA5uAK919Zg5DlDyxcye8+aZqktvj2GPhZz+LOgrJNCXJRSDZvWLkyMaF\nA2s378f8B2sZqX4+IqXqXuAu4P7mTprZeKC/uw80s5OBXwCjcxif5IkVK6Br1/CQtjn2WPjXv6KO\nQjJNNclFYs4cWLq08fmf+13PwNOVIIuUKnf/G/D+HoZcTCKBdvdXgAPN7NBcxFaQKiujjiBrNIvc\nfn36QH09rFsXdSSSSUqSi8SwYWF1bdILL2hxnojsUS9gecrzFYnXpDlF3AJO7d/az0yL94qRkuQi\nMX06jB8fjseMgepq9UAWEZHWqbNFZqjkovioJrkIvPmfP2fWjrMpP7g/48Y1dqxQ9woR2YMVwJEp\nzysSrzVrQspMaiwWI6abS9GYPx/OOivqKArfscfCrFlRRyHxeJx4hmYItZlIoZs0CWpqwpRxL31S\nKpKPotpMxMz6Ak+4+7HNnDsPuM7dzzez0cCd7t7swj3dg4vbUUfBs8+Gbaml7V58Eb7xDXj55agj\nkR1XY2IAAB6ASURBVFTtuf8qSS5kSpBFCkIUSbKZPQTEgO7AKqAS6Ai4u1cnxvwUOJfQAu4qd5/R\nwnvpHlykPvgAunWDjRuhTJ8tt8v770Pv3rBhA3RQMWveaM/9V/9JFColyCKyB+7+b2mMuT4XsRSF\nCROKcvHewoVhJlkJcvt16wYHHRQ6TfXrF3U0kgn6W6cQvfIK/OpXSpBFRHJl4sSoI8iKefNgyJCo\noyge6nBRXJQkF6KTT4YZM5Qgi4hIu8ydu2v7UGkfJcnFRUlyoTrggKgjEBGRAjd3rmaSM+m449QG\nrpgoSRYRESlR8+ZpJjmThg9XG7hiou4WhWD9+rAaQEQKUlQt4DKl5O/BELZUK7KfQUMDHHhg6MZQ\nXh51NMVh+/bwM125Urve5ov23H81k5zvqqrg4ouL7uYsIlJQKiujjiDj3nwT+vZVgpxJZWUwbJhm\nk4uFkuR8VlUV9pd+8MEwiyEiItEowvZvKrXIjpEjw9p6KXxKkvNVMkGurYWKiqijERGRIqPOFtkx\nYoSS5GKh9uERisfDo6EBpkyBc88NH3tdsbKK/s8pQRYRkeyZNw8uvTTqKIrPyJFw991RRyGZoIV7\nEauvh7FjQ/3S8OEwdSp0+f2v4eyzlSCLFAkt3JN8dMwx8PDDobevZM6WLXDwwbBuHXTqFHU0ooV7\nBWzOnPCRF4S/6ufOBT73OSXIIiKSNQ0NsGQJDBoUdSTFp1MnGDhQm4oUAyXJERs2rLEmbMgQ1YeJ\niOSlIlu4t2CBOltk08iR8NprUUch7aUkOWJduoQSC0iUWqivoohI/pk4MeoIMkqL9rJLHS6Kgxbu\nRSgeh82TfsqSo85g3LghTJ4cXo/FwkNERCQblCRn14gR8MADUUch7RVZkmxmk/5/e/ceLmVd7338\n/QHFIx4Qj6CgIhoq4CEgRV2loVailaZ7p+3sylIyrZ56rK0+SNnTTqyt6ZWFWyszNz6GpKWYWSwP\nmIocRWCDkme0AEVMRQ7f54/7XjoO63yY39yzPq/rmmvN4Z5Zn5m17rW+85vv/fsBJwFrgaeBsyPi\n9VR5UqibORGWToJJp/AVtyCbmVmFeGaLrjVsWPZGZN062Hzz1GmsvVK2W9wLHBgRw4GlwHcSZqk8\nz4NsZmaJeCS5a/XuDXvuCYsXp05iHZGsSI6I+yJiY37xEaD7VIoukM3MLJG33oJnnoH990+dpLZ5\nUZHiq5YD974ATEsdoiLmzoXrr3eBbGZWJOPHp07QaRYsgAMOgF69UiepbYcfDjNnpk5hHdGlPcmS\n/gTsWnoVEMDFEfH7fJuLgXURcUtXZkmhYUW9hvPZwXjD+cg18zi6/1apYpmZWVvV0BRwc+fC8OGp\nU9S+kSNh8uTUKawjkq64J+nzwDnARyJibTPbxfiSd/F1dXXUFWj6hzVrYLvt4PXXPcWbWXdQX19P\nfcM7ZGDChAlecc+qxle+ki128bWvpU5S2958E3beGVasgK08LpZMR1bcS1YkSzoB+BFwdESsbGHb\nwv6BvvvubAG9V17JVuE577ysYPY0b2bdh5eltmoyejRcfrn/B1XCYYfBNdfAEUekTtJ9deTvb8p5\nkq8BegF/kgTwSESMS5inS+zMP1ixYmcgW8/997/PGvk9omxmZpW2cSPMn59NUWZdb+RIePRRF8lF\nlaxIjoj9Un3vipk4kUNvu529Bz7MU09nb2KefTabemfUqMTZzMys21m2DHbcMTtZ1xs5EqZ1j2kJ\nalK1zG5ReyZO5M2rJ3H16Ns4aazYcsvs6j59YNWqtNHMzKyNauTAPR+0V1mjRsEjj6ROYe2V9MC9\n1ipcP1wj8yD74D2z7ss9yTVAghp4DS69FHr0gAkTUifpHjZuhJ12yhYV2XXXlre3zlfUnuTadOWV\n7yuQS6eBO+YY+NGPsvM+cM/MzCpt7lw4++zUKbqPHj1gxIisL3ns2NRprK1cJHe2/v3fN4LsYtjM\nzKqF2y0qr+HgPRfJxeOe5M52xhleSc/MzKrOypVZy9/AgamTdC/uSy4uF8lmZmbdwLx5MHRo1gJg\nlTNiBDz+OGzYkDqJtZV3FTMzs5aUrPpaVG61SKNvX9hlF1i4MHUSaysXyR1x9dUwa1bqFGZm1tVq\nYAq42bPhkENSp+iejjoKHnoodQprKxfJ7XXFFXDttbDbbqmTmJmZteixx7KP/q3yjj4a7r8/dQpr\nK8+T3B5XXAHXX5/N7davX+o0ZlblPE+ypbZqFQwYAK+9Bj17pk7T/SxbBqNHw4svZlNuW+V05O+v\nR5LbygWymZkVzOOPw2GHuUBOZe+9swMmn346dRJrCxfJbbFwIdx4owtkMzMrFLdapCVlLRcPPJA6\nibWFi+S2GDIE5s93gWxm1t0U/MA9F8npuUguHvckm5l1Mfck1wAJCvoaRMDuu2eF8l57pU7TfS1c\nCJ/4RNafbJXjnmQzMzNr1PPPZ1/33DNtju7uAx+ANWve+3lY9XOR3Jzly1MnMDMz65CGVgvPqpCW\nlM2X/OCDqZNYa7lIbsrEidnnIgX9eM3MzAzcj1xN3JdcLC6SGzNxIkyaBHfc4bfeZmZWaI89BiNH\npk5h4EVFisZFcrmGAnn6dOjfP3UaMzOrBuPHp07QLhs2wKxZcPjhqZMYwLBh8I9/uC+5KFwkl7ry\nShfIZma2qYJOAbdwIeyxB+y4Y+okBtliLscdB/femzqJtYaL5FKDB7tANjOzmvHwwzBqVOoUVur4\n410kF4WL5FJjx7pANrOaIOkESYslLZF0USO3HyPpNUmz89MlKXJa17r/fjjmmNQprNRHPwr33Ze1\nwlh1c5FsZlZjJPUArgWOBw4E/kXSAY1s+kBEHJqfLq9oSOtyEdlMCkcfnTqJlerfP1vcZdas1Ems\nJS6SzcxqzwhgaUQ8GxHrgMnAyY1s5+l7atiyZVmhvO++qZNYuTFj4I9/TJ3CWtJ9i+T//E+YMSN1\nCjOzrtAPKD1+/oX8unIfkjRX0l2ShlQmWkEV8MC9hlFkz2RafdyXXAybpQ6QRMM0b6edljqJmVkq\ns4C9IuJNSScCvwMGN7XxZSVFYl1dHXV1dV2dr7pMmFC4Qtn9yNXrqKNg7lxYvRq23z51mtpSX19P\nfX19pzyWogArykmKTsvpeZDNrMIkEREVG8+TNAq4LCJOyC9/G4iI+GEz9/kbcFhErGrkts77G1xU\nUuFWYN1nH/jDH2CIPyOoSmPGwLhxcMopqZPUto78/e1e7RYukM2se5gJDJI0QFIv4AzgztINJO1a\ncn4E2aDJJgWyFdPzz8OaNfCBD6ROYk05/niYNi11CmtO9ymSn3oKfvlLF8hmVvMiYgNwPnAv8CQw\nOSIWSfqypC/lm50qaYGkOcBVwOmJ4loXcD9y9Rs7Fu68EzZuTJ3EmlLz7Rb19dkJ4MHp6znqw1kb\ndl1ddjIz62qVbrfobG63oHDtFl/6Ehx0EFxwQeok1pyDDso+4D7iiNRJaldH/v7WfJHcYM0a2G47\neP116N27k4KZmbWCi+QacNllhTpw74ADYPJkGD48dRJrzqWXwttvZ92g1jXck9yCNWuyI0kh+7pm\nTdo8ZmZWMAUqkF94AVasgIMPTp3EWvLJT8LUqYX6kKJb6RZF8oIF8OST2fmFC987b2ZmVmumTctm\nTujZM3USa8khh8D69VmdYtWnWxTJBx0EBx6YnR8y5L3zZmZmtWbaNDjxxNQprDWkbAq4qVNTJ7HG\n1HxPcsOBe2vXZktAHn88bLGFD9wzs8pxT7JVyjvvwC67wJIl2Verfg88ABdeCHPmpE5Sm3zgnplZ\nFXORbJUyfTpcdBE89ljqJNZaGzbA7rvDo4/C3nunTlN7fOCemZlZVyrIgXtutSienj2zA/huvTV1\nEivnkWQzsy7mkeQaUJB5kg8+GK6/HkaNSp3E2mLGDDjnnGxiAS8A07k8kmxmZtbNPf88LF8OH/xg\n6iTWVkcckR07NWtW6iRWykWymZlZDZg2LTs43VO/FY8En/sc3HRT6iRWykWymZlZDbjrLvcjF9lZ\nZ2WrJL7zTuok1sBFspmZWcGtXp1Nd3rSSamTWHvtsw/svz/cc0/qJNYgWZEs6buS5kmaK+k+Sf1T\nZTEzM2vW+PGpEzTrjjuyuf+33z51EusIt1xUl5QjyVdExLCIGA7cAVyWMEvF1NfXp47QaWrludTK\n8wA/F7MuU+VTwE2eDKefnjqFddRpp8F998Hf/546iUHCIjki3ii5uA2wIlWWSqqlf/y18lxq5XmA\nn4tZd7RyZTaF2NixqZNYR+2wA5x6ajaNn6WXtCdZ0uWSngM+D/wgZRYzM7MimjoVxoyBbbdNncQ6\nwwUXwE9/CuvWpU5iXVokS/qTpPklpyfyrycBRMQlEbEX8Avgqq7MYmZmVovcalFbhg6FwYNhypTU\nSawqVtyTtCdwd0Qc3MTt6UOamXWAV9yzrvDKK9mMCC+9BFtvnTqNdZapU2HiRHj44dRJiq8jK+5t\n1tlhWkvSoIh4Kr94CjC3qW2L/M/FzMxqwGWXVeXBe7feCh//uAvkWnPSSfD1r8PMmV5BMaVkI8mS\nfgsMBjYAy4DzIsLHc5qZVRmPJJMtiVZlr0EEHHhg1r9aV5c6jXW2iRNh7lz4zW9SJym2jowkV0W7\nhZmZVS8XyVRlkVxfD+PGwZNPZvGstqxeDYMGwYMPwgEHpE5TXB0pkguz4l6tLD4i6QpJi/LnMUXS\ndqkztZekUyUtkLRB0qGp87SHpBMkLZa0RNJFqfO0l6QbJL0iaX7qLB0hqb+kv0h6Mj/Q94LUmdpL\n0haSHpU0J38+/zd1Jqst112XFckukGvT9ttnLRff/W7qJN1XYUaSJW3bMLeypK8CwyLii4ljtZmk\n44C/RMRGSf8BRER8J3Wu9pC0P7AR+DnwzYiYnThSm0jqASwBjgVeAmYCZ0TE4qTB2kHSaOAN4KaI\nGJo6T3tJ2g3YLSLmStoWmAWcXMSfCYCkrSPiTUk9gRnA/4qIGalztZVHkqm6keTly2HIEHjmGa+y\nV8vWrMlGk//8ZzjooNRpiqlbjCTXyuIjEXFfRGzMLz4CFHJEHCAi/icilgJFHccYASyNiGcjYh0w\nGTg5caZ2iYiHgFdT5+ioiHg5Iubm598AFgH90qZqv4h4Mz+7Bdnf28L/jKw6/Nd/ZdO+uUCubb17\nwze/CRMmpE7SPRWmSIaaXHzkC8C01CG6sX7A8yWXX6DABVmtkTQQGA48mjZJ+0nqIWkO8DJQHxEL\nU2eydho/PnWCd61fD5MmwXnnpU5ilTBuHDz0UHYQn1VWVRXJtbL4SEvPI9/mYmBdRNySMGqLWvNc\nzDpb3mrxW+DCsk+RCiUiNkbEIWSfGB0t6ZjUmaydqmj6t1tugX33hWHDUiexSthmG7jkErjwwqrq\n+OkWks2T3JiI+GgrN70FuLsrs3RES89D0ueBjwEfqUigDmjDz6SIXgT2KrncP7/OEpK0GVmB/OuI\nuCN1ns4QEa9Lugs4HLg/dR4rrnXrsgO5brghdRKrpHPPzX7mv/kNnHlm6jTdR1WNJDdH0qCSi80u\nPlLNJJ0AfAsYGxFrU+fpREXsS54JDJI0QFIv4AzgzsSZOkIU8+dQ7kZgYURcnTpIR0jqK2n7/PxW\nwEcp6N8tqx6//jUMHAjH+DOJbqVnz2w2k299C157LXWa7qNIs1vUxOIjkpYCvYCV+VWPRMS4hJHa\nTdIpwDVAX+A1YG5EnJg2Vdvkb1quJnvDeENE/EfiSO0i6RagDtgJeAUYHxG/SBqqHSQdCTwAPAFE\nfvr3iLgnabB2kHQw8CuyNy49yEbGr0ybqn08u0V1eOedbAnqm2+GI49MncZSOPdc2GwzuPba1EmK\nw4uJmJlZl3GRXB1+/nOYOhXuKdxbRussq1ZlqyzedhuMHp06TTF0iyngzMzMkkl84N7q1fC973lh\nie6uT59s+r/PfjYrmK1reSTZzMya5ZFkki8mMm7ce1O/mX3969lCMrff7hUXW+J2CzMz6zIukkla\nJM+YAZ/5DDz5JOywQ5IIVmXWroUjjoCzz4bzz0+dprp1pEiuqingzMzM7D1r18I558BVV7lAtvds\nsQVMnpz1JQ8eDGPGpE5Um9yTbGZmVqUmTIBBg+DUU1MnsWqz334wZUrWnzx7duo0tclFspmZWRW6\n/fZsurfrr3ffqTVu9Ohs1pNPfAKWLUudpva43cKqiqQ+wJ/J5sfdnWxe7H/kl0dExPpEuY4Fzo+I\nT6b4/maW2PjxFf12CxbAl78M06bBrrtW9FtbwXzqU/D3v0NdXTY94JAhqRPVDhfJVlUiYhVwCICk\n/wO8ERE/Lt9OaY4k6uZHLpl1YxWcAm7FCjjlFPjxj+Hwwyv2ba3Azj0Xtt0WPvxh+N3v4EMfSp2o\nNrjdwqrZux8wStpX0pOSbpa0ANhT0qslt58u6fr8/C6Spkh6TNIjkkZs8sDSTEn7lVx+UNJQSSMl\nPSxpVn7dvo3c93uSLii5vEjSHvn5z0l6VNJsSdfm1/WUdJOkeZLmS/KxyGbWqJdfzkYEzzgDzjor\ndRorkjPPhF/+EsaOhf/+79RpaoOLZCuS/YEfRcRBwItsOrLbcPknwA8jYgRwOnBDI481Ob8NSf2A\nHSNiPrAQGB0RhwGXA99vRa7IH+dA4JPAhyLiUGBzSWcAhwF9I2JYRAwFbmrtEzaz7uP55+Hoo7MC\n+XvfS53GiujEE+Hee7PuoC9+Ef75z9SJis1FshXJ0xExpxXbHQf8TNIc4HfA9pK2KNvmNqDhePHT\n88sAOwK3S3oCuBJoS3fXccDhwOP59z4a2Bd4Chgs6SpJYyLi9TY8ppl1AzNnwlFHZX3Il1ziA/Ws\n/Q45BGbNgnfeydp1pk9Pnai43JNsRVL6nngj73+Tt2XZth+MiA1NPVBEPCfpDUkfICuS/y2/6fvA\nPRHxs7zVYlojd19f9r23yr8KuDEiNjnCR9JQ4ERgnKRPR8SXm8pmZt3Hxo1Z7/EVV8B118GnP506\nkdWC3r3hppuyGVLOPhtGjoSJE2GvvVInKxaPJFuRvDu2kh+0tyrvVe5B1ubQ4D7gq+/eSRrWxOPd\nCnwH6BURi/PrtiNr5QA4u4n7PUPWQkHe77xnyff9jKSd8tv6SNpTUl+gR0RMAcaTH5hoZgXSBQfu\nzZ4Nxx6bFTIzZ7pAts73qU/BwoVwwAEwfDh84QuwaFHqVMXhItmKpLwH+dvAvcBDwPMl158PHJkf\nKLcA+GITj/db4F/IiuUGVwBXSnq8ke/X4DZgN0nz88d+GiAiFgATgPskzQP+COxCVkQ/kLdg3EhW\nmJt1KUknSFosaYmki5rY5ieSlkqaK2l4pTMWyoQJnfZQixdnB1l9/ONw+unwwAMwYECnPbzZ+2y9\ndfbr+9RTsPfecMwx2Qp9N9/snuWWqPKzaJmZWVfKP11ZAhwLvATMBM4o+cQESSeSzf39cUkjgasj\nYlQTj5dgxsX2qa+vp66urvMfWIIOvAZr1sBdd2ULPyxaBOedB9/4Bsya1UV5u0CXvbZdpEh5K5n1\nrbfgjjuydoyHH86mjfvYx+D441vfjlGk11YSEdGuLn+PJJuZ1Z4RwNKIeDYi1pHN5nJy2TYnk8+0\nEhGPkh3gWvhlK+rr61NHAODNN2HGjKzfeMwY2GMP+MUvYNw4eO65bPaB3r2rJ29rFCkrFCtvJbNu\ntVU2g8rdd8PSpVmbz/Tp2UF+/fvDaafBD3+Yval79tmsbz5l3pR84J6ZWe3px/tbkF4gK5yb2+bF\n/LpXujZasUVkswa89hq8+iqsXAnLl8NLL2VTuC1dCkuWwDPPwIEHwogR2ajxlClZUWxWTXbeOWv9\nOfPM7Hd72TL4619hzhz4yU+ylR9XrsxGmAcOhH79YPfds376m2+Gvn1hhx1gu+2y3+9ttslOvXrV\nxgwtLpLNzKxqPP00fO1rrd++vANiyZLsILjy20u/lp82bnzvtH49bNiQfV23Lju9/Tas5UXe3CHr\n4ZSywqBPH9hxx6xo2GOPbBRu1CgYPBj22w+2LJ9zx6yKSbDvvtnpzDPfu/6tt7I3fc88k70ZXL4c\nVq3KRqJXroTVq+H117Ov//xn9inKhg3Z7/+WW2YFc69esPnm75169nz/qUeP7CS997WpU0PWhvOD\nB8NVV3XRa1KUPjMzM2sdSaOAyyLihPzyt8kmhflhyTY/A6ZHxK355cXAMRGxyUiyJP+jMLPCam9P\nskeSzcxqz0xgkKQBwHLgDLKZXErdCXwFuDUvql9rrECG9v+DMTMrMhfJZmY1JiI2SDqfbIrEHsAN\nEbFI0pezm2NSRNwt6WOSniJbqKepecHNzLolt1uYmZmZmZXxFHBmZoakGyS9ki+S09jt/5ov0DNP\n0kOSDq50xpIszWYt2e6DktZJ+lSlsjWRo8W8kuokzZG0QNL0SuYry9HS78FOkqblC9A8IenzFY5Y\nnqe/pL9IejLPc0ET2yVfOKc1WatsP2vVa5tvm3Rfa8PvQZv2M48km5kZkkYDbwA3RcTQRm4fBSyK\niNWSTiA7MLDRxUe6WktZ8216AH8C3gJujIjbKxixPEtLr+32wMPAmIh4UVLfiFhR6Zx5lpayjge2\njIjvSOoL/A+wa0Ssr3DUhjy7AbtFxFxJ2wKzgJPbu3BOFWStpv2sxbz5dsn3tVa+tm3ezzySbGZm\nRMRDwKvN3P5IRKzOLz5CNqdyEi1lzX2VbOn5v3d9oua1Iu+/AlMi4sV8+yQFcv69W8r6MtAw43Nv\nYGWqAhkgIl6OiLn5+TeARWz6u1kVC+e0JmuV7WeteW2hCva1VmZt837mItnMzNrqi8C01CGaImkP\n4JSIuA4owswcg4E+kqZLminprNSBmnE9cKCkl4B5wIWJ87xL0kBgOPBo2U1NLZyTTDNZS1XNftZU\n3mrc15p5bdu8n3l2CzMzazVJHyabCWN06izNuAq4qORyVfzzbsZmwKHAR4BtgL9K+mtEPJU2VqO+\nA8yLiA9L2hf4k6Sh+ehdMvlH7L8FLkydpSWtyVpN+1kLeatqX2sha5v3MxfJZmbWKpKGApOAEyKi\npXaHlA4HJksS0Bc4UdK6iLgzca6mvACsiIi3gbclPQAMA6qxSD4S+D5ARDwt6W/AAcDjqQJJ2oys\nMPp1RNzRyCYvAnuWXO6fX1dxrchaVftZK/JWzb7Wiqxt3s/cbmFmZg1EEyNBkvYCpgBnRcTTFU3V\nuCazRsQ++Wlvsn+a46qgQG4yL3AHMFpST0lbAyPJeipTaS7rIuA4gLyvdzCwrEK5mnIjsDAirm7i\n9juBz8G7B8Y1uXBOBTSbtQr3s2bzVtm+1tLvQZv3M48km5kZkm4B6oCdJD0HjAd6kS8+AlwK9AF+\nmo8arYuIEVWatVTyKZxayhsRiyX9EZgPbAAmRcTCaswK/AD4haR5ZIX0/46IVSmy5nmPBD4LPCFp\nDtnP+9+BAVTZwjmtyUp17WetyVsq2b7Wyt+DNu9nngLOzMzMzKyM2y3MzMzMzMq4SDYzMzMzK+Mi\n2czMzMysjItkMzMzM7MyLpLNzMzMzMq4SDYzMzMzK+Mi2czMzMysjItkMzMzM7MyLpLNzMxsE5I2\nSJot6QlJt0rashMe86GO3G5WSV5xz8zMzDYh6fWI2C4/fzPweERcVbaNwoWE1SiPJJuZmVlLHgQG\nSRogabGkX0l6Augv6bOSHs1Hna+TJABJn5M0T9IcSb/Kr1uTf91a0h/y2+ZLOq309vz8N/JR7PmS\nLsyvGyBpoaRJkhZIukfSFpV+Max7cJFsZmZmjWkodjcDTgSeyK/fD7g2Ig4GtgFOB46IiEOBjcBn\nJQ0BLgbqIuIQ4ML8vg2jzicAL0bEIRExFLin9HZJhwH/BnwQ+BBwjqRh+TaDgGsi4iBgNfDpTYJL\nA/MWkccl3Slpiotpa6vNUgcwMzOzqrSVpNn5+QeBG4B+wDMRMTO//ljgUGBmPoK8JfAKsAPw/yLi\nVYCIeK3ssZ8ArpT0A+CuiCjvRT4SmBoRbwNIuh04Cvg98LeIaCjYZwEDG8neLyJOl/SliJjUjudu\n5iLZzMzMGvVmPjr8rryT4p+lVwG/ioiLy7Y7v7kHjoilkg4FPgZcLum+iLi8lbnWlpzfQFaYlz/+\njPzsHq18TLNNuN3CzMzMGqNWXP9n4FRJOwNI2lHSXsBfgNMk9Wm4vvS+knYH3oqIW4CJZKPRpY/9\nIHCKpC0lbQN8Mr+uuVzvDyntw/sLerM28UiymZmZNaapWSvevT4iFkm6BLhXUg/gHeArEfGYpO8D\n90taD8wBvlBy34OBiZI25vc5t/SxI2KOpF8CM/PrJkXEPEkDmslVblR+f7N28RRwZmZmZmZl3G5h\nZmZmZlbGRbKZmZmZWRkXyWZmZmZmZVwkm5mZmZmVcZFsZmZmZlbGRbKZmZmZWRkXyWZmZmZmZVwk\nm5mZmZmV+f9BTAST0TjbSgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x110113ed0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Now consider a setting where n \\gtrsim p\n",
"n, p, precision = 500, 100, 1.5\n",
"X, theta, y = simulate(n, p, precision)\n",
"estimate_theta, cov_theta, shape_precision, scale_precision = infer(X, y)\n",
"plot(estimate_theta, cov_theta, shape_precision, scale_precision)\n",
"\n",
"print \"One-sided p-value: {}\".format(stats.gamma.cdf(precision, shape_precision, scale=1.0/scale_precision))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment