Skip to content

Instantly share code, notes, and snippets.

@junpenglao
Last active April 25, 2017 05:34
Show Gist options
  • Save junpenglao/e2b272b877068eba7525d1d21a6e618d to your computer and use it in GitHub Desktop.
Save junpenglao/e2b272b877068eba7525d1d21a6e618d to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline\n",
"import numpy as np\n",
"import pymc3 as pm\n",
"import theano.tensor as tt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A simplistic PyMC3 port of Bob Carpenter's [The Impact of Reparameterization on Point Estimates](http://mc-stan.org/documentation/case-studies/mle-params.html) in STAN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Model 1: Chance-of-Success Parameterization"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using advi...\n",
" 0%| | 0/200000 [00:00<?, ?it/s]/usr/local/lib/python3.5/site-packages/numpy/lib/function_base.py:3858: RuntimeWarning: Invalid value encountered in median\n",
" r = func(a, **kwargs)\n",
"Average ELBO = -7.7679: 100%|██████████| 200000/200000 [00:13<00:00, 14483.89it/s]\n",
"Finished [100%]: Average ELBO = -7.7705\n",
" 97%|█████████▋| 2912/3000.0 [00:02<00:00, 1215.66it/s]/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/hmc/nuts.py:244: UserWarning: Chain 0 contains diverging samples after tuning. If increasing `target_accept` doesn't help, try to reparameterize.\n",
" \"try to reparameterize.\" % chain)\n",
"100%|██████████| 3000/3000.0 [00:02<00:00, 1114.93it/s]\n",
"/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/hmc/nuts.py:244: UserWarning: Chain 1 contains diverging samples after tuning. If increasing `target_accept` doesn't help, try to reparameterize.\n",
" \"try to reparameterize.\" % chain)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAACsCAYAAABmdA06AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW9x/HvM5NJZsi+sROigiA7QiWKFRBcWVxpxb24\n1HtrW20vSKkFtLWiVluKXKXuiBXBXpGlgrhQKgVtcWNRkSUQ9gSSQAiT9bl/TBwzZJvgJJPl8369\nzouZc55z5nfOaybkm+c5zxhrrQAAAAAA350j3AUAAAAAQEtBwAIAAACAECFgAQAAAECIELAAAAAA\nIEQIWAAAAAAQIgQsAAAAAAgRAhZwEmNMujHGGmMiwl0LAAAAmhcCFiDJGJNpjBkVguMQzgAAAFox\nAhYAAAAAhAgBC62eMeZlSWmSlhpjCiT9oGLTDcaY3caYHGPMryu1dxhjphhjthtjDhtjFhpjkio2\nr6n4N88YU2CMOdcYc4Yx5r2KtjnGmFeMMQmNd4YAAABoLAQstHrW2psk7ZY01lobI2lhxabzJfWQ\nNFLSNGPMWRXrfyrpSknDJHWUlCtpTsW2Cyr+TbDWxlhr10kykh6uaHuWpC6SZjTkOQEAACA8CFhA\nzR6w1p6w1n4m6TNJ/SvW3yXp19baPdbaIvnC0rU13Xdlrd1mrV1lrS2y1mZLekK+cAYAAIAWhhvx\ngZodqPS4UFJMxeOukt4wxpRX2l4mqV11BzHGtJM0S9L3JcXK94eN3JBXCwAAgLCjBwvwsfVomyXp\nMmttQqXFba3dW8Nxfl+xvq+1Nk7SjfINGwQAAEALQ8ACfA5KOj3Itk9LesgY01WSjDGpxpgrKrZl\nSyo/6Vixkgok5RtjOkmaFJqSAQAA0NQQsACfhyXdb4zJk3RtHW1nSVoi6W1jzDFJ6yUNkSRrbaGk\nhyStNcbkGWMyJD0g6WxJ+ZKWS/q/hjkFAAAAhJuxtj4jowAAAAAANaEHCwAAAABChIAFAAAAACFC\nwAIAAACAECFgAQAAAECI1OeLhpkNAwAQLL7rDQDQKtGDBQAAAAAhQsACAAAAgBAhYAEAAABAiBCw\nAAAAACBECFgAAAAAECIELDQr3pKyBm0PAAAAfBfG2qBnX2eadjQJ6VOWB902c+boBqwEQC2Yph0A\n0CrRgwUAAAAAIULAAgAAAIAQIWABAAAAQIgQsNCiMSkGAAAAGlNEuAsAGpLb5WRSDAAAADQaerAA\nAAAAIEQIWGiVyrwFOvzWn5X15+u1+4lrdHDBr1WcnVnv4yxYsEDGGHXu3LnKtpdeeknXXHONunbt\nKmOMbr311mqPYYypcZk5c2a9awIAAED4MEQQrY61VtmvP6jS/INKGvVjOdwxyl+/SAdfnao9D1xR\nbViqTl5enu655x61b9++2u3z589Xdna2LrroIi1atKjG46xbt67Kujlz5mj+/PkaN25ccCcFAACA\nJoGAhVbnxLYPVbR3i9pd93u5u/aTJEV16qm9T9+mRx99VH/+85+DOs7kyZPVv39/dejQQe+8806V\n7StXrpTD4eskXrFiRY3HycjIqLLuxhtv1ODBg9WrV6+gagEAAEDTwBBB+M2YMUPGGH355Ze65JJL\nFB0drbS0NL3wwguSpJdfflk9e/ZUTEyMRowYoe3btwfs/5e//EX9+/eX2+1WSkqKbrvtNh05ciSg\nzZNPPqlzzz1XSUlJSkhIUEZGhpYvD5yEIjMzU8YYzZ07V9OmTVOHDh2UkJCgsWPHas+ePd/5PAu/\n/lDOmCR/uJIkR1S0PN3O0ZtvvhnUMdauXav58+drzpw5Nbb5JlzV1wcffKDt27frlltuOaX9AQAA\nED4ELFQxfvx4jR49WosXL9agQYM0ceJETZ06VU899ZRmzpypF154QV999ZWuv/56/z5TpkzRT37y\nE40aNUpLlizRY489phUrVuiyyy5TWdm3U59nZmbq9ttv16JFi/Taa69p8ODBGjNmTLU9PA8//LC2\nbdum559/XrNmzdK6det04403BrSxtly2vKzGpbS01PfYWv8+JTm75UrpWuX1XClp2r17twoKCmq9\nPiUlJbrzzjs1adIkdevWLejrGqyXXnpJkZGRmjBhQsiPDQAAgIbFEEFUMWnSJN18882SpMGDB2vp\n0qWaO3eudu7cqbi4OEnS/v379fOf/1y7du2StVaPPfaYpk+frmnTpvmPc+aZZ+r888/X0qVLdeWV\nV0qS/vCHP/i3l5eXa+TIkdq6daueeuopXXrppQF1pKen669//av/eXZ2tiZNmqROvScqIjZZknT4\n77N0fNO7NZ6L6zHfv8mX36OYvqN8r+s9poj4tlXaOtyxkqTc3FzFxMTUeMxHHnlERUVF+tWvflVj\nm1Pl9Xq1aNEijR49WsnJySE/PgAAABoWAQtVXHbZZf7HiYmJatu2rQYOHOgPV5LUs2dPSVJWVpa+\n+OILlZeX64YbblBpaam/zZAhQxQbG6s1a9b4A9aGDRs0ffp0/fvf/1Z2dra/Z6lHjx5V6rj88ssD\nnvft21eSVHY02x+wEs6/XrFnj6nxXJb+9HyNnf2BIhLa1esa1GTbtm166KGH9MYbb8jtdp/SFxm7\nXc4aty9evFj5+fk1zjgIAACApo2AhSoSExMDnkdGRla7TvL1uBw6dEiSahwud/jwYUm+MDZy5Ej1\n6tVLs2fPVlpamiIiIvSb3/xGX3zxRZX9kpKSAp5HRUVJkmxZsX+dMy5VztiUGs9lwIABimy3TzLf\njoZ1uGNU7q06DLDce0xS1fOv7Gc/+5kuvPBCZWRkKC8vT5K06KNMFeWfUNo9r0lOlxyuqCr7Hcj3\n6vUNe/RiLeFKkubNm6fU1NSAkAsAAIDmg4CFWgXTQ/PNULa3335bbWLjFBXhrHb7ihUrlJ+fr4UL\nFwZMhV5YWHjK9Z3KEEFXSpq8Oz+p0rYkJ0tpaWnq87t/1Hi8PWs3qOzooWpDWNas6xQ7aJySRt1Z\nz7PwOXDggN5++23dfffdcrlcp3QMAAAAhBcBC7Vyu5zak3tCiz/Zq/Qp38725939uSTphmc/VER8\nO8k4dN0Ty3T4rVk1HuubIFU5PGzdulVr164N+runTnYqQwTbdBui4xvfkXf3RrnTfMMOy4sKdWLb\nRxp3+y1aWsvrpY6bHNCDtuDOc3XZxF+q+MA2pV45pdbetLrMnz9fZWVlzB4IAADQjBGw8J25Ejso\nbsg1yn3naU2eHKVhw4bJ7XYrKytLq1at0u23364RI0Zo1KhRioiI0M0336xf/vKX2r9/v6ZPn660\ntDSVl5ef0mtHxLfzBbwaDB48WFEdDgas83QfoqiOPZWz7HElDv+R/4uGJavJkydr6ZzP/W13PTpO\n0X1GKuXyn0vyfV9WZcOHD5czOlHG6ZI7rV/AtuKc3SrJ2S1JsqXFKj16SK+//rokadiwYUpNTQ1o\nP2/ePPXt21cDBw6s30UAAABAk0HAQkgkDrtFruQuWrNmjebMmSNjjLp06aKRI0eqe/fukqTevXvr\nlVde0bRp0zRu3DidccYZevB3v9d777yt1atXN1qtxjiUeu105b7/nI6sekq2tERRnXqo3YTfq0uX\nLpK+DViy5b7lFBR++U/lr33V/7xo90aNHz9ekvT+++9r+PDh/m2ffPKJNm7cGDDLIgAAAJofU/n7\ngeoQdEO0LJWHBtbly99eWussed/1+JkzRzf79kArYcJdAAAA4UAPFkLK7XISOAAAANBqOepuAgAA\nAAAIBgELAAAAAEKEgAUAAAAAIULAAgAAAIAQIWABAAAAQIgQsAAAAAAgRAhYAAAAABAiBCwAAAAA\nCBECFgAAAACECAELAAAAAEKEgAU0Im9JWYO2BwAAQHhFhLsAoDVxu5xKn7I86PaZM0c3YDUAAAAI\nNXqwAAAAACBECFgAAAAAECIELAAAAAAIEQIWAAAAAIQIAauZY1Y6AAAAoOlgFsFmrr6z0n3520sb\nsBoAAACgdSNgtTJMEw4AAAA0HIYIAgAAAECIELAAAAAAIEQIWAAAAAAQIgQsAAAAAAgRAhYAAAAA\nhAgBCwAAAABChIAFAAAAACFCwAIAAACAECFgAQAAAECIELAAAAAAIEQIWAAAAAAQIgQsoAnzlpQ1\naHsAAACEVkS4CwBQM7fLqfQpy4NunzlzdANWAwAAgLrQgwUAAAAAIULAAloxhiACAACEFkMEgVaM\nIYgAAAChRQ8W0ILQwwQAABBe9GA1Md6SMrldznCXgWaKHikAAIDwImA1MfyCDAAAADRfDBEEAAAA\ngBAhYAEAAABAiBCwAAAAACBECFgAAAAAECIELAAAAAAIEQIWAAAAAIQIAQsAAAAAQoSABQAAAAAh\nQsACAAAAgBAhYAEAAABAiBCwAAAAACBECFgAguYtKWuUfQAAAJqriHAXAKD5cLucSp+yvF77ZM4c\n3UDVAAAAND30YAEAAABAiBCwADSo+g4RZEghAABozhgiCKBB1XdYIUMKAQBAc0YPFgAAAACECAEL\nAAAAAEKEgAUAAAAAIULAAgAAAIAQIWABAAAAQIgQsBoYU04DAAAArQfTtDcwpqgGAAAAWg96sAAA\nAAAgROjBQqPw7vpcef+cr+KD22QiIuU543tKHDExoE1p/kHtffq2avfv8vMFcrhjJEnlJV4dWTVX\nJ75eJ4c7RgkX3Kzosy4IaJ//4es6vnm1Sn93SVD1HfjrFKm8XKqmB/HYZyt1ZMVsdbrrOUXEt5Mk\n5Sz/o45velfmEV8bhydOruQuij/3B/KcPsi/765Hxnx7IIdTqc8l6pi7rdzpAxU74DI5oxOCqg8A\nAADNAwELDc6btUkHF/5GntPOVuqVv1LZiWPK++fLOrjg1yr6w/gq7eMyxqtNtyEB60ykx//46PrX\n5c38VMmX36vi7J3KWfaEItudIVdSJ0lS6dEc5f/rNbUd/6AiIhruLe5oE6+1767QVXPWqux4ro7+\n+w0dWjRDbX/4W3nSB/jbRfcZpdgBl8paq2d+2EM3Pjxfxz5epmMblir16t/I3fmsBqsRAAAAjYuA\nhQaXv/ZVRcS1VerV98s4nJIkV3IXHZh3r5577jlJXQPaRyS0V1SnnjUe78SODYodNFptug9Rm+5D\ndHzzanl3feYPWLnv/kVtepzf4MHFOCKUkZGhqMWHJUnurv2156kf6diGJQEByxmb7D+fsWNHK3Gt\nQ3GDxurAK/cpe/FD6vTjZ+VwuRu0VgAAADQO7sGqJ2YFrL+ifV/JnT7AH64kKapDdzk8cXrjjTfq\nfTxbViITEeV/7nBFyZYWS/KFL2/WJiWO+NF3L7yeHFFt5ErqqNLc/XW2dUYnKnHERJUfz9PxLWsa\noToAAAA0Bnqw6olZAU+Bccg4q77VjDNCmzZtUtSgwPV5/3hJR1bOkXG55U7ro4QLblZkarp/e1TH\nHjq+6V216TFUJYd2qvjQTiVddJdsaYmOvPO0EobdIqcn7pRKLS0tlS0/KURbG9S+trxMZUdz/Pdp\n1cWdPlByOFW0d4ti+19c31IBAADQBBGw0OBcyZ1UtO+rgHWl+YdUVpCrI8UudfhmpdOlmAGXypN+\nthxt4lRyeI/y1y/SgfmT1OGmJ+RK6SJJih86QYcWzdDeOTdLkuLOuVpRnc5S3tpX5fTEK6bfqYWV\nor1b5HK56rXPN4GsrCBX+f9aoLLjuYobcm1Q+zpcUXJ64lRWkHsq5QIAAKAJImChwcUOGqfDyx5X\n7pqXFTdorMq9x3R4xZOSMXI4vh2lGhGTpORL7vY/d3fpI8/pg7Tvuf9W/rrXlDL2f3ztYlPU4Uez\nVZp3QA53tJyeOJXkHdDRj/5P7W94RLa0SEfee1aFW9frtNcSdLTbJYobNLbOOl1tT9O/li/U2Nkf\nBKw/8fV65a97rUr7soLDAYHMRHoUf/4Nih1c92t9w8pKJujmrYK3pExul7PuhqfYHgAAoCERsNDg\nYnqPUOnhPTr67zd0dN1rkozanPV9ec4YrPZl2SqvZd+IuFRFde6logNfB6w3xsiV6O/7Uu6qpxXT\n72JFtj1duWvmqfjANnW8bY7emNhLZw85T67kLgETT1TH4fJo8ODBiupwMGB98aEd1bdvk6AP/7FK\nY59cK6cnVs7YlID7zOpSXlKk8sKjckYnBb1Pa8AwXAAA0JwRsNAoEi64SXEZ41Waf0DONvFyRidq\n7zN36fzRIxTcFA81d/MUbl2n4kM7lDJusiTJu2ODovuOkrNNvAYMGCBP+kB5d35cZ8CqL+Nw+gJZ\n+4N1N66Gd+fHki2Xu3OvkNYFAACA8GEWQTQaR6RbkanpckYn6sSODSo9skd33XVXrfuUHj2koj1b\nFNWhe7Xby0u8OvLuX5R44R1yRLXxr7cl3m/bFJ8IeqKKxlJ2PE+5q1+QMyZJbU76kmQAAAA0X/Rg\nocEVH9yuEzs2KLLdGZIk754tOvrR3xQ35Bqdd9550hLfcLAj7z0rWauojj3lbBOvkiO+SS5kHIo/\n94fVHjv/XwvkSuqs6LO+71/nTh+gYxuWyZXUWf/7v7vk3fWZ4s65quFPtAZlxw6raO+XstZq2TKr\n3DXzVfDZSklWba+ZJocrqs5jAAAAoHkgYKHhOSJ0Yvt/lP/h36SyErmSuyj54p8opt9FAc0iU9J0\n7JO3dHzjOyov8crhiZU7rb8Shk6QK7lzlcOWHM7SsY+Xq8MtfwpYH3/edSo7nq/Db83SQx/GKWHY\nrfKcdnaDnmJtjm96R8c3vSM5nLp1RYKK3O0UO2iMYgdcJmeb+LDVBQAAgNAjYKHBRaZ2VfsbH62z\nXUy/i+s1xboruYvS7l1UZb0j0qOU0fdI8k2AEMyECe2vn1njttj+lyi2/yUB61JG3xtUjV3vWxbw\nPNh6AAAA0DxxDxYAAAAAhAgBCwAAAABChIAFAAAAACHSYgJWVlaWrr32WsXHxysuLk5XX321du/e\nHdS+U6dO1cUXX6zk5GQZY/Tiiy9WabN69WoZY2SM0a5HxlRZivZ+6W+bu/pF7Xv+bu3+0w/l8Xi0\n95m7lLf2VZVXmjocgFR6NFvZb/xeu//4A+3+43gdeuOhoD+3klSSk6Wrr7lWKSkp8ng86tGjh2bN\nmhXQJicnRxMnTlRqaqo8Ho++d845WrlyZZVjDR8+3P8Zr7z86U9/qtIWAACgJi1ikovCwkJdeOGF\nioqK0ksvvSRjjO6//36NGDFCn3/+uaKjo2vdf/bs2RowYIDGjBmjefPmVdvm7LPP1rp16yRJV81Z\n619/+K0/q9x7TJGVvqepvLhQMX1HyZXUSX+96/u6ZsYLyl+3UMUHtqntNb8JwRkDzV95iVcHF0yV\ncbp8k4YYo7w1L2vEiBEqv+JROSLdte5ftP9rHVwwVaWXjJLjgrsUHxWtnNx9euDNjfrjft9EIra0\nRPvn3avyE0eV8P2bFBedqK52o8aMGaNVq1Zp+PDhAcfs16+f5s6dG7AuPT09lKcNAABauBYRsJ55\n5hnt2LFDX331lbp16ybJ94tS9+7dNXfuXP3iF7+ocV9vSZny8/PlcDi0bdu2GgNWXFycMjIyJElR\niw9LkkrzD6nkcJbizrlKxuH0t02++L/9j0eOHKmEVV7Z0iIdXf+6ygrzmZobkFTw2UqV5h1Uxzue\nliuxoyTJlZquXc/+WHGfvlXrd5dZW67Dy5+Qu2t/LVmyxD8zo7trv4B2x7/6QCXZmWo34fdyp/m2\nLXp4mvr376/Jkyfro48+CmgfGxvr/5wDAACcihYRsJYsWaKMjAx/uJKk0047TUOHDtWbb75Za8By\nu5z+X85KcvdJkv5n0Wea8WX1U2lnzhztf1yw+T1JVtF9LqyzRoc7TpICghjQmp3Y9qGiOvbwhytJ\nciW019ChQ/Xhtg9rDVje3RtVcjhLSZf8pNbXKN73pUxElD9cSZIxRhdffLEef/xx7d27V506dfru\nJwMAAFChRdyDtXnzZvXp06fK+t69e2vLli0N9rrHN72nyHZnKDI1vdrttrxMBQUFOpH5qY79Z7Gi\n+14khzumweoBmpPinN1ypXStsr53794qyan9PqyiPb7PtS0tUUZGhnY9doWyZt+gI+/MVXlJ0bcN\njUOq5o8aUVFRkqRNmzYFrP/kk08UHx8vl8ulfv366bnnnqvvaQEAgFauRfRgHTlyRImJiVXWJyUl\nKTc3t0Fes2jvFyrN3afEkXdWu704O1P7n79bsY/5nkf3uVDJl97dILUAzVH5iYJq/+CQlJSkcm9B\nrfuWFfiG6eYseUQ/nnSvMk+/QkX7tyn/g1dUejRbba++X5LkSuosW1yokpwsuVK6SPINC/7mfsoj\nR474j3nBBRfohhtu0Jlnnqm8vDzNmzdPt99+u7L27NWM6dNCcs4AAKDlaxE9WOFQsOldyRGh6F7D\nqt3uSuyo9jf/UatXr1bCBTercOs65Sx/opGrBFooayVJ0b1H6MEHH5Q7rZ/ih1yt+KETdOLr9SrJ\nyfJt7zVMDk+ccv7+hIqzM1VWmK8nHntE76/+hyTpZws+VfqU5Uqfslzziofooe0ddctbBfr5ught\n6D5Rnu4ZemTmwyooqD3wAQAAfKNFBKzExER/T5W3pMy/vqaere/Klpao8MsP5DljcI0TVpiISEV1\n6K5hw4Yp/twfKGnUj1W45R8B07kDrZnDHVNtT9WRI0fqHErr8MRKktzpAwLWe04bKEkqPrTd/xqp\nV01VeeFR7X/+bu2ZfYOef/55xZ9/vSTJGZNU6+tEnzVMXq9XGzduDO6kAABAq9cihgj27t1bmzdv\nlhQ4acWBZf+UPO39z6tTedKKYBVu+1Dl3gLF9BkZ9D6R7X0TcJTk7VdUp571fk2gpXGlpFV7r9WW\nLVvkSkmrY9+q924FMv5H7i591PHHz6o0d59ky7V17p1KvvA2mYgoRbbrVssxKh3NmLobAQAAqIX0\nYI0bN07r16/Xjh07/OtK8w+qaO8XatPtnJC/3vFN78rhiZPnjO8FvU9Rlu9meldC+5DXAzRHbboN\nUdG+L1WSd8C/rjT/oNauXVvn59Zz+iDJ6ZJ358cB60/s2CBJAd9LJ/kCkiupk1zJXVRYWKiCz1Yq\nuveIOr9r6/iW1fJ4POrbt299Tg0AALRiLaIH64477tCTTz6pK664Qr/73e9U+PUG5f3zZUXEpihm\nwGX+dqX5h7R37u2KHzpBCUMn+Nd7d29UWWG+yo77hhkWH/hax12+X7yie54f8FqHDh3SiZ0fK3bg\n5TLOqpev+NBO5b7/vNr0GKqIhPb6+9+Ncle/qGMblsh9+iBFdTqrIS4B0OzE9L9Exz5epuy//VYJ\nF9wkySjvny+rS5cuKq3jc+v0xCk+Y7zy/7VAU6dO1YnMaBUf+Fr5/1qg6D4jA6Z+z/3Hi4ps103O\nNnEqyd2vQYN+ITmcShh2i7+NN2uTjq5/XZ4zz1NEfFvZokIVbHpXJ7Z9qJkzZ9b5ZeUAAADfaBEB\nKzo6Wu+9957uvfde3XTTTSooKpW7a38ljbxDjkhPpZZWsuW+pZK8D17x9zBJ0rGPl+vYx75hhdE9\nlwW0feWVV6TyMkXXMDzQGZ0ghydO+esXqfx4rm5aFSNvZLISR0xUTL9LQnPCQAvgiHSr3YSHlPvu\ns8pZ9rgkyd21v957+1UNf3pzpZbVf27jh06QI9KjhQsX6tDOXXLGJCrunKsVf951Ae3Kjucp991n\n/F/yff2tE7TYMVTOivu4JN+9WNZa5X8wX2Unjso4IhSZmq6UsZN03333Ndg1AAAALU+TD1jekjK5\nXXV/OW9aWpr+9re/SVKN91xFxLdT1/uWVVnf/vqZQddz7733atbBM2vc7oxOVOq4Sf7nmTNH13oP\nGNCaRcS1VepVUwPWpaenS/o2YNX0uTXGKO6cq7TtvWdr/YylXH5PwPPZM0dr6UntXYkd1e4HD9T/\nBAAAAE7S5ANW5UkrgnEqk1YAAAAAQCi0iEkuAKChVP7qh6bQHgAANG3GVnxhZxCCbhhq9e3B+i7t\ndz0ypl61AWiaqhtWKH33nxEN0b6FYm57AECrRA8WADQj9JABANC0Nfo9WMFOWgEArUF9fyZyXyoA\nAE1bowcsfjkAgG/xMxEAgJalyc8i2Nhqum/jG03x/g3a076ptm+qNbUm9e0hY5QBAADfTdCTXDzw\nwAMrJKU0bDl+HSXta6TXCiXqbjzNsWaJuhtTc6xZajl150yfPv3ScBUDAEDYWGub3DJjxgwb7hqo\nu2kvzbFm6qZm6mZhYWFhYWn5C7MIAgAAAECINNWA9UC4CzhF1N14mmPNEnU3puZYs0TdAAA0a/X5\nomEAAAAAQC2aag8WAAAAADQ7BCwAAAAACJGwBixjzKXGmK+MMduMMVOq2X6BMeZjY0ypMebacNRY\nnSDq/oUxZosx5nNjzLvGmK7hqPOkmuqq+S5jzEZjzKfGmA+MMb3CUefJ6qq7UrtrjDHWGDO4Meur\nSRDX+1ZjTHbF9f7UGHN7OOo8qaY6r7Ux5gcV7+3Nxpi/NnaN1QniWv+x0nXeaozJC0edJwui7jRj\nzPvGmE8qfpZcHo46T6qprpq7VvzM+9wYs9oY0zkcdQIAEFbhmr5QklPSdkmnS4qU9JmkXie1SZfU\nT9I8SdeGe8rFetQ9QlKbisf/Jem1ZlBzXKXH4yStaA7XuqJdrKQ1ktZLGtwc6pZ0q6Qnw11rPWvu\nLukTSYkVz9s2h7pPav9TSc83h7ol/UXSf1U87iUpsxnUvEjSLRWPL5T0crivNQsLCwsLS2Mv4ezB\nOkfSNmvtDmttsaQFkq6o3MBam2mt/VxSeTgKrEEwdb9vrS2seLpeUrj/ihtMzUcrPY2W1BRmP6mz\n7gq/lfSIJG9jFleLYOtuSoKp+Q5Jc6y1uZJkrT3UyDVWp77XeoKkVxulstoFU7eVFFfxOF7h//Lh\nYGruJem9isfvV7MdAIAWL5wBq5OkrErP91Ssa+rqW/dtkt5q0IrqFlTNxpifGGO2S3pU0s8aqbba\n1Fm3MeZsSV2stcsbs7A6BPseuaZiKNXrxpgujVNajYKp+UxJZxpj1hpj1htjLm206moW9OexYqju\nafo2AIQX0tDFAAACx0lEQVRTMHXPkHSjMWaPpL/L1/sWTsHU/JmkqyseXyUp1hiT3Ai1AQDQZDDJ\nRQMyxtwoabCkx8JdSzCstXOstWdIuk/S/eGupy7GGIekJyT9Mty1nIKlktKttf0krZL0UpjrCUaE\nfMMEh8vXE/SMMSYhrBXVz3WSXrfWloW7kCBNkPSitbazpMslvVzxnm/K/kfSMGPMJ5KGSdorqblc\nbwAAQiKc/1nvlVT5r/adK9Y1dUHVbYwZJenXksZZa4saqbaa1PdaL5B0ZYNWFJy66o6V1EfSamNM\npqQMSUuawEQXdV5va+3hSu+LZyUNaqTaahLMe2SPpCXW2hJr7U5JW+ULXOFUn/f2dWoawwOl4Oq+\nTdJCSbLWrpPklpTSKNVVL5j39T5r7dXW2oHy/fyTtbZJTCoCAEBjCWfA+rek7saY04wxkfL98rMk\njPUEq866jTEDJc2VL1w1hftUgqm58i/KoyV93Yj11aTWuq21+dbaFGtturU2Xb773cZZa/8TnnL9\ngrneHSo9HSfpi0asrzrBfB4Xy9d7JWNMinxDBnc0ZpHVCOrniDGmp6RESesaub6aBFP3bkkjJckY\nc5Z8ASu7UasMFMz7OqVSL9uvJD3fyDUCABB2YQtY1tpSSXdLWinfL5cLrbWbjTEPGmPGSZIx5nsV\n9x+MlzTXGLM5XPV+I5i65RsSGCNpUcXU0GENjkHWfHfF1NufSvqFpFvCVK5fkHU3OUHW/bOK6/2Z\nfPe73Rqean2CrHmlpMPGmC3yTWAwyVp7ODwV+9TjPXKdpAXW2qYweUuwdf9S0h0V75FXJd0azvqD\nrHm4pK+MMVsltZP0UFiKBQAgjEwT+X0DAAAAAJq9pn7DNAAAAAA0GwQsAAAAAAgRAhYAAAAAhAgB\nCwAAAABChIAFAAAAACFCwAIAAACAECFgAQAAAECIELAAAAAAIET+H43MNOhjF8O4AAAAAElFTkSu\nQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x111461b00>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"y = np.asarray([0, 1, 0, 1, 1, 0, 0, 1, 0, 0])\n",
"N = len(y)\n",
"with pm.Model():\n",
" # Set Transform to None so it won't use theta_log_\n",
" theta = pm.Uniform('theta', lower=0, upper=1, transform=None)\n",
" obs = pm.Bernoulli('obs', p=theta, observed=y)\n",
" trace = pm.sample(3e3, tune=1000, njobs=2)\n",
"pm.posteriorplot(trace[1000:]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Model 2: Log Odds Parameterization without Jacobian"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using advi...\n",
"Average ELBO = -6.2268: 100%|██████████| 200000/200000 [00:14<00:00, 14103.89it/s]\n",
"Finished [100%]: Average ELBO = -6.224\n",
"100%|██████████| 3000/3000.0 [00:02<00:00, 1171.33it/s]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAACsCAYAAABmdA06AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FHX+x/HXN7ubQiokoRNAmgoiCgqod4J0RBAUPbsn\nevo79exdBL1TUWwcomI7AbtgoXjoCYKCoNgREEQgCUVSSSF9d35/bLJkSduFTX8/H488yM58Z+Y7\nw2Y+85n5fr9jLMtCREREREREjl5QfVdARERERESkqVCCJSIiIiIiEiBKsERERERERAJECZaIiIiI\niEiAKMESEREREREJECVYIiIiIiIiAaIES+QIGGOuNMasCXRZERGRo2GM6WKMsYwx9vqui0hzpQRL\nREREpBEzxuwyxgwPwHqUnIkEgBIsERERERGRAFGCJVINY8zdxpjfjTE5xpjNxpiJVZSzjDH/MMbs\nMMakGWNmGmOCDivzhDEm0xiz0xgzptz0vxpjtpRuY4cx5tra3i8REWkajDELgARgiTEmF7igdNYl\nxpik0ph0X7nyQeViW7ox5l1jTKvS2V+U/nvAGJNrjBlsjOlmjFlZWjbNGPOGMSam7vZQpPFRgiVS\nvd+BPwHRwIPA68aYdlWUnQgMAE4GJgBXlZs3ENgKxAGPA68YY0zpvBRgHBAF/BV42hhzcoD3Q0RE\nmiDLsi4DkoBzLMuKAN4tnXUG0AsYBjxgjDmudPqNwLnAmUB7IBOYUzrvz6X/xliWFWFZ1jrAAI+W\nlj0O6ARMr819EmnslGCJVMOyrPcsy9prWZbLsqx3gN+AU6so/phlWRmWZSUBzwAXlZuXaFnWS5Zl\nOYF5QDugTek2llmW9bvlthr4FHdSJyIicqQetCwr37Ksn4CfgBNLp18H3GdZ1m7LsgpxJ0vnV9Xv\nyrKs7ZZl/c+yrELLslKBp3AnZyJSBXViFKmGMeZy4FagS+mkCNxPoZyVFE8u93si7rt9Zf4o+8Wy\nrLzSh1cRpdsYA0wDeuK+6dEC2BiQHRARkebqj3K/51Eac4DOwAfGGFe5+U5Kb/odzhjTBpiF+8Zf\nJO44lRnw2oo0IXqCJVIFY0xn4CXgBiDWsqwY4BfczSUq06nc7wnAXh+2EQIsAp4A2pRu4+NqtiEi\nInI4y4+yycAYy7Jiyv2EWpa1p4r1PFI6/QTLsqKAS1GMEqmWEiyRqoXjDiqp4B6MAuhTTfk7jDEt\njTGdgJuAd3zYRjAQUrqNktKnWSOPqtYiItLc7AeO8bHsC8DDpTcRMcbEG2MmlM5LBVyHrSsSyAWy\njDEdgDsCU2WRpksJlkgVLMvaDDwJrMMdvE4A1lazyEfAd8CPwDLgFR+2kQP8A3en5EzgYmDxUVVc\nRESam0eB+40xB4Dzayg7C3ec+dQYkwOsxz0QE5Zl5QEPA2uNMQeMMYNwD/B0MpCFO7a9Xzu7INJ0\nGMvy56myiFTGGGMBPSzL2l7fdRERERGR+qMnWCIiIiIiIgGiBEtERERERCRA1ERQREREREQkQPQE\nS0REREREJED8edGwHnWJiEhtCsS7dRSrRESkNtUYq/QES0REREREJECUYImIiIiIiASIEiwRERER\nEZEAUYIlIiIiIiISIEqwREREREREAkQJloiPCoqdtVr+SJcREZGGTzFBpPnw50XDGvpWmr0udy/z\nueyuGWf7Vb5sGZFmTMO0S5OmmCDSJGiYdhERERERkbqiBEtERESkCVAzRJGGwV7fFRARERGRoxfq\nsKkZokgDoCdYIiIi0uzVxUBGItI86AmWSANSUOwk1GGrtfIiIlI5f5/+6MmPiFRFCZZIA6IALyLS\nOOgGl4hURU0ERRqIotRERo4cSdJT55M86yLSlj2DMz/H7/Vcd911GGO49NJLqy03Y8YMjDGcccYZ\nXtNzcnK44IIL6N69O+Hh4cTExHDqqafy+uuv+10XEZGmquyGmK8//nIW5HL11VcTFxdHeHg4w4cP\nZ+PGjX6v5+Dm1SQ+No7dc66odP5LL73EscceS0hICL169eKFF17wmp+dnc1DDz3EaaedRmxsLDEx\nMZx22ml8+OGHftdFpLlQgiXSAJTkpLP/rXvIz88n/tx7aDXiOgoSfyR14YNYlsvn9axdu5bXX3+d\nqKioasvt2LGDf/3rX7Ru3brCvKKiIux2O/fccw+LFy/mzTff5LjjjuOyyy7j6aef9nvfRETEP5Zl\nkbrwIZYvX87s2bNZtGgRxcXFDB06lN27d/u8HldBLhkrX8IW3rLS+Tk/Lufaa6/lvPPOY/ny5Uye\nPJm///3vPP/8854ySUlJPPfcc5x55pm8/vrrvPPOO/Ts2ZOJEycyZ86co95XkaZITQRFGoDsb97H\ncjlZsmQJ/WasBcAWGcv+N+8mf9t6WvQ6rcZ1FBcXc+2113Lfffcxd+7casv+3//9H5dccglbt26l\npKTEa15sbCxvvvmm17SxY8eybds2Xn31VW655RY/905ERPyRv/1rCvdsZsHKlQwdOhSAwYMH07Vr\nVx5//HH+/e9/+7SezFX/ITi+K7aIVhTs+tFrnuVycuDLBVx22WU8/PDDAAwdOpS9e/cydepUrr76\nahwOB127dmXHjh20aNHCs+yoUaNITk7mscce4/rrrw/QXos0HXqCJQ3C9OnTMcbw66+/MmrUKMLD\nw0lISOA///kPAAsWLODYY48lIiKCoUOH8vvvv3st/+KLL3LiiScSGhpKXFwcU6ZMISMjw6vMs88+\ny+DBg2nVqhUxMTEMGjSIZcu8m23s2rULYwxz587lgQceoF27dsTExHDOOef4ddfQX/nbvyas2wBi\nYmI800I79cEWFU/e9vVVLld+FKuZM2fidDq5/fbbqy3/5ptv8v333/Poo4/6VcfY2Fjsdt2TEZH6\nUxexIvu7JexbcBvJs/5C0jMXsm/+beT9vsGrTEnWfowx5Pz4Xw58+Tq7n72MpGcuJGXhg5Rkpx31\nfub99jW2iFae5AogOjqac845h48++sindRTs3szBTatoNfL/Kp1fuGcLrrysCs3JL7vsMtLT01mz\nZg0A4eHhXslVmQEDBrB3715fd0mkWVGCJQ3K5MmTOfvss/nwww/p378/V111Fffeey/PP/88M2bM\n4D//+Q9bt27l4osv9ixz9913c/311zN8+HAWL17MzJkzWb58OWPGjMHpPJSA7Nq1i6uvvpr33nuP\nd955hwEDBjBu3DiWL19eoR6PPvoo27dv59VXX2XWrFmsW7euQhCyLBeWy1nlT0lJift3y6p2n13F\nhZQc2E9wXOcK8xxxCRSnJVe5bFkfgA7XvsT90x4is9/l9Jj6Kbsz8/nwhz0V+gDk52Zzyy238Pjj\nj9OqVatq62VZFiUlJaSnp/Piiy/yySef6OmViDQItRkrnFkpRPQdRdyEu4kffych7bqTuvBB8nd8\nV6EeWeveozhzH7FjbqLVsL9RuHcraUuf8CpTFis8MaGqn3KxojgtCUclMaF3794kJSWRm5tb7fGx\nnCVkLH+WqFMn4WjZvtIyxWlJAPTp06fCNgA2b95c7Ta++OILjj322GrLiDRXuh0tDcodd9zB5Zdf\nDrjvji1ZsoS5c+eyc+dOT7+iffv2cdNNN5GYmIhlWcycOZNp06bxwAMPeNbTs2dPzjjjDJYsWcK5\n554LwBNPHAp6LpeLYcOGsW3bNp5//nlGjx7tVY8uXbp4NZNLTU3ljjvuoEPvq7BHxgKQ/vEsDv6y\nosp9ccx0/xs79mYiThheZTlXQS5gERQaUWGeLTSSkow9VS5bJuOT52jRczChnftWW+6OO+6gZ8+e\nXHnllTWuc86cOdx4440AOBwOZs2a5fm/ERGpT7UZK1qeNcUz37JchHbpR3HGXnJ++JiwY/p71cMe\n3Zr48Xd4Pjvzsjiw6lVKctIrxIqymFCV8rHCVZCDPbpiH9myG2OZmZlERFSMGWWyvl6I5SwmevDk\nKsu4Yw+0bOndP6tsG4c/2SvvxRdfZP369Rr8SKQKSrCkQRkzZozn95YtW9K6dWtOOukkr0Ebyu6Y\nJScns2XLFlwuF5dccolXX6KBAwcSGRnJF1984Qma3333HdOmTWPDhg2kpqZ67hb26tWrQj3Gjh3r\n9fmEE04AwJmd6gmaMWdcTOTJ46rclyU3nsE5s9dgj2njmWa5vF9MaYKOfojf3E2fU/jHb3S45oVq\nyxUk/8L89+bz/fffY4ypcb0XXnghgwYNIi0tjcWLF3PjjTdis9m49tprj7rOIiJHozZjReEf28la\n8waF+37DlZcFuGOFvVXHCvUI6zbA63NwvPupU2WxoiwmVKV8rDgaxZl7yV73LvET78PYgwOyzvJW\nrVrFP/7xDy6//HIuueSSgK9fpClQgiUNyuF30oKDgyudBlBQUEBKSgoA3bt3r3R96enpgDvADhs2\njOOPP57Zs2eTkJCA3W5n6tSpbNmypcJyhzefCwkJAcByFnmm2aLisUXGVbkv/fr1I7jNXjDulrgl\nWfvZ88IUrzIdrnuFoBbRgPHcTSzPWZBDUGhkldvIzc0lc+XLRA88D2NzHFqHZWG5nLgKcjGOUIzN\nTsYnc5gyZQodO3bkwIED7jqVlOB0Ojlw4ABhYWGe/QSIj48nPj4egNGjR5OXl8ftt9/OVVddhcPh\nqLJOIiK1rbZiRUl2Kvvfvo/g2E60Gn4t9qh4CLJx4MvXKU6v2Fz78POzsbvPjZXFCk9MqIo51Gsj\nKDSi0phQ9lTp8H31KvPZXEIT+hLSvpdnHZazBLDcn20OghwhnlYTmZmZtGvXrsI2KmtGvmHDBsaP\nH89ZZ53Fyy+/XPW+iDRzSrCkUYuNdd8h/PTTTysNOGXzly9fTlZWFu+++y4dOx66C5mXl3fE2/a3\niaAtohVtL/ce5twW0Qpjc2CPbu1pD19ecVoyoQl9Kkwvk5aWhisviwNfzOfAF/O95uX9mkrer18S\nP/E+WvQcTHF6Mi+88EKFd5yAO1g//fTT3HzzzVVua8CAAcybN4/9+/d7HUMRkYbO11iRv+M7rMKD\nxE24G3vUoRtoVknhEW/7SJoIOuISKNj5Q4UymzdvJiEhodrmgcVpyTizU0ie9ZcK85Jn/YXI/uNp\nNfxvOOISANi0aZNXglXW9+r444/3Wnbjxo2MGjWKfv36sWjRIt1oE6mGEixp1EaMGEFQUBBJSUmM\nGDGiynJliVT5gLBt2zbWrl17xMmCv00Ejc1BSLselZYN6z6Qg7+sICsryzOtYPcmnNkphHUfWOU2\n2rZtS5uLHqkwPXXx4wTHdSH6tAs8HaXbXPQIb/9tsFe5m2++GafTyezZs6u8s1tm9erVREREVPru\nLBGRhszXWFGWSBnboebbxRl7KNy9udoWC9U5kiaCLboP5ODGz/jfipWMGHYW4H7h75IlS7wG7qhM\n/Pg7vZ6gAWStX0jRH9uJP/duz36EtD+WoLAo3njjDYYPP9RP+PXXX6dVq1acfvrpnmm//fYbI0aM\n4JhjjmHp0qWEhYX5fgBEmiElWNKodevWjbvuuosbbriBrVu3cuaZZxIaGkpycjL/+9//uPrqqxk6\ndCjDhw/Hbrdz+eWXc9ttt7Fv3z6mTZtGQkICLpfvL/Itzx7dBnt01W3mBwwYQEi7/T6tK2rgJA5u\nXsX48ePJbzMMV+FB9/tL2vWiRc9DSVFB0kb2v30fsWNvIqLPMEJDQwlNqDiwhbEFYwuP8ZoXmtCX\nIUOGeJWLiYmhpKTEa/rcuXNZv349w4cPp2PHjqSnp/Puu++ycOFCZsyY4Wl2IyLSWPgaK0I794Mg\nG2lLnyLq1Ik4czM5sOYN7FHxNY4IW5WyWOFPTAjrMZCQ9sdy1ZVXUHDSRQSFRpC1/j2K84r4wNmf\nJXcfesVI4uPjCe8zjLixN7FrxtmEdKg4sl/uxhUYm8MrJhibnZg/Xcq8ec/ToUMHhg8fzsqVK3n1\n1VeZPXu251yfkpLCiBEjKCoq4sEHH6wwuuBJJ53k1bxcRJRgSRPwyCOPcNxxxzFnzhzmzJmDMYZO\nnToxbNgwevRwPzHq3bs3b7zxBg888ADjx4+nW7duzJgxg+XLl7Nq1ar63QHAHhlHm4seITjpQ1I/\nfAQTZCesxyBanjUFY8q/TcECywVHGOh9ccIJJ/DRRx9x++23k5GRQVxcHMcddxxLly7l7LPPrrXt\niogESvl3BJbxJVYEx3cmbtztHFjzBimL/okjph0tz7yC/J3fU5C0sc7qb0wQ8edPY0TOp8x763ms\nkmJCOvSizUWPuPuFlWe53D9HIPKksTwyqS9PPvkkM2fOJCEhgWeffZa///3vnjKbN28mMTERgHHj\nKrba2LlzJ126dDmi7Ys0VcaPOzK1d0Un0kh0uXtZzYVK7Zpxtl/lj2SZIykv0oDVPLxlzRSrBPDv\nfA11c/5taHUqW0ZE/FJjrNKLhkVERERERAJECZaIiIiIiEiAKMESEREREREJECVY0ixV1gFaRERE\nRORoaRRBaZZCHTZ1BBYRaSQKip2EOmw1FxQRaQCUYImIiEiD5u9NMd0QE5H6pCaCIiIiIiIiAaIE\nS0REREREJECUYImIiIiIiASIEiwREREREZEAUYIlIiIi0kz5+9oSveZEpGYaRVBERESkmdIIjSKB\npydYIiIiIiIiAaIES0REREREJECUYImIiIiIT46kD5b6bUlzoz5YIs1IQbGTUIet1pcREZGmyd8+\nW6B+W9L8KMESaUYUGEVERERql5oIioiIiIiIBIgSLBEREakz6o8jIk2dmgiKiIhInVFTZRFp6vQE\nS0Sq5e/dZt2dFhERkeZMT7BEpFr+3m3WnWYRERFpzvQES5oEPTURERERkYZAT7CkSdBTFhERERFp\nCPQES0REREREJECUYImIiIhIg6HBlaSxUxNBEREREWkw1OxfGjs9wRIREREREQkQJVgiIiIiIiIB\nogRLREREREQkQJRgiYiIiIiIBIgSLBERERGpNRrlT5objSIoIiIiIrVGowJKc6MnWCIiIiIiIgGi\nBEtERESOmJp/iYh4UxNBaXAKip2EOmz1XQ0REfGBmn+JiHhTgiUNjr/BGhSwRURERKRhUBNBERER\nERGRAFGCJSIiIiIiEiBKsERERERERAJECZaIiIh4aFRAEZGjo0EuRERExEOjAoqIHB09wRIRERER\nEQkQJVgiIiIiIiIBogRLRERERBqtI+k3qL6GUpvUB0tEAqqg2Emow1bry4iIiID//QZBfQeldinB\nEpGAUqATERGR5kxNBEVERERERAJECZaIiIiIiEiAKMESEREREREJECVYIiIiIiIiAaIES0RERESa\nFX+Hadew7uIPjSIoIiIiIs2KvyPearRb8YeeYEmt010fEREREWku9ARLap3uEomIiIhIc6EnWCIi\nIiIiIgGiBEtERERERCRAlGCJiIiIiIgEiBIsERERERGRAFGCJSIiIiIiEiBKsERERJoovSZDRKTu\naZh2ERGRJsrf12SAXpUhInK09ARLREREREQkQPQESxq1gsSfOfDl64TNOp8i7IR1O4WWQ6/CFt7S\nU6Ykaz97XpjitZx5zP1vp5veJig0AgBXcQEZ/5tL/m/rCAqNIObPlxN+3J+9lsv6eiEHN62i3ZWz\nMEG2Guv3x5t3g8tF20sfrzAv56dPyFg+mw7XvYI9ug0AV155JYnz5nnKBIVF4YjtRPTgCwg7pr9n\neuJj4w7tR5CNoJBwHLEdCe1yEpH9xmALj6mxbiLS+BQUOwl11HzukYrK4kXR/u0Ye7DP8aLM4fFi\nypQpJL/5br3Fi7RlT3PwlxWeMkFhUfz5qxPJbzu80njhLlQxXohI4CnBkkarIPkX9r87lbCuJ7Nw\n0SIuf24FB75cwP6376PdFbMwdodX+ahBk2nRfSAAH1x/OhPnrMUEh3nmZ69fSMGuH4kdewtFqTtJ\nW/oUwW264WjVAYDdu3eT9dU7tJ78kE/B8kgFtYim9aSpADgPZpK94QNS3ptO6wv/SViXfp5y4X2G\n89lL/+TcZ9fgKsimcO9Wcr5fSs53S4ifNJXQjsfVWh0Dzd+LRl1kSnPlb5M/NfdzKx8v4s+9B2d+\njs/xoszh8eJ/yWurjBcl2Wn1Ei+s1NVVxovIfqOxLKtCvPjqoq61Vj+R5koJljRaWWvfwh7VmvhJ\n9zN27FgivrBwxHbij/m3kPvzp0Se7H1hYY9pS0iHYwEYNGgQIR+me83P3/Edkf3PpkWPgbToMZCD\nm1ZRkPiTJ2DefPPNtOh1Rq0nLibI7qknQGjnE9n9/F/J+W6xV8C0RcYyaNAgQkv3o0X3gUT1P4c/\n3riL1A8fpsO1LxPkCK3VugaKLhpFpDaVjxdlCY+v8aIy+Tu+495bbuC5jN6VxovMFS/WS7xY9tqd\nxLRuX2m8KF+ufLyYNGkSwZfMaTTxQqQxUB8sabQK924ltEs/r7uDIe16EBQWRd62dX6vz3IWY+wh\nns9BjhCskiLAHUxXrVpFy6F/PfqK+ykopAWOVu0pydxXY1lbeEtaDr0K18EDHNz8RR3UTkSk4auN\neBEWduiJ1uHxoiD5l3qJF1FRUX7Hi/379yteiASYEixpvEwQxlbxIayx2SlOS6ww/cDqeSQ+Pp6k\npy9g/PjxFKXu8pof0r4XB39ZQUluBvk7vqMoZSch7XthlRST8dkLzJgxA1tY1BFV1XI5K/xgWT4v\n68xOIygk3KfyoV1OgiAbhXs2H1FdRUSanKOIFymLHqo0XsybN6/KeBFz5hUBiRclJSV+xYuSkhK/\n44Xdble8EAkwNREUvzWU/i+O2A4U7t3qNa0kKwVnbibYytXP5iCi32jCupxMUIsoitN3s3HjUv7Y\n+xntLnsKR1wnAKJPv4iU96azZ87lAESdOomQDsdxYO1b2MKimTJlCv+652O/61m4ZzNJMyf4tYzl\ncr+7xpmbSdZXb+M8mEnUwPN9WjbIEYItLMp9HERE5KjiRdb69/jj9TsqxIuiL5+sMl5E9B15RPU8\nPF44Zta8TPl4cf311/sdL+Li4jigeCESUEqwxG8Npb9MZP/xpC99kswvFpCScgrF6cmkL38WjMGY\nQw9n7RGtiB11g+dzaKc+fPHSnXTu3ousde8Qd87t7nKRcbT762xKDvxBUGg4trAoig/8QfY379P2\nksfIz88n/ZNnydu2HuMIIeqUc4nqf06N9XS07krs6H9UmJ7/23qy1r1TYbozN90rwJrgMKLPuITI\nATVvq4yFBcbn4iIiTVr5eBHV/xxcBTk+x4uwY/qz95W/V4gXP/30Ex2ve6XSeGGVFJKx8mWveAE1\nx8LD48WSG8/gnNlrfI4Xb0ZE+B8vLMULXxzJzeWGckNa6p4SLGm0InoPpSR9N9kbPqBNm3cAQ4vj\n/kRYtwEUp1Zs8lFep06dCOl4PIV//OY13RiDo2U7z+fM/71ARN+RBLc+hocffpiiP7bTfsocnDnp\n/PHmXThiO3l1JK5MkCOMkHY9KkwvStlRefkWMbQ+fxoYgy0sEltknF+jULmKC3HlZWMLb+XzMiIi\nTVn5eJG9zr94YY+K9zteZH4xv0K8WLFico31PDxeDBgwgJB2+32OF4nPXkG3+5bXuJ0yruJC0tLS\nCG3d1+dlmiu9tFv8oQRLGrWYP19G1KDJLPtrT8a++BO28Jbseek6Qjoe7+Maqr5tl7dtHUUpO4gb\nfycAy5cvJ7zPMGwtorG1iCasy0kU7Py+xgTLXybIVmlC5quCnd+D5SLU52MgItL0lcWLkqw/3Ofx\nWowXBTu+I/yE4V7xYvny5WAbcvQ7Ur5Gh8ULm83PJyw7v8fpdCpeiASYBrmQRi8oOJQTTjgBW3hL\n8nd8R0nGbiL6ja12maSkJAp3b64ykXEVF5Cx4kVannUNQSEtPNOt4oJDZYryfe54XFecBw+Queo/\n2CJa0eKwl16KiDR3QcGhBMd38StelGSnBCReWA00XrRr107xopYUFDtrtbw0XHqCJY1W0f7fyd/x\nHcFturF8uY3MLxaQ/c0iogae5/XukYyVL4NlEdL+WGwtoinO2M2f/vR3MEFED76w0nVnffU2jlYd\nCT/uT55pw4cP58nnXsHRqiPO3AwKEn8i6tSJtb6fVXHmpLN+/XoKdm/BVZBD4d6t5P70CWDR+rwH\nCHKE1LgOEZHmoHy8ACjYvdnneJG1/j2/40Vol37kfLfUK16MHDmDhSuLa3dHq+DMSadwz6+lLxr2\njhcf/PcTLvwgrV7q1dQ1lD7rUveUYEnjFWQn//dvyfp6EROXOHFFdyB25PVE9B3hVSw4LoGcH/7L\nwY2f4SouICgsklHjx7A6/EwcsR0rrLY4PZmc75fR7opnvKZPnTqVOR9/T/p/Z2HswcSceSVhXU+u\n1V2szsFfPmPw4M8gyEZQSDiOVh2J7D+OyH5jsLWIrrd6iYg0OOXiBc5iHLGdfI4XoQknEnP6RX7F\ni+jT/oLzYJZXvBg5ciSs9K8PT6Ac/OUzDv5SebwYOHAgfFA/9RJpqpRgSaMVHN+Ztpc+Drjv+lR1\nlyii78gKQ+a+WU15R2wnEm55r+J6IiKIO/tmv+rY9uIZVc6LPHEUkSeO8pr22muvsaptzYGu811L\nger3uynzd2QmjeQk0ryVjxfVqSxeVKeqeBEUHFbr8SLu7Ft8Wm9ZvBCRuqMEq5nThac0Rmp2ISIi\nIg2VEqxmTsOOioiIiIgEjkYRFBERERERCRAlWHVszZo1XHnllfTp0we73U6XLl38Wn7Tpk1MmjSJ\n9u3bEx4eTu/evXniiScoKSkBDg3xmZ6ezk033cQxxxxDWFgYXbt25YYbbiA1NdVrfU6nkwNr32L3\nC1NIfOJc9rx4DdkbPgrIvor4oyQ7ldQPHiHp6QtIenoyKR88TEl2io/LppC27Cl2P/dXkp6cxJ4X\n/0bmFwtwFR0aJnnVqlUYY6r8Wb9+PQC7du2qttzbb79dK/svjUNycjLnn38+0dHRREVFMWnSJJKS\nknxa9t5772XkyJHExsZijOG1116rtNy8efM477zz6Ny5M8YYrrzyyirXabmcZG/4iL2v/J3EJyaS\nPOsi9r8qBa4lAAAWCUlEQVR9HyW5GUewdyJSXvm4FBUV5VdcAtiyZQupHz5K8r8vdseml64l+1vv\nayxnfjYZn81lzwtTqr1eA8jPz2f69On06NGDkJAQ2rRpw7hx4ygqKjrqfZXAUhPBOrZixQq+/PJL\nBgwYgDGGnJwcn5fdu3cvQ4YMoUOHDjzzzDPExcWxYsUK7rzzTlJTU3nssccIddjofNdS9r9xJ8UZ\ne4j506VED+pIdloyz7+6gJfe/4y2lz2BMe4XJo7KXEzWuneIPu0vhLTrRUHSz2R+/gqu4nxiTvtL\nbR0GES95eXnsf/tejM3h7rhtDAe+WMD+t+6l3V+fJSg4tMplXUUF7H/7fiyXk5g/XYo9Kp7CP7aR\nteZNSjL3Ej/hLgBOPvlk1q1bV2H5KVOmkJGRwSmnnAJAu3btKi13//33s2bNGkaNGlVhnjQPeXl5\nnHXWWYSEhDBv3jyMMdx///0MHTqUn3/+mfDw8GqXnz17Nv369WPcuHHMnz+fYqer0nKvv/46qamp\njBgxgvfeqziAQnlpS5+iYOf3RA+eTHDbHrgK8yhI3ohVogsukaPhKi7wiksvXXEKk6++yae4BFC4\n7zcGDrwIq81xxI6+kaCQcIoz92IV5XvKWJZF6qJ/eq7X3rpjEtu3beWBBx7g22+/Zd26dZ7rteLi\nYsaMGcPOnTu55557OP7449mzbz+rP1+B06n3ZzU0SrDq2NSpU5k2bRoAl156KWvWrPF52aVLl5KW\nlsbatWvp2bMnAGeddRa///478+fP57HHHgOgJHMvhXu20GrUDUT2Gw1AaEJfMIaMT5+jJGMPjtiO\nlGSn8PLLLxM9+EJPMhXW9SSsojyy171L5ElnYwuLDOTui1TqpZdeouTAftpf8wKOlu0BcMR3Ye+L\nfyP3x/9W+76xwj2bKcncS+sLHvIMmx/auS+u/Fyyv3kfV+nLPqOiohg0aJDXsomJiWzZsoXbbrsN\nm8092EtISEiFcnl5eXzzzTecc845tGzZMmD7LY3LSy+9xI4dO9i6dSvdu3cHoG/fvvTo0YO5c+dy\n6623Vrt8VlYWQUFBbN++nfnz5+OwBVXaB9Y6+R8YE8RnwEHXhyz8bjerSsuV7wN7cPNq8n79kraX\nP0VI2+6e6S16DAzA3oo0b7k/feIVlyZMOJv4/6b5FJcsy0X6sqc4e9gwfuh1tWd6aOe+XuUOv14b\nOXwYf/usAKv/hXz96XN0/NuLntcDZK1/j6x139B+ynPM2BUPu3LYNWMyF104uXYOgBwVNRGsY0FB\nR37Iyx4BR0VFeU2PiYnB5Tp0J9Ryul9kGBQc5r3tkNK7q6Vvky/cuw2Xy0XYMQO8yoV27Y9VUkT+\njm+PuK4i/li8eDEh7Xt5kisAR0xbQjoeT972r6td1nK6m8cGBbfwmh4UGu7+rltVL7tgwQIsy+KK\nK66odhvvv/8+OTk5NZaTpm3x4sUMGjTIk1wBdO3aldNPP52PPqq5abWv539jfCuX88PHhCb08Uqu\nRCQw8rd/fcRxqSBpI8XpyTXedPH1eg0g5/uPadHrDOxR8f7shtQTJViNyOTJk4mLi+OGG25g586d\nZGdn88EHH7BgwQJuu+02TzlHXGdCOvUh66t3KNz3G66ifAr3biXrq7cJPaY/jrhOAJjSYG9s3g8y\njc0BQHFaYh3tmTR3mzZtwhHXucJ0R1wCxWnV928J69IPe8v2ZK5+jaK0JFxF+eQn/kTOt4uJOGlM\ntc045s+fz8knn0yfPn2q3ca8efNo3bo1o0eP9m2HpEnatGlTpd+V3r17s3nz5jqti+UsoXDfVhxx\nncn8/FWS/30xiTMnsG/+reQn/lSndRFpiorSko44LhXudp8PCgoK2Df/NhJnTiB59iVkfDYXV3Fh\nuXV5X6/l5uZWer1Wkp2CMycVR0xb0v/7b5KenkziExMZNmwYP/74YwD3WgJFTQQbkTZt2rBu3Tom\nTJjAMcccA4AxhunTp3PnnXd6yhljaH3+dNKWPckf8w+9iDCs2ynETbjb89neyv3YuXDvVoLbdPNM\nL9z7KwCu/Nxa3R+RMhkZGYR1iagwPSg0EldB9d9DYw+m7SWPk/rhI+x75e+e6RF9R9JqxHVA5e97\nW7duHb/99huzZs2qdL1ly+zZs4eVK1dy0003YbfrlNmcZWRkVNpEtFWrVmRmZtZpXVz5OeAsIXfj\nZ9hj2hI7+gaMzUHWN++T8u402l46k5B2Peq0TiJNiSs/l6DQI4tLztx0AC688ELCjhtFyyFXULhv\nO1lr3qAkO5XWk+4HKl6vRZZesx1+vebMcQ9ak/X1QkLa9SB+/J1YzmJSty9hyJAh/PzzzyQkJARk\nvyUwdLVQSyzLqtDp8GgvzlJTU5k0aRLh4eEsXLiQ2NhYVq5cyb/+9S9CQkK46667PGXTl8+maO9W\nWo26HkdsJ4rTkzmw5g3SPnyU+PMfwJggguMSGD58OCvXvIE9ug0h7d2DXOSUjXBT2rFSpCGzSopI\nXfwYzrwsYsfd5h7kYu82sr56C4JsxI66vtL3vaV/8iwE2Zm5PY6nKukHU9bXZcGCBbhcrmpHchOp\na5ZV2izc5aT1+dOxR8YCENKpD3vmXk32N+97BngRkTpW2rTv0ksvZXEL98BIoQl9wXJxYPVrFKcl\ne55Olb9e++DeC5j4yLsVrtes0vUZRwjx5z1AkMPdMmPZrOvo3r07c+bM8fTDl4ZBCVYtWb16NUOH\nDvWaZlnVdAbxweOPP86uXbtITEz03EUdMmQITqeTqVOnMmXKFOLi4sj7fQN5W1bT+sJ/EdalHwCh\nnfpgj25LyrtTyd/+DS16uDvxv/baa3QbPJqU99wDb5jgFrQc+lcyPpmDLaLVUdVXxFctW7Ykr5I7\ngq6CnErvIJaX+/OnFCZtpP3fXsLRsh3g/r4HhbQg45NniTxpTIVlrJJi8n5dQ1i3AdhaRFe7/vnz\n59OvXz/69u1bbTlp+lq2bFnpk6qqnmzVJvffhcER28mTXIG7L0dI+2Mp2v97ndZHpKkJCo2o9EmV\nL3EpqHSAsBEjRrB47aHpYV1P4sDq1yhK+R1HXKcK12t//vOfifw4p8L1WtmAY6EdjvckVwCdOnXi\n2GOP5YcffgjAHksgKcGqJf3792fDhg0BXefGjRvp3r17hUB+6qmnUlxczPbt24mLi6M4dRcAIe16\nepULae/+XJyeDKUJVocOHWh78QxKctJxFeRgj2l3aPmOxwe0/iJV6d27N19tq9imvTgtCUdc9c0e\nilJ3ERQa4Umuyni+72nJFZbJ2/41roJcIvoMq3bdGzZsYMuWLTz99NM17YI0A71792bTpk0Vpm/e\nvJnjj6/b82WQIwR7TNuqC/g4UIaIVK6qvla+xKXK+m55Kx163cfrNXtMW4w9pMJayt596s8AapU1\nmZfAU4JVSyIjIxkwYEDNBf3Qtm1bvvrqKzIzM72SrK+/do9m06FDBwBs4e55hfu2eZ5ggbuvFYAt\n4tDdzjL2yFiIjMWyLLK//Qh7q46EJpwQ0PqLVGX8+PF8futtFB/4A0fpRWNJ1n4K92yh5ZnVj9xn\nC2+JqyCX4sy9XqM9Fe7d5p4fWfH7fvCXFQSFRRHW7ZRq1z1v3jzsdjsXX3yxv7skTdD48eO5/fbb\n2bFjh6cf7K5du1i7di0zZsyo8/q06DmY7O+WUJKThj0yDgBXYR6Fe7YQdkz/Oq+PSFPSovtAMj9/\n5YjiUtgx/cHm4JNPPoGIQ60o8nd8B0Bwaf9IX6/XjM1OWLcBFCRvwlVU4Bm8KWXfHn7etJnE8OMq\nfeVDZcq/6kFqj25x1bHU1FQWLlzIwoULSUpKIi8vz/O5/ChUq1evxm63M3/+fM+06667jvz8fEaO\nHMm7777LihUrmDp1Kk888QQTJ06kUyd3e94WPU/DFtGK9GVPkfPDxxQk/kzODx+TtuwpbJHxtOg5\n2LPO559/ntyf/0dB0s8c3LyalPemk//7BuLOvsXnoYJFjtY111yDPboNqYv+Sd5v68n77WtSFv0T\ne2QcEf0OBaeSrBQSHx/PQw895JkWccJwTHAYKe9NJ3fjCgoSfybr60Vkfv4KwW27V3gS6zx4gPyd\n3xN+/JkVRtAsr6ioiLfffpsxY8bQunXrwO+0NDrXXHMNXbp0YcKECXz00UcsXryYCRMm0KlTJ669\n9lpPucTEROx2u9f3FNzn9YULF7J8+XIAvv32Ww7+uoaDv3q/D7EoLckz3SopoiQ7xfM5NTXVUy7q\n1IkEhYaT8t50Dv66xv13s/BBrJJCogaeX4tHQqTpizhxlFdcWrx4cbVx6cDatzzTbGFRRA+azAsv\nvEDm6nnk7/rR/R6rr94mvM8wz83Aw6/XPv/88yqv16LPuASruICUhdPJ2/41B39dw9ixYwkKiSCy\n/zif96vsqVdtlRc3PcGqY5s2bWLyZO+XwpV9njZtGtOnTwcODZJx6zs/8MDmQ3fg4/4yg1++eouL\nr7oOV1Ee9qg2tBh0Id92nUiXu5exa8bZBIW0oO1lT3JgzZtkfb0IZ24GtohWtOh2KtFnXOz1vgWn\n00nW1wspyUohyBFCSMIJtL30CYLja3q8LRI44eHhtLnoYTJXvEza0icBCO18Iq2GXXPY+0EssFxe\n732zR7eh7WVPkrXmTQ58uQBXfja2yDgiTxxN1GkXVrhRcHDzKnA5Ca+heeCyZctIT0/Xu6/EIzw8\nnJUrV3LLLbdw2WWXYVkWw4YN45lnniEiIsLT9Kbs/F3+ewruc/zq1as9n+fMmXNo3ccu9fye9+uX\nZJW7WCtM2khh0kYANt08yjPdFt6Sthc/RubKl0n/7yywXIS0P5Y2F83QOVzkKAUFh3rFpUs+tWFv\n16fKuITl/fceffpFPDCpP3f98wmyv/kAW0RLok6dRPRpfzm0jcOu18aMeYWS0JhKr9eC4xJo85eH\nyVz9GmkfPQ5BNk4fM4LM02/yPAnzRWWDPlVHT7yOjBKsOjZkyBCfBrsoK3f4H0FIh2NpM/nBGpe3\nR8UTN/amGsvdcMMNPLG7a43lRGqbPao18RPvrb5MdBs637WU6dPP5rVyfxvBcQnEn3t3NUseEnXK\nuUSdcm6N5caMG+/XwDRq1948JCQksGjRokrnlb9w6XzXUl4rwOt7yqA76DzoDs/HXTPOrvRCJ+aM\nS4g545JKtzFkyBBYfmgZR6sOtD5/2hHsiYjUpHxcqurvtSwuHc4Yw6233sq/U3rVsI1D12tVbaNM\nSPtetL3oUc/nD2soL/VHCVYd0gWYSOOhu3wiItLcHcm1q653lWDVKV2wiYiIiEhj4e+1K+j6FTTI\nhYhIQBxJR2B1HhYREWl6jB99DI7uLbkC4PcTrMrKJz7m+2gxIkKl/RlrUlNb+KMtX7aMeDEBWEdA\nYtWRNnHx5zugc7lI81VZv626ijuKbUetxlilJoJHQW1MReRo+HsO0Tmn7qhZjIjIkVFsU4LlcaRf\nBvWpEpEjpXOIiIg0Nf7Gtl//OdrvbTT0pKzJJlj+Hnhd6IhIQ1cXozk19KAlIiJNy5G0GPA3Kavr\n2OZzH6wHH3xwORAX4O23B/YGeJ1NiY5PzXSMqqfjUz0dn5rV5TFKmzZtmv+3MsuppVjVHOlvo3bp\n+NY+HePa11yPcc2xyrKsevuZPn26VZ/bb+g/Oj46Rjo+Oj71/aNj1Dx/9P+u49vYf3SMdYzr80fD\ntIuIiIiIiARIfSdYD9bz9hs6HZ+a6RhVT8enejo+NdMxap70/167dHxrn45x7dMxroI/78ESERER\nERGRatT3EywREREREZEmQwmWiIiIiIhIgCjBEhERERERCZB6T7CMMf80xvxsjPnRGPOpMaZ9fdep\nITHGzDTG/Fp6jD4wxsTUd50aEmPMZGPMJmOMyxgzoL7r05AYY0YbY7YaY7YbY+6u7/o0JMaYV40x\nKcaYX+q7Lg2RMaaTMeZzY8zm0r+vm+q7TlI7ajpPGGNuLf0e/GyMWWGM6Vwf9WysfD0PG2POM8ZY\nimP+8+UYG2MuKHc+e7Ou69iY+XCOSCiNFz+UnifG1kc9G5p6H+TCGBNlWVZ26e//AI63LOu6eq1U\nA2KMGQmstCyrxBjzGIBlWXfVc7UaDGPMcYALmAvcblnWt/VcpQbBGGMDtgEjgN3ABuAiy7I212vF\nGghjzJ+BXGC+ZVl96rs+DY0xph3QzrKs740xkcB3wLn6/jQtvpwnjDFDga8ty8ozxvwfMMSyrAvr\npcKNjK/n4dK/sWVAMHCD4pjvfPwO9wDeBc6yLCvTGNPasqyUeqlwI+Pj8X0R+MGyrOeNMccDH1uW\n1aU+6tuQ1PsTrLLkqlQ4oGENy7Es61PLskpKP64HOtZnfRoay7K2WJa1tb7r0QCdCmy3LGuHZVlF\nwNvAhHquU4NhWdYXQEZ916Ohsixrn2VZ35f+ngNsATrUb62kFtR4nrAs63PLsvJKPyoG+cfX8/A/\ngceAgrqsXBPhyzG+BphjWVYmgJIrv/hyfC0gqvT3aGBvHdavwar3BAvAGPOwMSYZuAR4oL7r04Bd\nBfy3vishjUIHILnc593oAlmOgDGmC3AS8HX91kRqgb/niSkoBvmjxuNrjDkZ6GRZ1rK6rFgT4st3\nuCfQ0xiz1hiz3hgzus5q1/j5cnynA5caY3YDHwM31k3VGjZ7XWzEGPMZ0LaSWfdZlvWRZVn3AfcZ\nY+4BbgCm1UW9Goqajk9pmfuAEuCNuqxbQ+DL8RGRwDPGRACLgJsPa20gzYwx5lJgAHBmfdelqTDG\nBAFPAVfWc1WaOjvQAxiC+wnsF8aYEyzLOlCvtWo6LgJesyzrSWPMYGCBMaaPZVmu+q5YfaqTBMuy\nrOE+Fn0Dd/bbrBKsmo6PMeZKYBwwzKrvTnP1wI/vjxyyB+hU7nPH0mkiPjHGOHAnV29YlvV+fddH\naoVP5wljzHDgPuBMy7IK66huTUFNxzcS6AOsMsaA+0biYmPMePXD8pkv3+HduPsRFgM7jTHbcCdc\nG+qmio2aL8d3CjAawLKsdcaYUCAOaNZNMeu9iWBp58MyE4Bf66suDVHpo+w7gfHl2sGL1GQD0MMY\n09UYEwz8BVhcz3WSRsK4r/ZeAbZYlvVUfddHak2N5wljzEm4BxEar74rfqv2+FqWlWVZVpxlWV1K\nBwVYj/s4K7nynS+x7kPcT68wxsThbjK4oy4r2Yj5cnyTgGHgGXgsFEit01o2QPWeYAEzjDG/GGN+\nBkYCGg7Y27O473L9r3Qo+xfqu0INiTFmYmm738HAMmPMJ/Vdp4agdGCUG4BPcA9Q8K5lWZvqt1YN\nhzHmLWAd0MsYs9sYM6W+69TAnA5cBpxVet75UUPvNj1VnSeMMQ8ZY8aXFpsJRADvlX4PdKPGRz4e\nXzkKPh7jT4B0Y8xm4HPgDsuy0uunxo2Lj8f3NuAaY8xPwFvAlc2xtdXh6n2YdhERERERkaaiITzB\nEhERERERaRKUYImIiIiIiASIEiwREREREZEAUYIlIiIiIiISIEqwREREREREAkQJloiIiIiISIAo\nwRIREREREQmQ/wdiXjrAcOApAgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x112222a90>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"with pm.Model():\n",
" alpha = pm.Flat('alpha')\n",
" theta = pm.Deterministic('theta', pm.math.invlogit(alpha))\n",
" obs = pm.Bernoulli('obs', p=theta, observed=y)\n",
" trace = pm.sample(3e3, tune=1000, njobs=2)\n",
"pm.posteriorplot(trace[1000:]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Model 3: Log Odds with Jacobian Adjustment"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using advi...\n",
" 0%| | 0/200000 [00:00<?, ?it/s]Mean ELBO converged.\n",
"Finished [100%]: Average ELBO = -8.0688\n",
"\n",
"100%|██████████| 3000/3000.0 [00:02<00:00, 1135.67it/s]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAACsCAYAAABmdA06AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FNX+x/H32U3Z9E4nRBClhCKggHgFBJUmdq+iKGDD\n+/Nerw0LXdSLHS6iYlcsKGLHCwoiCAIKioJBFJEOSYAkJIQNSXZ+f2yyZEmHDWmf1/PkYXfmzMyZ\nYXe++51z5oyxLAsRERERERE5cbbqroCIiIiIiEhdoQRLRERERETER5RgiYiIiIiI+IgSLBERERER\nER9RgiUiIiIiIuIjSrBERERERER8RAmWyHEwxowwxiz3dVkREZETYYxJMMZYxhi/6q6LSH2lBEtE\nRESkFjPGbDXG9PfBepScifiAEiwREREREREfUYIlUgZjzP3GmD+NMZnGmCRjzKWllLOMMf8yxmwx\nxuwzxjxhjLEdU+ZJY0yaMeYvY8zAItNHGmM2FmxjizHm1qreLxERqRuMMbOBeOAzY0wWcFXBrGuN\nMdsLYtLYIuVtRWLbfmPM+8aY6ILZywr+TTfGZBljehpjWhljvi4ou88Y87YxJvLk7aFI7aMES6Rs\nfwJ/AyKAycBbxpjGpZS9FOgGdAEuBkYVmdcd2ATEAo8DrxhjTMG8FGAIEA6MBJ4xxnTx8X6IiEgd\nZFnWcGA7cJFlWaHA+wWzzgFOB/oBE4wxbQum/xO4BOgNNAHSgJkF884t+DfSsqxQy7JWAgb4T0HZ\ntkBzYFJV7pNIbacES6QMlmXNtSxrt2VZLsuy3gP+AM4qpfhjlmUdsCxrOzANuKbIvG2WZb1kWVY+\n8AbQGGhYsI35lmX9abktBb7EndSJiIgcr8mWZR22LOtn4GegU8H00cBYy7J2WpaVgztZuqK0+64s\ny9psWdZXlmXlWJaVCjyNOzkTkVLoJkaRMhhjrgfuAhIKJoXiboXKL6H4jiKvt+G+2ldob+ELy7Ky\nCxqvQgu2MRCYCJyG+6JHMLDeJzsgIiL11d4ir7MpiDlAC+AjY4yryPx8Ci76HcsY0xCYjvvCXxju\nOJXm89qK1CFqwRIphTGmBfAScDsQY1lWJLABd3eJkjQv8joe2F2BbQQC84AngYYF2/iijG2IiIgc\ny6pE2R3AQMuyIov8OSzL2lXKeh4tmN7Bsqxw4DoUo0TKpARLpHQhuINKKrgHowASyyh/rzEmyhjT\nHLgDeK8C2wgAAgu2kVfQmnXBCdVaRETqm2SgZQXLvgA8UnAREWNMnDHm4oJ5qYDrmHWFAVlAhjGm\nKXCvb6osUncpwRIphWVZScBTwErcwasDsKKMRT4B1gLrgPnAKxXYRibwL9w3JacBw4BPT6jiIiJS\n3/wHGGeMSQeuKKfsdNxx5ktjTCawCvdATFiWlQ08AqwwxqQbY3rgHuCpC5CBO7Z9WDW7IFJ3GMuq\nTKuyiJTEGGMBrS3L2lzddRERERGR6qMWLBERERERER9RgiUiIiIiIuIj6iIoIiIiIiLiI2rBEhER\nERER8ZHKPGhYTV0iIlKVfPFsHcUqERGpSuXGKrVgiYiIiIiI+IgSLBERERERER9RgiUiIiIiIuIj\nSrBERERERER8RAmWiIiIiIiIjyjBEqlBnLn5VVpeRESqx/Gcr3WOF6mdKvOgYQ19K3ISJNw/v8Jl\nt04dXIU1ETnpNEy71GmVOb8D/DZlAA5/e4XLO3PzK1VeRI5LubGqMs/BEhEREZGTxOFv10U3kVpI\nXQRFRERERER8RAmWiIiIiIiIjyjBEhERERER8RElWCIiIiIiIj6iBEtERERERMRHlGCJ1DJHUreR\n/N54tj99BTExMYwcOZIDBw6Uu9zatWsZMGAATZs2xeFw0KhRIwYNGsTKlStLLL9q1SoGDBhAZGQk\nISEhdOjQgTlz5niV2b59OzfccAPx8fEEBQVx2mmnMW7cOA4dOuSTfRUROVlq23MI851Z7P/ff9nx\n32Fsf/pykueMZf369ZVez5w5czDG0KxZs2Lz3njjDS6//HJatGiBMYYRI0aUXp/8fKZNm0ZiYiIO\nh4OYmBj69+/Pnj17Kl0nkdpOw7SL1CJ5mftJfvcB/GOaEXfJA0y/7DTuvfdehgwZwvLly7HZSr9m\nkp6ezqmnnsqIESNo3LgxKSkpPPPMM/Tu3Zvly5dz1llnecrOnz+fSy+9lGHDhvHOO+8QEBBAUlIS\nTqfTU+bQoUP079+f3NxcpkyZQnx8PD/88AMTJ07kjz/+4L333qvSYyEi4ku1aUh0y7JI/eAh8jKS\nie5/KzZHKBmr5tK3b1/WrVtXYrJUkvT0dP7973/TqFGjEue/9dZbpKamcv755zN37twy1zV8+HAW\nLlzIgw8+SLdu3cjIyGDp0qVecUOkvlCCJVKLHPz+QyxXPg0un4DNEcrVVw+mSZMm9O7dm48//pjL\nLrus1GX79etHv379vKYNGDCA2NhYZs+e7UmwMjMzGTlyJP/4xz+YNm2ap2z//v29ll2xYgV//PEH\nCxcu5IILLgCgb9++HDhwgCeffJLs7GyCg4N9tesiIjVKdT7U9/Dm1eTsSqLh1Y/iaNERgMCmbch6\nYzSPP/44//3vfyu0njFjxtChY0eaNmnCokWLis1fuHCh58LdggULgJL3e86cObz//vusXr2arl27\neqYPHTr0uPZPpLZTF0GpVSZNmoQxht9++40LL7yQkJAQ4uPjee211wCYPXs2bdq0ITQ0lL59+/Ln\nn396Lf/iiy/SqVMnHA4HsbGx3HjjjcW61z377LP07NmT6OhoIiMj6dGjB/Pne1/V3Lp1K8YYZs2a\nxYQJE2jcuDGRkZFcdNFF7Ny5s8r2//Dm1QS16obNEeqZdu655xIfH88nn3xS6fWFhIQQGBiIn9/R\nay1z584lNTWVu+++u8xljxw5AkB4eLjX9MjISFwuF5ZlVbo+IiIn4kRjROa6Bex+9Xa2PXkpO/47\njH1fTCf/cKZXmYNrP2PP7Ltp0jAOmyOUwCZtaHDlJBLun+/5a3bbqxhjiBlwO5G9rsEv1B1PUj6Y\nTN7BfSe8n9l/rMYeGu1JrgBsgSFcdNFFFY4FK1as4K233uL5557jg7U72Zvh9NqHhPvn0/LB/3le\n781w8sHanSUmlc899xy9e/f2Sq5E6jMlWFIrXXnllQwePJiPP/6Yrl27MmrUKB588EGef/55pk6d\nymuvvcamTZsYNmyYZ5n777+f//u//6N///58+umnPPHEEyxYsICBAweSn3+0L/3WrVu56aabmDt3\nLu+99x7dunVjyJAhnqt3Rf3nP/9h8+bNvPrqq0yfPp2VK1dy3XXXeZVxuVzk5eWV+1deQuLKzSEv\nPZmA2BbF5rVv356kpKQKHTuXy0Vubi7bt2/n9ttvB+Dmm2/2zF++fDnR0dGsX7+eDh064OfnR/Pm\nzZk8ebLXcerfvz+tW7fmvvvuIykpiaysLL7++mumT5/O6NGjCQkJqVB9RER87XhjxIGvnseR0JkG\nl48nqs9InH+tJWXuRCzX0XNffkYKoR0vZO7cucQNHUNg41NJ/WAyh7esLVaPjJVzyU3bQ8zAO5g+\nfTo5uzex7/MnvcpYlgvLlV/iX15e3tH3RWJE7r7t+JcSC7Zv305WVlaZxyc3N5dbbrmFe++9l1NP\nPbXCx7W0da1evZr27dszZswYYmNj8ff3p3v37nz99dcntG6R2kpdBKVWuvfee7n++usB6NatG599\n9hmzZs3ir7/+8rSo7NmzhzvuuINt27ZhWRZPPPEEEydOZMKECZ71nHbaaZxzzjl89tlnXHLJJQA8\n+eTR4OdyuejXrx+///47zz//PAMGDPCqR0JCAu+8847nfWpqKvfeey+7d++mSZMmAIwaNYo33nij\n3H167bXXyryB2OXMAiyv1qtC0dHRbNq0qdxtAFx11VXMmzcPgAYNGvDFF1/Qrl07z/zdu3eTnZ3N\nsGHDGD9+PF27dmXRokVMmTKF9PR0nnnmGQAcDgfLly/n8ssvp3379p7lb7rpJp599tkK1UVEpCoc\nb4yIOPsaIntd41mPX3RTkt8ew+HN3xN8Wk8Aos67EXB3u3Z8mY0joTO5B3aT+dMXBLX0bsHxi2hA\n3NB7AbjhhsH8+/VlpH/zKnmZ+/ELiwFg/xfTObRhcYn74f/E0dcxg/5NaAd3V22XMxO/iAbFykdH\nRwOQlpZGaGjxWFHoscceIycnhwceeKDUMhW1f/9+jhw5wuuvv07Lli156aWXCAwM5IknnmDAgAF8\n9913dOvW7YS3I1KbKMGSWmngwIGe11FRUTRo0IAzzjjDq7tamzZtANixYwcbN27E5XJx7bXXkpeX\n5ynTvXt3wsLCWLZsmSfBWrt2LRMnTuSHH34gNTXVc9Xw9NNPL1aPQYMGeb3v0KED4B5drzDBmjRp\nkqelqCynnHKK53XRq6UAxua7fv6PP/449913Hzt27GDmzJkMGTKERYsWeQKgy+XC6XTyyCOPcNdd\ndwHQp08f9u/fz8yZM5k0aRIRERE4nU7+/ve/k5KSwuzZs4mPj+f777/noYcews/Pj+eff95ndRYR\nqYzjjREh7fp4nX8Dm5yOCQjCuWODJ8HK2buZjOVv0/C1UaSkpALuGOEXXXxgiaBW3olFQJy71Sn/\nYKonwYo8ZxhhXYaUuB+f/fMcLpqx3L3+yIbl7nduvqvcMps3b+aRRx7ho48+wuFwlFu+PC6Xe5u5\nubl88cUXnth37rnn0rJlS5544gkNeiT1jhIsqZWioqK83gcEBJQ4DcDpdJKSkgJQaleI/fv3A+5A\n269fP9q1a8eMGTOIj4/Hz8+P8ePHs3HjxmLLFV4tLBQYGOjZZqH4+PgKjehkt7uTqLyMZHa9cKPX\nvKajX8EWHAGYgpYsbwcOHChWl9K0bNmSli1bcuaZZzJkyBASExMZN26cpwtkTIw76J9//vley11w\nwQW88MIL/Prrr5x99tm88sorfPPNN2zevJlWrVoB7oAaERHBLbfcwujRo+nUqVOF6iQi4kvHGyN2\nv3gzJXE53fdh5R1MJXnOWAJimjNjxgzu+GwH2Oykf/sWuft3FFvO5gjzem/8/AGw8o94ptnD47CH\nxZa43c6dOxPQcHfBwkfv6rA5QkuMBZkZ6QD0emY1toBfSlxn8tyJmCaJ3LIgAxa8xy+TLsTKzwMs\n9zrt/tj8A0tctiRRUVEYY2jXrp0nuQIIDQ2lZ8+e/PTTTxVel0hdoQRL6oXCpOHLL78sFmSLzl+w\nYAEZGRm8//77XklRdnb2cW+7sl0E7aHRNLr+Ga959tBojN0fv4gG5O7bXmzZpKQkevfuXem6BQQE\n0LFjR9atW+eZVrS7X0kKR5Rav349UVFRnuSqUOFohBs3blSCJSK1QmEMaHDVlBK7YduC3InS4S1r\nsXIOEXvx/Vx11VWM+dE9AJKVl3Pc2z6eLoL+sfE4/yqeuCQlJWEPj8MWEFTq9nL37SD/YAo7pl8N\nQNT0o/N2TL+asK5Die5/S4XrHxQURMuWLUudX9bjQ0TqKiVYUi+cf/752Gw2tm/fXqxlpqjCRMrf\n398z7ffff2fFihUVfq7IsSrbRdDY/Qls3LrEMkGndufQhsW4cg5hCwzBmZvPmtUr2bZtW4WGwz12\neN3s7GzWrFnj1f3xkksuYfz48SxcuNDT5RHcyafD4SAxMRGARo0akZaWxubNm71aBlevXg1A06ZN\ny62PiEhNUBgj8g6mEnbKGaWWK0ykjP3oeTT3wC5ydiaV2gpVnuPpIhh8ancOrV+Ec/t6HPHu87Qr\nJ5vPPvuM4FPPLnN7cUPHeLWgzbmlJwNH3c2RvZuJu+T+49qPSy+9lBkzZrBr1y7PuT8zM5Pvvvuu\n2L3LIvWBEiypF1q1asV9993H7bffzqZNm+jduzcOh4MdO3bw1VdfcdNNN9G3b1/69++Pn58f119/\nPXfffTd79uxh4sSJxMfHe/qZV1RhMpOQkEBCQoJP9iO8+2UcSvqGlHlTiOhxJZ98mMW1t/yTgMan\nc9f3Adz9g/tqqnP7epLnjCVm0B2EJrqffbV/wbPcPrAz3bp1IzY2lm3btvHss8+yZ88eZs+e7dlG\nYmIiI0aMYMKECbhcLrp06cKiRYt4+eWXGT9+vOfG6REjRvD0008zaNAgxo4dS3x8PGvWrGHKlCl0\n7dqVXr16+WSfRUSqWmGMmPrEU+Qd2ElgfAeM3Z/8zH0c3voTYR0vxNGiI44WncFmZ9/nT/PlRQ3I\nWr+Y9OVv4xced9yPpvCLaIhfRMn3V3Xr1o3AxsnFpge17k5gkzbs+/wpovqM9DxoONiyCO9+uVfZ\nbY8PJSSxH7GD7gDcz8sqqk+fPthDojB2fxzxHb3mHdm33dNrwso7Qt7BFD744AMAevfuTVxcHAD3\n3HMPs2fPZuDAgUyYMIGAgADP8xDvv//+4zgqIrWbEiypNx599FHatm3LzJkzmTlzJsYYmjdvTr9+\n/Wjd2t1i1L59e95++20mTJjA0KFDadWqFVOnTmXBggV88803ldqew99Owv3zyy9YxNapg8uc7xcW\nS8NrHiVt8cukfvwoo78KwhHflajzbsSYot0wLLBcUCTgBzY5jaVLl/Liiy9y6NAhmjZtSvfu3Xnl\nlVe8WqoAZs2aRdOmTZkxYwbJyckkJCTw9NNPc8cdd3jKJCQksGrVKiZNmsS4cePYt28fzZs355Zb\nbmHs2LHqFiIitcqjjz7Ki7/kkPnjfDJ/mg8Y7GGxBCV0wi/KfW9RQFwLYofcQ/rytxk6dCiu0IZE\n9b6Bw3/9iHP7+pNWV2NsxF0xkbQlr3Dgq+ex8nIJbHo6S5Ys4eJ3j3kWo+Vy/x2H7N++JWPFu573\nOdvXc+WVVwKwZMkS+vTpA0DDhg1ZtmwZd999NyNHjsTlctGzZ0+WLl1abrdzkbrIVOKKi54aKlJJ\nx5NgVWaZ4ykvUoMZH6xDsUqOW1Wff2taTDjebYjUc+XGKl1iFhERERER8RElWCIiIiIiIj6iBEtE\nRERERMRHlGCJiIiIiIj4iBIsERERERERH1GCJSIiIiIi4iNKsERERERERHxECZZIBTlz86u7CiIi\nIiJSw/lVdwVEaguHv73WP4zRmZuPw99e5cuIiIiI1FdKsETqkcomiVAzE0URERGRmkpdBEVERERE\nRHxECZaIiIiIiIiPKMESERERERHxESVYIiIiUqdo1NeqczzHVv8fUt9okAsRERGpUzSgT9XRsRUp\nn1qwREREpEZTC4iI1CZqwRIREZEarS48h1BE6g+1YImIiIiIiPiIEiwREREREREfUYIlIiIiIiLi\nI0qwREREREREfEQJloiIiIiIiI8owRIREREREfERJVgiIiIiIiI+ogRLRERERETER5RgiYiIiIiI\n+IgSLBERERGpMs7c/CotL1LT+FV3BURERESk7nL420m4f36Fy2+dOrgKayNS9dSCJfWSro6JiIiI\nSFVQC5bUS5W9mga6oiYiIiIi5VMLloiIiIiIiI8owRIREREREfERJVgiIiIiIiI+ogRLRMqk4XVF\nREREKk6DXIhImTS8roiIiEjFqQVLREREThq1cotIXacWLBERETlp9JgMEanr1IIlIiIiIiLiI0qw\nRERERKTG0OBKUtupi6CIiIiI1BiV7Ub625QBld6GMzcfh7+90suJVIQSLBERERGptXRfn9Q06iIo\nIiIiIiLiI0qwpE5Q/2sRERERqQnURVDqBD0MV0RERERqArVgiYhPHU9rologRUREpK5QC5aI+JRu\nNhYREZH6TC1YIiIiIiIiPqIES0RERERExEeUYImIiIiIiPiIEiwREREREREfUYIlIiIiIiLiI0qw\nREREREREfEQJloiIiIiIiI8owRIREREREfERJVgiIiIiIiI+ogRLRERERETER5RgiYiIyHFz5uZX\ndxVERGoUv+qugIiIiNReDn87CffPr3D5rVMHV2FtRESqn1qwREREREREfEQJloiIiIiIiI8owRIR\nEREREfERJVgiIiLioUErpD6o7Odc3wupDA1yISIiIh4atELqA33OpSopwZJ6wbntF9K/fYsjyZsx\nfgEM33UJ+ZEXYg+J8pTJy0hm1ws3lrh88zvmeF67cp0c+GoWh/9Yic0RSuS51xPS9lyv8hmrP6BT\npwexLpyCsdnLrd/ed+4HlwtKOIFn/ryQAwtm0HT0K/hFNARg3/xnOLRhMeYxdxlbUDj+Mc2J6HkV\nQS27epbd9tiQoyuy2Yl7JYpMRwMcCWcQ1nkg9pDIcusmIlIfLFmyhL1vjfHEiaBWZxLVd1SpcaLw\n/Fuo+R1zsDlCgZLjBHif3zNWf8ChX7+h8YjplYoTja57vNi8wjixdfRfnmmFcaJQaXHCGHN0RTY7\ntsAQ/GOaKU6InAAlWFLnOXdsIPn98QSd0oW4Sx4g/3Amy5Z9QHL2MhrfMB3j5+9VPrzHlQSf2t1r\nmgkI8rw+uOoDnFvXETPoTo6k/sW+z58moGEr/KObApB3cB8Z373H80sWMezTtCrbL1twBCsWL+DS\nmSvIP5TGwR8+ImXuJBr8fQpBCZ095UIS+xPWeQCWZfHS30/nuv+8ReaPn5O59jPiLhuPo1nbKquj\niEht4NyxgQsuGId/izM8cSL929kkzxlbapxY+MxdXDpzhWdaeXHijz9u8swvjBMNrnyoQsnV8bIF\nR9DgsvEAFY4TLudBcnZvUpwQOQG6B0vqvIwV7+IX3oC4y8YR1OpMQhPPY968eeTu207WL18WK+8X\n2YjApm28/ooGwMNb1hLWdTDBrbsTefbV+EU2wrntZ8/8tMUvEnz6OZx99tlVul/G5kePHj0IbNqG\n4NN60uCKSZjAYDLXfupVzh4WQ2DTNjiateWiiy4i6tzraTLqWWyOUFI/fgRXrrNK61kR6gsvItUp\nY8W7tGjRwitOxF3yYJlxovD8W9E4sWjRIs/8wjhR1YmLsfl56lfROBF8avcaGSdEahO1YEmN48zN\nx+Hvuyt6Obs3EdK+r1fw69atG7agcLJ/X0lYl8r1q7byczF+gZ73Nv9ArLwjgDuoOndsoMnNL/im\n8pVgCwzGP7oJeWl7yi1rD4kiqu8oUj98mENJywjrdMFJqGHp1BdeRKpTzu5NnH/jCP5XJE4ENm7t\n0zjhdLqTFMWJ2ul4fpv4+veM1B5KsKTGqeyPbSjnB7exYezFP+rG7kfuvm3FpqcvfYMDC2di/B04\n4hOJPPd6AuISPPMDm5zOoQ2LCT69F7kpf3Ek5S+izx+NlZfLgUUvENn7BuxB4ZWqf6G8vDws1zGt\nM5ZVoWUtVz75B/d57tMqjyPhDLDZydmVVO8Dp4jUc8ZGQEBA8cllxAk/v+dw2QMrHCd69OjBtHnJ\nJxwnisUIUJw4CXz+20TqNCVYUuf5xzQlZ/cmr2nbtm0jPysN7EWuLNn9Ce08gKCELtiCw8ndv5OM\nVXPZ+9a9NB7+tKdYRK9rSJk7iV0zrwcg/KzLCGzalvQV72IPiiC04/EFoZxdSfj7+5dfsIjChCw/\nK42M7+aQfyiN8O5XVGhZm38g9qBw93EQEanH/GOasmrVKuh79Pydl5FSZpyYd9cALnnkfa844R/b\nHCg5TvTs2ZOMMdedcJzY/sTFlVqmMCFTnBA5eZRgSZ0X1nUo+z9/irRlswnvehEuZybDhz8GxmDM\n0dsQ/UKjibnwds97R/NEglp2Zfcr/yBj5XvAaHe5sFgaj5xBXvpebI4Q7EHh5Kbv5eD3H9Lo2sew\n8nI48PXLNHxtFPtzDOFnXkJ414vKrad/g1P4bv77XDRjudf0w3+sKti+t/ys/V4JmQkIIuKcawnr\nVv62CllYYMovJyJSl4V1Hcr3nz9FuP1onNi/4Nky48Tf/vY3wjof9IoTsRfd4y5XQpzYsmVLsTiR\n/fsqjH9gpeJEzIB/FZteVpwompApToicHEqwpM4Lbd+XvP07OfjDRxxc+R5gOPfqvxPUqhu5qcW7\nfhTlFx5HYLN25Oz9w2u6MQb/qMae92lfvUBoxwsIaNCStGVvcmTvZrZs2ECnMXPY+859+Mc09xqx\nqSQ2/yC6detGYONkr+lHUraUXD44ktVLv+KiZ1dgDwrDHhZbqdGoXLk5uLIPYg+JrvAyIiJ1UWj7\nvtzWOYhHpj7uiRPBbf/m0zjxz3/+s1icaHLjTPIz91cqTgQ2bl1sellxosEVE8EYxQmRk0gJltQL\nkecOJ7zHleRl7MUeHMG7M67DP6Y5gc3aVXANpV++y/59JUdSthA7dAwAzi1rCenQn7i4OAIatiQo\n4Qycf/1YbuCsLGOzuxOyRsnlFy6B868fwXLhqPAxEBGpu6ZMmcIbhzp74oQ9JIpdL432WZxYt24d\nkVdNA47GCXtwBPbgiCqNEyUlZBWlOCFyfDRMu9QbtgAHAXEJ2EOiWLBgAXkHdhLaeVCZy+QdTCFn\nZ1KpAcqV6+TA4heJOu9mbIHBnulWkSFtXUcOV/gG5JMl/1A6ad+8hj00muBjHpIsInWHHmlQOUXj\nxOEta30aJ5555hnFCZF6Qi1YUucdSf6Tw1vWEtCwFQDOnUlc/MxHhHe/3OsZJAe+fhksi8AmbbAH\nR5B7wD3IBcZGRM+/l7jujO/m4B/djJC2f/NMcyR0JnPt53z88cdk/vglzm0/E37WpVW7k2XIz9xP\nzq7fsCyLzz+3SFv2Flk/LwQsGlw+AZt/YLnrEJHaSSOfVcyR5D959NFHObwlD3DHiYPfzyszTixZ\nEkzmT19UOE5cddVVjPnR/X9RGCf8o5uRn3WgRsUJlzOTnN2bFCdEToASLKly1f4cCJsfh/9cQ8bq\neZCfi39Mc1544QUmb2rgVSwgNp7Mn/7HofWLcOU6sQWF4YjvRGSva/CPaVZstbn7d5D543wa3zDN\na3rE2VeTfyiDUaNGcTDXENl7BEGndKnSXSzLoQ2LOLRhEdjsjFgQSY6jIWFdhxDWeSD24Ihqq5eI\nSI1h8+OLL74gdc06T5yIueD/CO14vlexonHigi+ewhUYekJxYv//pmP8AmpUnLAFhuAf3Uxxwgcq\n+/un2n8FOb6KAAAXZklEQVQvic8owZIqV90PkQ2Ia0Gj6x73mjZy5GAmH1On0I4XVGroXP+Y5sTf\nObfYdFtAELGD/83WqYMrvN+Nhk0tdV5YpwsJ63Sh17TYwXdWaL0t7vvc631l6iQiUl8ExLVg+SvL\nyz0/Fo0TFTmflhcnKqMicSIhIQH4Fah4nLAsS3GhilT37x+pProHS0RERERExEeUYIlIrVPZG/d1\no7+IiIicLOoiKCK1jrpdiIiISE2lBKuOWb58OS+//DJr1qzht99+o1mzZmzdurVS61i1ahWTJk1i\n1apV5Obm0rJlS8aOHcvVV19dYvmpU6fywAMP0KtXL5YvX15s/q5du9j3xTQOb1mDy5mFPTSGkLZ/\nI6r3iOPYQ6kr8g6mkrb4JQ5vXUf4c3byGicS3e9m/MIblLlc+vK3MWZIyTPt/rS45yPP2/zDB8lY\n8S4tW/6TPXv20KhRIwYPHszEiROJi4srcRVbtmwhMTGRw4cP88cff3Dqqace9z6KlGbHjh3ceeed\nfPXVV1iWRf/+/Zk2bRrx8fHlLvvggw+yZs0a1q5dy4EDB3jttdcYMWJEsXJ9+vRh6dKlxaZHnXcz\n4WdeDIArJ5uDaz7BuWUtuQd2YVkuzl7WkexG5xF8Ws8T3k+R2qZobAILR0JnovvdXOHlc/ftIH35\nWzi3rydo+hHygmMIO2MQ4d3c37ms9YvY/8W0Epc1j+GJVYXS0tKYPHkyH374IcnJycTFxdG/f39e\nf/31E9lNqWJKsOqYxYsX8+2339KtWzeMMWRmZlZq+fnz53PppZcybNgw3nnnHQICAkhKSsLpdJZY\nfsuWLTz88MM0aFDyj+KtW7fSq1cv8kwk0f1uxR4SSV5GCnlpuyu9b1J3uHKdJM95EGP3J3bwnbx0\nw5lcedMdJL/7II1HPostwFHqsqEdL2TBk//i0pkrvNaXMnciwad290yzLIvUeVPIPbCLx56aStu2\nbUlKSmLChAmsWbOGlStXYkzxB4P+4x//ICIigsOHD/t2p0UKZGdnc9555xEYGMgbb7yBMYZx48bR\nt29ffvnlF0JCQspcfsaMGXTu3JkhQ4bw5ptvllm2Y8eOpHS83muaX0RDz+u8gylk/jSf0A79iTj7\najCG04K3sPKNR4g+fzRhXUq5mCFSBx0bmzCG9GWzSX73QQ49cVW5y+fs+YPkOQ/iiO9AzIB/8t4/\n+3HFYx9iHTkaT4JanUmj6548ZkmLlHlT6Jp4erHk6pxzzsEYw8MPP0xCQgK7d+9mxYoVSM2mBKuO\nGT9+PBMnTgTguuuuK7FFqTSZmZmMHDmSf/zjH0ybdvTqSv/+/Utd5rbbbuPaa69l06ZN5OXlFZs/\nevRomjZtil/vsRi7Pm7ilvXzQvLSk2ly8wv4RzXh4osHE/e/fex+8Ray1v2vzOfB+IXH0qNHDwI/\n3n90fRu+Blc+IR36eablpe0mZ9dGoi+8ndtuuw1wX9G32Wzcdttt/P7775x++ule637nnXf46aef\neOCBB7jzzoqNwCVSWS+99BJbtmxh06ZNnhbSjh070rp1a2bNmsVdd91V5vIZGRnYbDY2b95cboIV\nFhZGRtM2pc73i2hE01tfxuZ/9KLG61Mn8N43P5Gxap4SLKlXjo1NAP5xCex+8RZmzZoFnF7qspbl\nYv/8p3G06ESDy8YB0LdvX8IWZnuVswdHFBv63rljA67DBxl23XCv6Q888ABZWVmsX7+e8PBwz/TC\nHkUa1r3m0i/eOsZmO/5xS+bOnUtqaip33313meUKv9DvvPMOP/74I++++y6XXXZZsXJ//vknCxcu\n5M0332T8r/qoyVGHN68msMnpngAG4B/ZiMBm7cjevLrSD9w8tGExtpBIr+fIWPm5gHs45KIiIyMB\ncLlcXtPT0tK46667ePLJJ8nP16AYUnU+/fRTevTo4dX99JRTTqFXr1588skn5SZYJ3KeL7auUlqL\nAxq1xrl9g8+2I1IblBWbPvnkE+g5ptRlndvXk7t/B9EX/l+lt3tow9dg9+P666713F/sOuJk5yuv\nE97jCjo++m2Jy+n+4ppLv3rFY/ny5URHR7N+/XoGDRrExo0bady4MTfddBPjxo3DbndfJXH422n+\n7/fY/dL/EdX7Bro8vpK9W/aDy+U18EDWhsUABAUFkTxnHM6dG7D5BRJ06llE9bsZe1B4ifWQuu/I\nvu0En9qj2HT/2Hiyf6t4qyu4+8s7t68nrNtQjO3olTz/2BYENk8k47v3WLPmKtq0aUNSUhIPPfQQ\nAwcOpG3btl7rGTNmDG3atGH48OHq2y5V6tdff+Xiiy8uNr19+/bMnVv8mUkn4qeffuLwD1dh5eXg\nH9OcsK5DCetU/vP+nDs2lPjgXJG6rKzYlJT0PSFl3JaYszMJACsvlz1v3s2R5M00eDUaZ3wPInuP\nwOYfWOJyrtwcDv22nOBWZxEdHX20LsmbsfJysIdEkfrRoxzesgaMDUdCZ6LOuwn/yEYlrk9qBg3T\nLh67d+8mOzubYcOGMWLECBYtWsQNN9zAlClTuOeee7zKpi95Ff/oJoR0KL37YH7WAQBGjRqFX3RT\nGlw5mcg+Izn85xpS3p+AZblKXVbqNtfhLGyO0GLTbY4wXM6sSq3r0K9LwHIRmtjPa7oxhgZXTMIv\npilnnnkmYWFhdO/enZYtWzJv3jyvst9++y1vvvkmzz33XOV3RqSSDhw4QFRUVLHp0dHRpKWl+Ww7\n5557LtOmTaPB5eOJu+QB/KOacGDBf0n/bk6Zy7344osc2b2J8B5X+qwuIrVBWbGpvO9mfpa72/q+\nTx8j6JQzaPj3KYwZM4asn79k32dPlLrc4T9WYR3JJuSYGJaf6V5f2pJXwGYn7rLxxAy4nSPJf5L8\n7gO4crJLWp3UEGrBqqUsyyrWjcnP78T+O10uF06nk0ceecTTRaVPnz7s37+fmTNnMmnSJCIiIvj2\n22/J2vA1jUdMK3GQgCKV9Kzjl3a3uqe16IQtMJh9nz6Oc8uPBLXqdkJ1Fsna8DUBDVsR0OCUYvP2\nL5jBkd2beOGFF5i4LJ3c/Tv43zdvE5N4LnFXTMAYG1Z+Lnte+xdBZwxl0Jt/AX+Rtf7nk78jIj72\n0EMPAfDIn+6eBcGte5Dy4cMcXPk+4d0uLtZ9FsC5/Rf+NW0yIYnnEdq+70mtr0itVvCbJ6R9XyL/\ndh0A99wzmEc+/5X0pa+Tu28H/rHNiy2WtWExtuDIEn4PudfnF9mI2KFjPL+3/CIbs3f23RxK+gbQ\nRZCaSi1YtdTSpUvx9/f3+jtRMTExAJx//vle0y+44AJyc3P59ddfAbj11lsJ7Xg+fmGxuJxZ7hYH\nlwvLysflzMLKK7j3JSisxPU5Cu6TOZKy5YTrLLWTzRFaYkuVy5lZ4tXD0uTs3kTegZ2EJJ5XbF72\nnz+QvXEpMYPv4tZbb8XRPJGwzgOJHXw3h7es4fDm7wE4uOYTXM4swrpe5Pk8W7k5gHvgl8qOxClS\nnqioqBKvhpfWslXUiT40O6Rtb6y8I+Smbi02L2fP76TMm8J5551HzIB/ndB2RGqjsmJTed/Nwt88\njoTOXtODTjkDgCMpfxZbJi/rAM6t6whp19uri7u7LgXra9HJ62J2YJPTMQHBHEn+87jOByd6DpGK\nUQtWLdW1a1d++OEHn66zffv2Zc4vvLF648aNwEay1v2vWJkd06/2PGPFP7ac57mU1foldZp/bDy5\n+7YXm567b3v5n5sisjZ8DTY/Qtr1Kb6ugh+QgY1P85oe2MT9Pnf/Dmjdg9x9O8g/lMau524oto4u\nXbrQqVMn1q1bV+E6iZSnffv2ngtWRSUlJdGuXbsyly36kO3cgsdd3DP3Zyb9VvKDt0u/Cd77/Hsk\ndSsp708goIG7C23byV+XsxcidU9ZsalLu3ZsK3PZFuWsvfhvntK6uFdofcZ4nQ8qSgNjnBxKsGqp\nsLAwunXzbfe6Sy65hPHjx7Nw4UI6dOjgmb5gwQIcDgeJiYkALFmyhKtfXOm17IHFL4HLRfT5t+IX\n6R59J7BJG+whUSxcuBASb/OUdW5ZC7hHqZL6KfjU7qQteYXc9L2eG3XzMpLJ2bWRqN7FE52SWPm5\nZG9cRlDLrsWGvAWwh7ivNubs+d1res7uTe75oe4W24geVxDawTu4Hd7yIwdXf8Bbb71VbCh3kRM1\ndOhQ7rnnHrZs2ULLli0B9zMDV6xYwdSpU6t024eSvsH4BeIfl+CZlntgF8nvjcMvshENrphIUFDx\nroMi9UFZsWnoHY8xI7X0ZYNadgW7P86/fvR6JuPhwt88jYv/5jm04Wv84xIIaNiy2Dy/8Fj3aJ5b\n12FZlqcVK2fXRqwj2QQ0Oq3YMlJzKMGqY1JTU1m6dCkA27dvJzs7mw8++ACAdu3aea6OLl26lH79\n+vHqq69y/fXuh1AmJiYyYsQIJkyYgMvlokuXLixatIiXX36Z8ePHExrq7rrVp08fHAsOeW3XFhgC\nLheO+I6eacZmJ7L3CObPf4bQXfkEn3Y2eWl7SF/2JoHxHXC06FTlx0NqptBOF5L54+ekzptC5LnD\n+fTTfFLmTcEvLJbQzgM95fIyUtg16yYiel1DZK9rvNZxePMPuJyZJV75Awg+7WzSl73J/vlP8/zz\n0Ti3pZN7YCfpK97BHhZH8Gnu4aD8Y5rjH+PdLz4vIwWA7t27ew2lLeILN998M88++ywXX3wxDz/8\nMMYYxo8fT/Pmzbn11ls95bZt20arVq2YMGECEyZM8Ex3bl9PfnYG+Yfc3QyP7P2DQwXPsQppc467\nzI4NHFz1Aa+03svhrXuxcrLJ2rCYw5tXu0c0KxiePf9QOsnvjYf8PCLOuZbcfdtZtWoVObt+AyCg\nYSuM34l3QRepDY6NTWBI/3Y2fmGx3Hrrrcx42P37qqTYZA8KJ6LHlWR8NwcTEIyjRSemTl1Pxndz\nCEns5zX0O0DO3s3k7ttGVN8bS61PZO8bSHl/Avs+/g+hHS8g/3AG6ctm4xfdjJB2vavsOMiJU4JV\nx/z6669ceaX3TY+F7ydOnMikSZOAo4NkHPssoFmzZtG0aVNmzJhBcnIyCQkJPP3009xxxx3HVZ/Q\nDv145uozuPGucWStX4TdEea+AbT3DWUPkCF1mi3AQcNrHiFt8cvs+/wprv3Sjl/jRKL73XzMjfcW\nWC733zGyNizG5ggj6NQzS95GYDCNhj9F+vJ3ePzxx0nesQt7aDTBrc4i4pxhJd7gL3IyhISE8PXX\nX3PnnXcyfPhwLMuiX79+TJs2zXMhC0o/T6cvf5ucHUefUZX543wyf3R3Ewpp8zkA9tBoLMtiwoQJ\npCSnYGx+BMQlEHvRvV4/zHL3byf/oPuCQuoHkwHo+dbRbTUd/Qp+EQ19ewBEaqhjYxO474GK7nez\n13eztNgU0esabAFBZP70BQe//4jnVzYh/KzLiDj76mLbOrRhMdjshLTvU2p9ghI60+DyCaQvf5uU\njx7B5u8gqFU3ovqOKnXYd6kZlGDVMX369MEqGMmmIuWOvdkxICCAhx9+mIcffrhS2200rPRuLcOH\nD2f8r9Glzpf6yS+8AXGXPgi4+4SX1I/cL6IhLe77vMTlG1w+vgLbiCN20B38Vcr6SxPaoT/75j9T\n4fIilRUfH1/scQHHSkhIKPF8Xtb5tpB/VBMaXjW51O9WIUd8x2LfsfKWEanLisamUsuUEpuMMYSf\ndSnhZ10KlP1diu5/K9H9by1xXlFBrbppxOVaSAlWPacbJEVEREREfEfDtIuIlKCyQ9lq6FsREanp\nFNtODrVgiYiUoLKtu2rZlcpy5ubj8LeXX1BExEcU204OU5H7dQpUuKBUn+MJ2Md+0bY9NsSXVRI5\n6XxxT0lll/ltyoBKf/f0A7sYX4x8U6ti1Yl+LnW+Fqk9SrunuFBl445iW7UpN1apBauO0ZUJkeqh\n+xlFRKSuUWw7ProHqwZTv1cRkdpD52wREd3nBWrBOqkq22SqqwYidVtlzwnH0+2iHnTVqDHUg0BE\nROdCqMQ9WJMnT14AxFZtdYppAuw+ydusLXRsSqdjUzodm7Lp+JTuZBybfRMnThxwIiuoplhVl+k7\nUbV0fKuWjm/Vq4/HuPxYZVlWjf2bNGmSVd11qKl/OjY6Njo2Oj46NvrT/3vt/tPx1fGt7X86xiX/\n6R4sERERERERH6npCdbk6q5ADaZjUzodm9Lp2JRNx6d0Ojb1k/7fq5aOb9XS8a16OsYlqMxzsERE\nRERERKQMNb0FS0REREREpNZQgiUiIiIiIuIjSrBERERERER8pEYnWMaYKcaYX4wx64wxXxpjmlR3\nnWoSY8wTxpjfCo7RR8aYyOquU01hjLnSGPOrMcZljOlW3fWpCYwxA4wxm4wxm40x91d3fWoSY8yr\nxpgUY8yG6q5LTWKMaW6MWWKMSSr4Pt1R3XUS3yvv3GCMuavgM/CLMWaxMaZFddSzNqvo+dcYc7kx\nxlLcqpyKHF9jzFVFzmXvnOw61nYVOE/EF8SLnwrOFYOqo541RY0e5MIYE25Z1sGC1/8C2lmWNbqa\nq1VjGGMuAL62LCvPGPMYgGVZ91VztWoEY0xbwAXMAu6xLGtNNVepWhlj7MDvwPnATuAH4BrLspKq\ntWI1hDHmXCALeNOyrMTqrk9NYYxpDDS2LOtHY0wYsBa4RJ+buqMi5wZjTF9gtWVZ2caY24A+lmX9\nvVoqXAtV9Pxb8B2bDwQAt9f3uFVRFfwMtwbeB86zLCvNGNPAsqyUaqlwLVTBY/wi8JNlWc8bY9oB\nX1iWlVAd9a0JanQLVmFyVSAEqLnZYDWwLOtLy7LyCt6uAppVZ31qEsuyNlqWtam661GDnAVstixr\ni2VZR4A5wMXVXKcaw7KsZcCB6q5HTWNZ1h7Lsn4seJ0JbASaVm+txMfKPTdYlrXEsqzsgreKNZVX\n0fPvFOAxwHkyK1cHVOT43gzMtCwrDUDJVaVV5BhbQHjB6whg90msX41ToxMsAGPMI8aYHcC1wITq\nrk8NNgr4X3VXQmqspsCOIu93oh/KUgnGmATgDGB19dZEfKyy54YbUayprHKPsTGmC9Dcsqz5J7Ni\ndURFPsOnAacZY1YYY1YZYwactNrVDRU5xpOA64wxO4EvgH+enKrVTH7VXQFjzCKgUQmzxlqW9Yll\nWWOBscaYB4DbgYkntYLVrLzjU1BmLJAHvH0y61bdKnJsROTEGWNCgXnAv4/pWSD1iDHmOqAb0Lu6\n61KXGGNswNPAiGquSl3mB7QG+uBugV1mjOlgWVZ6tdaqbrkGeN2yrKeMMT2B2caYRMuyXNVdsepQ\n7QmWZVn9K1j0bdwZcb1KsMo7PsaYEcAQoJ9Vk2+oqwKV+OwI7AKaF3nfrGCaSJmMMf64k6u3Lcv6\nsLrrIz5XoXODMaY/MBbobVlWzkmqW11R3jEOAxKBb4wx4L5w+KkxZqjuw6qQinyGd+K+jzAX+MsY\n8zvuhOuHk1PFWq8ix/hGYACAZVkrjTEOIBaol90xa3QXwYKbEgtdDPxWXXWpiQqauMcAQ4v0jxcp\nyQ9Aa2PMKcaYAOBq4NNqrpPUcMb9a+8VYKNlWU9Xd32kSpR7bjDGnIF7wKChunfluJR5jC3LyrAs\nK9ayrISCQQFW4T7WSq4qpiLx7WPcrVcYY2JxdxnccjIrWctV5BhvB/qBZ6AxB5B6UmtZg9ToBAuY\naozZYIz5BbgA0BDB3p7FfeXrq4Kh7F+o7grVFMaYSwv6AfcE5htjFlZ3napTwWAotwMLcQ9U8L5l\nWb9Wb61qDmPMu8BK4HRjzE5jzI3VXacaohcwHDiv4Byzrr4PvVvXlHZuMMY8ZIwZWlDsCSAUmFvw\nGdDFmUqo4DGW41TB47sQ2G+MSQKWAPdalrW/empc+1TwGN8N3GyM+Rl4FxhR33pWFVWjh2kXERER\nERGpTWp6C5aIiIiIiEitoQRLRERERETER5RgiYiIiIiI+IgSLBERERERER9RgiUiIiIiIuIjSrBE\nRERERER8RAmWiIiIiIiIj/w/EgCE2Duxs+sAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x112b91b70>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"with pm.Model():\n",
" alpha = pm.Flat('alpha')\n",
" theta = pm.Deterministic('theta', pm.math.invlogit(alpha))\n",
" \n",
" def logp(value):\n",
" return tt.sum(tt.sum(pm.Bernoulli.dist(p=theta).logp(value)) + \\\n",
" tt.log(theta) + \n",
" tt.log(1-theta))\n",
" obs = pm.DensityDist('obs', logp, observed=y)\n",
" \n",
" trace = pm.sample(3e3, tune=1000, njobs=2)\n",
"pm.posteriorplot(trace[1000:]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Model 3.1: Jacobian Adjustment using `pm.Potential`\n",
"`pm.Potential` is the `increment_log_prob` equivalent in PyMC3"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using advi...\n",
" 0%| | 0/200000 [00:00<?, ?it/s]Mean ELBO converged.\n",
"Finished [100%]: Average ELBO = -8.0688\n",
"\n",
"100%|██████████| 3000/3000.0 [00:03<00:00, 883.05it/s] \n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAACsCAYAAABmdA06AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0FFXax/Hv7U4nnX0HAiQgiCj7CAoqIyKIKIig4ojb\ngMvoO+q4oriwqOOI4g6OoiCjKILIOKIoKC4gCC64g4KIECAQQnbI3l3vHx1amqwdmmz8Puf0obvq\nVtWtIl1PP3Vv3TKWZSEiIiIiIiKHz9bQFRAREREREWkulGCJiIiIiIgEiBIsERERERGRAFGCJSIi\nIiIiEiBKsERERERERAJECZaIiIiIiEiAKMESqQNjzFhjzKpAlxURETkcxpj2xhjLGBPU0HUROVop\nwRIRERFpwowxW40xgwOwHiVnIgGgBEtERERERCRAlGCJVMMYM8EY85sxJt8Ys8EYM6qKcpYx5h/G\nmC3GmL3GmGnGGNshZR4zxmQbY343xpxz0PRxxpify7exxRhz3ZHeLxERaR6MMXOBFOAdY8w+4OLy\nWZcZY1LLY9K9B5W3HRTbMo0xbxhj4spnryz/N8cYs88Yc4oxpqMx5uPysnuNMa8ZY2Lqbw9Fmh4l\nWCLV+w34MxAN3A+8aoxJqqLsKKAPcCJwPnDVQfP6AhuBBOBRYLYxxpTP2wMMB6KAccCTxpgTA7wf\nIiLSDFmWdQWQCpxnWVYE8Eb5rP5AZ2AQMMkYc0L59JuAkcAAoDWQDTxbPu/08n9jLMuKsCxrDWCA\nh8vLngAkA1OO5D6JNHVKsESqYVnWQsuy0izLcluWtQD4FTi5iuKPWJaVZVlWKvAUMOagedssy3rR\nsiwX8DKQBLQs38YSy7J+szxWAB/gSepERETq6n7Lsgoty/oe+B7oWT79euBey7J2WJZVjCdZuqiq\n+64sy9psWdaHlmUVW5aVATyBJzkTkSroJkaRahhjrgRuA9qXT4rA0wrlqqT49oPeb8Nzte+A3Qfe\nWJZVUN54FVG+jXOAycBxeC56hAE/BmQHRETkaLX7oPcFlMccoB3wljHGfdB8F+UX/Q5ljGkJPI3n\nwl8knjiVHfDaijQjasESqYIxph3wInAjEG9ZVgzwE57uEpVJPuh9CpBWi22EAIuAx4CW5dt4r5pt\niIiIHMryo+x24BzLsmIOejkty9pZxXr+VT69u2VZUcDlKEaJVEsJlkjVwvEElQzwDEYBdKum/Hhj\nTKwxJhm4GVhQi20EAyHl2ygrb80acli1FhGRo0060KGWZZ8HHiq/iIgxJtEYc375vAzAfci6IoF9\nQK4xpg0wPjBVFmm+lGCJVMGyrA3A48AaPMGrO7C6mkXeBtYB3wFLgNm12EY+8A88NyVnA5cCiw+r\n4iIicrR5GLjPGJMDXFRD2afxxJkPjDH5wFo8AzFhWVYB8BCw2hiTY4zph2eApxOBXDyx7b9HZhdE\nmg9jWf60KotIZYwxFtDJsqzNDV0XEREREWk4asESEREREREJECVYIiIiIiIiAaIugiIiIiIiIgGi\nFiwREREREZEA8edBw2rqEhGRIykQz9ZRrBIRkSOpxlilFiwREREREZEAUYIlIiIiIiISIEqwRERE\nREREAkQJloiIiIiISIAowRIREREREQkQJVgiR0hRqatelhERkebJ35igGCLSOPjzoGENfSvip/YT\nlvhVfuvUYUeoJiJNgoZpFzmEP3FEMUSkXmiYdhERERERkfqiBEtERERERCRAlGCJiIiIiIgEiBIs\nkVrSzcMiIiIiUpOghq6ASFPhdNh1s7GIiNRJUakLp8Pe0NUQkXqgBEtERETET/4mTP5epANdqBNp\nqtRFUKSelGRsI33BRFKfuIjtT49h75KncBXm17jcunXrGDp0KG3atMHpdNKqVSvOPfdc1qxZU2n5\ntWvXMnToUGJiYggPD6d79+7Mnz/fp8zvv//ORRdd5C0zcOBAvv7664Dsp4jI0eBAwlTbV125ivaR\n+f4zbH/mUlKfuJD0+fdSkrHV7/XMnz8fYwxt27atMO/ll1/mwgsvpF27dhhjGDt2bNX1cbl46qmn\n6NatG06nk/j4eAYPHsyuXbv8rpNIc6UES6QelOVnkv763VhlxSSOvJu4s66naNt3ZLx5P5blrnbZ\nnJwcjj32WB5//HGWLVvG9OnTycnJYcCAAXz55Zc+ZZcsWcLpp59Oq1atmDdvHm+//TbXXnstRUVF\n3jKZmZn079+fn376iZkzZ3qTr4EDB/Lzzz8HfudFRKROLMsi480HKNyyjrjB15E48h4st4v01++h\nLG9vrdeTk5PDLbfcQqtWrSqd/+qrr/Lbb79x1llnERUVVe26rrjiCh588EHGjRvHsmXLmDNnDj17\n9vSJMyJHO3URFKkHeV/+F8vtosWFk7A5IwCwR8aTPm8ChZvWEtb51CqXHTRoEIMGDfKZNnToUBIS\nEpg7dy4nn3wyAPn5+YwbN46///3vPPXUU96ygwcP9ln2ueeeIz09nZUrV9KxY0cAzjzzTDp06MDk\nyZN54403ArLPIiJyeBYvXkzxzg20vORfONv1ACCkzfHsfP5q8r5cRNzg62q1njvvvJOePXuSlJTE\n8uXLK8xftmwZNpvnmvvSpUurXM/8+fN54403+OKLL+jdu7d3+ogRI/zZLZFmTy1YUi+mTJmCMYZf\nfvmFs88+m/DwcFJSUpgzZw4Ac+fO5fjjjyciIoKBAwfy22+/+Sz/wgsv0LNnT5xOJwkJCVx99dVk\nZWX5lJkxYwannHIKcXFxxMTE0K9fP5Ys8e2WsXXrVowxzJw5k0mTJpGUlERMTAznnXceO3bsOGL7\nX7j5C0I79vEmVwDO5G7YoxIp2LzW7/WFh4cTEhJCUNAf10gWLlxIRkYGt99+e7XLrl27lk6dOnmT\nqwPr+/Of/8y7775LWVmZ3/UREalPRyKmpKVn+JSpTUwpy01n2yPDyf/ufXI+e5UdM64g9am/sOfN\n+/1qYarK4sWLsUfEeZMrAFtIOKHHnkzBr7WLHatXr+bVV1/l2WefrbLMgeSqJv/+978ZMGCAT3Il\nIhUpwZJ6NXr0aIYNG8b//vc/evfuzVVXXcU999zDc889x9SpU5kzZw4bN27k0ksv9S4zYcIEbrjh\nBgYPHszixYuZNm0aS5cu5ZxzzsHl+mPo9K1bt3LNNdewcOFCFixYQJ8+fRg+fHilV+MefvhhNm/e\nzEsvvcTTTz/NmjVruPzyy33KuN1uysrKfF6W21XxZVnV7rO7tJiynHSCE9pVmOdISKF07/ZaHTu3\n201paSmpqanceOONAFx77bXe+atWrSIuLo4ff/yR7t27ExQURHJyMvfff7/PcbLb7QQHB1dYf0hI\nCIWFhRV+iIiINFaBjCmjRgyn3Z2LvfdM3fvKx2yKORnHkNsJGXIbG0oTGD58OC0vfqDCPVW5axZS\nmr2L+HNuJm7Q3yhO28jedx/zKeN2uyuPIdXElPXr1+OoIna48jJwlxRWe3xKS0v529/+xvjx4zn2\n2GPrcoh91vXFF1/QtWtX7rzzThISEnA4HPTt25ePP/74sNYt0tyoi6DUq/Hjx3PllVcC0KdPH955\n5x1mzpzJ77//7u33vWvXLm6++Wa2bduGZVlMmzaNyZMnM2nSJO96jjvuOPr3788777zDyJEjAXjs\nsT+CmdvtZtCgQWzatInnnnuOoUOH+tSjffv2zJs3z/s5IyOD8ePHk5aWRuvWrQG46qqrePnll2vc\np/hzbyGi++Aq57uL9gGWT+vVAXZnJGVZO2vcBsDFF1/MokWLAGjRogXvvfceXbp08c5PS0ujoKCA\nSy+9lIkTJ9K7d2+WL1/Ogw8+SE5ODk8++SQAnTt35sMPPyQzM5P4+HhPHd1u7/1ch7YMiog0VoGO\nKYltviTsuFMAiD3zau98y3LjbN+L0qw08r99j9AOvi04QdEtSBwx3vvZVZBLzqcvUZafSVCk5zx7\n1VVXkepnTMnKysLmrHjflM0ZCXjiiy04tMp1PfLIIxQXF3P33XfXuN2aZGZmUlJSwn/+8x86dOjA\niy++SEhICNOmTWPo0KF8/vnn9OnT57C3I9IcKMGSenXOOed438fGxtKiRQv+9Kc/+dxUe/zxxwOw\nfft2fv75Z9xuN5dddplP17W+ffsSGRnJypUrvQnWunXrmDx5Ml999RUZGRneq4CdO3euUI9zzz3X\n53P37t0BSE1N9SZYU6ZM8bYUHXDe9FUV1hUU09L73nL/0VIU6K52jz76KHfddRfbt2/n2WefZfjw\n4Sxfvtwb0NxuN0VFRTz00EPcdtttAJxxxhlkZmby7LPPMmXKFKKjo7n++ut55plnuPLKK3nmmWcI\nCwvjoYce4vfffwdq31VERKShBTqmFG3/yZtgFe/eTO6q1yje9SvuglzAE1OC4iqOwhfa0TexCE70\ntDq58jK8CdaUKVNYZvWqcZ8OjimHY/PmzTz00EO89dZbOJ3Ow16f2+0ZkKm0tJT33nvPGytPP/10\nOnTowLRp01iwYMFhb0ekOVCCJfUqNjbW53NwcHCl0wCKiorYs2cPQJVdGzIzMwFP4Bw0aBBdunRh\n+vTppKSkEBQUxMSJEysdGS8uLs7nc0hIiHebB6SkpFQYzja4ZVrFShhPQlKWm87O5/+44umYBm2u\nn40tLBow5S1ZvlxF+d4rkTXp0KEDHTp04KSTTmL48OF069aN++67z9sF8kBr1FlnneWz3JAhQ3j+\n+edZv349p556Kh06dOC1117jhhtu8B7XE088kVtvvZXHHnuMpKSkWtVHRKShBTqmhBd5Hp1RlpdB\n+vx7CY5PJm7wdQRFJYLNTs5nr1KaWbFb96HncRPkAMBylXinpaSkENyyQ807Zf64yBUbG8v2PRVj\nh7u8npX1jDjgH//4B2eeeSb9+vUjJycHgJKSEizLIicnh5CQEEJDq279OlRsbCzGGLp06eJNrgAi\nIiI45ZRT+Pbbb2u9LpHmTgmWNGoHkoYPPvigQtA8eP7SpUvJzc3ljTfe8EmKCgoK6rxtf7sI2iPi\naHXlk97p79zUnwsX7MDYHQRFt6B0b2qFZUv3bseZ0s3vugUHB9OjRw++++4777SuXbtWu8zBLVMX\nXnghI0eOZNOmTQQHB9OxY0f+7//+j+TkZFJSUvyuj4hIU1BTTDl/1vcAFG5Zh1W8n4TzJxAUleCd\nb5UV13nbdeki2LVrV9ZtWFyhTOne7dijEqvtHrhhwwa2bdtW6X7GxsZy8803+4w4W5PQ0FA6dKg6\nQVTvB5E/KMGSRu2ss87CZrORmppaoWXmYAcSKYfD4Z22adMmVq9eXelDFUtd1T97CvzvImjsDkKS\nOnmn9+nTB/NmOgChx/Zl/08f4S7ejy0kHICiHetx5e0h9Ni+NdblUAUFBXz99dc+3R9HjhzJxIkT\nWbZsmbfLI3iST6fTSbduvomc3W7nhBNOADz3by1YsIDx48cjItJc1RRTHOXn7AOJlLHbvfNKs3ZS\nvGMD9siECsvVRl26CI4YMYI5c+ZQlPojzhTPed1dXEDh5i8J7zKg2vXMnz+/wrOppk6dyrp161i4\ncGGlsbEmo0aNYvr06ezcuZM2bdoAnkeEfP755xXudRY5minBkkatY8eO3HXXXdx4441s3LiRAQMG\n4HQ62b59Ox9++CHXXHMNAwcOZPDgwQQFBXHllVdy++23s2vXLiZPnkxKSoq33/jBHHabzyhQRak/\nAHDJC2twLt1faV22Th1GSFJ6nfYjqu8F7N/wKXsWPUh0v9G4i/eT/ekcgpM6e/v7A6xYsYJBgwbx\n0ksveW/cvu6664iLi6NPnz4kJCSwbds2ZsyYwa5du5g7d6532W7dujF27FgmTZqE2+3mxBNPZPny\n5cyaNYuJEycSEeHpSlJaWsqdd97JgAEDiIqKYv369Tz88MN07dq1xiHeRUSasppiSpG7G852PXC2\n6wU2O3vffYKok0fh2pdNzqrXCIpKrHHk2Kq0b9/e5yJcbYwYMYKQ1sez993HiT1jHDZnBLlrFwIW\nUX0v9Cm77dERXJ0xltmzZwPQr1+/Cuv7z3/+Q0hICGeccYbP9A0bNrBhwwYACgsL2bZtG2+++SYA\nAwYMIDExEYA77riDuXPncs455zBp0iSCg4N57LHHKCgoYMKECX7tm0hzpgRLGr1//etfnHDCCTz7\n7LM8++yzGGNITk5m0KBBdOrkCVZdu3bltddeY9KkSYwYMYKOHTsydepUli5dyqefftqwOwAERSbQ\ncsy/yP5oFhn/+xfGFkRop37Ennk15qD+9pZl4XK5fJLCvn37MmvWLF544QX2799PmzZt6Nu3L7Nn\nz/ZpqQKYOXMmbdq0Yfr06aSnp9O+fXueeOIJbr75Zm8ZYwy//vor8+bNIycnh7Zt23qHNq5s+HYR\nkeakupgS5PbcWxSc2I6E4XeQs+o19ix6EEdMErED/krh799QlPpjvdXVZrOReNFksj+ZTdaHz2GV\nlRLSpjMtx/zLc1/YwSw3JaX+D65UVOrijTfe4P777/dO+/TTT72x85NPPvEmZC1btuTDjz7hngl3\nMm7cONxuN6eccgorVqyosZu6yNHE+HElpm6XbEQaqUOfY1KTrVOH+bWMv+UPLCNyFDMBWIdilVBU\n6sLpsNdc8BBH8hxf15jQGLchcpSrMVapBUukEfH3R0Fdf0SIiDRnToddF7hEpMEowRJpRPz9UaAf\nBCIiIiKNi8bUFBERERERCRAlWCIiIiIiIgGiBEtERERERCRAlGBJs1BU6mroKoiIiIiIaJALaR40\nOISISPOlEVNFpClRgiUiIiKNmi6iiUhToi6CIiIiIiIiAaIES0REREREJECUYImIiIiIiASIEiwR\nEREREZEAUYIlIiIiIkeMv49S0aNXpKnTKIIiIiIicsRoFEg52qgFS0REREREJECUYIk0Yep2ISIi\nItK4qIugSBOmbhciIiIijYtasERERERERAJECZaIiIiIiEiAKMESERGReqN7QUWkudM9WCIiIlJv\n/L13FHT/qIg0LWrBEjmK1OXKsa42i4iIiNSeWrBEjiK6ciwiIoejqNSF02Fv6GqINGpKsERERESk\nVnShTqRm6iIoIiIiIiISIEqwREREREREAkQJloiIiIiISIAowZJGR6PWiYiIiEhTpUEupNHRDbQi\nIiIi0lSpBUtERERERCRAlGCJiIiIiIgEiBIsERERERGRAFGCJSIiIiIiEiBKsERERERERAJECZaI\niIiIiEiAKMESERERkUbD3+dh6vmZ0tjoOVgiIiIi0mj4+zxMPQtTGhu1YImIiIiIiASIEiwRERER\nEZEAUYIlIiIiIiISIEqwREREREREAkQJloiIiIiISIAowRIREREREQkQJVgiIiIiIiIBogRLRKql\nBz6KSHX0nRcR8aUHDcsRV1TqwumwN3Q1pI70wEcRqY7OESIivpRgyRGn4CsiIiIiRwt1ERQRERER\nEQkQJVgiIiIiIiIBogRLREREREQkQJRgiYiIiEiTVZeRLDX6pRxJGuRCRERERJosfwfTAg2oJUeW\nWrBEREREREQCRAmWiIiIiIhIgCjBEhERERERCRAlWCIiIuKlm/9FRA6PBrkQERERL38HDNBgASIi\nvtSCJSIiIiIiEiBKsERERERERAJECZaIiIiIiEiAKMESEREREREJECVYIiIiIiIiAaIES0QCqi5D\nPGtYaBERacwU28QfGqZdmpSibT+Q89mrlKRvxgQFE9rxJGIHXuVTpiw3nZ3PX13p8sk3z8fmjADA\nXVpE1oczKfx1DTZnBDGnX0n4Caf7lM/94k32r/+UpLFP16p+u+dNALebVpc/WmFe/vfLyFo6nTbX\nzyYouiUAe5c8yf6fPvKWSZydQL6zJdGnXExoh97e6dseGf7Himx2bCHhOOLb8mD4N7j2d8AeHlOr\n+tUHf4d4Bg3zLCL179B4csXOkbhizsYeHustU108yZmQ7X1fm3jy6KOPkvbScySNfRpjs9dYv93z\nJtB/1cPQ/+4K86qKJ+agWGELjcIRn1xjPEmcHUu+swXO9n8istc5jSqeHElFpS6cjpr/Hw5QbBN/\nKMGSJqNo+0+kvzGR0GNOJHHk3bgK88n5bC7p8++l+LHRFcpH9RtN2LF9faaZ4FDv+7y1b1K09Tvi\nz72Vkozf2fvuEwS37Igjrg0AZXl7yf18AS1GP1CrYFhXtrBoWlwwEYDnLuzIJTfdy56FU2jxlwcJ\nbd/LWy6822Aiew3FsizcRXkUp21k+vTpZO4rJvGCiTjbnnDE6igiTZO/PyKPFpXFk5Ur3yS9YCVJ\nf30aE+TwKV9ZPImMjPS+r008+ee//0nciElHNJ4kJiZiP/suAFz7s8n76q0a48mLf+nM5Q+/Sv43\n75K/7p2jJp7oeW9yJCnBEr81VMDOXf06QVEtSLzgPm+AcsQns/uVW5k9ezbQzqd8UEwrQtocX+X6\nCresI7L3MMI69SWsU1/2r/+Uom3fewNi9kcvENa5/xEPNMYW5K3nyJHDaLGilB3PjSN/3WKfgGiP\njPfZn7Bj+/LFW8/Q9oQTyfjfQ7S5bhY2h/OI1lVEmhZdda9cZfHkzSljOOmkk9j3wwdEnuh7DCqL\nJ3b7H3GwNvHkktGj+SjxyMaT4OBggg6qp7NdzxrjyXnnDSN2tY2o3uex+7W7FE9EAkAJlvitoa76\nFKdtJLzrQJ+rfyFJnbCFRvHWW29B71v8Wp/lKsUEhXg/2xwhWGUlgCdYFm3/idbXPh+QuvvDFhKG\nI641Zdm7aizbsmVLYgdeRcZ//8n+DSuJ7DmkHmooItK0VRZP+vTpgy00ioJNayokWDWpTTx59NHF\n9J62NjA7UEv+xBN7eKziiUiAaJALaTqMDWOveE3A2IP46aefKkzPWfEy2x4dQeqTF7Nn0QOUZGz1\nmR/SujP7f/qIsn1ZFG5ZR8me3wlp3RmrrJSs5c8TM+Cv2EOj6lRVy+2irKwMy+3yvrCsWi/rytuL\nLSS8VuWd7f8ENjvFOzfUqa4iIkedauJJ6d5tFaYHIp7Ex8fXqaoHxxHFE5GmQS1Y0mQ44ttQnLbR\nZ1pZ7h5c+7LJKnGQdGCi3UFEr6GEtj8RW1gUpZk7yF27kN2vjifpiidwJCQDEH3aGPYsnMLOZ68E\nIOrkCwhpcwI5q1/HHhpNRI+6Xb0r3rmB1Gnn45hW+2Ust2ekoR07dpD1wXO49mcT1feiWi1rc4Rg\nD43CtS+75sIiIlJpPNm2bZvnPHpQ17/q4snPd5zpLXak4snq1ath9fl+LXMgnrj2ZZP7+XzFE5EG\noARLmozI3iPIfPdxslfOJar3ebiL8slcOgOMwWb7ozE2KCKO+LNv9H52JncjtENv0mb/ndw1C0g4\n7w5PucgEksZNpyxnNzZnOPbQKEpzdpP35X9pddkjWGXFZH08i4JNazGOEKa3uRfoUGM9HS2OIX7o\nP3jnpv6cN32Vd3rhr2vJXbOgQnnXvkxSp3kCaPI0z0Ac0f0vI7LPebU+NhYWmFoXFxE5qlUWT664\n4hEwBmNqF08eeughaDvGU64W8eT6669n+9wFGEcIUSeNJKp3zef4nj17kt5zbIXpVcWTnTt3wrQ/\nEjLFE5GGoQRLmoyIrgMpy9xB3ldvkbdmAWAIO+HPhHbsQytXBu5qlg2KSiSkbReKd//qM90YgyPW\n2/ZF9ofPE9FjCMEtOpC98hVKdm+m9dXP4srP5J577iFs2ASfG4UrY3OEEpLUiT59+hCSlO6dXrJn\nS+Xlw2JocdFkMIY1U85nwL9/8GuUKXdpMe6CPOzhcbVeRkTkaFZZPDn9kr8Q2rEPpRkVuwge7EA8\n+eqrr7wJFtQcT74u+c0bT3bPuwtHfHKN8SQiIoKcpE4VplcVT1q0aIFt6N1gDPbQSOyRCYonIg1A\nCZY0KTGnX0FUv9GU5e7GHhaNPTyWnS9eT/9hA1lZqzVUfVmuYNMaSvZsIWHEnQAUbVlHePfBnu2E\nRTNkyBCW//5NjQHRX8ZmJ6Q8gLZv3x5jW+/X8kW/fwOWG2fbLgGtl4hIc3ZoPHl9+uU44pMJqeW5\n1Bj/4slf772Zx3d64klo+z9RdATiicPhIKiShKy2FE9EAkODXBzlmuJTxm3BToIT22MPj6VwyzrK\nsnZw/fXXV7tMWd4einds8CYyh3KXFpH10QvEnnkttpAw73SrtMj7ft++fbW+sbi+7Nmzh+xP52CP\niCPskIdaiohI9Q6OJ0uXLqUsawcRvc6tdpkD8eTkk0+udH5V8WT//v1/lCkpbHTxxLU/R/FEJEDU\ngnWUa0rPSClJ/43CLesIbtkRgKIdG8j7chFRfS/k1FNPhcWe/cj6eBZYFiGtj8ceFk1pluemZIyN\n6FP+Uum6cz+fjyOuLeEn/Nk7zdm+F/nr3sUR1xbXviw++vgj4i+cfOR3tAqu/EyKd/5S/qDhfIrT\nNtJ9zlW4i4tpceEkbI6QmlciIiKVxpPzn3yLqL4X+jz7sLp4cu+997JyzuYK664qnsyYMYOSk8fi\n2pdF0bbviTp51JHf0SocHE/efdcie+Wr7Pt+GWApnogEgBIsaTpsQRT+9jW5XywCVymO+GTih9xA\nRI+zfIoFJ6SQ/+377P9xOe7SImyhkThTehJz2hgc8W0rrLY0czv53ywh6a9P+UyPPvUSXPtzyXz/\naUxQMFOnTmXG3oZ7uv3+n5az/6flYLNjCwnHEdeW8TfdxMw9HbCHRTdYvUREmpxK4snzzz/P/Rtb\n+BSrLp507twZ8E2wqosng7IjeWW+J57EDBhL6DEnHum9rNLB8WTs0hiKnS2J7D2cyF7nKJ6IBIAS\nLGkyghPb0eryR2ssF9FjiF9D4jrik0m5dWGF6bbgUBKG/fHw4jvuGMaMGlr7Wl06tcp5kT3PJrLn\n2T7TEobdWqs6trvr3Uqn33ffMGb52QLZGBWVunA6an8jtr/lRUQOVlk8GTduGPcfcj4NZDyZM2cO\nn7Ss3XDp4Iknq6YOq7SXSVXxZGsV5Q92aDypzTIi4h8lWCLS4PztqtpQ3VRFREREaqJBLkRERERE\nRAJECZZtiFnNAAAV2UlEQVSIiEgT0RRHfhU5Wvn7fdX3u/lQF8FGZNWqVcyaNYuvv/6aX375hbZt\n27J169ZaLbt161aOOeaYSudlZ2cTExPj/fz7778zfvx4li9fTmlpKa6EjsSecVWVQ5iL1FVZXgbZ\nH71I4dbvAAtn+17EDbqWoKgWtVh2DzmfvUrRth9xF+Zij0wg7Pg/E91vtE+5cePGsXbtWnbu3Inb\n7aZjx45cc801/P3vf8du/+M+rTPOOIMVK1ZU2M6TTz7JLbfcUmG6yOHavn07t956Kx9++CGWZTF4\n8GCeeuopUlJSalz2nnvu4euvv2bdunVkZWUxZ84cxo4dW6E77e55Eyje/lOF5WPPvJaok8736U7r\nLi0m74s32b/hU8ryMrCFRBCS1InEUfdg7I7A7LRIE3Y4MQugdO92cla9SlHqj1ilRfR8qz0ZyQOI\n6nO+t4yrMI/c1a9TuPlLXPuzsYXHEtaxD9GnXcr2Zy4FIC8vj6eeeoqlS5eyceNGXC4XXbp04c47\n72TkyJFHYtclwJRgNSIfffQRn332GX369MEYQ35+vt/ruPvuuxkxYoTPtMjISO/7zMxM+vfvT2Rk\nJDNnziQsLIzRN95L+vx7SLriCRwJyYe9HyLgeRZM+nzPD7eEYbeCMeSsnEv66/eQNG4GtmBn1cuW\nFJE+/z4st4uYP19OUFQixbs3kbtqHmXZafDEhd6yhYWF3HTTTXTs2BFjDMuWLePmm29m8+bNPP30\n0z7r7dGjBzNnzvSZ1r59+4DutwhAQUEBZ555JiEhIbz88ssYY7jvvvsYOHAgP/zwA+Hh4dUuP336\ndHr16sXw4cN55ZVXqi3rSGxP/Nk3+kwLim7p89lylbFn4WTKctOJ7jcaR0IKroJcirZ+i+V2YzRm\njBzlaopZNSne9Svp8+/BmdKd+KE3YQsJ5/aB8UyY/6W3jGVZZCx6kNKsncT8+XIc8W3Lk7LXKN61\nGevpMRhjSE1N5d///jfjxo1j4sSJ2Gw2Xn/9dUaNGsWMGTO44YYbjuShkABQgtWITJw4kcmTPc9Z\nuvzyy1m1apXf6+jQoQP9+vWrcv5zzz1Heno6K1eupGNHz/M/WnxSws6Z15Cz6jUSR06oW+VFDrHv\n+2WU5aTT+trnccS2Bjw/BNNe+Bv7vnu/2mfAFO/cQFl2Gi0ufsA7lLGzXQ/chfvI+/K/FBQUEBbm\neYDn/PnzfZYdMmQIaWlpvPTSSxUSrMjIyGq/HyKB8uKLL7JlyxY2btzIscceC3gS/E6dOjFz5kxu\nu+22apfPzc3FZrOxefPmGhMsW3AYIW2Or7ZM3ldvUZL+G62v/jdBUYne6eGdT6vlHok0bzXFLLiw\nymUty03mkidwtutJiwvu807/29+G8a8tbbyfy7LTKN75M3Fn30hkr6EAOFN6gDFkffBvNm3aROfO\nnTnmmGPYsmWLN84BnH322Wzfvp1HHnlECVYToHuwGhGb7cj/d6xdu5ZOnTp5kyvwPMne2bYLhb99\nheVW/18JjMLNXxDSurM3UAE4YloR0rYLBZu/qHZZy1UGeH44HszmDAfLwrKsapePj48nKEjXj6Th\nLF68mH79+nmTK4BjjjmG0047jbfffrvG5QMdD/K/eY+wzv19kisR+cPhxKyi1B8pzdxO1EnVd9+z\nXKWAZ9j+g9lCPC3abrcbgPDwcJ/k6oA+ffqQlpZW885Ig1OC1czcfffdBAUFER0dzYgRI/jxxx99\n5tvtdoKDgysuGOTAKiumLHtXPdVUmruSvak4EtpVmO5ISKF0b2q1y4a270VQbGuyV/yHkr2puEsK\nKdz2PflfLybiT+dgP6R7oWVZlJWVkZOTw6JFi3j55Zd9WggO3Dj87bffEh0djcPhoEePHsyePTsA\neypS0fr16+nWrVuF6V27dmXDhg0B3VbJnt9IffJitk07n7SXbiT/+w985pfl7cGVn4EjphWZ7z9D\n6pOj2fbYKNLn30NJ+paA1kWkqTqcmFW8w/OdtspK2fXK7Wybdj7bp1/GP/7xD9ylxQetqx0hyd3I\n/XwBxbt+xV1SSHHaRnI/n4+zQ29OOOGEarezcuVKjj+++tZqaRx0ibeZCAkJ4brrrmPIkCH8fdGv\nlGbu4L2VC3m3T1+fe6uydweRv2EjyTe/jj00iq1Th2FZbkp2bQLAXeT/fV8ilXEX7sPmjKgw3eaM\nxF20r9plTVAwrS57lIz//Ytds//unR7RYwhxZ11f4Ub/gs1fkrHogQNLE9XvImbvP5HZ5WW2Th3G\n6aefzmWXXcZxxx1HTk4Or7zyCtdccw27du3ivvvuQySQsrKyiI2NrTA9Li6O7OzsgG3HmdyN8C5n\n4Ihrg7t4P/t/+pispc/g2p9FzKmXAODKzwIg94s3PYNajLgTy1VKzqp57H79blpfNb3WN/GLNFeH\nE7Nc+zIB2Lv4ESJPHE7sGX+leNdmZs2aBW2/8nYbNMbQ4qIp7F3yOLtfudW7fGjHk0g4fwJFpS6c\njspviHzhhRdYu3Ytr776qndadeWlYSnBagCWZeFy+XbFO9zuTElJSTz//PMA3PblEpzJ3Qjt0Ju0\n2X8nd80CEs67A4CIP51L3rp3yHz3CWIHX8euXbvIXj6Tspx0z4qMGjWl4VllJWQsfgRXQS7xw2/3\nDHKRtoncz18Hmx0eOc+nvDO5K62ufBJ38X6Ktn1P3pdvgTHEnn6lt8wDDzzgs8z555/PqFGjeOih\nh7jllluIiKgYWEUau5g/X+7zOaxTP/b895/krXnDO3LZgS61xhFC4oWTsDk8LcDBrTqx84Vryf9m\nCbFnjKvfios0J+XfsfCuA73fSWdKD/5v6HFMmDCB0r3bvRe6M5dOpyRtI3Fn34AjPpnSTM8gF3v/\n9zDB00b6XDw8oCj1B9LfmEx4tzO598cY7j3o4qE0Tvo13QBWrFiBw+HweR0JQVGJhLTtQvHuX73T\nHDGtSDjvDorTN5P2wrW0bt2a4p2/ePsN2yMqXnEVqQubM6LSq37uovxKrxIebN8PH1Cc+iMtLppC\nRNeBOJO7Ed33AmIHXs2+797n+++/991WSDghSZ0Ibd+L2AF/JfqU0eStfZOy/L3VbmfMmDEUFRVV\n6EorcrhiY2MrbamqqmUrkMJPGIBVVkJpxlYA7KGekWSdbbp4kyvwxAhHXFt1ExTh8GKW7cB3rH0v\nn+lDhgwBPN14AQp++4qCn1cQP+w2InudgzO5G5G9ziFh2O0Ubvmad955p8K6i3dtYs+iB3G260H8\n0H/Uad+k/qkFqwH07t2br776qh63aHw+hXc+jbBO/SjLSmPl3Wcx6MVfyFz2LPbIRHUTkYCpqt96\n6d5UHAnVPweoJGMrNmcEjtgkn+khrY8D4OeffwYiK1nSI7hVJ7DclOWmExSZUGNdjTE1lhHxR9eu\nXVm/fn2F6Rs2bKBLly5AfXTv8fxdB8W0wgSFVFNMf/8ihxOzKrt3y5fnO3bgokdI0nE+c31jW3fv\n9JKMrex5YxLBLTqQOPIejF0/25sK/U81gMjISPr06VNjuerHSatZWd4eindsIKxTxWGpjc2OIyGZ\njh07Upa/moJfPiPq5AsOc4sifwg7ti/Zn8ymNGc3jphWAJTlplO882diB/y12mXt4bG4i/ZRmp3m\nM6JTcZrnXsE2bdrAd3lVLu958KohKLpVtdt57bXXCA0NpXv37tWWE/HXiBEjuOOOO9iyZQsdOnQA\nPA+EX716NVOnTgWocC9hZUqzPSOG3bHwe6b8sqRWXYL2b/gUExSCI7E9AMYeRGjHPhRtX4+7pMj7\nDLqyvD2UZu0g7Ni+dd1NkWbjcGJWaIfeYHdQ9Ps3Pt+npUuXAhCc1AnwxDbwtEqFHtTaVZy2ESiP\nbTmeaaVZO0lfcB9BMa1ocdFkbI5qLpJIo6MEqxHJyMhgxYoVAKSmplJYUEDiyLsBz5WV4PIrKEWp\nP5I+/17iz72ZiG6DAMj6eBZYFrPGX0rRtt8ozdpB7tqFYGxEn/IX7zYsVxnZn87BmdwNExLG9Olb\n2P3KFBwJKdU+l0jEXxE9zyb/m3fJWPQgMadfARhyPptLUGQCEb3O8ZYry91DUFAQEadcQsxpYzzL\ndh9M3lf/Y8/CKUSf8pfyBw3/Su7n8wludSynnXYaLHmfgt++Yv8PHxJ6bF+CohI9ow1u+Zp93y8j\notdQgiLjAfjss8+YOnUqF1xwAe3btyc3N5eXX36ZxYsXM3Xq1Bof+irir2uvvZYZM2Zw/vnn889/\n/hNjDBMnTiQ5OZnrrrvOW64sdw87Z15D9GljvH//4DnPuwpyce33dDMs2f0r+x1O3nyzEPAM8Vy0\n/Sfy1r5J6HGnEhTdAqu4gH0/fUTh5i+IGTDW52He0f0vo/CV29jz5hSiTh6FVVZK7up52EIiiOw9\nvH4Oikgj5k/MOvQ7aw+NIrrfaHI/n48JDsPZriclu3/lgS8WEN5tkPdCYdhxp5Kz8hUylzxB9KmX\n4IhrS2nWDnJWz8MemcioUaOYuH4Frv05pC+YCK4yovtfVqFlLbhlR0yQo06t4BoYo34owWpE1q9f\nz+jRo30nvu250hl92hiC+19WPtECy+29qRIgOCGF/G/f57rrriMnLx9baCTOlJ7EnDYGR3zbP9Zn\nDGXZaWRuWIG7eB9PfZlCePeziD5lNMZ+ZO4Fk6OTLdhJyzEPkf3RLPa++zgAznY9iRt07SHPACkf\n9MVye6cERbek1RWPk7tqHjmfzcVdmIc9MoHInkOJOvUv3mcEOWJaYWGR89lcXAU52EIicMS2Jn7Y\nrYR3GeBdX1JSEm63m0mTJrF3717vMO3z5s1jzJgxiARaeHg4H3/8MbfeeitXXHEFlmUxaNAgnnrq\nqUMGVDlwPnf7LJ+z6rXylliP/G+WkP/NEka/De3uehcAe0QclmWRu+pVXIV5GFsQwYntSThvvM/f\nP3hiRMtLHiJ7xX/Y+/ajYLPjbNeDxAvu815VFzma+ROzKvvORp82BltwKPnfvkfel29hj4jlnvHj\nebmw9x/bCAmj1RWPk7NqHrlfLMK1Lwt7RBxhHU8muv+l3nNDaWYqrrw9AGS8eX+Fura5fjZB0S1r\n1Qp+KA2MUT+UYDUiZ5xxRoUHqFb2xXGm9PAG2AMiegwhoscQtk4dVu2XzdjstLhosvfzbzWUFzkc\nQVEtSBx1T/VloltiWVaFv8PghBQSR06odllHfDItRt1bYz3atjuG999/v+YKH0RX+eRwpaSksGjR\nomrLBEW3rHA+B2h16dRKyx98jnfEtqblxRV/fFUlpHVnWo15uNblRY42tY1ZlX1njTFEnTzKpzfQ\nAw8M45VDYltQVCIJ595c7TYq+50nTYsSLBFp9nSVT+qDknIREQElWPVKwVdEpPnyN5FXEi8i9c3f\n36L67Vo3SrDqkYKviIiIiDQU/RatH+bQe36qcbijhjcrdc3o/f2jrq78tkc08pMcvfzpn17Td6ky\nvzw4VFf56l8gHsjUoLEqkOf4g+l8L9K01SZm+Rur6hLb/F3G31gIR0U8rDFWqQWrjnRPh0jzpqt8\nIiJytNPv3bqxNXQFGouiUldDV0FEmrC6nEN03hEREWl+mm0Llr/Nk7paLSKHoy5X+X55cKhf5Y+C\nbheNho61iEjd1MdAGo19sI5a34N1//33LwUSjmx1mrTWQFpDV6KJ0THzn45Z3ei4+a8hjtneyZMn\n+5d1HkKxym/6btQPHef6oeNcP47241xzrLIsS68AvKZMmWI1dB2a2kvHTMdMx63xvnTMjo6X/p91\nnJvTS8dZx7mxvHQPloiIiIiISIAowQqc+xu6Ak2Qjpn/dMzqRsfNfzpmRwf9P9cPHef6oeNcP3Sc\na+DPc7BERERERESkGmrBEhERERERCRAlWCIiIiIiIgGiBEtERERERCRAlGAFiDFmmjHmF2PMD8aY\nt4wxMQ1dp6bAGDPaGLPeGOM2xvRp6Po0ZsaYocaYjcaYzcaYCQ1dn6bAGPOSMWaPMeanhq5LU2GM\nSTbGfGKM2VD+3by5oeskh6emc4cx5rby/+8fjDEfGWPaNUQ9m7ranqONMRcaYyzFvLqrzbE2xlx8\n0HlsXn3XsTmoxbkjpTxefFt+/ji3IerZGGmQiwAxxgwBPrYsq8wY8wiAZVl3NXC1Gj1jzAmAG5gJ\n3GFZ1tcNXKVGyRhjBzYBZwE7gK+AMZZlbWjQijVyxpjTgX3AK5ZldWvo+jQFxpgkIMmyrG+MMZHA\nOmCk/taaptqcO4wxA4EvLMsqMMb8H3CGZVl/aZAKN1G1PUeXf6eWAMHAjYp5/qvl33Qn4A3gTMuy\nso0xLSzL2tMgFW6ianmcXwC+tSzrOWNMF+A9y7LaN0R9Gxu1YAWIZVkfWJZVVv5xLdC2IevTVFiW\n9bNlWRsbuh5NwMnAZsuytliWVQLMB85v4Do1epZlrQSyGroeTYllWbssy/qm/H0+8DPQpmFrJYeh\nxnOHZVmfWJZVUP5R8atuanuOfhB4BCiqz8o1M7U51tcCz1qWlQ2g5KpOanOcLSCq/H00kFaP9WvU\nlGAdGVcB7zd0JaRZaQNsP+jzDvSjV44wY0x74E/AFw1bEzkM/p47rkbxqy5qPM7GmBOBZMuyltRn\nxZqh2vxNHwccZ4xZbYxZa4wZWm+1az5qc5ynAJcbY3YA7wE31U/VGr+ghq5AU2KMWQ60qmTWvZZl\nvV1e5l6gDHitPuvWmNXmuIlI42KMiQAWAbdYlpXX0PWRI88YcznQBxjQ0HVpbowxNuAJYGwDV+Vo\nEQR0As7A0yK70hjT3bKsnAatVfMzBviPZVmPG2NOAeYaY7pZluVu6Io1NCVYfrAsa3B1840xY4Hh\nwCBLN7d51XTcpFZ2AskHfW5bPk0k4IwxDjzJ1WuWZf23oesjh6VW5w5jzGDgXmCAZVnF9VS35qSm\n4xwJdAM+NcaA56LjYmPMCN2H5bfa/E3vwHNfYSnwuzFmE56E66v6qWKzUJvjfDUwFMCyrDXGGCeQ\nABz1XTLVRTBAypuf7wRGHNSXXSRQvgI6GWOOMcYEA5cAixu4TtIMGc+vv9nAz5ZlPdHQ9ZHDVuO5\nwxjzJzwDDY3QvSp1Vu1xtiwr17KsBMuy2pcPArAWz/FWcuW/2sTD/+FpvcIYk4Cny+CW+qxkM1Cb\n45wKDALvoGVOIKNea9lIKcEKnBl4rlB9aIz5zhjzfENXqCkwxowq77t7CrDEGLOsoevUGJUPoHIj\nsAzPoANvWJa1vmFr1fgZY14H1gCdjTE7jDFXN3SdmoDTgCuAM8vPZd9p6N2mq6pzhzHmAWPMiPJi\n04AIYGH5/7cu3viplsdZAqCWx3oZkGmM2QB8Aoy3LCuzYWrcNNXyON8OXGuM+R54HRirHlweGqZd\nREREREQkQNSCJSIiIiIiEiBKsERERERERAJECZaIiIiIiEiAKMESEREREREJECVYIiIiIiIiAaIE\nS0REREREJECUYImIiIiIiATI/wNl+F2a1CkVoQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x113f99860>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"with pm.Model():\n",
" alpha = pm.Flat('alpha')\n",
" theta = pm.Deterministic('theta', pm.math.invlogit(alpha))\n",
" obs = pm.Bernoulli('obs', p=theta, observed=y)\n",
" constrain = pm.Potential('constrain', tt.log(theta)+tt.log(1-theta))\n",
" trace = pm.sample(3e3, tune=1000, njobs=2)\n",
"pm.posteriorplot(trace[1000:]);"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment