\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 2 chains for 10_000 tune and 10_000 draw iterations (20_000 + 20_000 draws total) took 32 seconds.\n"
]
}
],
"source": [
"with model3:\n",
" idata3 = pm.sample(10000, tune=10000)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"pm1 = idata3.posterior['pm1'].mean().item() # mean value of model indicator variable"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.67195"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pm1"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Posterior: p(model 1|data) = 0.67\n"
]
}
],
"source": [
"print(f'Posterior: p(model 1|data) = {pm1:.2f}')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Posterior: p(model 2|data) = 0.33\n"
]
}
],
"source": [
"print(f'Posterior: p(model 2|data) = {(1-pm1):.2f}')"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = 2.05\n"
]
}
],
"source": [
"print(f'Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = {pm1/(1-pm1) * (.5/.5):.2f}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So we have no good evidence that would allow us to choose between our 2 hypotheses. The data isn't particularly consistent with our \"null hypothesis\". A priori, the alternative hypothesis entails many credible values of theta that are much more consistent with the observed data (e.g., theta = .2). However, this alternative hypothesis also entails many values of theta that are **highly** inconsistent with the observed data (e.g., theta = .9999). So the \"null\" suffers because there is poor agreement with the data (i.e., likelihood) whereas the alternative hypothesis suffers because it is too agnostic about the possible values of theta.\n",
"\n",
"Let's compare this to a traditional, frequentist procedure to compare these models. We will first find the most likelihood of theta permitted each model, find the likelihood that each of these values yields, and then take the ratio of these likelihoods.\n",
"\n",
"To get a quick approximation of the maximum likelihood associated with our alternative hypothesis, we can plot the posterior and request a mode from the kernel density estimate."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAE5CAYAAACkkJZkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7IklEQVR4nO3dd3xUVf7/8ddJ75WEVBIIkAACAQQEKQIKIoiKBXRVWL4uX921rVh21wKy+7Og7td1195w1VV2QUVsKEonICtNIKGXkIT03mfm/P6YEElIQhImc2eSz/PxmAfMnXvvfO4h3HfOvefeq7TWCCGEEB3NxegChBBCdA0SOEIIIexCAkcIIYRdSOAIIYSwCwkcIYQQdiGBI4QQwi4kcIRoglLqWqXUA42mXaaU0kqpy230HfFKqUVKqV62WF8z33G1UupfSqmDSimLUmpdR32XEOcjgSNE064FHjjfTBcoHlgIdFjgYN2OZGArcKoDv0eI83IzugAhRIf6jdbaAqCU2mR0MaJrkx6OEI0opZYCc4DoukNoWil1/KxZfJRS/1BK5SmlcpVSHyilghqtw00p9UelVJpSqloplamUekEp5VX3+WXA2rrZvzvrey6r+3y2UuqHuvWXKaV2KqXmtHVbzoSNEI5AejhCnOvPQBgwHJhRN60aCKz7+9+AL4BbgERgCWDGGlJnfABcDTwLbAH61a03Hrge2AH8DngZuBfYXrfc/ro/ewHLgWcACzAOeEsp5a21fs1mWyqEHUngCNGI1vqIUioXqNFabz0z/UzvA9igtb6n7u/fKqUSgTuUUnO11lopNRaYBczRWv+zbr41SqkC4AOlVLLWepdS6ky4pJ79PXU1PHXW97oA64BI4C5AAkc4JTmkJkTbfdno/c+AJ9C97v2VQA2wou7QmptSyg34tu7zcef7AqVUH6XUR0qpDKC27nUH1h6VEE5JejhCtF1Bo/fVdX961f0ZDngAZc0sH9rSypVSfsB3QAXwB+AI1gC7C5jXjnqFcAgSOELYXj5QBYxt5vPM8yw/CogDxmqt60eW1fWShHBa8gMsRNOqAe92LvsN8AgQqLX+/jzfQRPf41P3Z+2ZCUqpYOCadtYjhEOQwBGiafuBEKXUXcB/sfZYWkVrvU4p9RGwXCn1V+BHrCPN4oGrgEe01geBg4AJmFc3oKAaOIB1VFsJ8LJSaiHgCzwG5PHLSLlWUUrFYR1tB9ZDeRal1A1177drrU+0ZX1CXAgJHCGa9hZwCfAUEAScAOa2YflbgXuwnnN5FGuYHAdWA9kAWut8pdTdWHtD6wFXYEJdYF0HvIB1aHQm1qHYIVjvTNAWE4B3G037T92fvwaWtnF9QrSbkkdMCyGEsAcZFi2EEMIu5JCaEE6o7mLQln5h1Fprs73qEaI1pIcjhHN6h18uCG3q1dLoOCEMIedwhHBCSql4oFsLs5RqrQ/YqRwhWqUzBI7Tb4AQQnQiqrkP5JCaEEIIu5DAEUIIYRcSOEIIIexCAkcIIYRdSOAIIYSwCwkcIYQQdiGBI4QQwi4kcIQQQtiFBI4QQgi7kMARQghhF3K3aNGizKJKvtl7mrJqE8PighmdEIpSzd65QgghmiWBI5pkMlt4bf0RXvr+MDVmS/30IT2CeP7GwSSE+RlYnRDCGcnNO8U5as0W7l+2iy/3ZDFtUCSPTEkiPMCTz3ZmsGT1AUxmC2/cfjGX9Ao1ulQhhONp9hCIBI5ooMZk4Z6PdrB6XzZ/uiqJ+eMSGnyeXlDBr5duJ7Ookn/OG8HF8SEGVSqEcFASOF2FxaJ5d8tx3tl0jLyyaqKCvLl+aDRzL+2Jn2fLR1CrTWZ+9+EO1qTmsPDq/vz60p5NzpdTWsWs17eSV1rNh78ZyaCYoA7YEiGEk5LA6Qq01jz62V7+te0kl/YO5aKoQH7OKGbLkXy6B3jy2LT+TB8U2eRJ/4oaE3d9sIP1B3P587UXcdslcS1+V2ZRJTe+lkJ5jYll80eRGOHfUZslhHAuEjhdwYqfTrHgP7u5c3wCj1yZWB8sO08W8vjKvezNKGF0QihPXN2fpIiA+uUO55Rx38c7Sc0q4emZA5k1vEervu9kfgU3vr4FswX+c+coenbz7ZDtEkI4FQmczq6kqpaJz68jNsSHFXeOxsWl4b+52aL517YTPP/tQYoraxnZM4R+kQGkF1Sw7mAuPh6uvHTzECYkhrfpew/nlDHr9RQ83VxY9r+jiA3xseVmCSGcjwROZ/fM12m8vuEIn/9uDANjApudr7C8hg+3neCLPVmkF1QQHuDFpKRw7rosgVA/z3Z99/7MEma/kUKQjwcfz7+EqCDv9m6GEML5SeB0ZhU1Ji556nvG9g3j5VuGGlLDrvQibntrG6F+Hnw8fxQRgV6G1CGEMFyzgSO3tukEPtuZSUmViV+PjjeshuTYIJbOG0FuaTW3vLmVgvIaw2oRQjgmCZxO4P2tJ+gfGcCwuGBD6xgWF8zSeSM4VVTJXR/8RI3Jcv6FhBBdhgSOkzuUXUpqVgmzhsc6xD3OhseH8NwNg9h2rIDnVqcZXY4QwoFI4Di5L3/OQimYelGE0aXUuyY5mlsv6cFbm46x5Uie0eUIIRyEBI6T+3JPFiPiQwgPcKyT9H+6qh/xob788ZOfqTaZjS5HCOEAJHCc2OGcMg7llDFtUKTRpZzDx8ONRTMGcCK/gqWbjxtdjhDCAUjgOLF1B3IAmJjUtos17WV83zAmJoXzjx8OU1xRa3Q5QgiDSeA4sfUHc+kd7kdMsONe3f/QlERKq028u+WY0aUIIQwmgeOkKmpMbDtawGV9w4wupUX9IgO4on933tl0jNIq6eUI0ZVJ4DiprUfzqTFbGJ/o2IED8LsJvSmpMrHip1NGlyKEMJAEjpNafyAXb3dXhjvBA9CSY4MYHBvE+1tP0AlupSSEaCcJHCe17mAuoxJC8XJ3NbqUVrntkjiO5JaTcjTf6FKEEAaRwHFC6QUVnMivYFyfbkaX0mrTB0US5OPOB1tPGF2KEMIgEjhO6EwvYXRv5wkcL3dXbro4ltX7sskpqTK6HCGEASRwnNDWo/mE+HrQJ9zP6FLaZNbwWMwWzee7M40uRQhhAAkcJ6O1ZtvRAi7pFeIQN+tsi4QwPwbHBPLZrgyjSxFCGEACx8mcKqwko6iSUb1CjS6lXa4dEs3ejBIOZZcaXYoQws4kcJxMyhHr+ZtLnDRwpg+KwtVFSS9HiC5IAsfJbD2aT6ivB73bcf7miy++QCnF8ePHbV9YI9XV1SxYsIDw8HB8fX2ZNm0ax48fJ8zfkzG9u/HZzkwslobX5JSUlLBw4UJGjBhBYGAgERERXHfddRw8eLDBfOvXr2fChAmEh4fj6elJr169WLBgASUlJR2+XUKI9pPAcSJaa7YezeeSXqEOf/7m3nvvZenSpTz//PMsX76cvLw8rrjiCqqqqrhuSDQZRZX890Rhg2VOnjzJm2++yZQpU1i+fDmvv/46WVlZjBw5kvT09Pr5CgoKGDJkCC+//DKrV69mwYIFvPfee9xyyy323kwhRFtorZ391WWcyCvXcY98of+55Vi7ll+1apUG9LFj7Vu+tdLT07Wrq6t+77336qedOnVKu7u76zfffFOXV9fqfo9/rf+wYk+D5crKynRFRUWDafn5+drX11cvWrSoxe984403NKDz8/NttyFCiPZodn8tPRwDzZ07l4svvpgvv/yS/v374+Pjw7Rp0ygoKODw4cNMmDABX19fLr74Yvbs2cPWuutvkiO9uffee4mIiMDLy4vhw4fz7bffNli31ppFixYRHh6Ov78/t99+e5OHnKqqqnj44YeJjY3F09OTwYMH89VXX13Qdp2pZebMmfXToqOjGTNmDF9//TU+Hm5MGRDBVz9nUWOy1M/j6+uLt7d3g3WFhIQQFxdHTk5Oi98ZGmo9p1VTU3NBtQshOo4EjsFOnjzJE088wV/+8hfeeOMNtmzZwvz585k9ezazZ89m+fLlmEwmZs+eTcqRPLr5efDMow/w7rvv8uijj/Lpp58SGxvLtGnT2LRpU/16X3rpJRYvXsz8+fNZvnw53t7ePPzww+d8/w033MDSpUv505/+xKpVqxg+fDgzZsxg165d9fNYLBZMJlOLL7P5l6d6pqWlERMTg59fw/NM/fr1Iy0tDYAZyVEUV9ay/mBui+2Tm5vL4cOH6d+//zmfmc1mqqur2bVrF3/5y1+YOXMmERGO86htIUQjLXV/nOTltObMmaNdXV314cOH66c99NBDGmhwOOrLL7/UgB58/9t69pLlWimlly5dWv+52WzWAwYM0JMnT9Zaa20ymXRkZKS+8847G3zf5Zdf3uCQ2po1azSg161b12C+sWPH6htuuKFBnUCLr/Hjx9fPf8cdd+jBgwefs72PPvqojoyM1FprXWMy66GLv9W//fCnFtvotttu0yEhITovL++czxITE+u/f8qUKbq8vLzFdQkh7KLZ/bWbESEnfhEfH09CQkL9+969ewMwceLEc6ZlZWUx0r8IrTU33nhj/ecuLi7ceOONLFmyBID09HSysrK45pprGnzXzJkzWbNmTf37NWvWEBERwaWXXorJZKqfPmnSJJYuXVr/ftGiRdx9990tboe/v3+D900NatBa1093d3Vh2qBIlm1Pp6zahJ/nuT+Kr776Kh988AErVqyoP2R2thUrVlBcXMzPP//M4sWLufHGG+tH4gkhHI8EjsGCgoIavPfw8Dhn+plp2lRDgDbj5+eHj0/Dp3x2796diooKqqurOX36NADh4Q0fPd34fV5eHqdPn8bd3f2culxdf7kLdY8ePYiJiWlxO87eyQcHB1NUVHTOPEVFRQ2265rkaP6ZcoLVe09z/bCG6//888+55557ePbZZ7nuuuua/M4BAwYAMHr0aPr168f48eNZu3Ztg7AWQjgOOYfjRAK83biodxxlZWVUVFQ0+Cw7OxsfHx88PT3rz2M0PtHe+H1ISAjR0dFs3779nNfWrVvr55s3bx7u7u4tviZNmlQ/f1JSEunp6ZSXlzf4vrS0NJKSkurfD+0RRGyI9zkXgW7ZsoXZs2dz55138tBDD7WqbYYOHQrA0aNHWzW/EML+pIfjBHTdQ8sSu/szYsQIlFIsX76c22+/vf7z5cuXM2bMGABiY2OJiIhg5cqVXHnllfXr+eSTTxqsd9KkSbzwwgv4+fk1CILG2npIbfLkyQB8+umn3HrrrQBkZmayceNGXnnllfr5lFJcMziaV9YdJre0mjB/T/bt28f06dO58soreemll87bNmds3rwZgJ49e7Z6GSGEfUngOIGMokoAEiMC6NevHzfffDN33303JSUl9O7dmzfffJO0tDReffVVwHo47OGHH+bBBx+kW7dujB07lhUrVpCamtpgvVdccQVTpkzhiiuu4JFHHmHAgAGUlJSwa9cuqqqqePrppwHreab4+PhW1xsTE8P//M//cP/996O1JiwsjEWLFhEXF1cfQACLFy9m8eLFxDy4ki/2ZDKtjy9XXnklfn5+3Hvvvfz444/18wYEBNSPVLvtttvo27cvycnJ+Pj4sGPHDpYsWcKoUaOYMGFCu9pYCGEHLY0ocJKX05ozZ44eNmxYg2nvvvuuBnRpaWn9tJc+26QB/fr7y7TWWpeXl+u7775bh4eHaw8PDz1s2DD9zTffNFiPxWLRjz32mO7WrZv28/PTt9xyi/7www/PufCzqqpKP/HEEzohIUG7u7vr7t276ylTpugvvvjigratqqpK//73v9fdunXTPj4+eurUqfro0aMN5lm4cKEG9NQXN+gZ/9ik165d26pRcC+99JIeOnSoDggI0L6+vvqiiy7SixcvbtBmQgjDNLu/Vlo7/TPmnX4Dzuf3y3ax8VAu2x+9vFOOwHpjwxGe+iqNdQ9eRnw3X6PLEUJcmGZ3UjJowMFpbb1/2kgnuH9ae109OAqlYOUueTCbEJ2ZBI6DO1lQQVZxldM+jqA1IgO9GdkzhJW7MugEPW4hRDMkcBzcmfunjeoVYnAlHeva5GiO5pWzN0MeMSBEZyWB4+BSjuTTzc+DhLC2P//GmUy9KBIPVxd5MJsQnZgEjgPTWrPlSOc+f3NGoI87lyWGsWp3JmaLHFYTojOSwHFgB7PLyCmtZmzvbkaXYhfXDokmp7S6/jCiEKJzkcBxYBsPWW/dP6ZP1wiciUnh+Hm68dlOOawmRGckgePANh3Oo2c3X2KCfc4/cyfg5e7KlRdF8M3e01TVms+/gBDCqUjgOKhqk5ltRwsY00UOp51xTXIUpdUm1qa1/IRPIYTzkcBxUDtOFFFZa+4yh9POGJ3QjW5+njJaTYhOSALHQW06nIuri2JUQue94LMpri6KqwdHsjYtl+LKWqPLEULYkASOg9p0KI/BMYEEeJ37cLTO7trkaGrMFr76OcvoUoQQNiSB44CKKmrYk1HMmD5hRpdiiEExgSSE+fLJjlNGlyKEsCEJHAe0+XA+WsPYLnb+5gylFDOHxrD9eCEn8svPv4AQwilI4Dig79OyCfR2Z0hskNGlGOa6IdEoBZ/skMEDQnQWEjgOxmzRrDuQy4TEMNxcu+4/T1SQN6MTQvlk5ym5g7QQnUTX3aM5qF3phRSU1zCxX3ejSzHc9UNjSC+oZPvxQqNLEULYgASOg1mTmoOri2J83645YOBsUwZE4OPhyoqfZPCAEJ2BBI6D+SE1h+HxwQR6d73h0I35erox9aJIvvw5S251I0QnIIHjQNILKjiQXcrlcjit3vVDoymrNrF632mjSxFCXCAJHAfyQ939wyYmhRtcieO4pFcoUYFeMlpNiE5AAseBrEnNplc3X3p18qd7toWLi+K6odFsPJRLdkmV0eUIIS6ABI6DKKs2se1oAZP6Se+msZlDY7BoWCk39BTCqUngOIhNh3KpMVuYmCTnbxpLCPMjOTaIFT9lyDU5QjgxCRwHsSY1hwAvNy6ODza6FId0/bAYDmSXsi+zxOhShBDtJIHjACwWzdq0HMYnhuPehe8u0JKrB0Xi4erCCrmhpxBOS/ZuDmDXqSLyy2u4XM7fNCvIx4NJ/cL5fFcmtWaL0eUIIdpBAscB/CB3F2iV64fGkF9ew/oDuUaXIoRoBwkcB7AmNZthccEE+XgYXYpDG58YRqivhxxWE8JJSeAYLKOokrTTpXI4rRXcXV2YkRzF96k5FFXUGF2OEKKNJHAM9kNqNoAMh26l64fGUGO28MUeefy0EM5GAsdga1JziA/1ISHM1+hSnMKAqAASwnz5fHem0aUIIdpIAsdA5dUmUo7kM6lfd5RSRpfjFJRSXJMczY/HCsgoqjS6HCFEG0jgGGjT4TxqzBYmyc062+Sa5CgAVkkvRwinIoFjoB9Sc/D3dGN4zxCjS3EqcaG+JMcGsXKXBI4QzkQCxyAWi+b7tBzGJYbJ3QXa4ZrkKFKzSjiYXWp0KUKIVpI9nUF+zigmr6xahkO307RBkbgo+Fx6OUI4DQkcg3yfloOLgvF9JXDaI9zfi0t7d2PlbrmDtBDOQgLHID+kZTO0RzAhvnJ3gfa6Jjma9IJKdpwsMroUIUQrSOAY4HRxFXszSpgoh9MuyJQB3fF0c+FzeTCbEE5BAscAaw/kADBJ7i5wQfy93Lm8X3e+2JOFSe4gLYTDk8AxwPepOUQHedO3u5/RpTi9GclR5JfXsPlIvtGlCCHOQwLHzqpqzWw+nMekfuFydwEbuCwxDH8vN1bulMNqQjg6CRw7SzmaT2WtmYlydwGb8HRz5aqLIlm97zSVNWajyxFCtEACx85+SM3B292VS3qFGl1Kp3HNkCjKa8ysqbvzthDCMUng2JHWmrUHcri0dze83F2NLqfTGNkzlO4BnnKrGyEcnASOHR3Pr+BUYSXjE+VR0rbk6qK4elAU6w/Kg9mEcGQSOHa06VAuAGN7dzO4ks7n2iHR1Jo1X/182uhShBDNkMCxow2H8ogN8SYu1MfoUjqdMw9mWykXgQrhsCRw7MRktrD1SD5jeofJcOgOcObBbNuOFZApD2YTwiFJ4NjJ7lNFlFabGNdHDqd1FHkwmxCOTQLHTjYczMNFwegECZyOcubBbJ/JaDUhHJIEjp1sPpzHwJggAn3cjS6lU7tWHswmhMOSwLGDyhozu08VMUou9uxw0wZF4eqiZPCAEA5IAscOdqYXUmvWjOwZYnQpnV6Yv6f1wWy7MuXBbEI4GAkcO9h2tAAXBcPig40upUu4ZnAUpwor+e+JQqNLEUKcRQLHDn48VkD/qAACvOT8jT1ceVEEvh6u/Ht7utGlCCHOIoHTwapNZnacLGREvJy/sRdfTzeuHhzFlz9nUVZtMrocIUQdCZwO9vOpYqpNFkb2kvM39nTT8Fgqasx8IdfkCOEwJHA62LZjBQAMj5fAsachsUH0Dvdj2X/lsJoQjkICp4NtP15An3A/Qnw9jC6lS1FKMeviWHaeLOKQXJMjhEOQwOlAWmt2nixiWJyMTjPCdUOjcXNRLJPBA0I4BAmcDnQsr5ziylqG9AgyupQuqZufJ5MHdOc/P52iokYGDwhhNAmcDrTzZBEAQ3pID8cov760J8WVtXyyQ+48IITRJHA60K70Ivw93egd5md0KV3WxXHBDIoJ5J3Nx7BY5M4DQhhJAqcD7UwvZFBsIC4u8vwboyilmHdpT47mlrO+7omrQghjSOB0kMoaM6lZpQyJlcNpRrtqYCTdAzx5Z9Mxo0sRokuTwOkgP2cUY7ZoGTDgADzcXLh9VDwbD+WRdrrE6HKE6LIkcDrIzpPWG0cmxwYZW4gA4Fcje+Dr4crLa48YXYoQXZYETgfZebKIuFAfQv08jS5FAEE+Htw+Op4v9mRyOEcuBBXCCBI4HUBrzY6ThQyR3o1DuWNMT7zcXPn7D4eNLkWILkkCpwNkFVeRU1ot1984mFA/T24fFceq3ZkcyS0zuhwhuhwJnA6w51QxAINiAg2uRDR2x9heeLi58LL0coSwOwmcDrA/sxhXF0W/yACjSxGNhPl7cuvIOD7blSG9HCHsTAKnA+zNLKF3mB9e7q5GlyKacOdlCXi6ufK3NYeMLkWILkUCpwPszShmQJT0bhxVNz9P5oyOZ9WeTA7KowuEsBsJHBvLKbUOGBgQLedvHNn8cb3wcXflxTUHjS5FiC5DAsfG9mVar2S/SHo4Di3E14N5Y3ry1c+n2ZdZbHQ5QnQJEjg2ti/DuvPqL4Hj8O4Y0wt/LzdelHM5QtiFBI6N7cssIT7UB38vd6NLEecR6OPOHWN68d3+bPacKjK6HCE6PQkcG9ubWSznb5zIvDHxBHq783/fybkcITqaBI4NFVfUkl5QKSPUnIi/lzvzx/Vi7YFcfjpRaHQ5QnRqEjg2tC/Lev7moijp4TiTuaPjCfX1kBFrQnQwCRwb2pdhHaEmPRzn4uvpxp3jE9h4KI8fjxUYXY4QnZYEjg3tzSwmKtBLHknghG69JI4wf09e+PYAWmujyxGiU5LAsaG9GcX0l8NpTsnbw5XfXpbAtmMFpBzNN7ocITolCRwbqagxcTSvnIui5XCas7p5RA/C/T35+/dyJ2khOoIEjo2kZpWgtQwYcGZe7q787/gEUo7my7kcITqABI6N7D0zYEB6OE7tlhE96ObnwUvfy90HhLA1CRwb2ZdZTKivBxEBXkaXIi6At4cr88f1YtPhPH46Ib0cIWxJAsdG9maUMCA6EKWU0aWIC3TrJXGE+HrwkpzLEcKmJHBsoNpk5mB2qVx/00n4eLjxm7G9WH8wl13pRUaXI0SnIYFjA4eyyzBZtAwY6ERuGxVHkI87f5dzOULYjASODeyteySBDInuPPw83bhjTE++T8up//cVQlwYCRwb2JtZjL+nG7HBPkaXImzo9tHxBHi58cK3B4wuRYhOQQLHBvZlltA/KgAXFxkw0JkEeLnzuwm9WXsgl82H84wuRwinJ4FzgUxmC6lZJVwkz8DplOaMjicm2Ju/fJmK2SL3WBPiQkjgXKCjeeVU1VpkhFon5eXuysNXJpGaVcKKHaeMLkcIpyaBc4H2ZVpPKA+QEWqd1tWDIkmODeL51QcorzYZXY4QTksC5wLtzSjBy92FhDBfo0sRHUQpxePT+5FTWs3fZJi0EO0mgXOB9mYU0y8yADdXacrObFhcCDePiOXtTcfqe7VCiLaRveQFsFg0+zNL5ILPLuKRK5MI9nHnT5/8LAMIhGgHCZwLcLKggtJqk1zw2UUE+Xjw+PT+7D5VzOsbjhhdjhBORwLnAuyVAQNdzozBUVw1MIL/++4gqVklRpcjhFORwLkAezNKcHdV9O3ub3Qpwk6UUvzl2oEEenvw+2W7qDaZjS5JCKchgXMB9mUWkxjhj4ebNGNXEuLrwbPXDyTtdCkvrpFRa0K0luwp20lrzd6MYhkw0EVN6ted2cNjeX39Ef57XB7UJkRrSOC0U2ZxFYUVtQyQW9p0WY9N709UkDcL/rObihq5IFSI85HAaaf6RxLILW26LD9PN56/cTAnCyp45us0o8sRwuFJ4LTTvoxiXF0U/SIlcLqyS3qFMu/Snvwz5QSbDskdpYVoiQROO+3NLKF3mB9e7q5GlyIM9tCURBLCfHlo+W5KqmqNLkcIhyWB0057M4oZIBd8Cqx3lP7rTcnklFbz5Of7jS5HCIclgdMOOSVV5JRWywi1DvbZZ58xaNAgPD096dmzJ3/9619bnP/+++9HKcWDDz7YYHpaWhojR44kMDCQ2bNnU1ZW1uDzDRs2EB0dfc70pixduhSl1DnzDo4Nonf6V7w4dxzrDuQAcPz4cZRS9S9fX18SEhL41a9+xcaNG89Z99y5c7n44ovPW4MQzkoCpx32ZVqvMJeHrnWczZs3M3PmTEaMGMGqVauYN28ejzzyCC+++GKT8+/fv5933nmHgIBze51z586ld+/e/Pvf/2b//v089dRT9Z9ZLBbuv/9+nn76afz8/C6o5uE9Q3BV8NhnexuMWnv++edJSUnhq6++4vHHHyc/P59x48bx5JNPXtD3CeFs3IwuwBntOVWMUtBfRqh1mMWLFzNmzBjeeustACZPnkxhYSGLFy/mt7/9LR4eHg3mv/fee7nvvvt4//33G0wvKytj27ZtrFq1irCwMIqKinj++efrQ+ftt9/G3d2d22677YJrdnNxwcfTjVOFlfxtzSFu7u8NQGJiIpdccgkA48ePZ+7cuTzxxBMsWrSI8ePHc9lll13wdwvhDKSH0w670gvpE+6Hn6fkdUfZtWsXl19+eYNpZ0InJSWlwfTly5eTmprKH/7wh3PWU1NTA4C3t3Xn7+PjUz+tpKSExx9/nL/97W8opWxSt5uLYvbwWN7adIyD2c3fa23hwoVERUXx2muv2eR7hXAGEjhtpLVmV3oRg2OCjC6lU6uqqjqnF+Pp6QlAampq/bTKykoWLFjAM888g6/vuQ/BCwkJoWfPnvz973+noKCAN954o/48yZ///Gcuv/zy+t5HW5jNZkwmU4OXxWIB4I9T+xHs486Sbw40u7yrqysTJ05k69atbf5uIZyV/IreRicLKiisqCW5R5DRpXRqvXv3Zvv27Q2m/fjjjwAUFPxyK5mnn36ayMhIbr311mbX9fLLL3PjjTfypz/9iT59+vDyyy9z+PBh3n77bfbs2dOu+oKCgpqcHhoaSqCPO49P78/vXl/d4jpiYmLIzs5u1/cL4Yykh9NGu9KLAEiODTK0js7uzjvvZOXKlbz55psUFhayevVqXnjhBcDaOwA4duwYzz//PC+++GKLh8SmTp1KTk4OBw4cIDU1lR49evDAAw/w+9//npiYGF5++WV69OhBjx49eOWVV1pV34YNG9i+fXuD129+85v6z2cMjmJI3c9IWXXT1+ZoLQ9xE12L9HDaaFd6Ed7uriTKIwk61Lx589i9ezd33XUX8+fPx8fHh2effZZ77rmH7t27A/CHP/yBqVOnkpSURFFREWAddVZdXU1RURGBgYH1QeTj40Pfvn0BWLNmDbt372bZsmXs3r2bxx9/nC1btgAwatQoxowZw6BBg1qsb8iQIeeMavviiy/q/66U4v4r+vDFo/DZzkxmX3/uOjIyMuq3RYiuQHo4bbQrvYiB0YG4uUrTdSRXV1f+8Y9/kJuby549e8jOzq4/13LmzwMHDvDJJ58QHBxc/0pPT+cf//gHwcHBZGRknLNek8nE/fffz5IlS/D29mbdunVMnDiRpKQkkpKSmDRpEuvXr7fJNiSEWX8pWXcgh7TTDQcQmEwmfvjhB0aNGmWT7xLCGUgPpw1qTBb2ZZYwZ1Sc0aV0GWeCBOCVV15h9OjRJCUlAfDWW2+dcwHm7NmzGT9+PHfddRdhYWHnrO+1114jODiYWbNm1U+rqKio/3t5ebnND3X5eLiy6PN9fPSbS+p7XIsXLyYzM5M777zTpt8lhCOTwGmD1KwSakwWkmODjS6l09u6dSubNm0iOTmZkpISPvroI1avXs2mTZvq52nqqnwvLy9iY2ObvLalsLCQJ598ktWrfzmZP27cOB5++GHeeecdAH744QeeeeYZm27LJaE1rN24hSXup4lwKeHjjz/mm2++qb8OR4iuQgKnDeoHDMgItQ7n7u7OsmXLWLRoES4uLowdO5bNmzczcODAdq9z4cKFzJgxg6FDh9ZPGzJkCEuWLOHRRx8FrHcFGDx48AXXf7aP/v7/APjjxx7ExUQzevQoNmzYwNixY236PUI4OtUJRsrYbQN+v2wXmw7n8eOfJtnsQkHRNWw7ms+sN7bywBV9uXdSH6PLEaIjNbtzlDPfbbArvYjk2CAJG9FmI3uFMvWiCF5dd4TTxVVGlyOEISRwWqmgvIZjeeVy/Y1otz9O7YfZolmyWp4OKromCZxW2n7cenX7iJ4hBlcinFWPUB/mjenJJzsy2F13PlCIrkQCp5W2HyvAw82FQTHySALRfr+bkEA3Pw8Wf7Ff7jQguhwJnFb68XgBybFBeLrJI6VF+/l7ufPg5ER+OlHIF3uyjC5HCLuSwGmF8moT+zJLGBEvh9PEhbvx4lj6RQbwzNdpVNWajS5HCLuRwGmFHScLMVs0w+X8jbABVxfF49P7kVFUydubjhldjhB2I4HTCtuPFeCiYFic3GFA2MbohG5MGdCdl9ceJqdEhkmLrkECpxW2HStgQFSgPOFT2NSfrupHrdnCc6ubf1CbEJ2JBM55VNaY2ZlexEg5nCZsLC7Ul3mX9mT5jlP8fKrY6HKE6HASOOex9Vg+NSYLY/uee+dhIS7U7yb2JsTHgz/LMGnRBUjgnMfGg3l4urlID0d0iAAvdxZMTuTH4wV8vfe00eUI0aEkcM5jw6FcRvQMwctdrr8RHWPW8FiSIvx56qtUGSYtOjUJnBZkFlVyOKeM8XI4TXQg6zDp/pwqrOTltYeNLkeIDiOB04INB3MBGCeBIzrYpb27MXNoNK+sO8LeDBlAIDonCZwWbDiUS0SAF33C/YwuRXQBT0zvT7CPBw8v30Ot2WJ0OULYnAROM6pqzaw/kMtliWHy/BthF0E+Hvzl2ovYn1XCa+uOGF2OEDYngdOMTYfyKK8xM3VgpNGliC7kyosimD4okr99f4ifThQYXY4QNiWB04yv9mYR4OXGqF6hRpciupj/d91AooK8uftfOykorzG6HCFsRgKnCTUmC2v2Z3NF/wg83KSJhH0Fervzyq+Gkl9Wwz0f7aDGJOdzROcge9MmfJ+aTUmViemD5XCaMMZF0YE8PXMgmw/ns+A/u7FY5C4EwvnJ3Sib8PH2dCIDvRjXR4ZDC+NcPyyG3LJqnvk6DV8PV/7fdQNxdZEBLMJ5SeA0kllUyYZDudw9obf85xaG+99xvSivNvH3Hw6TV1bDX2cNJsDL3eiyhGgXOaTWyHspx1HATRfHGl2KECilWDA5kSdnDGDtgRymvbSRH9Ky5UafwilJ4JylsLyGD1JOMH1QFLEhPkaXI0S9OaPj+ff/XoK7iwvzlv6X297+kbVpOZjl3I5wIqoT/KZksw1Y8k0ar6w7wur7x5EY4W+r1QphMzUmC+9vPcGr646QV1ZNVKAXM5KjuSY5iqQIf7lIWTiCZn8IJXDq7M8sYcY/NnH14Cj+b1ayLVYpRIepMVn4bn82y39KZ8OhPMwWTd/uflyTHM2MwdJDF4aSwGlJSVUtN76aQn55Dd/9fhzBvh62qEsIu8gvq+arvadZuTOD/54oxEXBjMFR3DOpDwlhch9AYXcSOE05mlvG5sN5vLvlOOkFFbwzdzhjZSi0cGLpBRW8v/UE76ecoNpk5toh0TxwRV9igqXHI+xGAqcp76cc5/GV+4gL9eHpmQMZndDNlnUJYZi8smpeX3+E91JOgIbbR8Xxuwm9pfcu7EECpylFFTWU15iJDPDCRa65EZ1QZlElL645yPKfTuHr4cadlyXw60vj8fGQS/BEh5HAEaIrO5hdypJvDrAmNRt/TzemD45kyoAIhseH4Osp4SNsSgJHCAE/nSjkX9tO8tXPWVTWmnFzUfTs5kvvcD8SwvxIjPBnQFQA8aG+0usX7SWBI4T4RUWNiZ9OFLLtaAFpp0s5mlvGiYKK+gtJ/T3dGNc3jMkDujNlQARe7q4GVyyciASOEKJl1SYzh3PK2JdRwo6ThaxJzSGvrJpgH3dmj+jBHWN6EurnaXSZwvFJ4Agh2sZi0Ww5ks/7W4/z3f5svNxdmTM6nvlje8loN9ESCRwhRPsdzinjpe8PsWpPJj7urtw6Ko7/GdOTcH8vo0sTjkcCRwhx4Q5ml/L3Hw7z5Z5M3FxduOniGG4e0YP+kQHtuo+bxaIxa43ZotEaPNxcOu1jQfLLqjmUU8aJ/HKO51eQWVRJebWJGrMm0NudiABPkmODGd4z2NmDXAJHCGE7x/PKeX3DEZb/dIpas6ZnN18ujgumX2QAoX4eeLu7UllrpqLGTFFFLQXl1eSX1ZBXXlP/9/zymiYfn+3j4Yqvpxv+Xm7EBPsQH+pDjxAf4kN9ie/mQ0ywj0MPYtBaczSvnN3pRaSdLiU1q4S006XkllbXz+PuqogM9Mbfyw03F0VJlYnMokqqTRaUgpE9Q5g+KIppAyOd8fClBI4QwvYKymv4Zu9pvtt/mj2niskvr2lyPi93F0J9PQn18yDU14NQP09CfT3w9nDFVSlcXBQuSlFVa6a82kRZtYniylrSCys4kV9BaZWpfl1KQWSAFz1CfQjy9sDDzQVPNxdclEJj7SlpwFL3Fw83F6KDvIkJ8SY22Ife4X4E+dh2J34yv4L1B3NIOZrPj8cKyCuztoOHmwt9u/uRFBFAUoQ/fbv707ObL5GBXri5Nnw6TK3ZQmpWCWvTcvl8dwZHcstxd1VMTArn+qExXJYYjoebUzxRRgJHCNGxtNbkl9dQVFFLVa0ZL3dXfDxcCfJxv6A7G2itKaqo5Xh+OScLKjieV1F3WKqcsmoTNSYL1SYLFq1RKJSy7vHOHOKrNpnrA+CMcH9PEiP86RPuT9/ufvSN8KdPuB/+rXyaalFFDTtOFrLhYB7rD+ZyLK8cgOggb0b2DGFEzxCGxgXTq5vvOcHS2m3en1XCpzsy+GxXZv1owRmDo5iQFM6wuOBmaz3TXsfyyzmZX0F6QQUnCypIL6ygoLyG8mozFTXWdlPK2l4uZ/05c0g0j03v3+aazyKBI4TouqpqzWQUVXIiv5xD2WUczC7jYHYph3JKqar95bBedJA3vcJ8CfX1IMjHg0Bvd1xdFCazhdyyGnJLq+rOw1QA1p7bqF6hjO8bxvjEcOJDfWz+TCKT2cLGQ3ms2HGKb/dn1wWFtZcXHeyNr6cbbi4ulFbVUlxZS2ZRJSVn9QjBGrCxIT508/PA19MNXw83PNxc6nqD1l6hRWssWjMsLpjrhsRcSMkSOEII0ZjFojlVWMmB7FJrAGWXciy/gsLyGgrLayittu64lYJQXw+6+XkSH+rLoNhAkmOCGBoXbNfzSeXVJnaeLGLHyUKO55eTUVhJZa2ZWrPG38uNIG93wgOsNfbs5ktcqCHnvCRwhBCirc7ceUGB3Oqn9ZptqFYdXFRK9VdKfa+UqlBKZSqlFiulzhuZSqlApdS7SqlCpVSxUupDpVRoE/Ndo5T6WSlVpZTar5Sa1ehzD6XUc0qpjUqpSqVUkyEzd+7cumOSDV9paWmt2UwhhGjA1UXh6qLOGzYrV65k4MCBeHl50b9/f5YtW9aq9X/88ccMHToUPz8/oqOjuf3228nMzKz/PCsri4ceeojBgwfj5+dHbGwsc+bMaTCPMzlv4CilgoE1WHsS1wCLgQXAk61Y/zLgMuAOYC4wHPis0frHACuAtcBU4EvgI6XU5LNm86lbRwWwpaUvTEpKIiUlpcErPj6+FaUKIUTbbdq0ieuvv54JEybw9ddfM23aNG6++Wa+/fbbFpf7/PPPufnmmxk9ejQrV67k2WefZcOGDUyfPh2LxXpe6aeffuLTTz/l5ptvZtWqVTz33HNs27aN0aNHU1ZWZo/Nsy2tdYsv4I9AIRBw1rSHse78A1pYbhTWkBp31rQRddMuP2vaauCHRst+BWxqNO3M4b+7rWXXf1Zvzpw5etiwYVoIIexl8uTJesKECQ2mTZ06VV966aUtLjdr1iw9dOjQBtNWrlypAb1//36ttdaFhYW6tra2wTwHDhzQgF66dKkNqu8QzeZJaw6pTQVWa61Lzpr2MeANjD/Pctla6w1nhduPwLG6z1BKeQITgH83WvZjYJRSKvCsZeVcjRDCoVRXV7N27VpuuummBtNnz55NSkoKxcXFzS5bW1tLYGBgg2lBQUEAZ37JJigoCDe3hkPK+/bti4+PDzk5OTbYAvtqTeAkAQ1OgmitT2Lt4SS1Zbk6qWctlwC4NzFfal1tfVtRXwP79+8nICAAT09PxowZw/r169u6CiGEaJUjR45QW1tLUlLDXWG/fv2wWCwcPHiw2WXnzZvHxo0b+ec//0lJSQkHDx7kscceY8KECfTv3/x1MHv27KGioqLFeRxVawInGChqYnph3WcXstyZPxvPV9jo81YZMmQIL7zwAqtWreLDDz/EbDZzxRVX8OOPP7ZlNUII0SqFhdZd1ZmeyRnBwcENPm/KtGnTWLp0KfPnzycwMJDExETMZjOffPJJs8tYLBbuu+8++vTpw+TJk5udz1G19vLfpg5nqWamt2e5xu9VM9NbdN999zV4P23aNPr3789TTz3FZ5991pZVCSFEqzW+2PPMIbGWLgJdu3Ytd955J/fddx9Tp04lOzubRYsWcd1117FmzRpcXc8dCPzHP/6RlJQU1q9fj7t76+6K4EhaEziFQFAT0wNpugdz9nJhTUwPOmu5wrOmNZ6H86z/vLy9vbnqqqtYtWrVhaxGCCGadKYnU1RU1GD6mfeNez5nW7BgATNmzODZZ5+tn5acnExSUhIrV65k5syZDeZ/5ZVXeO655/joo48YOXKkTeq3t9YcUkuj0bkapVQs4EvT52iaXa7O2ed2jgC1TcyXBFiA5g+AtoGtbzUhhBAACQkJuLu7n3OtX1paGi4uLvTt2/xp6LS0NJKTkxtMS0xMxNvbmyNHjjSYvmLFCu655x6WLFnCrFkNLlN0Kq0JnK+BKUop/7OmzQIqgZbOyH8NRNRdZwOAUupioFfdZ2itq7Fef3Njo2VnASla6+aHeLRCZWUlX3/9NcOGDbuQ1QghRJM8PT2ZMGEC//nPfxpMX7ZsGaNGjTpnFNrZ4uLi2LFjR4NpqampVFZWNrh2cN26dfzqV7/i7rvv5sEHH7Rp/XbX0pjpuuOQwUAW8B1wOTAfKAP+0mi+w8DbjaZ9AxwFZgLXAgeAjY3mGQOYgBexXiS6BGvvZnKj+aYCNwBvYT23cwNww/Hjx7XWWhcVFekxY8bo1157Ta9Zs0Z//PHHeuTIkdrDw0Nv3769A4ecCyG6so0bN2pXV1d933336bVr1+qHHnpIK6X06tWr6+c5fvy4dnV11e+99179tBdffFErpfQDDzygv/vuO/3BBx/ovn376vj4eF1WVqa11nr//v06MDBQDx48WG/evFmnpKTUvw4fPmz3bW2l5vOkpQ/1Lzv7/sAPWHs1WcCfAddG8xwHljaaFgS8i/VcTAnwL6BbE+u/FtgLVGM93Da7iXmO1wVNg9e7776rtda6srJSX3fddTomJkZ7eHjogIAAPWXKFJ2SktLhrSuE6No+/fRTPWDAAO3h4aETExP1Rx991ODzY8eONdhfaa21xWLRr7zyih44cKD28fHRUVFR+qabbtJHjhypn+fdd989Z5935jVnzhw7bV2bNZslcvNOIYQQtnRhN+8UQgghLpQEjhBCCLuQwBFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsAsJHCGEEHYhgSOEEMIuJHCEEELYRWufh9NpLVq0iCeffNLoMoQQDm7hwoUsWrTI6DKcmvRwhBBC2IUEjhBCCLuQm3cKIYSwJbl5pxBCCGNJ4AghhLALCRwhhBB2IYEjhBDCLiRwhBBC2IUEjhBCCLuQwBFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsAsJHCGEEHYhgSOEEMIuJHCEEELYhQSOEEIIu3D6xxM8+eST3wDdLmAVUUCmjcpxVtIG0gYgbQDSBmdcSDvkLVy48MomP9Fad+nXokWLtNE1GP2SNpA2kDaQNrBHO8ghNSGEEHYhgQNPGl2AA5A2kDYAaQOQNjijQ9rB6c/hCCGEcA7SwxFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsItOHThKqd8qpY4ppaqUUj8ppcaeZ/6BSqn1SqlKpVSGUuoJpZSyV70dpS3toJS6TCm1UimVpZSqUErtUUrNs2e9HaGtPwtnLddHKVWqlCrr6Bo7Wjv+Pyil1P1KqTSlVHXdz8Qz9qq3I7SjDaYopVLqfgby6v5v9LVXvbamlBqnlPq8bv+mlVJzW7GMzfaLnTZwlFKzgL8BTwFDgC3A10qpHs3MHwB8B2QDw4F7gYeAB+xScAdpazsAo4GfgRuAi4BXgTeUUrfYodwO0Y42OLOcB/AxsKHDi+xg7WyDF4DfAo8A/YCrcOK2aMc+oSewEthYN//lgDfwlV0K7hh+wF7gPqDyfDPbfL9o9BWtHfUCtgFvNpp2CHi6mfnvAkoA77OmPQZkUDd83BlfbW2HZtbxb2CF0dti7zYA/g94F5gLlBm9HfZsAyARqAX6GV27gW1wA2AGXM+aNgHQQDejt8cG7VEGzD3PPDbdL3bKHk7db6bDgG8bffQt1t/gmzIK2Ki1Pjv1V2O9p1C8rWu0h3a2Q1MCgEJb1WVP7W0DpdQ0YDrW3+icWjvb4BrgKHClUuqoUuq4Uuo9pVR4B5baYdrZBv/FGrp3KKVclVL+wBxgu9Y6r8OKdSw23S92ysDBejNPV6zdwLNlAxHNLBPRzPxnPnNG7WmHBpRS04FJwBu2Lc1u2twGSqlI4E3gNq11aceWZxft+TnoBcQBs7H28G4DkoBVSiln3G+0uQ201seBK7BedV8NFAMDsf4i0lXYdL/ojD84bdH4NgqqiWnnm7+p6c6mre1gnUmpS4F/AfdqrX/siMLsqC1t8AHwqtZ6a8eWZHdtaQMXwBNr6G7QWm/EGjojsB7Ld1atbgOlVATwNvBPrNt8GVAK/NtJQ7e9bLZf7KyNlof12GvjBA7n3LQ+43Qz89PCMo6uPe0AgFJqDPA18ITW+tWOKc8u2tMGE4GFSimTUsqEdafjW/d+fseV2mHa0wZZgElrffCsaYcAE9DiYAsH1Z42+B1QrrV+WGu9U2u9AbgVGE/bDkk7M5vuFztl4Gita4CfsHaHz3YF1pEpTUkBxiqlvBrNnwkct3WN9tDOdkApNQ5r2DyptX6xwwq0g3a2wUAg+azXE1hH9CQD/7F9lR2rnW2wGXBTSiWcNa0X4AacsHmRHaydbeCDNaTOduZ9p9x3NsG2+0WjR0p04AiMWUANcAfWIZ1/wzoqI67u86eB78+aPxBrmn+MdTjwTKyjMxYYvS12bofLgHLgOay/2Zx5hRm9LfZqgyaWn4vzj1Jr68+BC9Yd9HqsQ4KH1P19K+Bi9PbYqQ0mAhZgIdAHGAp8A5wEfI3enna2gR+//CJVgfWXqWSgRzNtYNP9ouEN0MGN+1usKVxd959n3FmfLQWON5p/INbrDKqwHlJYiBMPiW5PO9S91028jtu7biN/Fhot6/SB0542ACKx9uhKgRzgQ6C70dth5zaYDeyoC6ZcYBXQ3+jtuIDtv6yZ/99LW2gDm+0X5fEEQggh7KKrHIcUQghhMAkcIYQQdiGBI4QQwi4kcIQQQtiFBI4QQgi7kMARQghhFxI4Qggh7EICRwghhF38f7jSqv9raiZlAAAAAElFTkSuQmCC\n",
"text/plain": [
"