Skip to content

Instantly share code, notes, and snippets.

@fonnesbeck
Last active October 18, 2018 02:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fonnesbeck/129b16ca33591501b9c7c469a0ef647c to your computer and use it in GitHub Desktop.
Save fonnesbeck/129b16ca33591501b9c7c469a0ef647c 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": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import pymc3 as pm\n",
"import arviz as az\n",
"import warnings\n",
"warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n",
"warnings.filterwarnings(\"ignore\", category=FutureWarning)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some fake data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"MALE, FEMALE = 0, 1\n",
"y = np.array([34, 42])\n",
"gender = np.array([MALE, FEMALE])\n",
"N = np.array([104, 110])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the simplest model -- two independent probabilities."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as affection_status:\n",
" \n",
" π = pm.Beta('π', 1, 1, shape=2)\n",
" δ = pm.Deterministic('δ', π[1] - π[0])\n",
" like = pm.Binomial('like', N, π, observed=y)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"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: [π]\n",
"Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 1673.65draws/s]\n"
]
}
],
"source": [
"with affection_status:\n",
" \n",
" trace = pm.sample(1000, njobs=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the posterior distribution of the probability difference"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4FFXbwOHfpPdCeiAkQCCht9B7EUGxoQJWEBTl1c+Gr6iogIpiL4gFkSIqooCvAiJK772XAIEkkFCSkN6T3fn+mJBCQgokmd3kua9rruyeObvzzCTZZ8+ZM2cUVVURQgghTI2F3gEIIYQQZZEEJYQQwiRJghJCCGGSJEEJIYQwSZKghBBCmCRJUEIIIUySJCghqoGiKPcqirJLUZRkRVGyFEWJUBTlDb3jEsKcWekdgBDmTlGUxsCvwB/Am0A64A2E6BmXEOZOEpQQN88NrTfiDLBJVdVsneMRok6QLj4hbpKqqoeBGcBLwBVFUdYoivKUoijWOocmhFmTBCXETVIUpQEwGJgK9AZWA+8AOxVFcdUzNiHMmSJz8QlxcxRFmQ+0VFW1e7GyNsBBYKaqqq/rFpwQZkxaUELcBEVRbIHRaK2mQqqqHgW2AUP0iEuIukASlBA3xxOwAxLKWHcR8KrdcISoOyRBCXFzkgAj4F/GOv+C9UKIGyAJSoiboKpqJrAHGKkoiuXVckVRmgE9gQ16xSaEuZNBEkLcJEVRBgFrgHXAl4AL8BbgDHRQVfWCjuEJYbYkQQlRDRRFGYyWlDoCucB64GVVVU/rGpgQZkwSlBBCCJMk56CEEEKYJElQQgghTJIkKCGEECZJEpQQQgiTVJXbbchoCiGEENVFqaiCtKCEEEKYJElQQgghTJIkKCGEECZJEpQQQgiTJAlKCCGESZIEJYQQwiRJghJCCGGSJEEJIYQwSZKghBBCmKSqzCQhhChGVVX2RSfx56ELHDqfzJn4DLLzDFhZKgR5ONLK34UhrXzo18IbexvLit9QCFFCVe4HJVMdCVFg06l43l11gpOX07C3tqRDgBshvs442FiSk28kMiGD/eeSSM7Mw9Xemke6B/JYryA8nGz1Dl0IU1HhVEeSoISoguTMXF5ZdoS/j10iyMOB/wwI5va2fjjalu6MyDcY2RWZyA87ovjn+GUcbax4ekAw43oHYWslLSpR70mCEqK6HI5JZuKP+4lLy+aFW1owvneTSieaiLg0Zq4OZ+2JOEJ9nZn1QEea+zjXcMRCmDRJUEJUh62nE5iwaC/uDjbMfqgTHQLcbuh91p24zMtLD5Oek8+797Tl3s6NqjlSIcyGJCghbtba45f5z0/7aerlyA/juuLtYndT7xeflsOziw+w4+wVHu/dhFdva4mlRYX/q0LUNZKghLgZe6MSeWjuLkJ9nflhXDdcHayr5X3zDEZmrDrBgu1R3NHen09GtsfaUq76EPWKJCghbtSZ+HRGfLWdBo42LH2qR42MwPtm0xlmrg5nYKg3Xz3UCTtrGTwh6g25YaEQNyI9J58nF+3D0kJh4WNda2x4+FP9mvHO3W3YcDKOsfN3k56TXyPbEcIcSYIS4hqqqvLy0kOcjU/nywc60tjDoUa393D3QD4b1YE9UUk88v0uMiRJCQFIghKilB93RvPXkUu8PDSUnsGetbLNuzo0ZPaDnTh0PpmnftxHTr6hVrYrhCmTBCVEMacvp/HOqhP0a+HFk32b1uq2h7bxZea97dhyOoEXfz2EwSinfUX9JnPxCVEgN9/Is78cxNHWig/vb4ei1P7Q75FhASRn5vLuX+G42Vvzzt1tdIlDCFMgCUqIAnM2n+HExVTmPNIZb+ebu9bpZkzo24zEjDy+2XQGXxc7/m9Qc91iEUJPkqCEAM7Gp/PF+ghub+vHkNa+eofD5KEhXE7N5pO1p2jp58LgVj56hyRErZNzUKLeMxpVXl1+BFsrC6be2UrvcABQFIX3RrSljb8rLyw5yJn4dL1DEqLWSYIS9d5v+86zKzKR125rqWvX3rXsrC355pHO2FhZMOGHvaRl5+kdkhC1ShKUqNfi0rKZseoEXZs0YFRYgN7hlNLQzZ4vH+xE1JVMXlhyCKOM7BP1iCQoUa+9vfIE2XlG3hvRFgsTnbC1RzMPXr+9JWtPXGbu1rN6hyNErZEEJeqt3ZGJrDh0gaf6N6OZl5Pe4ZRrbM8gbm3tw4drTnI0NkXvcISoFZKgRL1kMKpM+/MY/q52TOzXTO9wKqQoCjNHtMPD0ZZnfzlAdp7MNCHqPklQol76de95jl9M5dXbWmJvYx4ziLs72vDxyPacjc/gi3Wn9Q5HiBonCUrUOylZeXy05iRdgtwZ3s5P73CqpFewJ/d1bsSczWcJv5SqdzhC1ChJUKLembXuNImZuUy9o7XJTiOUk5PDpEmT8Pb2xtHRkdtvv52oqCgAptzWEhd7a15dfqTUqL4lS5YwYsQI/Pz8UBSFBQsW1H7wQlQTSVCiXjkTn86C7VGM7hJAm4aueodzXc8++ywLFizgo48+YunSpSQkJHDLLbeQnZ2Nu6MNbwxvyYFzyfy0K7rE65YuXUpUVBTDhw/XKXIhqo/cUVfUKxN+2Mv2M1fY+N/+eNbQTQhvVkxMDEFBQcybN49HH30UgNjYWJo0acJXX33F448/jqqqPDpvNwfPJfPvi/3wddUuMDYajVhYWJCeno6zszPz589n7NixOu6NENcld9QVddfYsWMJCwtj1apVtGrVCgcHB26//XYSExOJiIhgwIABODo6EhYWxuHDh9kXncQ/xy/zRO8g5n75KcHBwdja2tKiRQsWLlxY4r1XrVrFLbfcgre3Ny4uLnTv3p1//vmnRJ1p06bh6enJgQMH6N69Ow4ODnTs2JEtW7bc1H5d3c6IESMKyxo2bEjv3r1ZvXo1oI3qe+fuNuQajLy18lhhPQsL+ZcWdYf8NQuzdu7cOd58803eeecd5syZw/bt25kwYQKjR49m9OjRLF26lPz8fEaPHs3M1SfwdLLl+LLPeOedd5gwYQKrVq3innvuYdy4caxcubLwfSMjI7njjjtYtGgRy5Yto2fPngwbNoxt27aV2H5mZiZjxozhySefZNmyZdja2nLPPfeQmZlZWMdoNJKfn1/uYjAUDRsPDw+nUaNGODmVvDarZcuWhIeHFz4P9HDkmQHB/HXkEjvOXKnuQyuE/lRVrewihEkZM2aMamlpqUZERBSW/fe//1UBdeHChYVlq1atUgHVf/zX6oe/blAVRVEXLFhQ4r0eeeQRNSwsrMztGAwGNS8vTx0yZIj62GOPFZZPnTpVBdR169YVlh04cEAF1NWrV5eIE62L/LpLv379Cus//vjjavv27UvFMWXKFNXPz69EWVZuvtrzvXXq0M82q/kGY2F5WlqaCqjz58+/ztETQncV5h253YYwa0FBQTRrVnShbXBwMAADBw4sLGvaVFvfQEnHPv4EFhYW3HPPPeTn5xfWGTRoEIsXL8ZgMGBpaUlMTAxTpkxh7dq1XLx4EbXgXG2vXr1KbN/a2pr+/fsXPm/VSpsNPSYmprBs2rRpPPPMM+Xuh7Ozc4nnZY0uVFW1VLmdtSWv3daSp3/ez5I953mwW+NytyOEOZEEJcyam5tbiec2NjalyjefSQLgrrZeJCfFYDAYcHUtewTfxYsX8ff358477yQtLY233nqL4OBgHB0defPNN4mLiytR38XFpcR5n6vbz87OLixr3LgxjRo1Knc/iiced3d3kpOTS9VJTk4utb8At7X1pUuQO5+uPcXdHf1xsJF/a1E3yF+yqNNy8418t0WbYLVrUANiYzOxsrJi27ZtZQ4o8Pb2JiIiggMHDrB69WqGDh1auC4rK+uGYhg3blypQRjX6tevHxs3bgQgNDSU8+fPk5GRgaOjY2Gd8PBwQkNDS71WURQmDw3lvm92MH9bFE8PCL6hOIUwNZKgRJ32865oLqZorRkLC4WBAwdiMBhISUnhlltuKfM1VxORrW3RMPTo6Gi2bdtGu3btqhxDVbv4hgwZAsDvv//Oww8/DMCFCxfYsmULX331VZmvDwtqwOCW3nyz6QwPdWss/9iiTpC/Y1FnpefkM2t9BB0D3IktKAsJCeGpp55i9OjRvPzyy4SFhZGdnc2xY8c4deoUc+fOJTQ0lEaNGjFp0iTefvtt0tLSmDp1Kg0bNryhOIKCgggKCqp0/UaNGjF+/Hief/55VFXFy8uLadOmERgYWJiwAN566y3eeuutwnNp/701lIFv/MgzM8IZ2tIDgL179+Lk5ISXlxf9+vW7ofiF0IskKFFnfb8lkisZubx1R1NWFiufPXs2LVq04LvvvuPNN9/ExcWFVq1aMX78eEBrOS1fvpynn36a++67j0aNGjFlyhQ2btzI0aNHayX2L774AkdHR1588UUyMzPp168fixcvxs6u6I6/RqOxxPD0EF9nGiYe5Ke5c/ip2L7Onj27RBeiEOZCZpIQddKV9Bz6frCBPs29+OaRznqHU2tOXExl2OdbeGFwC54b3FzvcIQoj8wkIeqnLzdEkJVn4KVbQ/QOpVa19HPhllY+zNsWSXpOfsUvEMKESYISdc75xEx+2nmOkWEBBHub9p1ya8IzA4JJycrjx53RFVcWwoRJghJ1zqdrT6Eo1NsurvYBbvRp7sncLWfJypU77wrzJQlK1CkRcWn870AsY3oG4edqr3c4uvm/gc1JSM/llz3n9A5FiBsmCUrUKZ+vi8DO2pIn+zbVOxRddW3SgK5NGvDtprPk5EsrSpgnSVCizjh5KY2Vhy/wWK8gPEz0Xk+16ekBwVxKzWbFoYt6hyLEDZEEJeqMz9aewtHGiif61O/W01V9m3vS3NuJeVsjqcLlJEKYDElQok44diGF1UcvMa53E9wcbPQOxyQoisK43k04fjGVXZGJeocjRJVJghJ1wqf/nsbFzorxvZvoHYpJuadjQ9wdrJm3NVLvUISoMklQwuwdjklm7YnLPNGnKa721nqHY1LsrC15qFsg/564zLkrmRW/QAgTIglKmLaDi7WlHJ/+ewo3B2vG9gqqnZjMzCM9ArFUFOZvl1aUMC+SoIRpSzmvLdexLzqJDSfjmdC3Kc520noqi4+LHcPb+fHb3hjSsvP0DkeISpMEJczaZ2tP4eFow5geQXqHYtLG9W5Cek4+v+6NqbiyECZCEpQwW3uiEtlyOoGn+jXD0VbuHFOedo3c6BLkzoLtkRiNMuRcmAdJUMJsfbHuNJ5ONjzcPVDvUMzCIz2COJ+YxZaIBL1DEaJSJEEJs3TgXBJbTifweJ+m2NtY6h2OWbi1tQ8ejjb8JLOcCzMhCUqYpVnrI3BzsJbWUxXYWllyf1gA68LjuJSSrXc4QlRIEpQwO0djU1gfHsf4Xk1wknNPVfJg18YYjCpL9lx/ZKQQpkISlDA7X66PwNnOijFy3VOVNfZwoE9zT37Zc458g1HvcIQolyQoYVZOXkrj72OXGNszCBe57umGPNQtkIsp2Ww4Ga93KEKUSxKUMCtfbojA0caScb1kzr0bNailNz4utvy0SwZLCNMmCUqYjTPx6aw8fIGHewTi7igzlt8oa0sLRnVpzKZT8ZxPlPn5hOmSBCXMxlcbzmBrZSH3e6oGo7sEoIDcEl6YNElQwiycu5LJ/w7G8mDXQDzlbrk3zd/NngEh3vy2N0YGSwiTJQlKmIWvN0VgqShM6Cutp+oyqksAcWk5MlhCmCxJUMLkpWXns3RfDPeHNcLX1U7vcOqMAaHeeDnbskS6+YSJkgQlTN7+c0kYjCpP9m2mdyh1irWlBfd1bsR6mVlCmChJUMKkZecZOByTzO3t/Gns4aB3OHXOqLAAjCos3SczSwjTIwlKmLRDMcnk5ht5Us491YggT0d6NPVgyd7zchsOYXIkQQmTlZ1nYP+5ZAI9HGjT0FXvcOqs0V0DOJ+YxY6zV/QORYgSZKZNUXkHF0P4Srh0GNLjwdYJvFtC3/9Ck74l6xoNsONL2P8DJJ8D+wYQejsMfB0cGlRqcwf++IJ7DavxSkuHD1ZBs0EweCq4NiqqdHI1rJ2ubcMrBIbOhMbditbH7oPvBmnb7ftSNRyE60g+D+vegjPrIDcDPIKh25PQ6dGa2d7lY7D5I4jZA+lx4OAB/h2gz0vQqHPJutE7YMMMuHBAe+7fEQZMgcAeANza2hdXe2t+2XOeXsGeZWzruLZv57ZDfi74toHeL2i/z6tSL8CqlyBqK9i5QOcxWiyKoq03GuHbPmDIg4nbwVI+ekTFFFWtdLNe2v/13awwuHK67HX3fg9t7yt6/sczcGBR6XrereGJdWBtX+6mjBs/wGLjjNIrXBrCExvA2QcSz8LsbuDdSktMy5/QksOzB8DeTas/byikxMIze8C6hkYApl2GOf0h7ULpdQNf1xJ4dUqKhq96QF5G6XWWNvD4OvBrpz0/uwl+vBeMeSXrWVjDI8sLv1hM+/MYP+86x87XBtGg+CwdcSfg+yGQk1p6W/fMgfajtMcL79CS08hFcOpv7XdffP2+BbDiOXhoGTQffHP7L+oKpaIK0sUnKs/OVfvAff4IvBoDfSYVrdv8YdHj87uLklPocPjvWRg0VXsedwx2fl3+dpLPwaYPAMiy94Gez2kfdgCpsbDxPe3xmfVgyIX2o7XWQOhwyEqEmL3a+mO/w7kdcMv06ktOqgqRW+DioaKyje8WJad75sCkU9AwrGDd+1rrqrJi92sf9OUJX1mUnLo8Dq9dgMHTtOeGXDjyW1HdVZO05GTnBk9t1RY7N61sVdHvb1SXAHINRn4/EFtyW2te05KTpQ08+gc8dxjcCu7BtfplyM3Ulsgt4NMaWg6HrhO09af+1n7mpMOGdyH4FklOokokQYnKe/QPrTXg1hhsnWHgG2Droq1LPFtUr/gHZO8XwdEDej0HNs4F65eWuxn12P+wULVv/LZNemhdie1HgWeIVuHocq3LyFDQKrAsmNXcquCbvyFH64paOw0a94A2I25ipwukXoQtH8MXHWHhcLh0VCs3GuHo79pjzxAtTmcf6PF0wfo8OP5H+e+dmQi7voWve8F3A7QP+/Ioxe4gHDocbByh5Z1FZfkFQ8YvHChq8ba5F3zbakube7WyhFNw4SAALf1caB/gxpI95yjsVclIgLMbtcdN+kHT/uAeCGHjtLLsZIhYW9A6U7UkBmBVMNOHIVf7ufUTyLwCt5bRIhaiHJKgROXZOpV8bsjVzjUBOPsVlRdvXXgUXLtkYQnuQdrj+HDIz7nuZuJP7yl8bOHgUey9grWfOSmQHAWNuwMKnFoD2SkQsR6s7MC/E+z6WusKu/XdKu1iCYZ8CF8FP4+CT1tr52FSzkOLodCw4DxPUqQWT/H4rn186XDp91ZV7cN/6Xj4OFRrjVw+BkF9tERQnlZ3gr279jh8pdatWTwJNhuo/Szxe6g4ttFdAjh1OZ0D55OL1qnGil9v5wo+bbT4E8/CiRXausCekBIDO76CsPHaOUIhqkDOVIobt/2Loq6mTo8UlWcWGw12tYUFWqsLQDVAVhI4+5b5tpcuxuJ99YlVsXn3rr4etG/3AV3hlrdg/dsws7HWQrtzlvZNfvPH0P4BaNhJq5+fW9TCqkhChNZFeWgxpF/WynzbQYcHoe394FhsIEGJfXUu+3FGsamEUmLh4M/a+ycX3O7Co7nW8mo3GtwCKo7PxR8e+xt+HAF75moLaImi/2sQMuyGYrujvT9vrzzOL7vP0amxu9ayq+zr7/4Kfh2jtTBBa6V1eQL++I/2O+z/ilZuNAKq9oVFiApIghI35uBi7bwCaN/6ez1ftO66A2+Kl5d9fvRITArJmblQ5udXGa/v9aw2Wi7tIrg00kaHrZoExnwY9KZ2nuqvl+FKhDbSrfcL0POZ8vfty4LWkbMf9HxWS0zeLa+zS1Xc1/lDtXNsdm5aq6L9AxDQpfx4rpV6EX55UDsfV1xuptaiyc3Quv2qGJuTrRV3tPPnz0MXeGN4K5yr8nq/9vDcQa3FZOuijeSL3ad15w6dqbXEfnkITv+rjewLHgx3fF4y2QtxDeniE1V38Gftm7Fq1M7xPLC46DwQlPzQKT76Kydd+6lYFo2yu8bcrWdJtSh2zVPxrsCrrwftvNZVVrZa96GlFcSf0kaM9X5Bi+nXMVr338iFWhfTP1O08ybluXqOJy9Tiz875fp1S+xr2nViLVbn6nsbcrX3zkkp6iatrG2fQ+IZ7fHtn8CUS/DAEq1levAn7dzbjcQGjOoaQFaegZWHL5Y8xpV8Pa6NtOQEsGYKeDbXBnKsflnrjuwzSTsvGb5SKxOiHJKgRNUc+An+eFpLTk36wsPLSnb5gPZt+qorBR+kRgMkRWmPvUJLdt0VuJSSzarDF7EJ6FBUmFWsm+lKhPbT1hXcgsqO758pBS2fZ7RrhHJSte6mVndB94lanTMbyt/HF49row4dPLVkN+9W+LyDNiLv6j5c5d5Ei6d4fNc+9m1X9HjCBhj2oXZu7shv2hDwT1vDP29o1xtVRsKposcdH9GG7IcMBceCjtHIzdrPEr+HSsQGdAxwo4WPE7/sOa+tUyyq9PpCV0dQDpmhfXE4s17rguz3srbYuWplQpRDEpSovAM/wp/PaMkpeDA8+JvWlXSttvcXPd76CWRc0b715xZ8Cy9+vdTvE2GaK0xz5ced0RhUlVaDH9Wu0wFtyHpOOhz+FRJOamVtRoBFGX+6ZzbA6X+0IdfW9kWtlasXhV59z4rOfzj7Qp8X4dn9MPYvrRsu/bI2nPzzDjBvGFwsGFxgYQFt7tEeJ5zU4kyPgx2zi7bZ+u6i97Z3h24TtOHeEzZprYu8TO183tc94Js+RSPnyovvqgOLIC8LTv4NGXFamV1BwvTvqJ3fAji6DC4d0Zajy7Qyzxbaxb2gjRyc5ooy3Y3XGx7g0PlkTqTaFA3YiNykxZUUDXvnFe1LcBnDxq+OoGw2CFoM0coUS7Cw0rr3FKXouRDlkL8QUXkb3y8a1RWxFmb4lFz/3GFtGHJAV+2b/YFFWldO+MqiOt6ti1oy1/h59zkGhfrQKCgE+k2GDe9o51m2f15Uydkf+r9a+sVGI/zzOgR0K0qAAV20GSxOrNCGYe9boJU3v7Xy+xzUS1uGfQBHl8L+RdqMCpeOFF0M2/81OPWPdi3U8idKvr7/5JIzXxTn30FbhsyAE39qs25EbdVmfmja//oxdRmvtb4MubDqRW0pLmx80ePbP9ZaadnJ8E3vonILa21dGToHumNzyIIle84z7dZ3iy7U/eGukhWHzgSbMibw3fW1di7qgSVFZSFDtS84+3/Q/oayEqHjw9ffRyGQBCVqyh2fa8ORDywqPdXRdWaRSMzIZVyvIO1Jv/9q1xOtexuyrmjdaMGDtK43Z5/SLz6wSBvm/MS6ojJ7d3jgF+1i0x/u0loed87SEk5V2blo1/+EjdO2U/xaJGcfGP8PrJsOEQVTHXk2r/xUR9Z20G6ktiRGluzWLEvDzjB2lXZdVsweyErWWrK+bbWLZIu32Jr2gzErtKmOYvcXvL5TiamOruVoY8WtbXxZvj+GV4YNxm7cGm2IffR2LSn6tNZamMWnOroq44o2grLzY+AdWlQ+ZIY2aOPq+bEOD93cJQCiXpCpjoTuVFVl2OdbUFX4+/k+KEqxUW8FM0rQT06o16ZtEQk8NHcXn4/uwF0dGuodjqibZKojYfp2nL1C+KU0xvUOKpmchG56NPWgcQMHftkt94kS+pEEJXQ3f1sU7g7W8k3dhFhYKIzqEsCOs1eISihjUlohaoEkKKGrc1cyWXviMg92a4ydtcwuYEru69wICwWW7JVWlNCHJCihqwXbo7BUFB7pHqR3KOIaPi52DAz1Zum+GPIMRr3DEfWQJCihm/ScfH7be57b2vrh61pD92oSN2VUl8bEp+WwITxO71BEPSQJSuhm6d7zpOXk89jVoeXC5AwI8cLb2VabWUKIWiYJSujCaFRZuCOaDgFudGzsrnc44jqsLC24P6wRG0/GcTElS+9wRD0jCUroYuOpOCITMqT1ZAZGhgVgVOG3vTF6hyLqGUlQQhfztkbh42LLbW39Kq4sdBXo4UjvYE8W7z5HvgyWELVIEpSodacup7E1IoFHewRhbSl/gubg0R6BXEzJ5t/jl/UORdQj8ukgat38bVHYWlnwQNfGeociKmlQSx8autmzcEeU3qGIekQSlKhVSRm5/H4ghrs7NKSBYyVvwS50Z2mh8EiPQHaeTST8UmrFLxCiGkiCErVq8Z5zZOcZeax3kN6hiCoaFRaArZUFP+yI1jsUUU9IghK1Js9gZNGOaHo28yDU10XvcEQVuTvacFcHf37fH0tKZp7e4Yh6QBKUqDVrjl3iYko2j/Vqonco4gY92iOIrDwDv+2TC3dFzZMEJWrN/G1RNG7gwMBQb71DETeoTUNXOge6s2hnNEaj3CJO1CxJUKJWHDqfzL7oJMb2DMLSQu75ZM7G9Awi+komG07K/HyiZkmCErVi/rZInGytuD+skd6hiJs0rI0vfq52zNl8Vu9QRB0nCUrUuLjUbFYduch9nRvhbGetdzjiJllbWjC+dxN2RSZy4FyS3uGIOkwSlKhxP+6MJt+oMrZnkN6hiGoyumtjXOys+HaTtKJEzZEEJWpUdp6Bn3adY1CoN0GejnqHI6qJk60Vj/YIYs3xS5yJT9c7HFFHSYISNerPQxe4kpErQ8vroDE9tbkU526RVpSoGZKgRI1RVZX526II8XGmZzMPvcMR1czL2Zb7Ojdi2b5Y4tKy9Q5H1EGSoESN2RWZyImLqTzWKwhFkaHlddGEPk3JMxqZvy1K71BEHSQJStSYeVsjcXew5u6ODfUORdSQIE9HhrXx5ced0aRkyfRHonpJghI14nxiJv+euMwDXRtjZ22pdziiBv2nfzBp2fl8vzVS71BEHSMJStSIhdujsFC0WzSIuq1NQ1eGtvZl3tZIkjJy9Q5H1CGSoES1y8jJZ8ne89zW1g8/V3u9wxG14IVbWpCRm88cGdEnqpEkKFHtftt7nrTsfB7rFaR3KKKWhPg6M7ydPwu3R5GQnqN3OKKOkAQlqpXBqDJ/exQdG7sgHMsuAAAgAElEQVTRqbG73uGIWvT84OZk5xn4dtMZvUMRdYQkKFGt1p64TPSVTB7v3VTvUEQta+blxN0dG/LDjmjiUuW6KHHzJEGJavX91kgautlza2sfvUMROnhuUHPyjSqzN0ToHYqoAyRBiWpzJCaF3ZGJPNYrCCtL+dOqjwI9HBnVJYCfdp2TOfrETZNPEVFtvt96FkcbS0Z2CdA7FKGjF29pgZ21Je+uOqF3KMLMSYIS1eJSSjYrD19kVJfGuMg9n+o1Tydb/m9gMOvC49h8Kl7vcIQZkwQlqsXCHVEYVVWGlgsAxvYKonEDB95ZdZx8g1HvcISZkgQlblpmbj4/7zrHra19CWjgoHc4wgTYWlny2m0tOXU5ncV7zusdjjBTkqDETVu2L4aUrDzG95Z7Pokit7b2oVuTBnz67ymZSFbcEElQ4qYYjSrztkXRPsCNzoFyYa4ooigKbwxvRVJmLh//c1LvcIQZkgQlbsr68DgiEzIY37uJ3PNJlNKmoStjegSxaGc0+88l6R2OMDOSoMRNmbv1LH6udgxr46t3KMJEvXRrCL4udry67Ai5+TJgQlSeJChxw47GprDzbCJjewZhLRfmiutwsrXi7bvacPJyGt/JbOeiCuRTRdyweVsjcbCxZHTXxnqHIkzc4FY+3NbWl8/XneaszDAhKkkSlLghl1OzWXH4AiPDAnC1lwtzRcWm3dEaWysLXvv9CKqq6h2OMAOSoMQN+WFHFPlGuTBXVJ63ix2vDmvJzrOJ/LgzWu9whBmQBCWqLC07j0U7ohnSyodAD0e9wxFm5IGuAfRt4cWMv04QmZChdzjCxEmCElW2ePc5UrPzmdg/WO9QhJlRFIUP7m2HrZUlL/56UKZBEuWSBCWqJCffwNwtkfRs5kGHADe9wxFmyNfVjrfvbsOBc8l8u1lG9YnrkwQlquT3/bHEpeUwsX8zvUMRZuzO9v4Mb+fHp/+e4mhsit7hCBMlCUpUmsGo8u3ms7Rp6ELvYE+9wxFm7p2729DA0YYXfz1Idp5B73CECZIEJSptzbFLRCZkMLFfsExrJG6am4MNH97fnlOX05m5OlzvcIQJkgQlKkVVVb7eeIYgDweGyrRGopr0a+HF2J5BLNgexSa5uaG4hiQoUSlbIxI4EpvCk/2aYWlhHq2n//3vf7Rr1w5bW1uaNGnCJ598Um79559/HkVReOmll0qUh4eH061bN1xdXRk9ejTp6SVnQti8eTMNGzYsVV6WBQsWoChKmXWnTZuGp2dR12lUVBSKohQuzs7OhIWF8euvv163jqOjI82aNeOhhx5iy5YtFcZjCl4ZFkoLHyde+u0QiRm5eocjTIgkKFEhVVX5fO1pfF3sGNGpod7hVMq2bdsYMWIEXbt2ZcWKFYwbN47Jkyfz2WeflVn/+PHjzJs3DxcXl1Lrxo4dS3BwML/++ivHjx/n3XffLVxnNBp5/vnnee+993BycqqRffnoo4/YsWMHy5Yto3nz5owaNYqVK1eWWeevv/7ijTfe4MqVK/Tt25fp06fXSEzVyc7aks9GdSQlM49Xlh2WWSZEEVVVK7uIemrr6Xg1cPJK9YftkbW/8Y3va0sVDRkyRO3Tp0+JshdeeEF1d3dXc3JyStUfNGiQ+vrrr6uBgYHqpEmTCsvT0tJUQI2Li1NVVVV/+eUXNSwsrHD9nDlz1K5du6pGo7FScc2fP18F1LS0tFLrpk6dqnp4eBQ+j4yMVAF1xYoVhWUGg0Ft0aKFetttt123zlVvvPGGCqgbNmyoVGx6m7PpjBo4eaW6eFe03qGI2lFh3pEWlCiXqqp8tvYUvi52jOwSoHc4lXbw4EEGDx5comzIkCEkJSWxY8eOEuVLly7lxIkTvPLKK6XeJzdX63Kyt7cHwMHBobAsNTWVN954g88//7zWBo1YWFjQoUMHoqKiKqw7depU/P39+eabb2o+sGowvncTejbzYPqK4zLLhACki09UYMeZK+yJSuI/A5pha2WpdziVlp2djY2NTYkyW1tbAE6cOFFYlpWVxaRJk5g5cyaOjqWnbWrQoAFNmjRh1qxZJCYmMmfOHMLCwgB4++23GTx4MN27d69yfAaDgfz8/BKL0Vi5WRWioqLw9a14oIqlpSUDBw5k586dVY5PDxYWCh+PbI+NlQXPLzlInswyUe9Z6R2AMF1a6+k0Pi62jAwzn9YTQHBwMHv27ClRtnv3bgASExMLy9577z38/Px4+OGHr/tes2fP5v777+e1116jefPmzJ49m4iICL7//nsOHz58Q/G5uZU9C4eHh0epMqPRSH5+PqmpqcydO5fdu3cza9asSm2nUaNGXL58+YZi1IOfqz3v3tOWp3/ez6x1p3lxSIjeIQkdSYIS17Xj7BV2RyUy/c7W2FmbT+sJ4KmnnmLixIl899133HfffezevZuPP/4Y0FoWAJGRkXz00UesX7++3C66YcOGERcXR0xMDM2aNcPS0pI777yTF154gUaNGjF79mzef/99AF555RX+85//VBjf5s2bC7sNr5ozZw7Lly8vVfeuu+4qfGxtbc2LL77IxIkTKz4IYJYDDm5v58f68EZ8uSGCfiFedA5soHdIQieSoESZireeRpnRuaerxo0bx6FDh5g4cSITJkzAwcGB999/n//7v//Dx8cH0JLJsGHDCA0NJTk5GdBaKzk5OSQnJ+Pq6lqYuBwcHGjRogUAa9eu5dChQyxZsoRDhw7xxhtvsH37dgB69OhB7969adeuXbnxdezYsdSov2tH5l316aef0rt3b5ydnWnSpEmprsvyxMbGFu6vOZl2Zyt2R13h+SUH+evZPjjbyT3H6iM5ByXKtPl0ArsjE5nYr5nZtZ5AayV9+eWXxMfHc/jwYS5fvlx4rujqz5MnT7J8+XLc3d0Ll/Pnz/Pll1/i7u5ObGxsqffNz8/n+eef54MPPsDe3p6NGzcycOBAQkNDCQ0NZdCgQWzatKla9yU4OJiwsDBCQkKqlJzy8/NZv349PXr0qNZ4aoOznTWfjuxAbFIW0/48rnc4QifSghKlGI0qM1eHE9DAnge6mfft3K8mHoCvvvqKnj17EhoaCsDcuXNLXTA7evRo+vXrx8SJE/Hy8ir1ft988w3u7u6MGjWqsCwzM7PwcUZGhsl0q7311ltcuHCBp556Su9QbkhYUAOeGRDMF+sjGBjqze3t/PQOSdQySVCilD8PXeDExVQ+H93BrEbuFbdz5062bt1Khw4dSE1NZfHixaxZs4atW7cW1rk6Gq84Ozs7AgIC6N+/f6l1SUlJTJ8+nTVr1hSW9e3bl5dffpl58+YBsH79embOnFn9O1SBkydP4unpSW5uLpGRkfzyyy/8/fffTJs2jX79+tV6PNXl/wY1Z9PpBF77/QidAt3wc7Wv+EWizpAEJUrIyTfw0T8nae3vwh3t/PUO54ZZW1uzZMkSpk2bhoWFBX369GHbtm20bdv2ht9z6tSp3HnnnXTq1KmwrGPHjnzwwQdMmTIF0GZ0aN++/U3HX1VXp2eys7PDz8+PHj16sHnzZvr06VPrsVQna0sLPhvVgds+38JLvx1i0bhuWJjJVFvi5ilV6I4wjX4LUaPmbY3krZXHWTS+K32al+7iqnWbPtB+9ntZ3ziErhbvPsery48w5baWPNG3qd7hiOpR4TcNGSQhCqVm5zFr/Wl6B3uaRnISosDoLgHc0sqHD/85SURcmt7hiFoiCUoUmrPpLEmZeUweGqp3KEKUoCgKM+5pg4ONJS8vPYzBKB069YEkKAHA+cRM5mw5y53t/WnbyFXvcIQoxdvZjjeHt2L/uWQWbo/SOxxRCyRBCQDeXnkcS0Xh1duk9SRM1z0dGzIgxIsP15zk3JXMil8gzJokKMHmU/H8c/wyzwwMlmG8wqRpXX1tsbRQeGW53DuqrpMEVc/l5huZtuIYgR4OPN6nid7hCFEhfzd7XrutJdvPXOGXPef1DkfUIElQ9dzC7VGcjc/gzeGtzPaiXFH/PNA1gJ7NPJix6gQXU7L0DkfUEElQ9VhcWjafrzvNgBAvBrU0vwlFRf2lKAozR7TDYFR5bfkR6eqroyRB1WPTVxwnN9/Im3e01jsUIaqssYcD/701hA0n4/nfwdIT+wrzJwmqnlpz7BKrDl/k2UHBNPEsfSdZIczBmJ5BdGrsxlsrjnMlPUfvcEQ1kwRVD6Vk5vH6/47S0s+FJ/s10zscIW6YpYXCzHvbkZ6Tz1sr5bYcdY0kqHpoxl/HSczI5cP72mFtKX8Cwry18HHm6QHB/HHwAuvDzef29qJi8ulUz2w5Hc+ve2OY0LcpbRrKjBGibpjYvxnNvZ14/fejpOfk6x2OqCaSoOqRtOw8Xl1+hKaejjw3qLne4QhRbWytLJl5bzsupmbz4d/heocjqokkqHpCVVVe/99RLiRn8cF97czyNu5ClKdzoDtjegTxw85o9kUn6h2OqAaSoOqJpfti+OPgBZ4f3IKwoAZ6hyNEjfjvrSH4u9ozedkRcvINeocjbpIkqHrgTHw6b/5xjO5NG/D0gGC9wxGixjjaWjHjnjZExKUze8MZvcMRN0kSVB2XnWfg/34+gJ21BZ+N6oil3C5b1HH9Q7y5p2NDvt4YwclLcnNDcyYJqo57768THL+Yykf3t8fX1U7vcISoFW8Mb4WznTWTl8nNDc2ZJKg6bOm+GBbuiGZ87yYy156oVxo42jD1jlYcPC83NzRnkqDqqAPnknjt9yP0bObBq8PkJoSi/rmzvX/hzQ3PJ8rNDc2RJKg6KC41m6d+3IePiy2zH+yElcwWIeohRVF45562WCjw2u8y47k5kk+uOiYn38CTP+4jNSufOY+E4e5oo3dIQuimoZs9Lw8NZcvpBH4/IDOemxtJUHWI0ajyyrIjHDiXzCcj29PSz0XvkITQ3SPdA7UZz1ceJ0FmPDcrkqDqkPf/Duf3A7H899YQhrX10zscIUyChYXC+/e2IzPHwPQVMuO5OZEEVUd8vzWSbzef5ZHugfynv9xCQ4jimhfMeL7i0AXWHpcZz82FJKg64M9DF3h75XGGtvZl2p2tURS5GFeIa03s34xQX2deWX5EuvrMhCQoM7fpVDyTfj1I16AGfDa6g8wUIcR12FhZ8NnoDqRm5/HKssMyqs8MSIIyY9sjEpjww16aezvz3aNhMkO5EBUI9XVh8tBQ1p6IY/Hu83qHIyogCcpM7YlKZPzCvQR6OPDj491wdbDWOyQhzMJjPYPoHezJ2yuPczY+Xe9wRDkkQZmhA+eSeGz+Hvzc7Pjp8e40kGudhKg0CwuFj0e2x9bagueXHCQ336h3SOI6JEGZmX3RSTw6bzceTjb8/Hh3vJxt9Q5JCLPj42LHzBFtORyTwrt/ndA7HHEdkqDMyK6zV3j0+114ONqw+InuMju5EDdhaBs/xvduwoLtUSzfH6N3OKIMkqDMxLaIBMbO34Ovqx1LnuyBv5u93iEJYfZeGRZKtyYNeHX5EY7GpugdjriGJCgzsOFkHOMW7KFxAwd+mdADHxdpOQlRHawtLZj9UCfcHWx46sd9JGXk6h2SKEYSlIn742AsTyzcS7C3E4snyDknIaqbp5MtXz/cibjUHP5v8QHyDDJowlRIgjJh87dF8twvB+kc6M7iCTJaT4ia0rGxO+/c04atEQm8vPQwRrkLr0mw0jsAUZqqqnzy7ylmrY9gSCsfvnigo1yEK0QNGxkWQFxqNh/9cwoPRxum3N5Spg3TmbSgTExOvoFJvx1i1voIRoUF8NVDnUokJ1VVeffddwkICMDe3p6+ffty8ODBCt/333//5YEHHiAoKAhFUZg2bVoN7oUQtev48eMMGjQIBwcH/P39efPNNzEYDOW+Zs+ePTz22GMEBwfj4OBASEgICZt/4qHOvswtmHxZ6EsSlAm5kp7DQ9/tYvn+WF4Y3IKZ97YtdTfcmTNn8vbbbzN58mRWrFiBk5MTgwcP5tKlS+W+999//83hw4cL/4mFqCuSkpIYPHgwiqLwxx9/8Oabb/Lxxx8zderUcl+3ZMkSzpw5w+TJk/nrr794+umn+fTTTzmx+F2Gt/Nj5upwFu8+V0t7IcqkqmplF1GDwi+mqr1mrlNbTPlLXXEotsw6WVlZqouLizp9+vTCsvT0dNXT01OdMmVKue9vMBgKH3t4eKhTp06tlrhr3Mb3tUWI63j33XdVNzc3NSUlpbDs/fffV+3t7UuUXSsuLq5U2bfffqsC6smIM+qj3+9SAyevVOduOVsjcYuK8460oEzAhpNx3Pv1dnLyjSx5sgfD2/mXWW/79u2kpqYycuTIwjJHR0fuuOMOVq9eXe42LCzkVy3qptWrV3Prrbfi4lJ0B+nRo0eTlZXFpk2brvs6Ly+vUmUdO3YEICXxCnMe7cywNr68vfI4n/x7SmY/14F8aulIVVXmbjnL+IJrnP58phcdAtyuWz88PBxLS0uaN29eorxly5aEh4fXdLhCmKTw8HBCQ0NLlDVu3BgHB4cq/19s374dCwsLQkJCsLWyZNYDHbm/cyO+WHea6SuOy+i+Wiaj+HSSmp3H5KWHWX30Ere29uGTkR1wtC3/15GUlISTkxOWliVH9Lm7u5OZmUlubi42NjIUXdQvSUlJuLmV/mLn7u5OUlJSpd/n0qVLzJgxg0ceeaSwNWZlacH797bDxd6a77dGcjEli49HdsCpgv9VUT3kKOvgaGwKT/+8n5ikLF67LZQn+jQtNZxVVdUSo5Curi9r2OvVrgcZEivqq+v9X1T2fyI3N5eRI0fi5OTEp59+WmKdhYXC67e3xN/NnhmrjnPvV9v58sGONPdxrpbYxfVJF18tMhpVFmyLZMTX28nJM7JkQncm9G1W5j/Rpk2bsLa2LlwGDRqEu7s7aWlppYbPJicn4+DggLW13BNK1D/u7u4kJyeXKk9JSSmzZXUtVVV59NFHOXbsGH/99Rfu7u6l6iiKwvjeTVg4rivx6TkMn7WVBdsi5bxUDZMWVC2JvpLBy0sPsysykf4hXnx8f3s8nK4/bVHnzp3Zs2dP4XNnZ2diY2MxGAxEREQQEhJSuK6sPngh6ovQ0NBS55rOnz9PRkZGpf4vXnjhBf744w/+/fffCuv3ae7F38/3YfLSw0xbcZz1J+P56L52eMv8mDVCWlA1zGBUWbg9iqGfbeH4hVQ+uLcd88d2KTc5gZaQwsLCCpeQkBB69uyJi4sLv/32W2G9zMxMVqxYwbBhw2p6V4QwScOGDWPNmjWkpaUVli1ZsgR7e3v69etX7mvfe+89Zs2axY8//kjv3r0rtT1vZzvmje3C23e3YXfkFQZ9son52yLJlzn8qp1lFWYUqHRFodl/LomnFu1jyd7z9Az2ZOG4rnRv5nHD54qsrLQG74wZM3BzcyMtLY0XX3yRmJgYFi5ciJOTEwA//PADnTt3ZsyYMYVdHNHR0axdu5bjx4+zfPlyXF1dURSF6OjoUqMCTUr0Nu1nUC994xAmq02bNnz77bds2LABf39/1q5dy6uvvsoLL7xQ4otbcHAwhw4d4q677gLg559/5umnn2bMmDEMGDCAmJiYwsXW1hZHR8frblNRFNo3cmNYWz/CL6Xxw45o1hy7TLC3EwEN5EL4SppeYY3KXCylyoW6VRKXmq2+9OtBNXDySrXrjH/V/x2IUY1GY7W8t9FoVN955x21YcOGqp2dndq7d291//79JerMnz9fBdTIyMhSZdcugYGB1RJXjZELdUUlHDt2TB0wYIBqZ2en+vr6qq+//rqan59fok5gYKA6ZsyYwudjxowp838CUOfPn1/pbRuNRnX1kYtqz/fWqYGTV6rjF+xRj8QkV9Oe1WkV5h1FrfxJPjkbWIHEjFzmbjnLwu1R5BqMjO/dlGcGBsuQ1Jux6QPtZ7+X9Y1DiApk5xn4bvNZvttyltTsfAa39Oa5QS1o28hV79BMVYVdSZKgqsGV9By+3xrJwu1RZOYZGN7OnxcGN6epl5PeoZk/SVDCzKRm57FgWxTfb40kJSuPbk0aMKZnELe08sHaUk77FyMJqiYdu5DCgm1R/HHoAnkGI8Pb+fPswGC5PqI6SYISZiotO4+fd51j0c5oYpKy8HWx46FujRndtbHceFQjCaq6ZeTk8/fRSyzZc57dUYnYW1syolNDHusVRLC3JKZqJwlKmDmDUWVDeBwLd0Sx5XQC1pYKQ1r7cld7f/qFeGFrVW/v9SYJqjrkGYzsjkxk+f5YVh+9SGaugcYNHHikeyAjwwJwdZALZGuMJChRh5yJT2fRjmj+PHSBxIxcnO2sGNralzva+9OzmUep2+vUcZKgblRqdh6bTsaz9sRlNoTHkZqdj7OtFbe382NEp0Z0CXKXqYVqgyQoUQflGYxsi0jgz0MX+OfYZdJz8vFwtKF/iDcDQr3oE+xVH774SoKqrHyDkWMXUtl59gpbTiew8+wV8o0qDRxtGBTqzeBWPvRt7oW9Tb1tjutDEpSo47LzDGw8GceqI5fYfCqelKw8LBTo1Nid/iFe9G3hRSs/l7rYupIEdT2ZufmcuJjKrshEdp1NZF90Euk5+QA083JkcCsfbmnpQ8fG7lhaSEtJN5KgRD2SbzByKCaZDeHxbDwVx9HYVAAcbCzp1NidLkEN6BLkTofGbjjYmP3lK5KgcvINnLuSyZn4DCITMjh5KZWjF1I5G5/O1Vu7NPd2olvTBnRr4kG3Jg1kXi1TIglK1GNxadnsOpvI3qhEdkclEX4pFVUFKwuFYG8n2jR0pY2/C60butLSz8XcrrmsmwnKYFTJyM0nM8dAek4+yZm5xKXlEF+wxKVlcyk1h8iEdGKTsih+jzFfFzvaNHShtb8rrf1d6BTojmcF8+KZi48//php06aRnp6udyjV5vW+2v2t3tmcq3Mkoi5wcnJi2rRpTJo0Se9Qbkhqdh77o5PYG5XE0QspHI1NISFd+99QFAjycKSZlxPNvB0J9nKimbcTzbyccLU3yfNZppWgoq9kMPHH/dpcIqqKUVVRVUr+RPtpLJh30Vis3tXElJ13/UkZLS0UPJ1s8Ha2I8jTkSaejjTz0n4GeTriYmeSv6hq4e/vz8WLF/UOo1pJghLVzc/PjwsXLugdRrVQVZW4tByOxqZw7EIqJy6mciY+nciEDPIMRR/ZXs62NPNyJMDdAX83exq62ePvZo+/mx2+rnbYW1vqMeirwg3WanvQytICfzc7FEXBQgEFBQsLbeJFBbC4Wq4oKEqx5wX1LBQFR1srHG2scLS1xKHgp7uDDV7Otng529LAwQaLenrOaNKkSXWuBSVEdXJycjLb1lNZFEXBx8UOHxc7BrX0KSzPNxg5n5RFRFw6Z+LTiYhL52x8OptPxxOXlsO17RJbKwvcHWxwd7TB3cEadwcbHGwssbO2xM7aouCnJfbWlozr3aT29s8cu/hEPSLnoISoVrn5Ri6nZhObnMWF5Cwup+aQnJlLUmYuiRl5JBU8zso1kJ1nIDvPSHa+AVUFO2sLwt+utlv7mFYLSogqcw3QOwIh6hQbKwsCGjhU6bYgqqqSk28kJ79273klLSghhBB6qLAFVeeu/BJCCFE3SIISQghhkiRBCSGEMEmSoIQQQpgkSVBCCCFMkiQoIYQQJkkSlBBCCJMkCUoIIYRJkgQlhBDCJFV6Jonp06f/DXjWbDh1hj9QN6ZLNg1yPKuXHM/qJcfzxiRMnTp1aLk1VFWVpZqXadOmqXrHUJcWOZ5yPE15keNZc4t08QkhhDBJkqBqxnS9A6hj5HhWLzme1UuOZw2pymzmQgghRK2RFpQQQgiTJAlKCCGESZIEJYQQwiRJghJCCGGSJEHdJEVRbBVFmaUoSoKiKBmKovypKEqjCl7Tt6BerKIoqqIoY2spXJOkKMp/FEWJVBQlW1GUfYqi9Kmgfr+CetmKopxVFOWp2orVHFTleCqK4qcoys+KooQrimJQFGVBLYZqFqp4PEcoivKPoijxiqKkKYqyS1GUO2sz3rpEEtTN+wy4F3gA6AO4ACsVRbEs5zVOwFHgOSCrxiM0YYqijAI+B94FOgLbgdWKojS+Tv0mwF8F9ToC7wGzFEW5t3YiNm1VPZ6ALZAAzAR21UqQZuQGjmc/YD1we0H9v4DfK/rSJcomw8xvgqIorkA88Jiqqj8VlAUA0cAwVVXXVOI90oFnVFVdUJOxmipFUXYBh1VVfaJY2Wlgqaqqr5ZR/31ghKqqzYuVzQVaq6raozZiNmVVPZ7XvHYlkKCq6tiajdJ83MzxLFZ/N7BFVdVJNRRmnSUtqJvTGbAG/rlaoKrqeeAE0FOvoMyFoig2aMfwn2tW/cP1j1+PMuqvAcIURbGu3gjNyw0eT3Ed1Xg8nYGk6oqrPpEEdXN8AQNaF0lxlwvWifJ5ApZox6u48o6f73XqWyGTGd/I8RTXd9PHU1GUp4FGwKLqDa1+kARVBkVR3ikYvFDe0r+8twCk77Tyrj1WFR2/suqXVV5fVfV4ivLd0PEsOC/6IfCQqqrRNRFYXWeldwAm6jPgxwrqnAO6o33D8kQ7F3WVN7C5ZkKrUxLQWqDXfhv1pvS31qsuXad+PnClWqMzPzdyPMX13fDxLEhOi4BHVVX9s2bCq/ukBVUGVVUTVFUNr2DJBPYBecAtV19bMMS8JdpoH1EOVVVz0Y7hLdesuoXrH78dwOAy6u9VVTWveiM0Lzd4PMV13OjxVBRlJNoX3LGqqi6tuQjrPmlB3QRVVVMURfke+FBRlDi0b/CfAIeBtVfrKYoSDnypquqXBc+dgOCC1RZAY0VROgCJqqqeq819MAGfAIsKRjptA55CuwHcNwCKovwAoKrqowX1vwGeURTlM+BboBcwFm2Yv6j68aTgbw+0SySMBc9zVVU9XpuBm6gqHU9FUUajtZxeAjYrinK19ZWrqmpiLcdu/vS+IZW5L4AdMAstOWUCK4CAa+qowIX5/10AAACMSURBVLRiz/sXlF27LNB7f3Q6hv8BooActG+sfYut2whsvKZ+P2B/Qf1I4Cm998GUlhs4nmX9LUbpvR+mslTleBY8L+t4bqztuOvCItdBCSGEMElyDkoIIYRJkgQlhBDCJEmCEkIIYZIkQQkhhDBJkqCEEEKYJElQQgghTJIkKCGEECZJEpQQQgiT9P9dXBbyvtTWTwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"az.plot_posterior(trace, var_names=['δ'], ref_val=0);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Paramter estimates"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>sd</th>\n",
" <th>mc_error</th>\n",
" <th>hpd_2.5</th>\n",
" <th>hpd_97.5</th>\n",
" <th>n_eff</th>\n",
" <th>Rhat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>π__0</th>\n",
" <td>0.33</td>\n",
" <td>0.04</td>\n",
" <td>0.0</td>\n",
" <td>0.25</td>\n",
" <td>0.42</td>\n",
" <td>2000.39</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>π__1</th>\n",
" <td>0.38</td>\n",
" <td>0.05</td>\n",
" <td>0.0</td>\n",
" <td>0.30</td>\n",
" <td>0.48</td>\n",
" <td>2099.90</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>δ</th>\n",
" <td>0.05</td>\n",
" <td>0.06</td>\n",
" <td>0.0</td>\n",
" <td>-0.08</td>\n",
" <td>0.17</td>\n",
" <td>2094.05</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat\n",
"π__0 0.33 0.04 0.0 0.25 0.42 2000.39 1.0\n",
"π__1 0.38 0.05 0.0 0.30 0.48 2099.90 1.0\n",
"δ 0.05 0.06 0.0 -0.08 0.17 2094.05 1.0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pm.summary(trace).round(2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another approach is a GLM:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"with pm.Model() as affection_status_glm:\n",
" \n",
" μ = pm.Normal('μ', 0, sd=5)\n",
" θ = pm.Normal('θ', 0, sd=1)\n",
" π = pm.math.invlogit(μ + θ*gender)\n",
" like = pm.Binomial('like', N, π, observed=y)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"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: [θ, μ]\n",
"Sampling 2 chains: 100%|██████████| 3000/3000 [00:03<00:00, 928.44draws/s]\n"
]
}
],
"source": [
"with affection_status_glm:\n",
" \n",
" trace_glm = pm.sample(1000, njobs=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here you get an effect size, on the logit scale (so perhaps harder to interpret)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4FFXbwOHfpPcKIYQ0IEDovUoXEFBUioIKihRFRV8VFf1QAfW1AmIBsSGgiCigvICAgPQivQkBAgRCaCG9kWSz8/0xSUhCOklms/vc1zVXdmfOzDwnk+yzZ+bMGUVVVYQQQghTY6V3AEIIIURhJEEJIYQwSZKghBBCmCRJUEIIIUySJCghhBAmSRKUEEIIkyQJSogKoChKPUVRViuKkqwoSrSiKF8qiuKsd1xCVGeSoIS4Q4qieAKbAXtgKPA8MAhYomdcQlR3NnoHIIQZeBlwBO5XVTUNQFGUM8B+RVF6q6r6t67RCVFNSQtKiDs3DPg9JzkBqKp6ADibvUwIUQ6SoIS4A4qiOAANgZOKotjknYBjQEt9IxSi+lJkLD4hyk9RFD8gqpgip1VVbVRV8QhhTuQalBAV4wNgRYF5/wWCqz4UIcyDJCgh7kwsYARiVFXdn3eBoihZwA1dohLCDMg1KCHugKqqN4HTaNehCmoMHKnaiIQwH5KghLhzvwH3K4rikjNDUZTuaKf3lukVlBDVnXSSEOIOZd+oewQ4j3YtyhP4BDikquogPWMTojqTBCVEBVAUpR7wOdALSAOWAq+pqpqia2BCVGOSoIQQQpgkuQYlhBDCJEmCEkIIYZIkQQkhhDBJkqCEEEKYpLKMJCG9KYQQQlQUpaQC0oISQghhkiRBCSGEMEmSoIQQQpgkSVBCCCFMkiQoIYQQJkkSlBBCCJMkCUoIIYRJkgQlhBDCJEmCEkIIYZLKMpKEEKIQGQYjG05cY/XRy5y4kkhkbCqKouDhaEsLf3e6NajJ0Db+uDvZ6h2qENVKWZ4HJUMdCZFHllHlt/2RzNxwmuikdHzdHGgb5EmQtxOKAtcT0zl4MY6z0Sk42FrxWMcg/tOnAW4OkqiEoBRDHUmCEqIcLsWl8tzPhzgSGU/bIE8m9g6he4OaWFvd/j93PCqBhbsiWHbwEt7Odrz3YHP6N/PVIWohTIokKCEq2u6zMTz380EyDUbefbAZD7TyQ1FK/F/jeFQCb6w4xrGoBMZ2rcvrA0KxtZbLwMJiSYISoiJtPR3N+EX7CfRy4ptRbalX06VM66cbsnh/zUkW7r5A71Af5jzaBkc760qKVgiTJglKiIqyM/wGYxbso15NF34e1xFPZ7tyb2vxPxd484/jtA7wYMGYDnJdSlgiSVBCVISwq4kM+2o3dTwcWfJUJ7zuIDnl+PPYFV5YcohWAR4sGtsBJzvpVCssijwPSog7dT3pJmMX7MfJzpoFY9pXSHICGNi8Np8/0pqDF+MYv2g/6YasCtmuEOZCEpQQxTBkGZm4+BCxKRl8/0R7ars7Vuj2BzavzSfDWrIzPIbXlx+jDGc0hDB7ck5BiGLM3HCavRGxfDq8Jc393StlH0Pb+nMlIY0Zf50m0MuJl/o2rJT9CFHdSIISoghbTl3nqy1neaRDIINb+1fqvp7rFcKFmFQ+23SGZnXc6dukVqXuT4jqQDpJCFGIuJQM7pm9DQ8nW/43sSsOtpXfFfxmZhbD5u3iQkwqq5/vSpC3c6XvUwgdSScJIcrjrZXHiUvN4NPhraokOQE42Frz1WNtsVIUnvv5IBkGY5XsVwhTJQlKiALWHb/K6qNX+M/dDWjqVznXnYoS4OXEx8NacDwqkZkbTlXpvoUwNZKghMgjJd3A9FX/0ri2GxN61Nclhnua+vJIh0C+2XaOXWdv6BKDEKZAEpQQeczeeJorCTd578Fm2Og4Tt5b9zUm2NuZycuPkpYh90cJyyQJSohsJy4nMn9nBI90CKRtkKeusTjZ2fDBkOZExqbx6cbTusYihF4kQQkBGI0qb/5xDA9HWyb3b6R3OAB0qufNIx0C+G77OY5HJegdjhBVThKUEMDS/ZEcvBjP/w1sjIdTxQxlVBFeH9AYbxd7Ji8/iiFLevUJyyIJSli8G8npfLg2jI51vRjSpo7e4eTj7mjLO/c35d/LiXy/47ze4QhRpSRBCYv3wZ9hpGYY+O/gZqV68GBV69/Ml75NavHpxtNExqbqHY4QVUYSlLBou8/GsPzgJZ7qXo8QH1e9wymUoii880BTrBSF99ac0DscIaqMJChhsTIMRt5aeRx/T0cm9mqgdzjFqu3uyHO9Qlj/7zV2nJF7o4RlkAQlLNa3288Rfj2Zdx9opstj19PT05k0aRI+Pj44Oztz7733EhERUWT5sV3r4uek8vhzr9C+fQfc3d3x9fVl8ODBnD4tXdGF+ZEEJSzShZgUPt90hoHNfekV6qNLDC+88AILFixgxowZLFu2jBs3btC3b19u3rxZaHkHW2tGt3Dh4u5V+DbtyLJly/j666+5cuUKHTt2JDIysoprIETlktHMhcVRVZXRP+zjwIU4Nr7cA193hyqP4dKlSwQHBzN//nwef/xxAKKioqhbty5z585l3Lhxha6XnJzM+B8PcPzaTba80hNvF3tiY2MJDAzk1VdfZerUqVVZDSHuhIxmLqq30aNH065dO9asWUOTJk1wcnLi3nvvJTY2lvDwcHr16oWzszPt2rXj6NGjuesZjUY+/PBDQkJCsLe3p2HDhixcuBCAP49dZevpaPq4XGLUsEH4+Pjg5uZGp06d+Ouvv/Ltf9q0adSoUYNDhw7RqVMnnJycaN26Ndu3b7+jeuXsZ8iQIbnz6tSpQ9euXVm7dm2R67m4uPDu0DakZWQx4y9tMFkvLy+CgoK4fv36HcUkhKmRBCVM3sWLF3n77bd57733+Oabb9i1axdPPfUUI0aMYMSIESxbtgyDwcCIESNyH5n+/PPP89577/HUU0+xZs0aBg8ezJgxY/h1xR9MX/Uvzeq4Udc+lUGDBvHjjz+yfPlyunTpwoABA9i5c2e+/aempvLEE0/w9NNPs3z5cuzt7Rk8eDCpqbe6fBuNRgwGQ7FTVtatMfXCwsLw9/fHxcUl374aN25MWFhYsb+PEB9XRnUOYum+SE5fSyI6Oprw8HCaNGlyp79qIUyLqqqlnYSock888YRqbW2thoeH58579dVXVUBduHBh7rw1a9aogHrixAn1zJkzqqIo6oIFC/Jta9SoUapfSFM1+PXV6uGLcfmWZWVlqZmZmWq/fv3UJ598Mnf+1KlTVUDdtGlT7rxDhw6pgLp27dp8caKdBi9y6tGjR275cePGqS1btrytvlOmTFFr165d4u8lNjldbfb2OnXsgn3qqFGjVC8vL/XGjRslrieECSkx78gj34XJCw4Opn79W4++CAkJAaB37963zYuKiuLs2bNYWVkxePBgDAZDbpmGrTvz4+KfmdLen5YBHly6dIkpU6awceNGrly5ktv6uuuuu/Lt39bWlp49e+a+z2mpXLp0KXfetGnTmDhxYrH1cHXNf59VYTcFq6paqpuFPZ3teLpHPaZ+NJu4DT+xfPlyvL29S1xPiOpEEpQweR4eHvne29nZ3TY/Z97Nmze5ceMGWVlZuLsX/rDBx1q4YTQauf/++0lKSuKdd94hJCQEZ2dn3n777duu5bi5uWFldetseN595QgMDMTf37/YeuRNPJ6ensTHx99WJj4+/rb6FqVW3HFiN86j+YPP8OCDD5ZqHSGqE0lQwux4eXlhY2PDzp07cxPLioOX+GFnBK/e04j6gXUIDw/n0KFDrF27lv79++eum5aWVq59jhkzJrcTRlF69OjBli1bAAgNDSUyMpKUlBScnZ1zy4SFhREaGlri/nbt2sUTox6jz+CRnAkZyKaT1+nTpFa5YhfCVEmCEmand+/eZGVlkZCQQN++fTl9LYllv0cz6O6uvPRIWxRFyU1E9vb2uetduHCBnTt30qJFizLvs6yn+Pr16wfA77//zsiRIwG4fPky27dvZ+7cucVu599//+W+++6jf//+LFnyPf0/28HH68PoFeqDtZXpjSUoRHlJghJmp1GjRkyYMIERI0Yw6ZVXWRXlQFZ8MrUdrBg/fh7fffcdoaGh+Pv7M2nSJN59912SkpKYOnUqdeqUbzTz4OBggoODS13e39+fsWPH8uKLL6KqKjVr1mTatGkEBQXlJiyAd955h3feeSf3Wtr169fp378/Li4uvPDCCxw6sJ97fZOYuf40ny/P4qWHehe1SyGqHUlQwizNmTOHhg0b8uHsOVy7FIGbmxtbIpsxduxYQGs5rVixgueee45hw4bh7+/PlClT2LJlC8ePH6+SGD///HOcnZ15+eWXSU1NpUePHixZsgQHh1s3DhuNxnzd00+cOJHbOaNXr175tjdtV0smDjmIrY6PqheiIslIEsJsHbuUwOC5O7mvRW1mj2itdziVatPJa4xduJ+PhjZnePtAvcMRojRkJAlhmVIzDLz862G8XeyYfn8zvcOpdL1DfWjp784Xf4eTKU/eFWZCEpQwO6qq8saKY4RHJzPjoZa4O9nqHVKlUxSFF/s05FJcGssPXCp5BSGqAUlQwuws2BXBysOXmdS3Id0a1NQ7nCrTs1FNWgZ48MXf4WQYpBUlqj9JUMKs7D0fy3/XnKRP41o82zNE73CqlNaKakBUfBrLD0orSlR/kqCE2biWeJNnFx8kwMuJWcNbYmWB9wT1bFiTVgEefCmtKGEGJEEJs5CSbuCpRftJSTcwb2Rb3BzM/7pTYaQVJcyJJChR7WVmGXl28UGORSXw+SOtaeTrWvJKZqxHw5o0r+PO11vPkmWUu0NE9SUJSlRrqqoyeflRtp6O5v3Bzekr49GhKArP9apPREwqfx67onc4QpSbJChRrX207hQrDkbxct+GjOggN6jm6NfEl/o1nZmzOZwy3IwvhEmRBCWqra+2nGXe1rM81jGQ53tbVo+9klhZKTzTM4Swq0lsPiWPghfVkyQooZ/DS7SpHOZsDuejdWHc39KPdx5oVqqH/FmaB1r5UcfDkTmbz0orSlRLkqCEfhIitamM5mwO55P1p3iwlR+zHm4pj5gogq21FU91r8eBC3HsPR+rdzhClJkkKFGtfLHpDJ+sP8Xg1nWY+XArbGTk7mINbx9ADRc75m45q3coQpSZ/HeLauOzjWeYueE0Q1rXYcZD0nIqDQdba568qy5bT0dzPCpB73CEKBNJUKJamL3xNJ9uPM3QNv58IsmpTEZ1DsLV3oa5W8L1DkWIMpEEJUzeZxvPMHvjGR5q68/Hw1pIciojNwdbRnUOYu3xq5yNTtY7HCFKTRKUMGmfbzrDpxtPM6ytPx8NleRUXmO61sXO2op5ci1KVCOSoITJ+mLTGWZt0E7rfTS0hUUO/lpRarjYM6J9AH8cjuJa4k29wxGiVCRBCZM0Z3O41iGiTR05rVdBxnSti8GosnBXhN6hCFEqkqCEycm5z2lI6zp8Mkw6RFSUIG9n7mniy+J/LpKaYdA7HCFKJAlKmJTvtp/Lvc9JeutVvHHd6pKQlskyeSy8qAYkQQmTsfJwFO+tOcnA5r5yn1MlaRvkSasAD77fcV4exSFMniQoYRJ2nLnBK78doVM9L2Y93EqSUyVRFIXx3epxISaVDSeu6R2OEMWSBCV0dzwqgad/3E/9mi5883g7HGyt9Q7JrN3TtBZ1PBz5fsc5vUMRoliSoISuEtIyGP3DPjyc7Fg4poPFPqq9KtlYWzGma132RcRxODJe73CEKJIkKKGbDIORlYcvk27IYuGY9tRyc9A7JIsxvH0Arg42fLtdWlHCdEmCErowGlXWHr/CjeR0vny0DSE+rnqHZFFc7G14tEMga49dITI2Ve9whCiUJCihi1kbThN+PZkeDX3o0bCm3uFYpCe6BGOlKCyQG3eFiZIEJarcysNRfLk5nOZ13GkT6KF3OBbLz8ORe1vUZum+SBJvZuodjhC3kQQlqlTY1UQmLz9Kh2AvejeuJY9q19m4rvVITjewdG/Zn2wsRGVTVLXUN+vJXX2W4PASCFsNV49CcjTYu4BPY+j+KtTtfqvc1eNw4Ae4sBsSoyArEzyDocVD0OlZsLG/bdOJNzN54MudpKQbWP1CV3yWPwTxFyAlGjKzr4MM/gZaDr+1kqrC9hlwYCHcTITgrnDvDHDzu1Vm52ew4W0Ysx4CO1XO7wW0um7+L1w+pL33aw29pkBQ54rf16fNIeFi0csfmAutH4NDi2Hls0WXcw+El44Vu6sR83bSOfo3XvDchRJ3HuxdoX5vuHsqeATcKhh9CtZM0urvXBO6vgRtn7i1PD0ZvmgDPk3g8T9KWVFhwUr8dmpTFVGIamT7TIg5c+u9IQ3OR8P5bTD0e2g+TJt/5i/Y913+da//Cxv/hYgdMHJ5vkWqqjJ52VEuxqayZHwnfFwd4NI/YEgvPp6jS+Hv96D1KGjYH34dBRnJ8MT/tOUpMbBtJjQdUrnJ6dxW+GkoGPOcCovYDgsHwagV+ZN3VbBzrrByMx3nUydrGdzInpEaA8d+g4id8NRmcPWFLAP88igkX4cRi2H3HFj1AtRoAEFdtPV2zIKUG3DP++WrkxAFyCk+kZ+DO/R+E148Bm9cgm6Tbi3b9smt14oCje+HMX/BlKsw+k+wd9eWhW+EqAP5Nvv9jvOsPX6Vyf0b0aGulzbTtyWE3gc93yg6ntPrtJ8dn4bG94FPUy1ZZmS3uLa8D1np0Hf6HVY8j8ybcGwZpMbemrdmkpacHDxgwg5tcvDQ5q2ZVPS2CjJmwam1EHu++HIvHYNpCfmn2q20Zfbu0KCf9rr1Y7eXe2Tpre3kfKEoypWj1Dm/DIA9dp1QJ1+Ax1eCYgVJl7UWI0DsWYgJ1xJxvZ7QJrvldHq99jPhEuyeC21HQ60mpf1tCFEsSVAiv8dXaqfzPAK1Uz293wJ7N21ZbJ57ZtqPh+E/QmBHsHWE4Lvyn5rL8wF84EIsH64No1+TWozvVu9WmQZ9oXZLcM9zGqmgrOwWi7Wd9tPGDlC1xBB9Gg4sgM4TtXjv1NVj8OdrMLMRLB8L6Una/MuHbrUqmw0F3+ba1GyoNu/Gabh8uPhtx56DjdNhVhNYMkL7QC+Ly4fgSvY+Wo4AO6eiy+6fr/20soU2jxe/3YgduS/nJXdj75UsLQHVaqbNPL5CS6pZGdp76+wbqXNO4ebM3zhNOza9ppS2RkKUSE7xifzsXfK/z8rQPqAAXGsXXQ7AkOdBeNllE1Izef7nQ/h5OPLJQy3L3ikiqIt2TezkKu3D8epx7cPTwR2WjwOnGtq1kPK6maC1lg79eOvakr279sHu5K29v3LkVnnvkMJfXz0Kfq3ybzszDU78T9t2xA5ABRtHLbF51y9bnDlJB6DdmKLLJVyC8A3a69B7wcWn+O0a0nJfutjbMn/neTrW8761PCNZS67eDbTrThd2a9cmw9Zoy4O6aK3lY8ug33vg7I0QFUUSlCjers8hM0V73WZU0eVizmofUqB9cAd2RlVV3vj9KNeT0ln+TBfcHcsxjFH78doH4N/vapNnXXhwLpzdrF0He2DurWRpyMhuYZVCxE4tcZxYqXXQUKwhpK/WOgm9D2zzjGqRGnPrtb1r4a9Tom+9vnIEDi7SruPcTAAUCOysbbvpYHBwK9vv4GYiHMu+phd0F/iEFl32wAJQjdrr9mNL3rZP09yX//HYydAT9bh2OIZa147fKpMaq11rGvINrHgaZoRopwA7PqP9rn4YAF71tNOwoF2vsrLWTgMLcQckQYmiHV4Cm7MveAd3g7teLLxcfCT8OFj7oLdzhWE/gJUVS/de5M9jV3l9QCgtA8p5v5ONHQybD/fNhvREcPcHoxG+7qZdk2n1qHbBfvssLZF4h8DAj7VeaEWJ2AEL7tVe12oGLR+BFg8X3doosqdr3vnZH8bxF+Hr7A4TnnWh03NaYvIMKkut8zu69NaXhOJaT1kGOPij9tq7Qek6bjToq/VGvHyIBnFbOWq/FQp2wLPO/pio3xsmnYKESK01ZecE//4OF3fDiCVa6+1/z8PFPdop2WaDYcAnxZ+OFKIYcg1KFO7wz1r3ZdWofft/ZMmt6w95xUdqH/bxF8DOBR77FWq3IPx6EtNXnaBrSA2eynvdqbwc3LTkBHBoEVw7Dv0/gLObYP3/QY2G8PAircXy6xNab7KiKHlGS09P0hJfzvWmwjjXyF8+93Xy7WUUK3KTVWZq9rYTS1XFIh1YoP10qqF1TCnKqT8h+ar2ut2Tpdu2lTWMXKH1knTyJl1xYK/ahMz6fW+VcauTp7yVlmztnLQW68ZpULcHhA6EFeO15D/gI2j1CBz6CbZ9XJaaCpGPJChxu0OLYeVzWnKq213rMm5fyFh58RdvJSd7dxj1BwR1Id2QxfNLDuNoZ82sh1tiVZHPdkpP1lp1TR7Urn+c3azN7/wcNLkfmg3REsKlfUVvI6gzTNgJHSdoZbd+pN2/830/7VpPWoERvmu3vPU6Jrzw174ttJ/u/vDCQa33o2INu7+EeV3hq7tg1xeQVMZnMF38R0vGAK1HFn8KM+c6lY2j1rIsLScveOBLeO0cp8ee5uH0KcRGZyc67xCtm3lh/vlK+4LS/wPtNOSlfVqLtP1YrScowNm/Sx+HEAVIghL5HfoJ/jdRS04hfeDR3wq/lyZvcnL0hCdWQkB7QBtn7+SVRGY81AIfNwf44V6Y5q7dfJpXZprWXTwj5da8jGTt3qabCYXHt2OWtqzvO9p7JftPOKd1Z2WT/2dRfJtp3/QnndJOIdbrpX3Arn4JZjTUWmE3s1s+fq21U2YAx5drvf2uHtNeg9Z6y9tBwqse3P02vHQcHv1Vu04TfQr+ehNmNYafhmkf7KVx4IfsF4rWhbsosefh3BbtddPB2jEp6NBi7ThMc4fz22/NP7IUrv0LmWk0d03ia68l1ErMvrm344TC95dz/1mbx6FWU60lhnLrdKBVgeMhRDnIX4/Ib8tHty6yh2+E/9bKv/w/R7VTPIcWa0kKIC0OvumZW+QNoH3gk/QOvbf4fe3//vZEtOZlbQrqCk+uyb8s516bTs/euqbTsL/WkePAQu1D+eQqcPQC/3alq6+NvdarrtlQrT6HFsPhxXDiDy0J5nRouHemdqPuzXitRZTDylZbVhgra2h4jzYlR8ORJdoXgPANEBeRf5SGwqTFadd4QLv+41W36LIHfiD3mlhx16kKs38+RO7JfXtP9s9rtXtTq10RHS22ZF+bzGkp2TlD3W5wYZd2n1dOt/uG/csWixB5SAtKVIrulTFC+cbpWsLIe/Nw8F1w/5cQHQaLHtCuBT3yS+EtiJJ4BEKvN7QkPHIFOObp2FGvBzyxSussYuusTcHdtHml6YzgUhPuegEm7tVubvYMLnmdw0tudd0vLukYMrTEClCreW5LttQa3gM1Q7VriDaOqL7NmW07luezXtauORWUc/9Z91fyX58b/LWWkH5/Wlve6bmiO9YIUQoyFp+oMNNX/csPOyP4eXxHutSvUfIKW7MvoPd4rXIDE2X2/Y7zvLv6BCufu6v8PTCFKF6JF6elBSUqxO6zMfywM4InOgeVLjkJk/ZwO39c7G2Yv7OEIZmEqESSoMQdS0438OqyIwR7OzF5QDE3kYpqw9XBlofbBbDm6BWuJtwseQUhKoEkKHHH3v/zJFHxacx4qCVOdtLvxlw8eVcwRlVl0e4IvUMRFkoSlLgjW09H8/M/FxnfrR7tgr30DkdUoAAvJ/o18eXnvRdJy8jSOxxhgSRBiXJLSTfwxvKjhPi48HLfhnqHIyrBmK51iU/NZPnBMo6+LkQFkAQlym3mX6e5nHCTj4Y2x8HWuuQVRLXTPtiT5nXcmb/zPEajdOQVVUsSlCiXY5cSWLDrPI91DKRtkJzaM1eKojC2a13ORaew9Ux0ySsIUYEkQYkyM2QZeeP3o3i72PNaf+m1Z+4GNq+Nj6s983dIl3NRtSRBiTJbsCuC41GJTB3UpHzPeBLVip2NFU90CWb7mRuculrMqO9CVDBJUKJMouLTmLXhNL0a1eTe5rVLXkGYhUc7BOJgayWtKFGlJEGJMvnvmhOoKrzzQLOyP75dVFueznYMaePP74ejiElO1zscYSEkQYlS2302hj+PXeXZnvUJ8JKnpFqaMXcFk2Ewsvifi3qHIiyEJChRKllGlXdWn6COhyPju1fAE3JFtRPi40qPhjVZtPsC6Qa5cVdUPklQolSW7ovk5JVE/m9gY7nnyYKN7VqXG8nprD5yRe9QhAWQBCVKlJCWyYy/TtGhrhcDmxfx+G9hEbo1qEEDHxe+33GeMjyqR4hykQQlSvTFpjPEpWYwdVAT6Rhh4RRFYUzXupy4ksiec7F6hyPMnCQoUazzN1JYsCuCEe0DaOrnrnc4wgQMbl0HL2c7vpcu56KSSYISxZrx1ynsbKx4uW8jvUMRJsLB1prHOgayKewaETdS9A5HmDFJUKJIRy/Fs+boFcZ1q0dNV3u9wxEmZFSnIGysFBbsitA7FGHGJEGJIn20LgwvZzvGd6urdyjCxPi4OTCohR+/7o8kIS1T73CEmZIEJQq1/Uw0O8NjmNgrBFcHGW9P3G5st7qkZmTx054LeocizJQkKHEbo1Hlo3Vh+Hs68linQL3DESaqqZ87PRrWZP6O8/LEXVEpJEGJ26w5doXjUYm83Lch9jZyU64o2rM96xOTksGv+yP1DkWYIUlQIh9DlpFZG04T6uvKA63q6B2OMHEd6nrRJtCDb7adIzPLqHc4wsxIghL5rDx8mfM3Unipb0OsreSmXFE8RVF4tmcIUfFp/O/wZb3DEWZGEpTIZcgy8sXfZ2hS241+TWrpHY6oJnqH+tColitfbT2L0SjDH4mKIwlK5Fp5+DIRMam82KeBDGkkSs3KSuGZnvUJv57MxpPX9A5HmBFJUALI33rqK60nUUb3tahNgJcjc7aclUFkRYWRBCUAaT2JO2NjbcWEHvU5EhnPtjM39A5HmAlJUEJaT6JCPNQ2AD93B2ZvPC2tKFEhJEEJaT2JCmFnY8WzvUI4dDGe7dKKEhVAEpRy2hHiAAAgAElEQVSFyzKqzNkSTmNpPYkK8FA7f2lFiQojCcrCbThxlXPRKTzbs760nsQds7ex5tleIRy8GM+OcGlFiTsjCcqCqarK3C1nCfZ2YmDz2nqHI8zEQ+38qe3uwOyNZ6QVJe6IJCgLtjM8hqOXEni6R30ZNUJUmJxW1IELcdKKEndEEpQFm7slHB9Xe4a0kTH3RMV6OPta1Iy/5FqUKD9JUBbqcGQ8u87GMK5bXRmxXFQ4extrXuzbkCOR8fx57Kre4YhqShKUhZq7ORx3R1se7RikdyjCTA1t40+jWq58sj5MRjoX5SIJygKduZbEXyeu8UTnIFzsbfQOR5gpayuFyQMaERGTypK9F/UOR1RDkqAs0Fdbz+Joa83ou+rqHYowc70a+dCxrhefbTxDcrpB73BENSMJysJcikvlf4cvM6JDAF7OdnqHI8ycoii8MbAxMSkZfLPtnN7hiGpGEpSF+W77eRQFxnerp3cowkK0CvBgYHNfvtt+juuJN/UOR1QjkqAsSExyOr/su8iDrerg5+GodzjCgrx6TyiZWUY+WndK71BENSIJyoIs2BVBusHI0z2k9SSqVt0azozpWpflBy9x4EKc3uGIakISlIVITjewaPcF+jWpRYiPq97hCAv0fO8G1HKzZ9r//iVLHg0vSkESlIX4Ze9FEtIymdCjvt6hCAvlYm/D/w1szLGoBH7bH6l3OKIakARlAdINWXy7/Ryd63nTOtBT73CEBbu/pR8dgr34eP0pElIz9Q5HmDhJUBZg5aHLXEtM55me0noS+lIUhWn3NyU+NYNPN57WOxxh4iRBmTmjUWXetrM09XOjW4MaeocjBE383HisYxCLdkdw7FKC3uEIEyYJysz9deIa56JTmNBDHkgoTMcr9zTC28We11ccxSDj9IkiSIIyY6qq8tXWswR5OzGgma/e4QiRy93Rlun3N+Xfy4n8sDNC73CEiZIEZcZ2n4vhSGQ8T3Wvh421HGphWgY086VP41rM2nCayNhUvcMRJkg+tczYvK3nqOFiz9A2/nqHUmX++OMPWrRogb29PXXr1mXWrFnFln/xxRdRFIVXXnkl3/ywsDA6duyIu7s7I0aMIDk5Od/ybdu2UadOndvmF2bBggUoilJo2WnTplGjxq1rgxERESiKkju5urrSrl07fv311yLLODs7U79+fR577DG2b99eYjymQlEU3nmgKVYKTPnjuDzYUNxGEpSZOh6VwLbT0YzpGoyDrWU8kHDnzp0MGTKEDh06sGrVKsaMGcPkyZOZPXt2oeVPnDjB/PnzcXNzu23Z6NGjCQkJ4ddff+XEiRO8//77ucuMRiMvvvgiH3zwAS4uLpVSlxkzZrB7926WL19OgwYNGD58OKtXry60zJ9//slbb71FTEwM3bt3Z/r06ZUSU2Xw83Dk1Xsase10NCsPX9Y7HGFqVFUt7SSqkWcXH1Cbvr1OTUjL0DuUom35SJsqSL9+/dRu3brlm/fSSy+pnp6eanp6+m3l7777bvXNN99Ug4KC1EmTJuXOT0pKUgH1+vXrqqqq6i+//KK2a9cud/k333yjdujQQTUajaWK64cfflABNSkp6bZlU6dOVb29vXPfnz9/XgXUVatW5c7LyspSGzZsqA4cOLDIMjneeustFVA3b95cqthMgSHLqD7w5Q615fT16rWENL3DEVWnxLwjLSgzFH49iT+PXeGJLkG4OdjqHU6VOXz4MH369Mk3r1+/fsTFxbF79+5885ctW8bJkyd5/fXXb9tORkYGAI6O2oC6Tk5OufMSExN56623+Oyzz6qsV6SVlRWtWrUiIiKixLJTp07Fz8+PefPmVX5gFcTaSmHmwy1Jy8hi8vKjcqpP5JIEZYbmbNYeSDi2q2UNCnvz5k3s7PI/48re3h6AkydP5s5LS0tj0qRJfPjhhzg7O9+2HS8vL+rWrcsXX3xBbGws33zzDe3atQPg3XffpU+fPnTq1KnM8WVlZWEwGPJNRmPpulhHRETg61tyT0xra2t69+7Nnj17yhyfnurXdOGNAaFsPhXNL/tkGCShked9m5nzN1JYeTiKcd3qWdwDCUNCQti3b1++eXv37gUgNjY2d94HH3xA7dq1GTlyZJHbmjNnDg899BD/93//R4MGDZgzZw7h4eF8//33HD16tFzxeXh4FDrf29v7tnlGoxGDwUBiYiLfffcde/fu5YsvvijVfvz9/bl27Vq5YtTT452D2XDyGu+uPsFd9WsQ6O2kd0hCZ9KCMjNzN4dja23FuG6W9zj3CRMmsHLlSr799lvi4uJYv349M2fOBLSWBcD58+eZMWMGs2fPLvYU3YABA7h+/TqnTp3i5MmTBAYG8vLLL/PSSy/h7+/PnDlzCAwMJDAwkLlz55Yqvm3btrFv37580/jx4wst+8ADD2Bra4u3tzdvvvkmL7/8Ms8880yp9lNdT5FZWSl8Mqwl1lYKk347LCOeC2lBmZPI2FR+PxTFqM5B+Lg66B1OlRszZgxHjhzhmWee4amnnsLJyYmPPvqI559/nlq1agHw+uuvM2DAAEJDQ4mPjwe01kp6ejrx8fG4u7vnJi4nJycaNmwIwMaNGzly5AhLly7lyJEjvPXWW+zatQuAzp0707VrV1q0aFFsfK1bt76t11/Bnnk5Pv30U7p27Yqrqyt169a97dRlcaKionLrW934eTgy/f6mvPzrEeZtPctzvUL0DknoSFpQZuSrrWexUhSe7m6Zg8JaW1vz5ZdfEh0dzdGjR7l27VrutaKcn6dOnWLFihV4enrmTpGRkXz55Zd4enoSFRV123YNBgMvvvgiH3/8MY6OjmzZsoXevXsTGhpKaGgod999N1u3bq3QuoSEhNCuXTsaNWpUpuRkMBj4+++/6dy5c4XGU5UGt67DfS1qM2vDaXm4oYWTFpSZuJKQxm/7IxnePgBfd8trPeWVk3gA5s6dS5cuXQgNDQXgu+++u+2G2REjRtCjRw+eeeYZatasedv25s2bh6enJ8OHD8+dl5p6a+SDlJQUkzmt9s4773D58mUmTJigdyjlpigK7w9pzpFL8byw5BB//qcb7o6W0xtV3CIJykx8vikcwKIfSLhnzx527NhBq1atSExMZMmSJaxfv54dO3bklsnpjZeXg4MDAQEB9OzZ87ZlcXFxTJ8+nfXr1+fO6969O6+99hrz588H4O+//+bDDz+s+AqV4NSpU9SoUYOMjAzOnz/PL7/8wrp165g2bRo9evSo8ngqkpuDLZ+PaM1D83bzxoqjzHm0jQx2bIEkQZmBiBsp/Lo/kpEdA/H3tNyeT7a2tixdupRp06ZhZWVFt27d2LlzJ82bNy/3NqdOncr9999PmzZtcue1bt2ajz/+mClTpgDaiA4tW7a84/jLKmd4JgcHB2rXrk3nzp3Ztm0b3bp1q/JYKkPrQE9euacRH64NY8neSB7tGKh3SKKKKWU4NWEa5zDEbV5YcogNJ66x9bWe1atzxNaPtZ89XtM3DmGyjEaVJ37Yy97zsax4tgtN/dz1DklUnBKbxNJJopo7cTmR/x25zJN3BVev5CREKVhZKXw6vBWeTnY889NBeUy8hZEEVc3N/OsUbg42FttzT5i/Gi72zB3ZhisJaby49BBGuT/KYkiCqsb2R8SyKew6T/eoj7uT9HIS5qtNoCdvD2rK5lPRfP73Gb3DEVVEElQ1paoqn6w/RQ0Xe568K1jvcISodCM7BjKkTR0+23SGzWHX9Q5HVAFJUNXUxpPX+ed8LC/cHYKTnXTGFOZPURTeH9ycxr5u/OeXQ1yMkafwmjtJUNVQhsHIf9ecIMTHhUc6SNdbYTkcbK2ZN7ItiqLw9E8HSMvI0jskUYkkQVVDi3ZHEBGTypv3NsbWWg6hsCyB3k7MHtGKsKuJvLFCnh9lzuTTrZqJSU7ns01n6NmoJj0b+egdjhC66NXIh1f6NeKPw5eZu+Ws3uGISiIXL6qZWRtOk5qRxZv3NtY7FCF09WzP+py5lsQn609Rv6YL/ZuV/EBHUb1IC6oaCbuayJK9FxnVKYgQH1e9wxFCV4qi8OHQFrQK8OClpYc5HpWgd0iigkmCqiZUVWX6/07g6mDLf+5uoHc4QpgEB1trvnm8LR5OtoxftJ/rSTf1DklUIElQ1cRvBy6x+1wMr/VvhKeFPcpdiOL4uDrw7ePtiE/N5KlFB7iZKT37zIUkqGogOimd/645SYdgLx5pL93KhSioWR13Ph3eksOR8UxeLj37zIUkqGrgndUnSMvI4v0hzbCykmfiCFGY/s1q80q/hqw8fJnZG2U4JHMgvfhM3Oaw66w6cpkX+zSQjhFClOC5XiFExKTy2aYzBHg5Maytv94hiTsgCcqEpaQbePOP44T4uPBMTxmtXIiS5AyHdCUhjdeXH6W2uwN3hdTQOyxRTnKKz4S9/+dJouLT+HBIc+xtrPUOR4hqwc7Giq9GtqVeTWcm/HiAU1eT9A5JlJMkKBO1/t+rLP7nIk93r0e7YC+9wxGiWnFzsOWHJzvgaGfNmAX7uJ4o3c+rI0lQJuha4k1eX36UZnXcmNSvkd7hCFEt1fFwZP7o9sSlZjBm4T5S0g16hyTKSBKUiTEaVV7+9TA3M418NqI1djZyiIQor2Z13JnzaBtOXE5k4s8HMWQZ9Q5JlIF8+pmY73acY2d4DFMHNaF+TRe9wxGi2usV6sO7DzZj86lo/u/3Y3KPVDUivfhMyIELcXyy/hT9m/oyvH2A3uEIYTYe6xjEtYSbfP53OLXcHOTUeTUhCcpEXEu8yYSfDuDn4ciHQ5ujKHJDrhAV6aW+DbmWmM4X2UlqZKcgvUMSJZAEZQLSDVk8/eMBUtIN/DS2Ix5OMtaeEBVNURT+O7gZN5LTeXvlcWq62nNPU3lEhymTa1A6U1WVt/44zuHIeGY93JJGvjJahBCVxcbaii8ebU0Lfw9eWHKI/RGxeockiiEJSmeLdl/g1/2XeKF3CP2b1dY7HCHMnpOdDfNHt6eOhyNjF+7nzDW5kddUSYLS0brjV5m+6l/6NPbhxT4N9Q5HCIvh5WzHwjEdsLOx4on5e7mSkKZ3SKIQkqB08s+5GF745RAtAzz4/JHWMkq5EFUswMuJH0a3J/GmgdHz95GQlql3SKIASVA6OHklkXGL9hPo5cT8J9rjZCd9VYTQQ7M67swb2ZZzN5J5atF+edihiZEEVcUiY1N5fP5eXOxtWDSmgzwdVwiddW1QgxkPteSf87G8/OthsoxyI6+pkK/uVehCTAqPfLOHDIORxRM64+fhqHdIQgjggVZ1iE5K5701J3FzOMYHQ+ReRFMgCaqKnItO5pFvs5PTuI40rCXdyYUwJeO61SM+NZMvN4fj7mjL6wNCJUnpTBJUFQi/nsQj3/6DqqoseaoTob5ueockhCjEpH4NSbyZydfbzuHmaMtzvUL0DsmiSYKqZMcuJfDkgr0oisIvT3WSx7YLYcIURWHaoKYkpmXyyfpTuDrY8HjnYL3DsliSoCrR32HXmPjzITyd7PhxbAfqyejkQpg8KyuFTx5qSUpGFm+v/BdrK4XHOsq4fXqQXnyV5Kc9Fxi3cD/1a7rw+3NdJDkJUY3YWlvx5aOtuTvUhym/H2fJ3ot6h2SRJEFVsCyjygd/nuTNP47Tq5EPvzzVCR9XB73DEkKUkb2NNXNHtqFXo5q8seIYv+6L1DskiyMJqgLdSE7n8fn/8PW2c4zqFMTXo9ribC9nUYWoruxtrPlqZFt6NKzJ5BVHWfzPBb1DsiiSoCrIgQux3Pf5DvZHxPHxsBa8+2AzbKyL//Wqqsr7779PQEAAjo6OdO/encOHD5e4r6lTp9K8eXPc3NxwdXWlXbt2LF26tKKqIoTZOXHiBHfffTdOTk74+fnx9ttvk5VV8qgR+/fv5/57B7D61QFc/uJRxg1/gMlfLa+CiAVIgrpjWUaVb7adZfjXe7C3tWLFs114uF3pnob74Ycf8u677zJ58mRWrVqFi4sLffr04erVq8Wul5iYyOjRo1m6dCnLly+nTZs2jBgxgmXLllVElYQwK3FxcfTp0wdFUVi5ciVvv/02M2fOZOrUqcWuFxkZSZ8+fTAYDCxatIilP/+El5MNn7z4OP+3aJM8Or4qqKpa2kkUcDEmRX1o3i41aPJq9alF+9T41IxSr5uWlqa6ubmp06dPz52XnJys1qhRQ50yZUqZY+nSpYs6aNCgMq+nqy0faZMQlej9999XPTw81ISEhNx5H330kero6JhvXkFfffWVamVlpcbFxeXOi74RoypWVqpXv2fVN38/phqyjJUau5krMe9IC6ocVFXl1/2RDPhsOycuJ/LJsBbMG9kWd0fbUm9j165dJCYm8vDDD+fOc3Z2ZtCgQaxdu7bMMXl7e5ORkVHm9YQwd2vXruWee+7Bze3WDfIjRowgLS2NrVu3FrleZmYmNjY2uLjc6oHr7uaKrY0NXUO8+XHPBcYt3EfSTRkFvbJIgiqjc9HJjPp+L68tO0qzOm6se7EbD7ULKPOQKGFhYVhbW9OgQYN88xs3bkxYWFiptmEwGIiPj2fx4sX89ddfTJgwoUwxCGEJwsLCCA0NzTcvMDAQJyenYv/Xhg4dipOTE5MmTeL69etcv36dl156CU9PT76b9jz/HdyM7WduMPSrXVyMSa3salgk6WJWSjczs5i75SzztpzF3saKdx5oysiOQeV+jlNcXBwuLi5YW1vnm+/p6UlqaioZGRnY2RU90vmePXvo3LkzADY2Nnz55Zc8+OCD5YpFCHMWFxeHh4fHbfM9PT2Ji4srcj0/Pz82b97Mfffdx+effw5A7dq1Wb9+PTVr1uSxmlDX25lnFh/kwbk7mftYGzrV8660elgiaUGVQFVV/jx2hXtmb+PzTWcY0NyXTa/04PHOwaVOTqqqYjAYcqec3kOFtbrU7AuvJbXImjdvzr59+9iwYQMTJ05k4sSJLFmypIy1E8IyFPW/Vtz/2ZUrVxg2bBht27Zl7dq1rF27lrZt23Lvvfdy8aJ2426XkBr8/mwXPBxtefTbPXy+6Yw8rqMCSQuqGAcuxPHfNSc4eDGeRrVcWTyuI3eF1CjzdrZu3UqvXr1y3/fo0YOHH36YpKQksrKy8rWi4uPjcXJywta2+OtZzs7OtGvXDoA+ffqQkJDA5MmTeeSRR8ocnxDmzNPTk/j4+NvmJyQkFNqyyvHJJ59gMBhYtmxZ7v9j7969adCgATNmzMhtVdWr6cLKiXfx5h/HmbXhNLvO3mD28Nb4ussN+ndKElQhTl1N4tMNp1n371Vqutrz0dDmDGsbgHU5T+e1bduWffv25b53dXUlKiqKrKwswsPDadSoUe6yws6Xl0abNm344YcfyMzMLDG5CWFJQkNDb7vWFBkZSUpKSrH/a2FhYTRt2jTf/5OdnR1Nmzbl7Nmz+cq6Otgye3gruobU4O2V/zLgs22892BzBjb3lUd23AE5xZdH+PUkJv58kP6fbWNH+A3+c3cDtrzSk+HtA8udnIDcm2lzpkaNGtGlSxfc3Nz47bffcsulpqayatUqBgwYUOZ97Ny5E39/f0lOQhQwYMAA1q9fT1JSUu68pUuX4ujoSI8ePYpcLygoiOPHj+frHZuens7x48cJDg6+rbyiKDzULoBVz3eljqcjz/18kHEL93M5Pq1C62NRStMXXTXz+6AOX4xTn/lpvxr8+mq18Vtr1Y/XnVTjUtIrfb/vv/++6ujoqH755Zfqxo0b1YEDB6re3t7q1atXc8ssXLhQtba2ViMiIlRVVdWIiAi1V69e6rfffqtu2rRJXblypTp69GgVUL/66qtKj7lCyX1QogrExsaqvr6+ap8+fdQNGzaoX3/9ters7Hzb/Yb169dXx4wZk/t+//79qo2NjTpw4EB19erV6qpVq9T+/furNjY26uHDh4vdZ6YhS/1m61k19M21apO31qrfbz+nZhiyKqV+1ViJecdiT/GpqsrW09HM23qWPedicXWwYUKP+ozvVg8v56J7z1Wk119/HaPRyAcffEBMTAzt2rVjw4YN1KpVK7eM0WgkKysrt/OEh4cHfn5+vPfee1y9ehUPDw+aNGnCmjVrGDhwYJXELUR14unpyaZNm5g4cSKDBg3Cw8ODl156iWnTpuUrl7cDE2in5tetW8f06dMZNWoUoHVO2rBhAy1btix2nzbWVozvXo/+zXx584/jvLP6BD/tucBr/RtxT1M57VdaSs4HXymYRdeUdEMWa45e4Ztt5wi7moSvmwPjutVlRIdAXGRg16q19WPtZ4/X9I1DiEqkqiqbTl7nw3VhhF9Ppk2gB6/1D5Uu6VBilraYBHUt8SaL91zg570XuZGcQQMfF57uUZ/7W/phZyOX4nQhCUpYEEOWkWUHLjFrw2muJ6XTPtiTZ3uF0LNhTUttUVl2glJVlQMX4liwK4J1x6+Spar0buTDE12C6RpSo9w32YoKIglKWKC0jCyW7rvIN9vOcTnhJk393BjbtS4Dm9fGwda65A2YD8tMUPGpGfx+KIql+yIJu5qEq4MND7cL4PHOQQR5O+sdnsghCUpYsAyDkT8ORzFv61nORafg7mjLsLb+PNIhkBAfi3gCt+UkKKNRZc+5GH7ZF8m6f6+SYTDSwt+d4e0DeLBVHXlwoCmSBCUEqqqy+1wMi/+5yF//XiUzS6V5HXfub+nHvS1q4+fhqHeIlcW8E5Sqqvx7OZHVR6+w+uhlLsWl4eZgw5A2/jzcLoAmfm4lb0ToRxKUEPlEJ6Xzx6EoVh29zNFLCQC0DvSgR8OadG9Yk5b+Hnd0T6aJMb8EZcgycigynk0nr7Pu+BUiYlKxsVK4K6QGQ9rU4Z6mvpZ2Hrf6kgQlRJEibqSw6shlNoZd5+ileFQV3B1taR/sSasAD1oFeNIiwB03h2p7c371T1BGo8q5G8nsj4hj97kYtp6OJj41ExsrhY71vLivhR/9m/riWUX3LokKJAlKiFKJS8lgR/gNtp2O5sDFOM5Fp+Quq+VmT7C3M8HezgR4OeLtYo+Xsx01XOxwtrfBztoKe1tr7KytsLOxws7aiowsIxkGIxlZRjKzf2YYjKQbjGRm3Zq0MmpuGVWFRzsGVlS1SkxQJnNhxmhUiUvN4EJsKuHXkzl7PZlT15I4HBlPfKr2QLAaLnbcHVqL3qE+dGtYozp/cxBCiFLzdLZjUEs/BrX0AyAhNZMjl+I5FpXAuegUImJS2BR2jRvJlfvQUjtrq4pMUCWq0hZUZGwq//f7MbKMKgajiiHLyM1MIzeS04lJycg3TL2djRX1ajjT0t+DtsGetAvypG4NZ4u5X2DmzJlMmzaN5ORkvUOpNG9211q9722TJwEL0+Di4sK0adOYNGmS3qGUy83MLGJTMohNySAmJYOUdIPWCjIYSTdkkW4wYjCq2Fgp2Ntkt6hsrLC11lpWtjZW2Fvfmmebp9Vla6Nga21FDRf7igrXtE7xXYxJ5T9LD2FjpWBtpWBjZYW9jVbhGq521HSxp46nEw18XAjwcjKni4Fl5ufnx5UrV/QOo1JJghKmqHbt2ly+fFnvMCyBaZ3iC/R24vdn76rKXVZbkyZNMvsWlBCmxsXFpdq2nsyRyXeSEGZMOkkIYclMqwUlRD7uAXpHIIQwYdKCEkIIoYcSW1AyjLcQQgiTJAlKCCGESZIEJYQQwiRJghJCCGGSJEEJIYQwSZKghBBCmCRJUEIIIUySJCghhBAmSRKUEEIIk1TqkSSmT5++DqhRueFUGj/AXIYnNqe6gHnVx5zqAuZVH3OqC5hHfW5MnTq1f7ElVFU1+2natGmq3jFIXcy/PuZUF3OrjznVxRzrU9Qkp/iEEEKYJEtJUNP1DqACmVNdwLzqY051AfOqjznVBcyvPoUqy2jmQgghRJWxlBaUEEKIakYSlBBCCJMkCUoIIYRJkgQlhBDCJFX7BKUoir2iKF8oinJDUZQURVH+pyiKfwnrTFMURS0wXS1QRskud1lRlDRFUbYoitLUBOvyhqIo+xRFSVQUJVpRlFWKojQrUGZBIfXdUwnxP6soynlFUW4qinJAUZRuJZTvkV3upqIo5xRFmXCn26xIZdm3oihDFEX5K/sYJCmK8o+iKPcXKDO6kOOgKoriYGJ16VlEnKEFyg1VFOWEoijp2T8HV3Y98uy7LPUp7O9fVRQlJU+ZUtW5EurRPfv/PCp7f6NLsU5zRVG2Zn8uRSmK8raiKEqBMrodmwql941YdzoBX6HdUd0XaANsAQ4D1sWsMw0IA3zzTDULlJkMJAFDgWbAr9n7cTWxuqwHnsyOsTnwO3AV8MpTZgGwoUB9vSo49uFAJjAeaAx8ASQDgUWUrwukZJdrnL1eJjC0vNvUuT6fAa8DHYAQYCqQBXTLU2Z0dp3zHgdfE6xLT0AFmhSI1TpPmc6AAZiSvc0p2e87mmB93Av+zoGzwA9lqXMl1WUg8D4wDEgFRpdQ3i37//vX7P/5oWifU5NM4dhU+O9H7wDu8OC6AxnAY3nmBQBG4J5i1psGHC9muQJcAabkmeeY/YfwtCnVpZDtuGR/MA7KM28BsLqSj8U/wLcF5p0BPiii/EfAmQLzvgN2l3ebetaniG3sBWbmeT8aSK7s2Cvg2OR8WNcoZptLgQ0F5m0ElphafQpZ/67s+nUpS52roF7JpUhQzwCJgGOeeW8CUdy6bUi3Y1PRU3U/xdcWsAX+ypmhqmokcBLoUsK69bKbx+cVRflFUZR6eZbVRfv2lHe7acC2Umy3vO6kLnm5op26jSswv6uiKNcVRTmtKMq3iqL43GnAORRFsUOL/68Ci/6i6Ng7F1J+PdBOURTbcm6zQlTgvl25/Tg4KopyQVGUS4qirFYUpfUdhFqiO6zLfkVRriiKsklRlF4FlhV1/KrDsRkP/Kuq6q5ClhVXZ1PQGdie/XmUYz3a2HzBecpU+bGpDNU9Qfbmph8AAASsSURBVPmitRZuFJh/LXtZUf5B+zY7AO2P1RfYpSiKd57t5mynLNu9E+WtS0GfoZ0W3J1n3jrgceBuYBLaaai/FUWxL3e0+dUArCnb78u3iPI22dsrzzYryh3vW1GU5wB/4Mc8s08BY4AHgEeAm8BORVEa3GnAxShPXa6gfVMfCgxBi3uToijd85Qp6viZ9LFRFMUdeAj4tsCi0tTZFBT1e89ZVlyZyj42Fc5G7wAKoyjKe2jnTYtT3LcbBa25XihVVdcW2N8e4BzwBDArb9GybLfQQCq5LgX2NQvoCnRVVTUrZ76qqr/kKXZMUZQDwAXgXmBFabZdSmX9fRVWPme+UkyZqhr+pFz7VhRlKPAJMEJV1Qu5G1PV3eT54qAoyi60LxPPAy9URMDFKHVdVFU9hfYBnWO3oijBwCtoZxHKvM1KUN59j0RLcHm/OJSlzqaguP+b4spUu2GDTDJBAbOBn0oocxHohPbHVgOIzrPMhzL8Uamqmqwoyr9AzjfZnB59vkBkge0W/GZSkiqpi6IonwIjgF6qqp4rrqyqqpcVRbnErfreqRtorb+C39CK+31dLaK8AYhB+4cq6zYrSnnqA+Qmpx+Bx1VV/V9xZVVVzVIUZT8VdxwKU+66FPAP2t9XjqKOn8kem2zjgeWqqsaWomzBOpuCon7vcKv+eh2bCmeSp/hUVb2hqmpYCVMqcACtN0/fnHUVrVt2Y6Cw88uFyu7mG4rWzAc4j3aQ+xYo060s262quiiK8hnwKNBbVdWwUtS3BlAnT33viKqqGWjx9y2wqC9Fx74b6FNI+f2qqmaWc5sVorz7VhTlYbQvI6NVVV1W0n6yuwa3oIKOQ2Eq8PfYivxx7q6AbZbZndRHUZQOQEtuP71XlIJ1NgW7gW4Fbk3oi9b7NyJPmSo/NpVC714adzqhdc2OQvuwaw1spkDXbLQu5RPzvJ8B9EDrDNERWI3WMyYoT5nJ2fOGoHXn/IWq6WZe1rrMyY6zN/m7x7pkL3fJrm9ntIuoPdH+gC9VZF3Quv5mAOPQkupnaL2SgrKXLwIW5Smf0818dnb5cdnrF+xmXuQ2K/nvqqz1GYH2BeM/BY5D3u7+U4F7gHpoH37zs9fpYGJ1eRF4EK1l1xT4AO300JA8ZbqgtXbfQPty90Z2Xaqqm3mp65Nnve+A02T3diuwrMQ6V1JdXLL/FlqhdTN/O/t1YPbyD4BNecq7o315/gXtc2kI2v9/3m7muh2bCv/96B1ABRxgB7T7IGKyD/AqIKBAGRWYlud9TrLJQEsIy4EmBdZR0LqjX0G7mL0VaGaCdVGLmKZlL3dE68FzPbu+F9C6nQdUQvzPon2LS0f7lts9z7ItwJYC5XsAB7PLnwcmlGWbVfC3Ver6ZL8v7DjkLfNp9u8/Pft4rAc6m2BdXgPCgTQgFtgODCxkm8PQvjBloPU2rdQP8zv8W3NFS2KvFbG9UtW5EurRs4i/mwXZyxcAEQXWaY522v8m2ufTVAokXT2PTUVO8rgNIYQQJskkr0EJIYQQkqCEEEKYJElQQgghTJIkKCGEECZJEpQQQgiTJAlKCCGESZIEJYQQwiRJghJCCGGS/h/dDSDWDfwj0gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"az.plot_posterior(trace_glm, var_names=['θ'], ref_val=0);"
]
}
],
"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.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment