Horseshoe prior fail

model_check_test_point.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import theano.tensor as tt\n",
    "import theano\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAElCAYAAAAhjw8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8VNX5/993MskkmWwkQ4AESMKeBEVBEVEiSrXu6Fe9VVvrVmnFrS2utS7V8sO2torVal1wqwu3rlgtiFvV1p1qFaICJhECIYSQkEzIJDP3/v6YmThJZpKZyazkeb9eeSVz77nnPnMvnM85z3POcxTDMBAEQRCEoWCKtwGCIAhC8iNiIgiCIAwZERNBEARhyIiYCIIgCENGxEQQBEEYMiImgiAIwpAxx9uAZOWTTz4xpaWlXZ2SklKOiLIgCPFHd7lc1V1dXb+fNWuWHuubi5iESVpa2tU5OTlqampqzF+aIAiCP7q7u/fbs2cPwG2xvrf0qMMkJSWlXIREEIREIjU1Vfd4S2KOiEn4yLMTBCERiUvbJA1iErNt2zbzeeedVzZr1qz95s+fX3700UdP+/vf/54Xaztmzpy5344dO/q5TJcuXTo6nPqeeeaZvP/973/p3s8nnHDC1Pfffz9zKDYmIjt37kypqqqqqKqqqpg+ffqM/ffff3/vZ4fDoQRTx0UXXVS6fv16S7g2zJgxY//m5uaUQOddLhfLli0L6z0mE5F4FwAPPPBAQX19/aDhgy+//NJSVVVVMVCZjRs3pj3xxBMjgr13vBExSVJ0Xefcc8+dNHv27PZPPvnk87feeqv6r3/96zfbtm1L61u2u7s7HiayYsWKMf6O67qOy+UKeN2aNWvyNmzYkBE1wxKEkSNHut5+++0Nb7/99obTTz995znnnLPD+9lisRgw+LN64IEHaisrKx3RslHXdR5++OF9XkyCeRfB8Mwzz9i2b9+eGgmbNm/ebFm1alV+JOqKBSImScratWuzU1NTjcWLF+/0HpswYULXFVdc0Qjw0EMPFaiqOunEE0+csnDhwqm6rnP11VePPeywwyoPP/zwiieffHIEwGuvvZZ9+umnT/LWccUVV4x/6KGHCsA94rjpppuK5s+fX3744YdXfPHFF+ng7sUtXLhw8ty5cyt/9rOflfhLFnrdddcVd3V1maqqqirOO++8ss2bN6cdcsgh0y+88MLSww8/vLKuri5t4sSJB3rLP/300yMuuuii0rffftv6zjvv5P3hD38YW1VVVfH1119bAF544YURCxYsKJ89e/b0N998MytKjzUh+PLLLy2HHnpo5XnnnVd22GGHVdbX16cuXry45MgjjyyfO3du5S233NIj0t///venfvzxxxnd3d1MmTLlgOuuu6543rx5FUcfffS07du39+shNzQ0mL3vbvHixb3e3RlnnDFp/vz55XPnzq3861//agO4/vrrx+7duzelqqqq4sILLywNVG5fZsWKFQULFiwor6qqqrjsssvGu1wuuru7ueCCC8oOP/zwisMOO6zyzjvvLHzyySdHbNq0KfOSSy6Z6G9E895772XOmzevoqqqqmLFihUjvce//vpry7HHHjv1iCOOqJg/f375v//9byvA7373u7Hr1q3LrqqqqrjzzjsLA5VLFGQ2VwTQH7lrnFFfF1E3jFJc0mE67/Itgc5XV1dnlJeXdwxUx1dffZX5+uuvr7fZbK6VK1fmffXVVxn/+te/1jc2NpqPP/748qqqqvbB7MjPz3e+9dZb1XfffffIu+++e9R9991Xt3Tp0qJZs2a133jjjdtffPHF3Jdeeqlfg7Js2bJ6TdMK33777Q0AmzdvTquvr7f88Y9/rDn88MNrA92vqqrKPm/evJYFCxa0nnnmmbu9x10ul/L6669Xr1q1KveOO+4oOvLII78ezPZkZsuWLel33HFHzZw5czoAbrnllq02m83V3d3NSSedNPV///vf7v3337/T9xq73Z4yd+7ctmXLltVfeeWVYx955BHbdddd1+BbZunSpUWzZ89uu/766xueffbZvBdffLHn3d177701NpvN1d7ebjr66KPLTz/99N1Lly7d+sILL9i87zFQuYKCgsDDpyTm008/TV+zZk3e6tWrq1NTU1m8eHHJk08+mT9hwgRHS0uL+d13390A0NzcnJKfn+969NFHC5cuXfrtQQcdtLdvXUuWLClbtmxZ7RFHHGFfsmTJOO/xoqKi7ueff/7rjIwM44svvki/4oorSl9//fUvr7nmmq0PP/xw4cqVKzcDtLe3m/yVi93TGBgZmewjXH755ePnzZtXcdRRR/XM5Jg9e/Yem83mAvjwww+zTzzxxGaz2UxRUZFz5syZ7R9++OGgAnjqqafuBjjwwAM7tm3bZgFYt25d9llnnbULYOHCha3Z2dlBNSSjRo3qOvzww+3hfL8TTzxxN8BBBx1k3759ez9XXqzQ7W0m/YtPrLq9Lar/d4qKihxeIQF46qmn8ufPn19+1FFHVdTV1aX7cwNaLBb9pJNO2gOw//77d2zdurXfc1q3bl322Wef3Qxw2mmntWRmZvbMSFy+fPmoefPmVRx33HHTmpqa0jZu3Og3FhNsuWjT5nCa1m1rt7Y5nFF7F2+88UZOdXW1dcGCBRVVVVUVn3zySXZtba1lypQpnXV1dek///nPx7388ss5eXl5A/4f2LFjh9nhcJiOOOIIO8APfvCDXd5zDodDufjii0sPO+ywysWLF0+oq6vz6+INtly8kJFJBBhoBBEtysvL97766qs9wbm77rrr2x07dpiPO+64HjHxbSgCYTabDV9XR9+huddfnJKSYrhcrqADkf5IT0/vZY+ifFfdYEFOHzsYqh1DoubrDKNtj1mp+TqD6bPCEsZg8H1W1dXVlscff3zU6tWrq/Pz813nn39+mb/nZTabe15kqO9r9erV2Z988kn26tWrq61Wq/H9739/amdnZ79GOthysWDjrs6MVofTvHFXZ8bMoqyovAvDMFi4cGHTrbfeuq3vuTfffHP9K6+8kvvoo48WvvzyyyP+8pe/1IVzjzvvvHPUmDFjulasWFHT3d2tTJky5cChlIsXMjJJUo4++ug2h8Oh3HPPPT2+146OjoDv85BDDml7+eWX851OJw0NDeb//ve/WYcccoi9tLTUUVNTk9HZ2ak0NzenfPTRRzmD3XvmzJltK1euLAB46aWXctra2vzOBjKbzUZXV1fABm3EiBHdn3/+ebrL5WLNmjU9wmi1Wl3t7e2J+W+zbMpeJTvHSdmUfm6MaLFnz56UzMxMV25urmvr1q2pH3zwwaDvKBAzZ85se/rpp/MBnn/++Vzvv5nW1taUnJwcp9VqNT777LP0r776ygqQmuqOJXsncQQqFw8mF6TvzbWYnZML0qP2Lo466qi2V199Nd87W3Hnzp0pNTU1aQ0NDWbDMDjzzDN3X3311fXV1dWZAFarVff3/2HUqFFOi8Wiv/POO1aAv//97z2B9ba2tpTCwsJuk8nEww8/XODt3GVnZ+sdHR0pg5VLFGRkkqSYTCYeffTRzb/61a/GrVixYnReXp4zPT3ddeWVV271V/70009v+fjjj7OOOOKISkVRjKuuumprcXGxE+Doo4/eXVVVVTlmzBjHlClTBozDAFx//fXbfvKTn0yYO3du5YwZM9oLCwu7/JU79dRTdx5xxBEVU6dO7bjhhhvq+55fsmRJ/fnnnz8pLy/PWVlZ2eFt2E499dTma6+9tvTxxx8f9cADD2wO7clEF5M1W4/miMQfBx98cMeECRM6Dz300Oljxoxx7LfffoPGugJx/fXXb1u0aNGEww47rGDGjBntBQUF3QAnn3xy61NPPTXy0EMPrSwpKeksLy/v+Y4LFy5sqqqqqqyoqLDffffddYHKxZpsi1mP1ojEy4EHHrj30ksv3XbGGWdMMQyDlJQUY9myZXUpKSlcddVVpYZhoCgK11xzzVaAM844o+naa68ttVgs+tq1a6t9Z4LdfvvtNddcc02poijMmTNnj/f4RRdd1PjTn/504nPPPWerqqpqTU1NNQAOOuigDl3XlXnz5lWcdtppTYHKJQpKoqlbsrBhw4bH8/Ly4rLSVBAEIRAtLS3VFRUV58T6vonpShAEQRCSChETQRAEYciImAiCIAhDRsQkfCRjsCAIiUhc2iYRkzBxuVzV3d3d8vwEQUgYuru7TS6Xqzoe95apwWHS1dX1+z179iA7LQqCkCD07LQYj5vL1GBBEARhyEiPWhAEQRgyIiaCIAjCkBExEQRBEIaMiIkgCIIwZBJmNpeqqnnAg8B0wAAuAL4CVgKlQC2gapq2W1VVBVgOHA90AOdpmrYuDmYLgiAIJNbIZDmwWtO0acAMoBq4Fnhd07TJwOuezwDHAZM9P4uAe2NvriAIguAlIcREVdVcoAp4CEDTtC5N01qAhcCjnmKPAqd4/l4IPKZpmqFp2vtAnqqqYxAEQRDiQqK4ucqAncDDqqrOAD4BrgBGaZq23VOmARjl+bsY8N3dcKvn2HYCIwtqBEEQwmPQXTsTRUzMwEzgMk3TPlBVdTnfubQA0DTNUFU1JEFQVXURbjcYmqbR1eV3D6fgDDSbcTqdYV8fLcSu0BC7QkPsCo190a60tLTg7hFW7ZFnK7BV07QPPJ+fwS0mO1RVHaNp2naPG6vRc74eGOdz/VjPsV5omnY/cL/no9HU1BS2gTabjaFcHy3ErtAQu0JD7AqNfdGuoqKioMolRMxE07QGYIuqqlM9hxYAG4BVwLmeY+cCL3r+XgX8WFVVRVXVOUCrjztMEARBiDGJMjIBuAx4QlXVNOAb4HzcYqepqnohUAeonrKv4J4WvAn31ODzY2+uIAiC4CVhxETTtE+Bg/ycWuCnrAFcEnWjBEEQhKBICDeXIAiCkNyImAiCIAhDRsREEARBGDIiJoIgCMKQETERBCGpsXe5WN/Ygb3LFW9ThjUiJoIgJDW1LQ7aHE5qWxzxNmVYI2IiCEJConfY0TeuR++wD1iuNM9CtsVMaZ4lRpYJ/hAxEQQhMamvhfY29+8+tDucPa4ta1oKlYWZWNNSYm6i8B0iJoIgxI0BRx/FpZCV7f7dh81NdnFtJRgiJoIgxI8BRh+mTCumyZWYMq39zk20Wfu5tiQQH19ETARBiB8DjD4GIsti7nFteUc3NQ0tMlqJIyImgiDEjYFGH0HjGd2U2BskEB9HEibRoyAIQlgUl0J9LdbiUiozM+NtzbBFRiaCICQ1ERndhEmgOE2w05r3JURMBEFIWuIddA+4YHKAiQX7KiImgiAkLfFe/R5wwWSYEwuSGRETQRCSlr6NeazdS4EWTMbT9RYvREwEQUha+jXmw9C9lCiImAiCsO8QY/fScAy0B0LERBCEqBFKYxuJYHrM3UsyEupBxEQQhOgRQmMb72B6WAQYCcV7llk8EDERBCF6hOB2CieVfLwb7UAjoaQUxiEiYiIIQtQIxe0UTCp5r9vMZW8HoKahhT01NdQ0tETM5qGid9gZv6uGLKN7WKV2ETERBCFh6Rdz8bjNnHWbANz5uJwdlNgbQq8rwuV7qK/F2tlGhaNhWO2xImIiCELi0jfm4nGbmUsmAWAtKaWiMBNrSWnodUW6vJdhuGAREizRo6qqKcDHQL2maSeqqloGPA0UAJ8A52ia1qWqqgV4DJgF7AJ+oGlabZzMFgRhEPQOu7tRLi716/Kyd7mobXFQmmfp3Zv3JHH0NsymTCtMriTFmgV7O3s+B0WfuiJe3kNINu1DJNrI5Aqg2ufz74A7NE2bBOwGLvQcvxDY7Tl+h6ecIAgJRk+AvK52wF5+oID1QDEXf26ogVxTfesaLHg/HFexD4WEERNVVccCJwAPej4rwFHAM54ijwKneP5e6PmM5/wCT3lBEBIIr0jUWUcP6PoJZyaXPzeUva6WDV7xCtK2RJlxlewLIBNGTIA7gasB3fO5AGjRNM3p+bwVKPb8XQxsAfCcb/WUFwQhgfCKRNnovAF7+cHM5OqHn9hEnXU0beZMt3gFaVvCzLhK8gWQCREzUVX1RKBR07RPVFWdH8F6FwGLADRNw2azhV2X2Wwe0vXRQuwKDbErNIZqlw0oKQp8vt3hZHOTnYk2K1mW4Jsjs9lM4fgSGF/S6/hB2XlsHjU6qPoGsy0chvK8XBmzcNZtwlwyyR0TShC7gr5HVGsPnsOAk1VVPR5IB3KA5UCeqqpmz+hjLFDvKV8PjAO2qqpqBnJxB+J7oWna/cD9no9GU1NT2AbabDaGcn20ELtCQ+wKjYHsGiyoHgzrGztoczhpaW2lsjD4XRL72mU4nbC7CZp2MGZnAx1NO+hoawVdh5QUz48Zw+UCx14omYSpuARGjoYRNhRTeE6avs9gyO+xcCzs7XT/RJCh2FVUFJziJoSYaJp2HXAdgGdkcqWmaT9UVfXvwOm4Z3SdC7zouWSV5/N7nvNvaJpmxNpuQRjW+Lplwpy9VJpnobaFfq4mvcOO8U01GArKxGn9xMrVvBP9/X/Bl//DqPkaGra6hcOLyQRZOWAYoCjucy4XOLvdP++92eNPJ9MKE8tRJleiTC6H0sko5tR+9vgVzgg8g32FhBCTAbgGeFpV1d8C/wUe8hx/CHhcVdVNQDNwZpzsE4R9grBGGWFOnfXFGyvpR30t1H0DCpCejjGpAuo2YXz0Dsb/PqKpweOkyMyCidNQDjgERo5GsY0C2ygYYcP45kt3Q5+VjWly5XcC5TIgfySm9laMnduhdhPGxg0Yn3+M4alTOXAOysHzYNr+KCkpgUXDzzOIxIgtGVEMY9h06I1t27aFfXEyuiHiidgVGvG2S9+4vlfDG45dAdeKBHm+lz3ehr+pCXZug/++DzsbMFJS0KfuT8bMQ/kqdxxjpk3Bmp4asA7fRj3Qd/RitLXCxvUY//0A49P3oXMv5I5AOeI4jIPnobS1DCoQNpuNxvf+NeB94kEE3FyDzpZN9JGJIAixoLjU3Xh3dqJ32MPqUX831Ra/ow3f8+VZRsDeu2EYKBs3YKx+BjZVg2KCafux/fCT2T5pJpl5eeTl5tLUuAvHnm4qA4hJv8WDg4yklOxcmDkXZeZcjO4u+GId+juvYqx6El7RYPYRKCeobrfYQERgxJaMiJgIguDuuVsyhuT/7xv/6DsyKLSaqd3dycT8dKjf2O9ehu7C+OQ/GK/8HbbWQkEhyhnno8yZj5IzghFdLlo9I5uCAistra0UWs2sb+zwO9rpOxIaaGV637JKahocOIeUA+dgNGzFeOMfGP9+DeODf6EccSzKCSpKTl7gZ+kRFN0rKMPA7SViIgiCmyH2qPvFP/rEGRrtTgoyze7fPvcynN0Y77+F8c9noXEbjB6Lcv7PUWZXoZjNfuvPcHZS3lrDhs7RtCupfkdDg42Ugi2rjB6LcvbPMI4/A+OllRhvvYLx79fg4Hlw6jmk+BOVvmtG2tuw19XybUFZUG6+ZETERBAEIPI5pYz8QthSA+MnAr1HLqa0FIzSKRjvrEFf8xw0N8H4CZh+di0cOGfQqbrOuk3Q3kaJk54Gui99R0oDxWwCzSrzRckrQDlnMcbRJ6M/dg+8uxY2fIpx4S9RpvR5bn2Fub6WOsto2oMUt2RExEQQhKigNDfCCBs0N0LByJ6RhWEYGOvew/n3hzE1NaBPLMd8ziVQORNFCS4rkrlkErS2Yi0upTIzuIZ5oNFHwFll/r7X6LEol/7aLYRvvIxx+69wff800hae3VOmnzBPrqTMR8z2RURMBEEIiaBnZflxmxk7tqE/+VfY8F+6C8dSe/ZVGJWzQu6pp1izBp0p1Vc8ghl9BMLfd1YmTOXriXMoePUpRq5+Bv3rz3Fe9Vsw+68/FMFKRkRMBEEIiWBjEb69c1fLbnj+UfjgbUhLQzlzEcbcYzDanFHpqdu7XOztdpFqUnrqH0pj3u87e2Ii49Ph29MvJvfAg0l76l6af3ku/GgxptlVEfw2yYGIiSAIIRFqD9/YVA33/8Gd7qR8BqYLfoGSl48VqMxIi4qNtS0OnLpBtsUckWB3v+/sGXX1uNkKqzCmTCPlkeV0P3A7+tYalFPOCTtNSzIiYiIIQkgE28O3t3dgf/Zx8v/9CsqIAlAvQDnsaJQYTI8dikvLH32/s7/JCkpBISNuuZuddy/D+OezGNu2YLroShRLekRsSHRETARBiBjetSWGAcrDf6agaRu7Z3+PgnN+gpIeu3hBvOITSmoqyo8uhuLxGE8/iP6nGzBdfiOKNTvmtsSa4TMGEwQh6hhbazDe/xf86QbSujv59vzryTj3kh4hSbYNoAbbjdEfiqJgOupETD+7Br79Bv1312I0+09lEk79iYqIiSAIEcHY2wGvvQRvr4YplaTctJyyuYf0jlmEuQFUrBrdvmI3lN0YlZmH4rj0RlzNTbhuuxpjV2O/Mom22+NQEDERBGHIGN9uRv/tL+CzD+Cks+BEFSPFjxfdz+6IwRCzRreP2A11N8aakZOpPvfX6Hs70Jf/BsPe3ut8wu32OARETARBCBvDMNDfegV92VXQ1QWX3uAWi+ZdfkcfpkzrgNv3BiJmjW4fseu7nXCgEZJ3ROPyIxbmkom4fnod7NyO/pelGN3dPfUAoW9XnKCImAiC0Itg4xrG3g6M+/+A8cR9UH4AphuXo6SngzUbOtojmjU3rD3ig6Dvdx1M7AKOkDwjGmfdJr92Z0yfgXLeFfD1eoyH76S2ea/fepI5hiJiIggC8F3Damz+ctC4hrF9K/rSJRjr/oNy2rmYLv01SnYOFJei5NtQ5hyJKdPqt3GMZxC+371DjOEEHCF5RjTmkkkBrzUdcgTK6edhfPQOk/6l+a0nmWMoIiaCMIzp1bh6G1bF6OXq6dsAG1+sQ192JXS0Y1qyFNOxp/Uszuvbs/fbOIYZhPe1pa87KWj63jvEGE6gEZL3e6dYswYcXSjHnIpy5PGYX3uB8s/f8JtwMlx3XrxHNSImgjCc8W1cPQ2rMqG8t6vHp4z+1j/R/3wLFIzCdP2f+mfL7YPfxjHMILyvLX3dSUHT596hxHCCbawHGl0oioJy5kVwwCEYKx/AWPder/OhuvN8bYr3qEbERBCGMz6Na8CGtbgUI9OK8dG7GE/cC6WT4bIbUApGDlq9v8ZxsAZ8wEY7CHfSQJgyrewtmUZ1uxJyDz7Yxnqw0UWHE75auBi9dDL6g390uxXDxNemeM8MEzERhH0EvcNOW/V6vvh2V9ANZTA9cyUlxb1+5M2XYeahcMypKLt2RMrsfgzUaPu6k6JR/0AE21gPNrqobXGwx0hh81lXwogC9LtvxWioD8kWfzZFa5JCsIiYCMK+Qn0tdc0dtG1viJirw9XchP6HX8FnH6KceRHKuZej5ORGdX/zaPeww60/mMa63eEc1BXmvf+44pGYrrgZFBP68psx9uwOyZ5gbYoVIiaCsK9QXEpJfibZY0ZHpCE2tm+l+dqLoGErpkt+hWnBSWGvEwmFaDeQ0ax/c5N90FGP7/2VwjGYLrsB9uxGv+tWDEdnxG2KFSImgrCPYMq0kl1eyfTxBQEbymCDyEbNRvTfXwNdXZiuWoYyY3Y0TA5IsuXw8lKaaZDVuIXxFmfQ1yhlUzAtutqdx2vFnRiGMWQ74jGzS8REEIYRwcQLjOrP0P/4a7BkMOL/3YdSMjGGFnoYwvTheJK+vZYKUxvWxi0hXddRPovG758F6/6D8a9/DtmOeMzsEjERhGHEYPEC45P/oN/1G7AVYrr2dyi5I+IzQhjK9OE4Yi6ZBFnZGPmFIT232hYHdQd/n/YpB2CsfAhja82Q7IjHzK6E2M9EVdVxwGPAKMAA7tc0bbmqqvnASqAUqAVUTdN2q6qqAMuB44EO4DxN09bFw3ZBSCYG2udDf3sNxt/uhQlTMF12I4o1y72ewztCGGTP9Ujib/OpZMC7N72+cX1Iz827mZf5gitg2RL0v/4B06//FPbGWvHYzyVRRiZOYImmaRXAHOASVVUrgGuB1zVNmwy87vkMcBww2fOzCLg39iYLQuzRO+y4vvgY1+efRHS0oK9+FuPxe6DyQEy/uBXFM/XW29NOthFC3Al3ZX1BAaYLfwE76jFWPjjodYkUW0oIMdE0bbt3ZKFpWhtQDRQDC4FHPcUeBU7x/L0QeEzTNEPTtPeBPFVVx8TYbEGIKn4bivpaqPsGtmz2G08IJ/Cqr34W49lHUQ6eh+mS61Es37lGvD1tf7O34p2+I5EZyqw3pXwGyvf/D+OdVzH+99HAhRMotpQQbi5fVFUtBQ4EPgBGaZq23XOqAbcbDNxC4xvh2uo5tt3nGKqqLsI9ckHTNGw2W9h2mc3mIV0fLcSu0EgmuxwbtmKYFJT23VjGlwDgyphFV1oaCpA6bb9+i/fq61tRLCm0GKmU2HIHva/9hSdpf/ZR0ucdTc4VN6D02YNkoOcVzL3aHU42N9mZaLOSZTEPejxY4vUeXfZ2nHWbMJdM8rtwsq9dg5UPhHHBZTR/+Rn64/eQv/xvmHLy/NuTMSuo+mPxvBJKTFRVzQKeBX6uadoeVVV7zmmaZqiqGtKcOU3T7gfu93w0mpr8b50ZDDabjaFcHy3ErtBIJrv0rBHQ2gpZIzD5niuZ7P69t9P940Oe4qLF4SAvw0JTU/eA99RffQHj7ytQDp5H1w8Xs2t3S1B2hXKv9Y0dtDmctLS29vLhBzoeLPF6jz2xkNZWTH5iIX3tGqx8Tzlvok1PWhsA49zL0H+7hKa7l2FadFVgowrH+v23MJBdoVBUVBRUuYRwcwGoqpqKW0ie0DTtOc/hHV73lee3d9/LemCcz+VjPccEYZ8hHFdJsAvy9NdedAvJrMNQLvylO2VKkHjdbxnOzkHvFWhWUbzzSIVNqLPMgi3vx12ljC1DOVHF+OgdjM8/DsvcWJIQYuKZnfUQUK1p2p98Tq0CzvX8fS7wos/xH6uqqqiqOgdo9XGHCYIwAPrrL2GsfAhmzUX5yZKQhAQIyU8fSNxikQYkGjGdYATe9759ywe0KYDoKMeeBmPGof/tXozOvRH7HtEgIcQEOAw4BzhKVdVPPT/HA7cBR6uquhH4nuczwCvAN8Am4AHBG8M4AAAgAElEQVRgcRxsFoSkwLcB09/4B8bTD8DMQzH95EoUs7lfmUFJkjUg4Szci4QA9b1vMGnifUXH3uXii2930Va9HqOrC9OPL4XmnRgvPhG2TbEgIWImmqa9CygBTi/wU94ALomqUYKQwHgbJm+22IHKdHa76NYNWlavIf2lh+HAOZguuqpHSMC3AWTQGEayrAHxrt0IxZUWynMY7L6FVjPrGzt6nr/XlsFsqm1x0La9gTpnBxXmWpTJlSjzj8N4/R8Ys49AKZscll3RJlFGJoIg+DDY+oFget3eMgYw7tM3GP3Sw3DAIZgW9RYSSOIYxgCE40qLxHPw3rfR7ux5/qGkiS/Ns5A9ZjQl+Zk9oz/l1B9Dbh76Y3djOAPn/YrnuhMRE0FIRAaJSwTT6HnLTP/yHUavWgEzZmP66dUo5tR+ZRMplXk8ieRz8D7/8pGZIdVpTUth+vgCssu/i7UomVZMZ/0UttZgrH2xV3m/Wy/HYd2JiIkghEHUe4CDxCX6Nnr+fP3WtBTKv3yH1Kfvg/0PxvTTa/wKyb5IIqwMj4Qw+X4PZeah7u1+X16J0bLru0J+tl6ORzxLxEQQwiFCPUCvCLQ7nL2OVbcr7C2ZFvS0YH9uL/2dV90pUvY7CNPPrkVJTY3YDKdEaKwHJIFWhodKr3fU53uYzrgAXE6M5x7vKbvBMhp7+iBbL8cAERNBCIcI9QC9IrC5yd7vWCizkPq6vfR317qFZPosTBe7hSTUuu1dLj6rb/UvPIneWCfJjDN/9HpHfb6HUjgGZcHJGO+9gVG7kdoWB+1KKt8WlMVFQHwRMRGEMIhUD9ArAhNt1n7Hgg0C6x12Muq+pDzLwJqWgv7v1zEeuxsqDsC0+DqU1LSw6q5tcdDa2e1feBK8sY5nD32o+L4jf99DOUGF7Fz0lQ9SmpuWMBMnREwEIY54/eq++alC9rX7jBL0/7yB8ehdUD4D0+Jf9RKSUOsuzbOQm57qt6FK5sY60RnsHSkZmSinngObqsn433v9YmfeNSqxdkGKmASB3mHHseHTxPUPC8MbzyjB2F6P8chymLa/O/tv2tB6q9a0FGYU57pHO4keIxkm7LR38fJXzew6YB6MLcN45hGMru9Gjj1rVJo7Yu6CFDEJhvpajLY9iesfFoY1pkwrir0NnvgLTJmO6ZJfD1lI+pHoMZJhwodb29nV0c2H2/ZiOvMn7pXxr77Qc97fGpVYIWISDMWlKNk5CesfFoY3RvVn6H/9PZRMwnRp7/1IQiXgCCTBYyTDhdljsyjITGX22CyUqfvBzEMx/vkMxm73VGF/a1RihYhJEJgyrVgqDhD/sBBz/DXuvlNHjZqN6PcshVHFmK64CSV9iFu1BhiB+IuRyOZYsWekNY0TpuYz0uqOhZlOPx9cLox/rIyzZSImgpDY+GncvVNHt22sQb/rN5CVQ+clN7LBnjL0hj2EEUg4U5iFoeMr4srI0SjzjsH491qMnQ1xtUvERBASGT+Ne2mehRGdeyh7dBkoCqZf3EKtYQ3YsIcygghllta+mM8rGegr4soJZ4ApBeOlp+Jql4iJICQgXvcW0LNDn9fdldndwaS/3Yaytx3TFTejjCoasGGP1ghC8nnFh77vWskrQDnyeIz3/0XHt7Vxcz2KmAhCItLHvWWvq2VDYwf2TV+j//lWaNwGP1mC0dWJa9fOnkWLQL/GREYQ+xb+RFw59jRIs+B44Um/e6n4puuJFiImghBFwg5S93Fv1VlH026y4Fr9HHzzNaaLrkTJtLoF57MPeoTH3yhksBGErCFJfpTsXJTvnUTO5+9ja67v6Tj4S9cTLURMBCGKBOtianc4e4mOKdPqFpL6WvQOO2WFOUxZt4bsjZ+h/OhnKDPnfic4Mw7pEZ5QRyH2Lhcbqmuxt8oakkQmmJXtyjGnQKaVsnee7ek4+EvXEy1ETAQhigTbuG9usvcXHY+ry9haQ8aqx8j99B2UU34EB82jrXo9G5o62VsyjZSCkT1B81DjGLUtDtpyC6lTZA1JIuNd2V7b0Irx/pt+BUXJzEI55lT47EOMb74C/KfriRYiJoIQRYJt3CfarP1Fxzvy+OpzjLUvoiw4CeX4M6C+lrrmDnfjMsSgemmehZwcK2X7BZ/uXog9PSvbnbshMyvgKFJZcBJk5aDHYb94ERNBSACyLOZ+omPKtMKObbDqKZTZVSjqhSiKAsWllORnkj1m9JCD6jIjKznwrmzPOmw+SsFIv6NIe5eLDXsMur93Cmz4FKNmY0xtFDERhCgQSlA7UCJR49MPMB67ByoPRDn/ChST+7+rKdNKdnkl08cXiAgMMwZaB+SNz9XsfyRkWtH/+XcgdolqRUwEIRoEkRixZ6ZXXf9EosbXX6Df/wconeTeJXGYbLcrhI83Pjd+VB7KkSfAf9/H2L4lZolqRUwihEyvFHoRRFoSb0+yzjq6VyJRY0sN+t1LoaAQ02U3oqRnxMRkIbnxdVkqC07GSLOw+/mV2AvHxSRRrYhJpJAU3YIPwaQl8fYky0bn9SQSNXY2oC+/GdIzMP38N+5GQBBCRMnOYfdBR5H92b/ZtqMlJolqoz9fLIqoqnossBxIAR7UNO22uBnjWRMg0yuFYPH2JL0Ye3aj33EjOJ2Yrr7VHWgVhBDQO+w97VDm8f+H8sGrlH74T5g7M+r3TtqRiaqqKcA9wHFABXCWqqoV8bJHtjEVhoJub0e/82Zo3Y3pshtQisbH2yQhGfH1kIyw0TpjHin/WYve0hz1WyetmACzgU2apn2jaVoX8DSwMM42CUJAAqVWMbq7aLntWtj2LVzwcwzdJbE3ITx8YnW1LQ62zDkenN10/EOL+q2TWUyKgS0+n7d6jkUFwxHe4jAJzAte/KVWMXQX+gO30/3FOpTzf46SkyexNyFsfD0kpXkWUovHo8+Yg3PblsEvHiJJHTMZDFVVFwGLADRNw2azhVWP438fs/P315P//+7FPH5Cz3GXvR1n3SbMJZNIsWb5v3bDVgyTgtK+G8v4krDuPxBmszns7xVNxK7+zMp2J9ybaLOSZTFjGAZt9/6Ovf99n9yLlpB+/GlB/Zvypd3Ru85II+8xNBLJLhtQUgTGr5aRmpGJ0xndzMHJLCb1wDifz2M9x3rQNO1+4H7PR6OpqSmsGxk5+eBy0vzEA5guWtJzXN+43t2LbG3t2XOiL3rWCGhthawRmMK8/0DYbDbC/V7RROzyT7EFOtta6GwD/fm/YaxdhXK8Svrxp31nV+FY2Nvp/hmE9Y0dtDmctLS29grmR4p4P69ADFe77F0ualsclOZZQlqwaktNC9uuoqKioMols5vrI2CyqqplqqqmAWcCq6JxIyUrh4zvn4rx0TsYjdu+OxHEWgIJzAv+0F9bhfGKhjLvGJRTftiTNXinvSuklPWyV8nwIlAW6rC3OoggSSsmmqY5gUuBNUC1+5C2Plr3yzz5TEhJwVj9XM8xEYp9j2BSfYdTp+9/dP39tzBWPggzD0X50cUoisLmJjtNdgcvVjez094VdAJHya01vAjUeYjWbpqhkMxuLjRNewV4JRb3Ssm3oRx+NMY7r2Kc+AOUfFkDsC/iTfVd5+ygwlwLAdyXIdfpcFLbAhWN1RiPLIep+2H6yRIUk1sEJtqsfFprMC43FXuXzsHFMtIQ+tN3bZKX0jwLtS3EdYQa9MhEVdU7VFU9IJrGJDrKsf8HGBhrno+3KUKU6En1nZ8ZsQWoPSvdd9ei33sbFJdguuR6lNS0njJZFjPzy3IZm5vB/LJcGWkIIZEII9RQRiYpwBpVVXcCjwNPaJq2NTpmJSZKQSHKnPnu0ckJZ6DkjIi3SUKE8ab6hgJ3sLOxI+Rgp786K7qbcP3lt3Rlj0C/+AasGf17l4F6nYIQCr6r4GPpgg96ZKJp2uVAEXAtcABQrarqa6qq/lhV1cHnMO4jKMeeDs5udG2FrB3Zx/Hnhw4n0Gk07UC/80acKalUn301tbokbhSiSJzyBIYUgNc0zaVp2j80TTsLmAOMBB4BGlRVfVBV1agtGkwUlNHFMHU/WPcexqYN8TZHiCL+gp2hBjqNPS3od9wEXQ70y27CMrooLL+2LH4VgiaIWabRICQxUVU1R1XVC1VVfRN4G/gAmAeUA+3APyNvYgJy8tnQ3QUbvxOTRJiaJ0QWf35oX4EZ7J0bezvQl/8GWpowXXYjmWUTwvdrS1ZqIUjiNcs0lAD8M7gXBf4fcB9QpGnaIk3T/q1p2hbgl0BZdMxMLFImV8CM2fD2GozODiAxpuYJ4RNsz99XYAZ650Z3F/o9S6G+1r251aTyoRkYp96mkNx4OzztjuiufofQRibvA5M1TTtB07SVmqb1+h+kaZoOjIqodQmM6QQVOtpxvLqK9Y0dFFrNsngsmQmj5x9ozr/hcqHffzt89TnKeVeg7HdQ0HUGGu3ImiYhHLwdns1N0XePBj2bS9O024Mo0zE0cxIf35kSzJiNae0LdEyropG8Ic3ECTdNghAc3vfmypjlv0AY+9H4m31lGAbG3/4Cn76PcuZFmObMD8lO3zUpMrNLGCre9ScTbVY621qieq+kXQEfN3x6sKZTfoTJsZeSD14Z8ohE3GRRxvPenHWb/J6OVM/feO4xjHfXopz4A0wLTgr5ekmPIkQSb4cnGklA+yJiEio+vmtlbCmm2VUUvL+azI7WIVUrjUiU8bw3c8mkiFTnL8air3keY/WzKEcci3Ly2WHVmwiLzwQhHERMQsS3B2vvcrHpsFMxnC6Ml1cOqV5pRKKL970Fk9Y9KPrEWPR/v47xzMMoBx2OcvZPURTF72Uy60/YVxExGQK1LQ6as2y0HHSUe1V84/Z4myQEINxZLQFneRWXYqSaMTo7cf3nDYxH/wwVB6Bc+IuefFv+8LozaxpaZN2IsE8hYhImeoed8btqyDK6yTjFk1H4xSfibZYQgLBntQSY5WXKtKJYMmBzNTx6F0ycimnxr1DMqQNW53VnltgbZN2IsE8hYhIu9bVgb4OdDZCbj/K9hRgfvo1R83W8LRP64Cv8E20hBtgHWN9hdNhhzXNQVILpshtRLOn9yvR1a/W4M0sC1ysIyYiISbgUl1KnZNOWW0htiwPluNMgdwT6ygdx2dvFhZFI1Ndi7WyjwtEQ8qyWQLO8jM1fwgO3w8gxmH55K4rnfF+3WKBZerJuRNjXEDEJE1OmlbL9ppGTY6U0z4KSnomy8Iew+Ut4fVVAF4bkWIoDEVg97vvejG+/Qb/rN5CT5xaS7JzvCvZxi5XmWcgyuhm/q0beubBPI2ISJN5tVX1n4WQ4OylvrSHD6d6rWzlsAYwrg3dfw7Ck+2+8JMdSzInIKMDz3oz/fYh+502QnoFpyW9R8vJ7l+sjXNa0FCocDVg75Z0L+zYiJkGyucne313RRxgUUwom9ULY3YRS87X/xktyLCUnxaUYzm545hFQFEy/uBWloLBfMb/CJe9cGAaImATJRJu1/6LC4lLs6dlssIzuGbEo0/aHA+dgvPJ3jF2N/eoRX3niEIrLUXF0wnOPgdOJ6Re3uLciCBJ558JwQMQkSLIs5n6LCk2ZVr4tKKNdSe01YjH94CIA9Kfuj7mdQggE6XI0Wnah//HX0L4H0xU3o4wtjYV1gpBUiJgMEX9pUJSCkSgnnwWffYjx6ftB1yWro2NMEO4no6XZLSQtzZiuuAmlbHLUzJH3LyQzIiZDJFAaFGXByVBcgv7U/Ride4OqS5I9xpa+7qe+bi+jdbdbSHbvcgvJpIqo2hPs+5cZgUIiImISJRSzGdOPLobmJrpefCqoHmegZI/SY40RPm4vY49HSJp3Yrr8RpTJkRGSgYQg6GSfMiNQSEBETKKIMqkC5fCjSXnjJZzf1vTqcfprVAKNcmTEEj4h9eI9bi8jZwT6H2+AXY2YLr8JZcr0yN1vACEIOtmnzA4TEhARkzAJdrSgnHYuZGQy4Z8rKM3xWX0dQu9S0tMPAT/PeaDdDJVRxXDvMmhqwHTZDShTgxOSnjrrarG3trGhutb/v40ICIHMDhMSkejvmDIIqqr+ATgJ6AI2A+drmtbiOXcdcCHgAi7XNG2N5/ixwHIgBXhQ07TbYm13oB3x+u6YqGTlYDrzIjIf+hPKmy/Bcae5C4aws5+/Hf2895mVHf29nRMJ350ug2pM/TznQO/OPWvrBmhuxHTpDe5p3kHa4a2zzjoaOuhJs9P3vZkyrTC5MqTvLAjJQCKMTNYC0zVN2x/4GrgOQFXVCuBMoBI4FviLqqopqqqmAPcAxwEVwFmesjEl0GjBn0tKOeQImDUX48UnMLbWAEPvXUZzb+eEjtGEGC/w95x9353XLeXaWof+++s8wfbfoJTPCMkOb51lo/N60uyMtzglUC4MG+IuJpqmvappmrd7/T4w1vP3QuBpTdMcmqbVAJuA2Z6fTZqmfaNpWhfwtKdsTAnk3/Y7VVhRMP1wMViz0B/8E0Z395Dv771PyFlwgyChYzQDuImCFcFe766+FqP+W7jjBmhvw/SL36BMcY8cBox/9LHDN7VOT/2NWyRQLgwb4u7m6sMFgHfLwmLc4uJlq+cYwJY+xw/xV5mqqouARQCapmGz2cI2zGw2B3W9DSgp8nfChuOy62lZehXpa58n+8eLw7bF9z5msxmnJfzv5Y9Z2e4Rz0SbNey9o4N9XqFjg/El/Y62O5y89eUOcizptBgWSmy5Qdnl2F1I632/A93FiFvvJnXi1O/ObdiKYVJQ2ndj6XfP3nb4K+vKmIWzbhPmkkmD7vAYvec1NMSu0BjOdsVETFRVfQ0Y7efU9Zqmvegpcz3gBCK2w5SmafcD3mXoRlNTU9h12Ww2wrm+l2+9dCrKvGPoeOEJOidPj8h003DtGoxiC3S2tdDZFt71odqld9ix19VSZx1N2ei8kLcvXt/YgdLlYHu7wdScXJqa/I/+fO0yvt2MfufNAJiWLKU1twB8bNazRkBrK2SNwDTIdwlYtnAs7O10/wxAtN7jUBG7QmNftKuoyF/vuD8xERNN07430HlVVc8DTgQWaJpmeA7XA+N8io31HGOA44mHr299ciWKegFG9WfoK+7AdNNylPTMwWoYHtTXUtfcQdueBmrTM/oFrgfD61Y8qNgSlBAZX32Ofs9SyMjE9PNbUMaM7Td5IpRguQTWheFO3GMmnplZVwMna5rW4XNqFXCmqqoWVVXLgMnAh8BHwGRVVctUVU3DHaRfFWu7+xLMXuF6hx0lPRPTBb+AXY0YT0rurh6KSynJzyR7zOiwpkAHvUYDMNa95x6R5IyAH/wEI3cEkOCxIkFIcOIuJsDdQDawVlXVT1VVvQ9A07T1gAZsAFYDl2ia5vIE6y8F1gDV7qLa+viY7sMge4V3OLp71h4YxSUw50iM995Af3dtXMyNFt4geLsjtCnLpkwr2eWVTB9fEFAQIjHLrGPtKvT7fgfjJ4B6Popi6jcjq9BqTtzZbIKQoCiGYQxeat/A2LZtW9gXD+ZzHGj9g95hZ0N1LW25heTkWClvrcFoaoSXNfeahutuRxlXFhW7/NHXnRMJvHV2drvo1g3GFRZQbBn6rDVf1jd20OZwku3J4BwKhmFgrH4W47nHoPJATBdfh+Fy+X1nQ7lPuOyLvvZoInaFRgRiJspg5RJhZLJPMNC6kb5b/FJcCnvtMO8YMKei37MU1xef9HORRSuhXzTcOd46DYjalOVwMwHYO7toeuQ+jOceI73qGEyX/hrFkh7wnUnGAUEIHRGTGOHr0zdlWlHmHAmpqXDUSdC8E154AmPLN71dOVFK6BeosRxMvAZyM3nrLB+ZSWVhZthTigcilLiIF8PhoPve35H/n3/SPPc4cq64EcWcGvH7CMJwR8QkAoTjy/cKilJ5AJx8NtRtgnXv9R41RCmhX8DGchDxGmhEE04D7PvcorHq3mjdjX77r8iq/pgdJ/yYzB8uQjH1/yef0Cv+BSFJSLRFi0lJoFxPg9EznXRyJXrrbow3/sHEojI2lx9GaZ4FU1pKbKebDpIvrDTPQm0LEXP/+D43IKxnGAij/lv0P98Cba2kLL6OogPmBGVHoHtHI84kCPsSMjKJAJHwsSs/+AmUzyD1qXuZ1rCejLovIxorCSb+Mli+sEi7f3yfm/fvSOSzMtb9B33ZldC5Fy67AWUAIelrRyBk2rAgDIyISQQIt5H1da8oZjOmn14DI0fDX3+PUbspsrGSEOMvsdjNz/e5RSKflaHr6C8+gX7vbZBvg9POQ0kZ/J0E8/4kKC8IAyNiEiP8+eX79nYVaxamX/wGMrPgHysx0tIHrSNoQo2/xGs3Pz92BvO9jQ47+j1LMf6xEuWwBbBkKcroIr/fNxyhlKC8IAyMiEmM8Ocm8ZthOH8kpl/eCiYF7rsNo6V5wDqCwd7lorpdYW/JtOBT3sdpNz9/rrbaFgd79tip+dy/689o2Op2a61fh3L2T1HOvZyU3BGBXXay7a0gRBwRkxjhTzgC9XaVMWMxXX4ztLWi/+mGHkEJd4V2OCIU7n4rPfuD2NtDum4gSvMsZLc2UmL0FwDj43dxLb0SZ1sbjstuxnTkCSjKIOurwhBKmfElCAMjYhIjfHdLHKhR8jbGxqgiTJfdAM1N6L+/FqNpR08djXZnSOLgT8ii1jh6ev3Ouk0Bi4TqZrKmpVBRXoo19zsBMBwO9MfuRv/r7+kcWcR/z7uZN83jgvo+vkIZrC0SgBeEgRExiTGDNko+Lhhl6n6YfnkL2NvRf38dHVu+ZX1jB4VWc0jBYH8joKG4zAYUIU+v31wyKXAlYbiZfAXA2FqLvvSXGO+uRTnuNJSrlrErIx9rming9wmYMyxIW0INwMtIRhhuiJjEmNI8C2aTQme3y29DY+QXYuxuwsgvBECZMBXTVUvB5STlj7/CWbeZRrszrGCwby883NlJg4mQt9H33QzK3uXi4/o2Pqlvc3/nMOMxhmGgv/kK+tIl0NGO6ec3Y/q/c7FmWJhflstIa1rA7xNom2N74Tg26NnYC8f5vc5LqAF4GckIww0RkxhjTUshIzWFbt3w29AozY0oI2wozY3fHRtbhumqZaSkpVH+2FLKmjaHd3OfXniGs5Pxu2qoaWjxK2qBetbhiFBti4NvdneyeXcntS0O9prTqc4tY685fcDrfG0wmnei3/UbjCfvg2n7YbpxOUrFgT1lB2vsfbc59hXVbx1m2gvH8a0jsut3ZSqxMNwQMYkDAzY0AXrtyuhiUq75HSm5I7DcdRP6O6+GfmPfur2bUW1v8CtqgXrWAzXageIPpXkWJoxIZ+KIdM8q+uB67bUtDto6u9n92j/Rb74Mvl6PctYiTJfdiJKT16vsYG4lr91ZFnMvUY1Woy9TiYXhhqRTiQO+wfi+DLRjn1IwEtN1f0C//w8Yj92NvrUW4+KrgeDSffjWrReXUuJ0b5PrryENK3WKb/zBZ390a1oKBxVnh1x3maMJ55N/xbr5c5gyHdN5l6OM7L37szf1f41lNO1KanDpWHzSxgz0LgRBCB4RkyRDsWZhuvxGjGcfwVj7Iru3fYtx/hXUOq0h5bbq2YwqwPmwGtlBcnv1rds7mugrgIbDgb7qSSyvr8KSmgZHnwzH/wAlK7tfXcY31VD3DeOL97KlqCIo8YvlFruS00sYLoibKwlRUlIwqRei/GQJztpN6DdfxsSPXiHHpMfVRx/q2pSahhb21LjjNuBOh+J69zX06xfBq8/DxHI49jSUydNRtn/rvxJDAQWsKUpCupUkEC8MF2RkkqToHXbIt5H323vY/eg9pL7wGFPf/ifKqT/CmH2E31Tr4RKt3nWJvYE6Zwcl7dsxvqhDf/4x+PYbyB8Jp/0YZWK5e1Zbc2PA0Y4ycRqkp4c8M2ygnTEjSaQzLQtCoiJikqx44hNGbi4pl9+IUf0Z+jMPYzx0B8baVZhOPw+lfEZEbhVuin1f/DXe1pJSyhvfhpf+hb5pA9hGwTmXwOgilLETvmvkC0YGrDdsl5VvfCeKLi+JyQjDBXFzJSne9SimkWMAUMpnYLr+TygX/gLa96D/6Qa6lv+GjV98PeSFcwPNeBpoBXmvcz6Nt6HrGJ99hHHfbfC3e2FHPcpZi+C621HGjO0tJFEi2PUlgiAEh4xMkhSluRFG2GjbVs/6VFuPC0qZcyTGzLkYr/8DXtYoW38VrbOOJPMHP0bJyw/rXgP2rgfq4fueKy7FqPkKvvka45E/Q+M2yCtA+b9zUY48HiU9A9fH/8HY8AlUzIKD5oZla7D4ri+J4fZjgrDPImISByLiry8uxfimms279rIn004t37mgOjBTO+s4Rs08AuNljREfrEX//D8ox5zi/kmPjNtF77BjOPaCORXFX8yiuBRjaw10d6M8/xh88DbstUPZFJSLrkSZORfF7PNPcFcDOLrcv6OMxDIEIbKImMSDCPjrTZlWdEsGE3QX9tZGSsdP6znnjXFgyaTygsUYJ56C/txjGC89jbHmeZQZs1Fmz4PKWSipqUP6Hkq3E7Kye4mi0d2F45P34J218NlH0LILIy0N5YBDUY46wR0098dB8yA1DWYcEr5NQSKxDEGILCIm8SCI9Rj2Lhc1DS2U2BuwlgQYwRSXktW+m4pxI9z7xXvo2+tWCotI+dm1GDVfY7z7Gsa6f2N89A6kpUHpZJSJ01AmlsOEaSjZOSF9D2NrDWTlYnyxDmPzlxjffAmbv6TF0QmWdKic6RavA+egZAzceKcUjISjTgz+/oIgJAwJIyaqqi4BbgdGaprWpKqqAiwHjgc6gPM0TVvnKXsu8GvPpb/VNO3ReNgcLsHMQKptcdC23T11tsJc67e8KdOKZXwJpqamXscD9bqVsikoZVNwnfJD+M/r0NQIdZswXn0Rw/Wsu1BBIRQUouTbIGeE+5ihg66DYbh/6zpG+x5o3A47t4OjEx1AMUFxCcrco8g9/HvsGTMeJVvEVJgAAA9RSURBVDUtjCcUPrGa8isIQm8SQkxUVR0HHAP4rkw7Dpjs+TkEuBc4RFXVfOAm4CDAAD5RVXWVpmm7Y2t1dCnNs1AzZjQl9oZeI5hINJZKw1YoLIIJUzH98GcYXQ6o24z+5Wc4v9pAV0cHGV+tR2lvBQVQUtw7Pyqm735nZkHhGJSp0zHyCsCcAjMPIyXfBoDFZkPpI3IxIUZTfgVB6E1CiAlwB3A18KLPsYXAY5qmGcD7qqrmqao6BpgPrNU0rRlAVdW1wLHAU7E1ObpY01KYPr4AKOh9IhKNZR83m5JmgckVKBhszJ9EmzmTnLKyoGMK+sb1bpt27YB8G/YuF/X1reQprkEXOUZ8QWSQKV0EQYgscV9noqrqQqBe07TP+pwqBrb4fN7qORbo+PDAk/nXXjjO/2ZPQRAw7UlxKSX5mWSP8Z/8cTCbvA14bYuD1s7uoFKIRDrdSLjbDQuCMDRiMjJRVfU1YLSfU9cDv8Lt4orGfRcBiwA0TcNms4Vdl9lsHtL1kcMG40v4rL4VhW5qd3cyffTQ7XLZ23E2bmXEtGkU7dyOObe41wZXwdjkZVa2k9rdnZSOzXOnfO9Du2eTqok2K7Oy83r+9lfWH77XB3uNl8R5j70Ru0JD7AqNWNgVEzHRNO17/o6rqrofUAZ8pqoqwFhgnaqqs4F6wHd58ljPsXrcri7f428FuO/9wP2ej0bTEHz4NpuNoVwfafIUFy0OB6Vj8yJil9dVZXzxKcoIG0ZjA4olI+zYzPTR7ufV2db/3PrGDtocTlpaW6kszKTYAp1tLX7L+qPv9aGQaO/Ri9gVGmJXaAzFrqKioqDKxTVmomna50Ch97OqqrXAQZ7ZXKuAS1VVfRp3AL5V07TtqqquAf6fqqqeqUYcA1wXY9OjSjBBdt/NnoJthAfEG2sYP9GdWLGzM2qB7KEuGJQFh4KQeMQ9ZjIArwDfAJuAB4DFAJ7A+63AR56fW7zB+GQmUB6rSDDYLoTgs3d7wUhMkyvdCwsD7NPeNx/XQPm5/DHUXQhlF0NBSDwSZTYXAJqmlfr8bQCXBCi3AlgRI7NiQ588VqHMSBpsRlQwWX/bWvdQt2kLJZPGkZ2bM/BamL4zykKYYSabRQnCvkkij0yGFz4zokKdkTTYjKhg9jmv27SFtva91G3a4vd8r9FN333qA+xbH46tgQhmdCUIQvwQMUkQhjKldTCxCMYtVDJpHNlZGZRM8p+S3VcE+trqz3aXvd2v6ysYYRvs/oIgJB4iJvsAwcYQBurdZ+fmMH1WJdm5/nNzhSoCzrpNfuM+4cY7whUhQRBig4jJMGIovfuBRMBfAN5cMilo19dQ7y8IQvwRMUlyArmT/BG13r2f2Wcp1qxhsxI91NlsgrAvImKS5ARyJ/kjar37AQLwwTa0Sd0gR3gqtyAkIyImSU6k3UnhMODkgWAb2mRukEOYzSYI+yoiJklOwruTgm1ok7hBluSSgiBiIvghki6nYBvaeDTISe1aE4QEQ8QkCYn6Ar4gXU5J3xgns2tNEBIMEZMkJOoL+Py4nPwKWLI3xknsWhOEREPEJAnxTvEdb3Hi2PBpxEcG/lxOfgUsyRtjiXUIQuQQMUlCeqb4Nm7BaNsTk5GBvzUq0WyMk96FJgjDDBGTZKa4FCU7JyYjg5ivQE92F5ogDDNETJIYU6YVS8UB+6abJsldaIIw3BAxGaYkekp3iWcIQnIhYjJMkZTugiBEEhGTYcpgSR+9AfC21j0JPYIRBCExEDEZRvi6tgYNqHsC4HWbtsgIRhCEQRExGUaE5NryBMBLJo2TTakEQRgUc7wNEGJHaZ6F2haCEgZTphUmV5INVEbfNEEQkhwZmQwjgl0rIgsGBUEIFREToT+yYFAQhBARMRH6IwsGBUEIEYmZCP3wxksEQRCCJSHERFXVy4BLABfwsqZpV3uOXwdc6Dl+uaZpazzHjwWWAynAg5qm3RYXwwVBEAQgAdxcqqoeCSwEZmiaVgnc7jleAZyJezLRscBfVFVNUVU1BbgHOA6oAM7ylBVI/DQpgiDsm8RdTICLgds0TXMAaJrW6Dm+EHha0zSHpmk1wCZgtudnk6Zp32ia1gU87SkrIGlSBEGID4ng5poCzFNVdSnQCVypadpHQDHwvk+5rZ5jAFv6HD/EX8Wqqi4CFgFomobNZgvbSLPZPKTro0Vfu2ZlO9ncZGeizUqWJX6vN1meV6IgdoWG2BUasbArJq2NqqqvAaP9nLreY0M+MAc4GNBUVZ0QiftqmnY/cL/no9HU1BR2XTabjaFcHy382VVsgc62Fjrb4mQUyfW8EgGxKzTErtAYil1FRUVBlYuJmGia9r1A51RVvRh4TtM0A/hQVVUdsAH1wDifomM9xxjguCAIghAHEsHN9QJwJPCmqqpTgDSgCVgFPKmq6p+AImAy8CGgAJNVVS3DLSJnAmfHw3BBEATBTSKIyQpghaqqXwBdwLmeUcp6VVU1YAPgBC7RNM0FoKrqpcAa3FODV2iatj4+pgv2Lhe1LQ5K8yyx29JXEISEQzEMI942xApj27ZtYV+8L/pCI8H6xg7aHE6yLWYqCzMTxq5AiF2hIXaFxr5olydmogxWLhGmBgtJzGCbbAmCMDxIBDeXkMR4MxELgjC8kZGJIAiCMGRETARBEIQhI2IiCIIgDBkRE0EQBGHIiJgIgiAIQ0bERBAEQRgyIiaCIAjCkBExEQRBEIbMsEqnEm8DBEH4/+3df+xVdR3H8SfT4g9raVoKYmlFrrSNkmGt2igQgZlkyxff1kwNLFyu9Wsl2sqBf2BlxrQiKZY2E18rSZaYYuXqH0whyRRbmLgggiZKNpvt2+iPz/myy5d7v1w433OvxOux3cH5nM/3nvf5nM+573M+59x74jCVn1NpMabOS9L6uu/RxCtxJa7E9dJ5/R/HdUBHUjKJiIiGJJlERERtSSbdu/nAVfoicR2cxHVwEtfBOWLjOpIuwEdERENyZhIREbXleSYtJF0IXAO8BZhi++GWeQuBecB/gU/bvrfN358GrASOB9YDF9n+zyjHeAdwejV5LPCc7Ult6m0Bnq/iHbQ9eTTjaLO8a4DLgH9URVfZXtOm3kxgKeWRy9+3vaThuL4OfIDySOgngUttP9em3hZ60F4HWn9JY4FbgbOAZ4C5trc0EUvLMk+plnki5Rb6m20vHVZnKnAX8FRVdKftRU3GVS13CyNsF0ljKO05G3gBuMT2hoZjOh24o6XoDcBXbH+rpc5UetBeklYA5wE7bZ9Zlb26iu9UYAsg28+2+duLgS9Xk9favqVOLEkm+/oj8CHge62Fkt4KDABnAOOB+yW9eeiZ9C2uA26wvVLSMkry+e5oBmh7bktc1wO7R6j+Ptu9fIboDba/0WmmpKOAbwPnAFuBhySttv14gzGtBRbaHpR0HbAQ+FKHuo22V5frPw941vabJA1Q+tTc/d9tVA0Cn7e9QdIrgfWS1rbZLr+1fV7DsbQz0naZBUysXmdT9rezmwzG9p+ASbB3m24DVrWp2ov2+iFwE+VgYMiVwC9tL5F0ZTW9T5+vEs5XgcmUA4j1VV/cL+l0K8NcLWxvqjrKcHOAlbZftP0UsBmY0lqhOkJ6P/CTqugW4INNxVotT8DtTS2jAVOAzbb/Up2xraS0bWNs32d7sJpcB0xocnkH0M36z6H0HSh9aVq1rRtje/vQ0bzt54FNwMlNLnMUzQFutb3H9jrgWEnjerj8acCTtp/u4TL3sv0bYNew4tY+1Olz6Fxgre1dVQJZC8ysE0uSSXdOBv7aMr2V/Xe24ylDToMj1BlN7wV22P5zh/l7gPskrZf0iQbjaHWFpD9IWiHpuDbzu2nHJn0cuKfDvF60Vzfrv7dO1Zd2U/pWT0g6FXg78GCb2e+StFHSPZLO6FFIB9ou/e5TA3Q+oOtHewGcaHt79f+/U4Yvhxv1djvihrkk3Q+c1GbW1bbv6nU87XQZ40cY+azkPba3SXotsFbSE9VRTCNxUYYXFlN2/sXA9ZQP78Z1016SrqYM59zW4W1Gvb0ON5JeAfwU+Iztfw6bvQF4ve1/SZoN/IwytNS0l+x2kfRy4HzK0Olw/WqvfdjeI6knt+weccnE9vRD+LNtwCkt0xOqslbPUE6xj66OKNvVGZUYJR1NubZz1gjvsa36d6ekVZQhllo7YbdtJ2k58PM2s7ppx1GPS9IllIuU02y33bGaaK82uln/oTpbq+38KkrfapSkl1ESyW227xw+vzW52F4j6TuSTmj6mlwX26WRPtWlWcAG2zuGz+hXe1V2SBpne3s15LezTZ1twNSW6QnAA3UWmmGu7qwGBiSNre7Ymgj8rrVC9SH1a+DDVdHFlLs5mjAdeML21nYzJR1TXUhF0jHADMrNBY0ZNk59QYflPQRMlHRadVQ3QGnbJuOaCXwRON/2Cx3q9Kq9uln/1ZS+A6Uv/apTAhwt1TWZHwCbbH+zQ52Thq7dSJpC+exoNMl1uV1WAx+TNEbSO4HdLUM8Tes4OtCP9mrR2oc6fQ7dC8yQdFw1JD2jKjtkR9yZyUgkXQDcCLwGuFvSI7bPtf2YJAOPU4ZKPjV0J5ekNcB823+j3DGxUtK1wO8pO2gT9hunlTSecqvpbMoY6SpJULbxj23/oqFYhnxN0iTKMNcW4JPD46ruqLqC0mmPAlbYfqzhuG4CxlKGSADW2V7Qj/bqtP6SFgEP215N6TM/krSZcmF1YLTjaOPdwEXAo5IeqcquAl5Xxb2MktgulzQI/BsYaDrJ0WG7SFrQEtcaym3Bmym3Bl/acEzA3uR2DlU/r8pa4+pJe0m6nXKGcYKkrZQ7tJYAljQPeJpyow6SJgMLbM+3vUvSYsoBDsAi28Mv5B+UfAM+IiJqyzBXRETUlmQSERG1JZlERERtSSYREVFbkklERNSWZBIREbUlmURERG1JJhERUVu+AR/RJ5LeSPkG8vTqWSLjgY3AhbYf6GtwEQcp34CP6CNJlwGfpTykaBXwqO0v9DeqiIOXYa6IPrK9nPK7Ug8C4yg/5x9x2Ekyiei/5cCZwI22X+x3MBGHIsNcEX1UPZBqI+XxBbOAt9X99daIfsiZSUR/LaX8/Px84G5gWZ/jiTgkSSYRfSJpDjATuLwq+hzwDkkf7V9UEYcmw1wREVFbzkwiIqK2JJOIiKgtySQiImpLMomIiNqSTCIiorYkk4iIqC3JJCIiaksyiYiI2pJMIiKitv8B8ZTPZF3FpC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "def gen_batch(n, w, beta):\n",
    "    d = len(w)\n",
    "    X = np.random.uniform(-10, 10, (n, 1))\n",
    "    X = np.sort(X, axis=0)\n",
    "    X = np.hstack([X ** i for i in range(d)])\n",
    "    t = X.dot(w) + np.random.normal(size=n) / beta ** 0.5\n",
    "    return X, t\n",
    "\n",
    "n = 200\n",
    "d = 21\n",
    "w_true = np.zeros(d) # lot's of ittelevant features\n",
    "w_true[1] = 100 # one coef is super large\n",
    "w_true[3] = -1 # one is small\n",
    "beta_true = 1e-4\n",
    "\n",
    "X_train, t_train = gen_batch(n, w_true, beta_true)\n",
    "X_test, t_test = gen_batch(n, w_true, beta_true)\n",
    "\n",
    "plt.scatter(X_train[:, 1], t_train, s=3, label='Train data', alpha=0.3)\n",
    "plt.scatter(X_test[:, 1], t_test, s=3, label='Test data', alpha=0.3)\n",
    "plt.plot(X_train[:, 1], X_train.dot(w_true), label='Ground truth')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend(ncol=3, loc=9, bbox_to_anchor=(0.5, 1.15))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "local_scale_log__      -24.040001\n",
       "global_scale_log__      -1.140000\n",
       "weight                 -19.299999\n",
       "sigma_log__             -5.060000\n",
       "obs                  -2015.530029\n",
       "Name: Log-probability of test_point, dtype: float64"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_storage = theano.shared(X_train.astype(dtype='float32'))\n",
    "t_storage = theano.shared(t_train.astype(dtype='float32'))\n",
    "with pm.Model() as model:\n",
    "    # Horseshoe prior\n",
    "    local_shrincage = pm.HalfCauchy('local_scale', 1., shape=(21,))\n",
    "    global_shrincage = pm.HalfCauchy('global_scale', 1., shape=())\n",
    "    scale = local_shrincage * global_shrincage\n",
    "    weight = pm.Normal('weight', mu=0., sd=scale, shape=(21, ))\n",
    "    sigma = pm.HalfCauchy('sigma', 1., shape=(), testval=100)\n",
    "    obs = pm.Normal('obs', x_storage.dot(weight), sigma, observed=t_storage)\n",
    "model.check_test_point()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, weight, global_scale, local_scale]\n",
      "Sampling 2 chains:  10%|█         | 206/2000 [00:05<00:48, 37.05draws/s]\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "Chain 1 failed.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRemoteTraceback\u001b[0m                           Traceback (most recent call last)",
      "\u001b[0;31mRemoteTraceback\u001b[0m: \n\"\"\"\nTraceback (most recent call last):\n  File \"/Users/ferres/dev/pymc3/pymc3/parallel_sampling.py\", line 82, in run\n    self._start_loop()\n  File \"/Users/ferres/dev/pymc3/pymc3/parallel_sampling.py\", line 123, in _start_loop\n    point, stats = self._compute_point()\n  File \"/Users/ferres/dev/pymc3/pymc3/parallel_sampling.py\", line 154, in _compute_point\n    point, stats = self._step_method.step(self._point)\n  File \"/Users/ferres/dev/pymc3/pymc3/step_methods/arraystep.py\", line 247, in step\n    apoint, stats = self.astep(array)\n  File \"/Users/ferres/dev/pymc3/pymc3/step_methods/hmc/base_hmc.py\", line 135, in astep\n    self.potential.raise_ok(self._logp_dlogp_func._ordering.vmap)\n  File \"/Users/ferres/dev/pymc3/pymc3/step_methods/hmc/quadpotential.py\", line 231, in raise_ok\n    raise ValueError('\\n'.join(errmsg))\nValueError: Mass matrix contains zeros on the diagonal. \nThe derivative of RV `sigma_log__`.ravel()[0] is zero.\nThe derivative of RV `weight`.ravel()[1] is zero.\nThe derivative of RV `weight`.ravel()[3] is zero.\n\"\"\"",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;31mValueError\u001b[0m: Mass matrix contains zeros on the diagonal. \nThe derivative of RV `sigma_log__`.ravel()[0] is zero.\nThe derivative of RV `weight`.ravel()[1] is zero.\nThe derivative of RV `weight`.ravel()[3] is zero.",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-47-e5af4f687374>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtest_point\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjitter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/dev/pymc3/pymc3/sampling.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, compute_convergence_checks, use_mmap, **kwargs)\u001b[0m\n\u001b[1;32m    438\u001b[0m             \u001b[0m_print_step_hierarchy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    439\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 440\u001b[0;31m                 \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_mp_sample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0msample_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    441\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPickleError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    442\u001b[0m                 \u001b[0m_log\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarning\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Could not pickle model, sampling singlethreaded.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/dev/pymc3/pymc3/sampling.py\u001b[0m in \u001b[0;36m_mp_sample\u001b[0;34m(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, use_mmap, **kwargs)\u001b[0m\n\u001b[1;32m    989\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    990\u001b[0m                 \u001b[0;32mwith\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 991\u001b[0;31m                     \u001b[0;32mfor\u001b[0m \u001b[0mdraw\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    992\u001b[0m                         \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtraces\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchain\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    993\u001b[0m                         if (trace.supports_sampler_stats\n",
      "\u001b[0;32m~/dev/pymc3/pymc3/parallel_sampling.py\u001b[0m in \u001b[0;36m__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    343\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m         \u001b[0;32mwhile\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m             \u001b[0mdraw\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mProcessAdapter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrecv_draw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    346\u001b[0m             \u001b[0mproc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_last\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuning\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstats\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwarns\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_progress\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/dev/pymc3/pymc3/parallel_sampling.py\u001b[0m in \u001b[0;36mrecv_draw\u001b[0;34m(processes, timeout)\u001b[0m\n\u001b[1;32m    247\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    248\u001b[0m                 \u001b[0merror\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Chain %s failed.\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mproc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 249\u001b[0;31m             \u001b[0msix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mraise_from\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merror\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mold_error\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    250\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"writing_done\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    251\u001b[0m             \u001b[0mproc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_readable\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.pyenv/versions/pymc4/lib/python3.6/site-packages/six.py\u001b[0m in \u001b[0;36mraise_from\u001b[0;34m(value, from_value)\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Chain 1 failed."
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    trace = pm.sample(start=model.test_point, jitter=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc4",
   "language": "python",
   "name": "pymc4"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}