From 672f82c548d0d1c1b20aff8e7d33716b9bdd8695 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sat, 21 Apr 2018 16:16:29 -0400 Subject: [PATCH 1/4] Problem 5.9.3 (in part) --- BDA 5.9.3.ipynb | 309 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 309 insertions(+) create mode 100644 BDA 5.9.3.ipynb diff --git a/BDA 5.9.3.ipynb b/BDA 5.9.3.ipynb new file mode 100644 index 0000000..d155ec7 --- /dev/null +++ b/BDA 5.9.3.ipynb @@ -0,0 +1,309 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#Problem 5.9.3\n", + "\n", + "Hierarchical models and multiple comparisons:\n", + "\n", + "1. Reproduce the computations in Section 5.5 for the educational testing example. Use the posterior simulations to estimate:\n", + " * for each school $j$, the probability that it's coaching program is the best of the eight;\n", + " * for each pair of schools $(j,k)$ the probability that the $j$th school is better than the $k$th \n", + "2. Reproduce (1) but for the simpler model where the population variance $\\tau$ is $\\infty$ so the eight schools are independent.\n", + "3. Discuss the differences between 1 and 2.\n", + "4. What happens when $\\tau=0$?" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " effect se\n", + "school \n", + "A 28 15\n", + "B 8 10\n", + "C -3 16\n", + "D 7 11\n", + "E -1 9\n", + "F 1 11\n", + "G 18 10\n", + "H 12 18\n" + ] + } + ], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import pystan\n", + "import matplotlib.pyplot as plt\n", + "\n", + "schools=['A','B','C','D','E','F','G','H']\n", + "effects=[28,8,-3,7,-1,1,18,12]\n", + "se=[15,10,16,11,9,11,10,18]\n", + "p55=pd.DataFrame(index=schools)\n", + "p55.index.name='school'\n", + "p55['effect']=np.array(effects)\n", + "p55['se']=np.array(se)\n", + "print(p55)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "pooled mean= 7.685616724956035\n", + "pooled variance= 16.580525632563663\n" + ] + } + ], + "source": [ + "print('pooled mean=',sum(p55['effect']*1/p55['se']**2)/(sum(1/p55['se']**2)))\n", + "print('pooled variance=',(1/sum(1/p55['se']**2)))" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1c0a010b4129370aa04f0b4b9f729b4d NOW.\n" + ] + } + ], + "source": [ + "stan_code='''\n", + "data {\n", + " real means[8];\n", + " real se[8];\n", + "\n", + "}\n", + "\n", + "parameters {\n", + " real theta[8] ; \n", + " real mu ; \n", + " real tau ; \n", + "}\n", + "\n", + "model {\n", + " \n", + " theta~normal(mu,tau) ; \n", + " means~normal(theta,se) ; \n", + " \n", + "}\n", + "\n", + "generated quantities {\n", + " real results[8] ; \n", + " \n", + " \n", + " for(i in 1:8) {\n", + " results[i]=normal_rng(theta[i],tau);\n", + " }\n", + "}\n", + "'''\n", + "sm=pystan.StanModel(model_code=stan_code)" + ] + }, + { + "cell_type": "code", + "execution_count": 264, + "metadata": {}, + "outputs": [], + "source": [ + "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)" + ] + }, + { + "cell_type": "code", + "execution_count": 265, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inference for Stan model: anon_model_1c0a010b4129370aa04f0b4b9f729b4d.\n", + "4 chains, each with iter=5000; warmup=2500; thin=1; \n", + "post-warmup draws per chain=2500, total post-warmup draws=10000.\n", + "\n", + " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", + "theta[0] 12.02 0.16 8.34 -2.36 6.78 11.04 16.71 31.15 2751 1.0\n", + "theta[1] 8.09 0.16 6.45 -4.95 4.02 8.05 12.37 20.76 1639 1.0\n", + "theta[2] 6.35 0.15 8.02 -11.67 2.0 6.96 11.29 20.8 2679 1.0\n", + "theta[3] 7.93 0.17 6.82 -6.17 3.85 7.98 12.22 21.32 1702 1.0\n", + "theta[4] 5.01 0.21 6.48 -8.91 1.17 5.27 9.39 16.85 934 1.0\n", + "theta[5] 6.19 0.17 6.8 -8.48 2.27 6.49 10.63 18.67 1624 1.0\n", + "theta[6] 11.18 0.14 6.73 -1.17 6.74 10.84 15.34 25.85 2277 1.0\n", + "theta[7] 8.77 0.16 8.16 -7.36 3.94 8.58 13.42 25.97 2654 1.0\n", + "mu 8.15 0.15 5.19 -2.15 4.84 8.09 11.43 18.47 1215 1.0\n", + "tau 7.11 0.2 5.17 1.37 3.35 5.84 9.47 20.71 700 1.01\n", + "results[0] 12.0 0.18 12.34 -9.06 5.18 10.65 17.58 40.77 4545 1.0\n", + "results[1] 8.13 0.16 10.95 -13.75 2.61 8.11 13.74 30.22 4412 1.0\n", + "results[2] 6.46 0.16 11.79 -19.82 0.76 7.19 12.86 28.9 5157 1.0\n", + "results[3] 7.82 0.15 11.12 -15.81 2.21 7.88 13.72 29.86 5431 1.0\n", + "results[4] 5.04 0.19 10.95 -20.09 -0.12 5.9 11.33 24.39 3484 1.0\n", + "results[5] 6.08 0.17 11.29 -18.73 0.73 6.67 12.38 27.52 4227 1.0\n", + "results[6] 10.98 0.16 11.13 -9.55 4.98 10.24 16.19 36.05 4864 1.0\n", + "results[7] 8.96 0.17 11.99 -15.44 2.81 8.58 14.76 34.5 4746 1.0\n", + "lp__ -18.29 0.26 4.93 -27.55 -21.75 -18.51 -14.86 -8.66 367 1.02\n", + "\n", + "Samples were drawn using NUTS at Sat Apr 21 15:56:14 2018.\n", + "For each parameter, n_eff is a crude measure of effective sample size,\n", + "and Rhat is the potential scale reduction factor on split chains (at \n", + "convergence, Rhat=1).\n" + ] + } + ], + "source": [ + "print(answers)" + ] + }, + { + "cell_type": "code", + "execution_count": 266, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF3BJREFUeJzt3X+0XWV95/H3pwGiFQckhFECGihojVPL1Ii2amuhtaGosQ6UoLa0wxRnjUzrUpeC0yrShULXjNRRZs2wBMWfQPFHMyWKjNTWWqUJiGKKlEijhFgJBFCqCIHv/LF39Hhyb+659x5yufd5v9a665797Gfv/TznnPs5z332OfukqpAkteGn5roBkqQ9x9CXpIYY+pLUEENfkhpi6EtSQwx9SWqIoa+xSnJfksP34PFekOTmMe7vU0lO6W//XpK/G+O+X5nkM+Pa3zSO+7wkt/SPzcv29PH16GLoz1NJNif5Qf+H/J0k70uy7yz2tzxJJdlrNu2qqn2r6tbZ7GOgTWcleTDJ9/qff0ryniRPGjje56vqaSPu60NT1auq46rqkjG0fZf7s6o+XFUvmu2+Z+Bs4D39Y/PJ4ZXjfi5NJsnqJDck+W6SO5N8Nsnyft2Uj7XGw9Cf315SVfsCvwA8G/jjuWrIbF8sdrP9ZVX1eOAA4LeAJwLXjTsM0lmofw9PATZOUecRfS4lOQL4APB6YD/gMOB/AQ8PVNsjj3XrFuqTvClVdTvwKeDfASQ5OMnaJNuTbEryBzvrJjk6yYZ+tPWdJO/sV/1t//uefsT3i339/5jkpiR3J7kqyVMG9lVJXpPkFuCWgbIj+tv7JflAkm1Jvpnkj3cGaz918oUk5yfZDpw1RR8frKqNwEnANrrwIMkLk2wZaNObktzejxZvTnJsklXAm4GT+r59pa/7uSTnJPkC8H3g8L7sPw0cOkneneTeJF9PcuzAis1Jfm1gefC/iV3uz+HpoiS/lGR9v+/1SX5pYN3nkvxpfx99L8lnkhw42f2T5A/6x3p7/9gf3Jd/Azgc+L99OxZPcT9P57l0VpIrklzWt/H6JD8/ya6PAv65qj5bne9V1ceq6lsTtGHCx1rjYegvAEkOBX4T+HJf9FFgC3AwcALw9oGwehfwrqr6N8DPAJf35b/c/96/nwb4Yrr53zcDLweWAp/v9z3oZcBzgBUTNO3ddKO6w4FfAX4X+P2B9c8BbgUOAs4Zpa9V9RDwl8ALhtcleRpwOvDsfsT4G8Dmqvo08Ha6keS+VTUYTL8DnAY8HvjmBIfc2cYDgbcCH09ywAhN3eX+HGrrAcCVwP8ElgDvBK5MsmSg2ivo7q+DgH2AN0x0oCTHAO8Afht4Ut+PSwGq6meAb9GP5Kvqh7tr9DSfSwCrgb+gG51/BPhkkr0n2PX1wM/2L/K/mhGmj3b3WGvmDP357ZNJ7gH+Dvgbuj/IQ4HnA2+qqvur6gbgvXThBvAgcESSA6vqvqr60m72/2rgHVV1U1XtoAvOowZH+/367VX1g8ENkyyiG6md2Y/qNgP/Y6AdAFur6t1VtWN4+ylspQuZYQ8Bi4EVSfauqs1V9Y0p9vX+qtrYt+HBCdbfAfx5P/q8DLgZOH4abZ3M8cAtVfXB/tgfBb4OvGSgzvuq6p/6++ZyutHyRF4JXFxV1/ehfibwi+nny0c0k+cSwHVVdUV/370TeAzw3OGd9+d5Xggs6/tyZ5L3jxD+kz3WmiFDf357WVXtX1VPqar/0ofDwcD2qvreQL1v0v2xAZwKPBX4ej+l8OLd7P8pwLuS3NMHwnYgA/sCuG2SbQ+kG50Ojp4H27G7baeyrG/LT6iqTcBr6aaK7khy6c5pjt2Yqg23109elfCbdPfxbB3Mrv9ZDN8//zJw+/vAZAH5E/uqqvuAu4b2NZWZPJdg4P6rqof58X8Fu6iqL1XVb1fVUrrR+y8D/22Kdk34WGvmDP2FZytwQJLHD5Q9GbgdoKpuqaqT6aYMzgOuSPI4YKLLrd4GvLoPg50/j62qvx+oM9llWu+k+69i8L+CH7Vjim0n1Z8TeAndVNMuquojVfX8/rhF18fdHWuqNixLkoHlJ9PdxwD/Cvz0wLonTmO/W/nJ+2bnvm+foO5UfmJf/eO5ZIb7Gt7vpM+l3qEDx/0p4BB+fP9MqqrWAx+nP3cwkakea82Mob/AVNVtwN8D70jymCTPpBvdfxggyauSLO1HZff0mz1Ed8LsYbr5953+N3Bmkmf02+6X5MQR2/EQ3b/x5yR5fD8l9DpgyrdNTiTJ3kmeTjfH/ES6qYThOk9Lckx/svJ+4Ad93wC+AyzP9N+hcxDwh/3xTwSeDqzr190ArOnXraSb895povtz0DrgqUlekWSvJCfRnRf5q2m2D7q59N9PclTf97cD1/ZTajM21XOp96wkL0/37qvXAj8EdpkyTPL8/mTzQf3yzwIvnaTulI+1Zs7QX5hOBpbTjbg+Aby1qq7u160CNia5j+6k7pp+vvb7dCdTv9BP5zy3qj5BN1K+NMl3ga8Bx02jHf+VbjR8K91c8UeAi6fZl5P6tt4DrKWbtnhWVU00mlwMnEv3X8a/0AX2m/t1f9H/vivJ9dM4/rXAkf0+zwFOqKq7+nV/Qncy/G7gbXT9A2Ci+3Nwp/0+Xkz3zpS7gDcCL66qO6fRtp37+mzflo8B3+7btGa6+5nE7p5L0J1oPYnuPvgd4OWTnBu5hy7kb+wfz0/3+/uzgTrTeaw1Q/FLVCTNRJKzgCOq6lVz3RaNzpG+JDXE0Jekhji9I0kNcaQvSQ2Z1UWyHgkHHnhgLV++fK6bIUnzynXXXXdn/8G33XrUhf7y5cvZsGHDXDdDkuaVJBNdO2oXTu9IUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDHnWfyH2kLD/jygnLN587ju+4lqT5wZG+JDVkwY30JxvRS5Ic6UtSUwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGLLi3bE6XH9qS1BJH+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaMlLoJ1mV5OYkm5KcMcH6X05yfZIdSU4YWndKklv6n1PG1XBJ0vRNGfpJFgEXAMcBK4CTk6wYqvYt4PeAjwxtewDwVuA5wNHAW5M8YfbNliTNxCgj/aOBTVV1a1U9AFwKrB6sUFWbq+qrwMND2/4GcHVVba+qu4GrgVVjaLckaQZGCf1lwG0Dy1v6slGMtG2S05JsSLJh27ZtI+5akjRdo4R+JiirEfc/0rZVdWFVrayqlUuXLh1x15Kk6Rol9LcAhw4sHwJsHXH/s9lWkjRmo4T+euDIJIcl2QdYA6wdcf9XAS9K8oT+BO6L+jJJ0hyYMvSragdwOl1Y3wRcXlUbk5yd5KUASZ6dZAtwIvB/kmzst90O/CndC8d64Oy+TJI0B0b6usSqWgesGyp7y8Dt9XRTNxNtezFw8SzaKEkaEz+RK0kNMfQlqSGGviQ1ZKQ5/RYtP+PKCcs3n3v8Hm6JJI2PI31JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ/xw1jT5oS1J85kjfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUkJFCP8mqJDcn2ZTkjAnWL05yWb/+2iTL+/K9k1yS5MYkNyU5c7zNlyRNx5Shn2QRcAFwHLACODnJiqFqpwJ3V9URwPnAeX35icDiqvo54FnAq3e+IEiS9rxRRvpHA5uq6taqegC4FFg9VGc1cEl/+wrg2CQBCnhckr2AxwIPAN8dS8slSdM2SugvA24bWN7Sl01Yp6p2APcCS+heAP4V+DbwLeC/V9X2WbZZkjRDo4R+JiirEescDTwEHAwcBrw+yeG7HCA5LcmGJBu2bds2QpMkSTMxSuhvAQ4dWD4E2DpZnX4qZz9gO/AK4NNV9WBV3QF8AVg5fICqurCqVlbVyqVLl06/F5KkkYwS+uuBI5MclmQfYA2wdqjOWuCU/vYJwDVVVXRTOsek8zjgucDXx9N0SdJ0TRn6/Rz96cBVwE3A5VW1McnZSV7aV7sIWJJkE/A6YOfbOi8A9gW+Rvfi8b6q+uqY+yBJGtFeo1SqqnXAuqGytwzcvp/u7ZnD2903UbkkaW74iVxJaoihL0kNMfQlqSGGviQ1ZKQTuZra8jOunLB887nH7+GWSNLkHOlLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhfjH6I8wvTJf0aOJIX5IaYuhLUkNGCv0kq5LcnGRTkjMmWL84yWX9+muTLB9Y98wkX0yyMcmNSR4zvuZLkqZjytBPsgi4ADgOWAGcnGTFULVTgbur6gjgfOC8ftu9gA8B/7mqngG8EHhwbK2XJE3LKCP9o4FNVXVrVT0AXAqsHqqzGrikv30FcGySAC8CvlpVXwGoqruq6qHxNF2SNF2jhP4y4LaB5S192YR1qmoHcC+wBHgqUEmuSnJ9kjdOdIAkpyXZkGTDtm3bptsHSdKIRgn9TFBWI9bZC3g+8Mr+928lOXaXilUXVtXKqlq5dOnSEZokSZqJUUJ/C3DowPIhwNbJ6vTz+PsB2/vyv6mqO6vq+8A64Bdm22hJ0syMEvrrgSOTHJZkH2ANsHaozlrglP72CcA1VVXAVcAzk/x0/2LwK8A/jqfpkqTpmvITuVW1I8npdAG+CLi4qjYmORvYUFVrgYuADybZRDfCX9Nve3eSd9K9cBSwrqom/oiqJOkRN9JlGKpqHd3UzGDZWwZu3w+cOMm2H6J726YkaY75iVxJaoihL0kNMfQlqSFeWnmOeMllSXPBkb4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGuIF1x5lvBCbpEeSI31JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhI4V+klVJbk6yKckZE6xfnOSyfv21SZYPrX9ykvuSvGE8zZYkzcSUoZ9kEXABcBywAjg5yYqhaqcCd1fVEcD5wHlD688HPjX75kqSZmOUkf7RwKaqurWqHgAuBVYP1VkNXNLfvgI4NkkAkrwMuBXYOJ4mS5JmapTQXwbcNrC8pS+bsE5V7QDuBZYkeRzwJuBtuztAktOSbEiyYdu2baO2XZI0TaOEfiYoqxHrvA04v6ru290BqurCqlpZVSuXLl06QpMkSTMxypeobAEOHVg+BNg6SZ0tSfYC9gO2A88BTkjyZ8D+wMNJ7q+q98y65ZKkaRsl9NcDRyY5DLgdWAO8YqjOWuAU4IvACcA1VVXAC3ZWSHIWcJ+BL0lzZ8rQr6odSU4HrgIWARdX1cYkZwMbqmotcBHwwSSb6Eb4ax7JRkuSZmak78itqnXAuqGytwzcvh84cYp9nDWD9kmSxshP5EpSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1JCR3qevubf8jCsnLN987vF7uCWS5jNH+pLUEENfkhpi6EtSQwx9SWqIJ3LnuclO8IIneSXtypG+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhriBdcWML9tS9IwR/qS1BBDX5IaMlLoJ1mV5OYkm5KcMcH6xUku69dfm2R5X/7rSa5LcmP/+5jxNl+SNB1Thn6SRcAFwHHACuDkJCuGqp0K3F1VRwDnA+f15XcCL6mqnwNOAT44roZLkqZvlJH+0cCmqrq1qh4ALgVWD9VZDVzS374CODZJqurLVbW1L98IPCbJ4nE0XJI0faOE/jLgtoHlLX3ZhHWqagdwL7BkqM5/AL5cVT8cPkCS05JsSLJh27Zto7ZdkjRNo4R+Jiir6dRJ8gy6KZ9XT3SAqrqwqlZW1cqlS5eO0CRJ0kyMEvpbgEMHlg8Btk5WJ8lewH7A9n75EOATwO9W1Tdm22BJ0syNEvrrgSOTHJZkH2ANsHaozlq6E7UAJwDXVFUl2R+4Ejizqr4wrkZLkmZmytDv5+hPB64CbgIur6qNSc5O8tK+2kXAkiSbgNcBO9/WeTpwBPAnSW7ofw4aey8kSSNJ1fD0/NxauXJlbdiwYcbbT3bpAU3NyzNI81eS66pq5VT1/ESuJDXE0Jekhhj6ktQQQ1+SGmLoS1JD/BIV/YhfuiItfI70Jakhhr4kNcTQl6SGOKevKTnXLy0cjvQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqI79PXjE33W8p8X7809xzpS1JDDH1JaoihL0kNMfQlqSGeyNWc84Ju0p5j6OtRyxcDafwMfe0x032Lp6Txc05fkhriSF/zjtM+0swZ+lowZjJ95AuFWjNS6CdZBbwLWAS8t6rOHVq/GPgA8CzgLuCkqtrcrzsTOBV4CPjDqrpqbK2XZslLSag1U4Z+kkXABcCvA1uA9UnWVtU/DlQ7Fbi7qo5IsgY4DzgpyQpgDfAM4GDg/yV5alU9NO6OSHvCdKeWxjUV5ZSWxmWUkf7RwKaquhUgyaXAamAw9FcDZ/W3rwDekyR9+aVV9UPgn5Ns6vf3xfE0X3p0mO5/DON6J9O4XoR2t40WllFCfxlw28DyFuA5k9Wpqh1J7gWW9OVfGtp22fABkpwGnNYv3pfk5inadCBw5whtn8/s4/w3Z/3LeXtmGxb+Ywjzp49PGaXSKKGfCcpqxDqjbEtVXQhcOEJbuoMlG6pq5aj15yP7OP8t9P6BfZyPRnmf/hbg0IHlQ4Ctk9VJshewH7B9xG0lSXvIKKG/HjgyyWFJ9qE7Mbt2qM5a4JT+9gnANVVVffmaJIuTHAYcCfzDeJouSZquKad3+jn604Gr6N6yeXFVbUxyNrChqtYCFwEf7E/Ubqd7YaCvdzndSd8dwGvG9M6dkaeC5jH7OP8t9P6BfZx30g3IJUkt8No7ktQQQ1+SGjLvQj/JqiQ3J9mU5Iy5bs84JLk4yR1JvjZQdkCSq5Pc0v9+wly2cTaSHJrkr5PclGRjkj/qyxdSHx+T5B+SfKXv49v68sOSXNv38bL+zRDzVpJFSb6c5K/65QXVP4Akm5PcmOSGJBv6sgXzXJ1XoT9wSYjjgBXAyf2lHua79wOrhsrOAD5bVUcCn+2X56sdwOur6unAc4HX9I/bQurjD4FjqurngaOAVUmeS3dJkvP7Pt5Nd8mS+eyPgJsGlhda/3b61ao6auD9+QvmuTqvQp+BS0JU1QPAzktCzGtV9bd073oatBq4pL99CfCyPdqoMaqqb1fV9f3t79GFxjIWVh+rqu7rF/fufwo4hu7SJDDP+5jkEOB44L39clhA/ZvCgnmuzrfQn+iSELtc1mGB+LdV9W3oQhM4aI7bMxZJlgP/HriWBdbHfurjBuAO4GrgG8A9VbWjrzLfn69/DrwReLhfXsLC6t9OBXwmyXX9JWJgAT1X59v19Ee6rIMenZLsC3wMeG1VfbcbKC4c/WdQjkqyP/AJ4OkTVduzrRqPJC8G7qiq65K8cGfxBFXnZf+GPK+qtiY5CLg6ydfnukHjNN9G+i1d1uE7SZ4E0P++Y47bMytJ9qYL/A9X1cf74gXVx52q6h7gc3TnL/bvL00C8/v5+jzgpUk2002rHkM38l8o/fuRqtra/76D7sX7aBbQc3W+hf4ol4RYKAYvbXEK8Jdz2JZZ6ed+LwJuqqp3DqxaSH1c2o/wSfJY4Nfozl38Nd2lSWAe97GqzqyqQ6pqOd3f3TVV9UoWSP92SvK4JI/feRt4EfA1FtJzdb59IjfJb9KNMHZeEuKcOW7SrCX5KPBCuku4fgd4K/BJ4HLgycC3gBOravhk77yQ5PnA54Eb+fF88Jvp5vUXSh+fSXeCbxHdYOryqjo7yeF0I+MDgC8Dr+q/X2Le6qd33lBVL15o/ev784l+cS/gI1V1TpIlLJTn6nwLfUnSzM236R1J0iwY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakh/x90grO/jaMYMQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig,ax=plt.subplots(1)\n", + "j=ax.hist(answers['tau'],bins=50,density=True)\n", + "j=ax.set_title('Posterior Distribution of Pop SD')" + ] + }, + { + "cell_type": "code", + "execution_count": 267, + "metadata": {}, + "outputs": [], + "source": [ + "predictions=answers.extract()['results']\n", + "def best_school(x,i):\n", + " return x[i]>=max(x)\n", + "def better_school(x,i,j):\n", + " return x[i]>=x[j]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 333, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Chance is empirical probability that given school is best\n", + " effect se Chance\n", + "school \n", + "A 28 15 0.2209\n", + "B 8 10 0.1119\n", + "C -3 16 0.0972\n", + "D 7 11 0.1116\n", + "E -1 9 0.0629\n", + "F 1 11 0.0829\n", + "G 18 10 0.1750\n", + "H 12 18 0.1376\n" + ] + } + ], + "source": [ + "print('Chance is empirical probability that given school is best')\n", + "p55['Chance']=[sum([best_school(x,i) for x in predictions])/len(predictions) for i in range(8)]\n", + "print(p55)" + ] + }, + { + "cell_type": "code", + "execution_count": 332, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Empirical Probability that school in given row is\n", + " as good as or better than corresponding column\n", + " \n", + " A B C D E F G H\n", + "A 1.00 0.61 0.64 0.60 0.68 0.66 0.53 0.58 \n", + "B 0.39 1.00 0.53 0.50 0.59 0.56 0.41 0.47 \n", + "C 0.36 0.47 1.00 0.47 0.55 0.52 0.37 0.44 \n", + "D 0.40 0.50 0.53 1.00 0.58 0.56 0.42 0.47 \n", + "E 0.32 0.41 0.45 0.42 1.00 0.47 0.33 0.38 \n", + "F 0.34 0.44 0.48 0.44 0.53 1.00 0.35 0.42 \n", + "G 0.47 0.59 0.63 0.58 0.67 0.65 1.00 0.57 \n", + "H 0.42 0.53 0.56 0.53 0.62 0.58 0.43 1.00 \n" + ] + } + ], + "source": [ + "compare=[[sum([better_school(x,i,j) for x in predictions])/len(predictions) for i in range(8)] for j in range(8)]\n", + "l=['A','B','C','D','E','F','G','H']\n", + "print('Empirical Probability that school in given row is\\n as good as or better than corresponding column')\n", + "print(' ',end='')\n", + "print('{0:10}'.format(''))\n", + "for i in range(8):\n", + " print('{0:>10}'.format(l[i]),end='')\n", + "print()\n", + "for j in range(8):\n", + " print('{0:<6}'.format(l[j]),end='')\n", + " for i in range(8):\n", + " print('{0:4.2f} '.format(round(compare[i][j],2)),end=\"\")\n", + " print()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 93bb64e3dc2a2c860f7ca6f550a6e119386871d1 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sat, 21 Apr 2018 16:29:41 -0400 Subject: [PATCH 2/4] Fixed the readme file. --- README.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/README.md b/README.md index f2310b7..85c4867 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,6 @@ # BDA Solutions to some problems from Bayesian Data Analysis, 3rd edition, by Gelman *et. al.* -This repository has selected solutions to Bayesian Data Analysis, 3rd edition, by Gelman *et. al.* - They are offered with no warranty and may be wrong. Comments welcome. From 54ba9a1a66e2acb7c8f7b21d526056e04f1e9457 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sat, 21 Apr 2018 16:40:37 -0400 Subject: [PATCH 3/4] finished off (sort of) 5.9.3 --- BDA 5.9.3.ipynb | 201 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 190 insertions(+), 11 deletions(-) diff --git a/BDA 5.9.3.ipynb b/BDA 5.9.3.ipynb index d155ec7..db62b5d 100644 --- a/BDA 5.9.3.ipynb +++ b/BDA 5.9.3.ipynb @@ -74,18 +74,15 @@ ] }, { - "cell_type": "code", - "execution_count": 69, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## First part" + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1c0a010b4129370aa04f0b4b9f729b4d NOW.\n" - ] - } - ], "source": [ "stan_code='''\n", "data {\n", @@ -277,6 +274,188 @@ " print()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Second part" + ] + }, + { + "cell_type": "code", + "execution_count": 340, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_3314be03c6a1d2db4b5293ee7101b10c NOW.\n" + ] + } + ], + "source": [ + "stan_code_2='''\n", + "data {\n", + " real means[8];\n", + " real se[8];\n", + "\n", + "}\n", + "\n", + "parameters {\n", + " real theta[8] ; \n", + "// real mu ; \n", + "// real tau ; \n", + "}\n", + "\n", + "model {\n", + " \n", + " // theta~normal(mu,tau) ; \n", + " means~normal(theta,se) ; \n", + " \n", + "}\n", + "\n", + "generated quantities {\n", + " real results[8] ; \n", + " \n", + " \n", + " for(i in 1:8) {\n", + " results[i]=normal_rng(theta[i],se[i]);\n", + " }\n", + "}\n", + "'''\n", + "sm=pystan.StanModel(model_code=stan_code_2)" + ] + }, + { + "cell_type": "code", + "execution_count": 341, + "metadata": {}, + "outputs": [], + "source": [ + "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 342, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Inference for Stan model: anon_model_3314be03c6a1d2db4b5293ee7101b10c.\n", + "4 chains, each with iter=5000; warmup=2500; thin=1; \n", + "post-warmup draws per chain=2500, total post-warmup draws=10000.\n", + "\n", + " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", + "theta[0] 27.97 0.15 14.96 -1.53 17.74 27.93 38.15 57.42 10000 1.0\n", + "theta[1] 8.13 0.1 9.94 -11.16 1.34 8.09 14.73 27.7 10000 1.0\n", + "theta[2] -3.24 0.16 16.22 -34.28 -14.16 -3.11 7.61 28.45 10000 1.0\n", + "theta[3] 7.01 0.11 10.91 -14.33 -0.39 7.0 14.49 28.16 10000 1.0\n", + "theta[4] -0.86 0.09 8.9 -18.45 -6.75 -0.86 4.95 16.99 10000 1.0\n", + "theta[5] 1.02 0.11 11.09 -20.34 -6.62 1.13 8.44 22.73 10000 1.0\n", + "theta[6] 17.89 0.1 10.11 -2.11 11.17 17.83 24.75 37.53 10000 1.0\n", + "theta[7] 11.87 0.18 17.87 -22.34 -0.36 11.86 23.8 46.91 10000 1.0\n", + "results[0] 28.25 0.21 21.3 -13.95 13.99 28.29 42.64 69.78 10000 1.0\n", + "results[1] 8.09 0.14 14.22 -19.52 -1.66 8.04 17.7 36.06 10000 1.0\n", + "results[2] -3.5 0.23 22.91 -47.94 -19.39 -3.27 12.05 40.73 10000 1.0\n", + "results[3] 7.01 0.16 15.6 -23.52 -3.69 7.07 17.64 37.45 10000 1.0\n", + "results[4] -0.83 0.13 12.62 -25.85 -9.24 -0.93 7.68 23.68 10000 1.0\n", + "results[5] 0.86 0.15 15.4 -28.99 -9.4 0.83 11.01 32.07 10000 1.0\n", + "results[6] 17.77 0.14 14.11 -10.2 8.31 17.82 27.31 45.08 10000 1.0\n", + "results[7] 11.77 0.25 25.4 -38.49 -5.21 11.54 28.6 61.73 10000 1.0\n", + "lp__ -4.0 0.03 1.99 -8.74 -5.13 -3.68 -2.53 -1.07 4933 1.0\n", + "\n", + "Samples were drawn using NUTS at Sat Apr 21 16:37:25 2018.\n", + "For each parameter, n_eff is a crude measure of effective sample size,\n", + "and Rhat is the potential scale reduction factor on split chains (at \n", + "convergence, Rhat=1)." + ] + }, + "execution_count": 342, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "answers" + ] + }, + { + "cell_type": "code", + "execution_count": 343, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Chance is empirical probability that given school is best\n", + " effect se Chance\n", + "school \n", + "A 28 15 0.4455\n", + "B 8 10 0.0575\n", + "C -3 16 0.0529\n", + "D 7 11 0.0588\n", + "E -1 9 0.0119\n", + "F 1 11 0.0276\n", + "G 18 10 0.1605\n", + "H 12 18 0.1853\n" + ] + } + ], + "source": [ + "predictions=answers.extract()['results']\n", + "def best_school(x,i):\n", + " return x[i]>=max(x)\n", + "def better_school(x,i,j):\n", + " return x[i]>=x[j]\n", + "print('Chance is empirical probability that given school is best')\n", + "p55['Chance']=[sum([best_school(x,i) for x in predictions])/len(predictions) for i in range(8)]\n", + "print(p55)" + ] + }, + { + "cell_type": "code", + "execution_count": 344, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Empirical Probability that school in given row is\n", + " as good as or better than corresponding column\n", + " \n", + " A B C D E F G H\n", + "A 1.00 0.78 0.84 0.79 0.88 0.85 0.66 0.70 \n", + "B 0.22 1.00 0.67 0.52 0.68 0.64 0.31 0.45 \n", + "C 0.16 0.33 1.00 0.36 0.46 0.45 0.21 0.33 \n", + "D 0.21 0.48 0.64 1.00 0.65 0.61 0.30 0.44 \n", + "E 0.12 0.32 0.54 0.35 1.00 0.47 0.16 0.33 \n", + "F 0.15 0.36 0.55 0.39 0.53 1.00 0.21 0.36 \n", + "G 0.34 0.69 0.79 0.70 0.84 0.79 1.00 0.58 \n", + "H 0.30 0.55 0.67 0.56 0.67 0.64 0.42 1.00 \n" + ] + } + ], + "source": [ + "compare=[[sum([better_school(x,i,j) for x in predictions])/len(predictions) for i in range(8)] for j in range(8)]\n", + "l=['A','B','C','D','E','F','G','H']\n", + "print('Empirical Probability that school in given row is\\n as good as or better than corresponding column')\n", + "print(' ',end='')\n", + "print('{0:10}'.format(''))\n", + "for i in range(8):\n", + " print('{0:>10}'.format(l[i]),end='')\n", + "print()\n", + "for j in range(8):\n", + " print('{0:<6}'.format(l[j]),end='')\n", + " for i in range(8):\n", + " print('{0:4.2f} '.format(round(compare[i][j],2)),end=\"\")\n", + " print()" + ] + }, { "cell_type": "code", "execution_count": null, From 006e57b3b93d174ca9a0a366973319d9ac4cd3e4 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sat, 21 Apr 2018 16:46:34 -0400 Subject: [PATCH 4/4] Pointed out a place where things look fishy --- BDA 5.9.3.ipynb | 90 ++++++++++++++++++++++++++++--------------------- 1 file changed, 51 insertions(+), 39 deletions(-) diff --git a/BDA 5.9.3.ipynb b/BDA 5.9.3.ipynb index db62b5d..abfbe93 100644 --- a/BDA 5.9.3.ipynb +++ b/BDA 5.9.3.ipynb @@ -81,8 +81,18 @@ ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 349, "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1c0a010b4129370aa04f0b4b9f729b4d NOW.\n" + ] + } + ], "source": [ "stan_code='''\n", "data {\n", @@ -118,16 +128,16 @@ }, { "cell_type": "code", - "execution_count": 264, + "execution_count": 350, "metadata": {}, "outputs": [], "source": [ - "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)" + "answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=00)" ] }, { "cell_type": "code", - "execution_count": 265, + "execution_count": 351, "metadata": {}, "outputs": [ { @@ -135,31 +145,31 @@ "output_type": "stream", "text": [ "Inference for Stan model: anon_model_1c0a010b4129370aa04f0b4b9f729b4d.\n", - "4 chains, each with iter=5000; warmup=2500; thin=1; \n", - "post-warmup draws per chain=2500, total post-warmup draws=10000.\n", + "4 chains, each with iter=500; warmup=250; thin=1; \n", + "post-warmup draws per chain=250, total post-warmup draws=1000.\n", "\n", " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", - "theta[0] 12.02 0.16 8.34 -2.36 6.78 11.04 16.71 31.15 2751 1.0\n", - "theta[1] 8.09 0.16 6.45 -4.95 4.02 8.05 12.37 20.76 1639 1.0\n", - "theta[2] 6.35 0.15 8.02 -11.67 2.0 6.96 11.29 20.8 2679 1.0\n", - "theta[3] 7.93 0.17 6.82 -6.17 3.85 7.98 12.22 21.32 1702 1.0\n", - "theta[4] 5.01 0.21 6.48 -8.91 1.17 5.27 9.39 16.85 934 1.0\n", - "theta[5] 6.19 0.17 6.8 -8.48 2.27 6.49 10.63 18.67 1624 1.0\n", - "theta[6] 11.18 0.14 6.73 -1.17 6.74 10.84 15.34 25.85 2277 1.0\n", - "theta[7] 8.77 0.16 8.16 -7.36 3.94 8.58 13.42 25.97 2654 1.0\n", - "mu 8.15 0.15 5.19 -2.15 4.84 8.09 11.43 18.47 1215 1.0\n", - "tau 7.11 0.2 5.17 1.37 3.35 5.84 9.47 20.71 700 1.01\n", - "results[0] 12.0 0.18 12.34 -9.06 5.18 10.65 17.58 40.77 4545 1.0\n", - "results[1] 8.13 0.16 10.95 -13.75 2.61 8.11 13.74 30.22 4412 1.0\n", - "results[2] 6.46 0.16 11.79 -19.82 0.76 7.19 12.86 28.9 5157 1.0\n", - "results[3] 7.82 0.15 11.12 -15.81 2.21 7.88 13.72 29.86 5431 1.0\n", - "results[4] 5.04 0.19 10.95 -20.09 -0.12 5.9 11.33 24.39 3484 1.0\n", - "results[5] 6.08 0.17 11.29 -18.73 0.73 6.67 12.38 27.52 4227 1.0\n", - "results[6] 10.98 0.16 11.13 -9.55 4.98 10.24 16.19 36.05 4864 1.0\n", - "results[7] 8.96 0.17 11.99 -15.44 2.81 8.58 14.76 34.5 4746 1.0\n", - "lp__ -18.29 0.26 4.93 -27.55 -21.75 -18.51 -14.86 -8.66 367 1.02\n", + "theta[0] 12.07 0.66 9.38 -3.89 5.93 11.21 17.23 34.12 200 1.01\n", + "theta[1] 7.54 0.39 6.89 -6.07 3.02 7.48 12.15 21.17 312 1.01\n", + "theta[2] 5.45 0.46 8.58 -12.34 0.2 6.07 10.92 20.96 342 1.0\n", + "theta[3] 7.16 0.37 7.21 -7.68 2.72 7.31 11.83 20.44 383 1.01\n", + "theta[4] 4.51 0.43 6.58 -10.0 0.23 5.07 9.18 15.99 235 1.02\n", + "theta[5] 5.53 0.41 7.27 -10.59 1.01 6.16 10.74 18.24 313 1.01\n", + "theta[6] 11.03 0.53 7.5 -2.83 5.64 11.03 15.77 27.15 199 1.01\n", + "theta[7] 8.3 0.39 7.93 -6.88 3.08 8.28 13.26 24.5 414 1.0\n", + "mu 7.66 0.41 5.62 -3.36 4.06 7.92 11.47 18.73 192 1.02\n", + "tau 7.74 0.71 6.11 0.4 3.91 6.29 9.9 22.49 75 1.05\n", + "results[0] 11.81 0.57 12.82 -9.18 3.71 10.54 17.23 43.38 513 1.01\n", + "results[1] 7.26 0.52 12.33 -16.66 0.94 7.04 13.57 33.52 563 1.01\n", + "results[2] 5.43 0.42 12.64 -22.88 -1.07 5.99 12.55 29.09 902 1.0\n", + "results[3] 7.01 0.41 11.92 -16.05 1.03 7.18 13.39 29.57 861 1.0\n", + "results[4] 4.62 0.46 11.55 -23.93 -1.48 5.45 11.65 26.17 638 1.0\n", + "results[5] 5.83 0.44 11.71 -20.83 0.05 6.61 11.82 28.91 724 1.0\n", + "results[6] 10.76 0.56 11.87 -11.69 3.67 10.23 17.43 36.29 443 1.0\n", + "results[7] 8.01 0.54 12.75 -16.27 1.62 7.82 14.62 31.63 561 1.0\n", + "lp__ -18.65 1.21 5.8 -27.82 -22.04 -19.13 -16.31 0.06 23 1.18\n", "\n", - "Samples were drawn using NUTS at Sat Apr 21 15:56:14 2018.\n", + "Samples were drawn using NUTS at Sat Apr 21 16:45:21 2018.\n", "For each parameter, n_eff is a crude measure of effective sample size,\n", "and Rhat is the potential scale reduction factor on split chains (at \n", "convergence, Rhat=1).\n" @@ -172,12 +182,12 @@ }, { "cell_type": "code", - "execution_count": 266, + "execution_count": 352, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF3BJREFUeJzt3X+0XWV95/H3pwGiFQckhFECGihojVPL1Ii2amuhtaGosQ6UoLa0wxRnjUzrUpeC0yrShULXjNRRZs2wBMWfQPFHMyWKjNTWWqUJiGKKlEijhFgJBFCqCIHv/LF39Hhyb+659x5yufd5v9a665797Gfv/TznnPs5z332OfukqpAkteGn5roBkqQ9x9CXpIYY+pLUEENfkhpi6EtSQwx9SWqIoa+xSnJfksP34PFekOTmMe7vU0lO6W//XpK/G+O+X5nkM+Pa3zSO+7wkt/SPzcv29PH16GLoz1NJNif5Qf+H/J0k70uy7yz2tzxJJdlrNu2qqn2r6tbZ7GOgTWcleTDJ9/qff0ryniRPGjje56vqaSPu60NT1auq46rqkjG0fZf7s6o+XFUvmu2+Z+Bs4D39Y/PJ4ZXjfi5NJsnqJDck+W6SO5N8Nsnyft2Uj7XGw9Cf315SVfsCvwA8G/jjuWrIbF8sdrP9ZVX1eOAA4LeAJwLXjTsM0lmofw9PATZOUecRfS4lOQL4APB6YD/gMOB/AQ8PVNsjj3XrFuqTvClVdTvwKeDfASQ5OMnaJNuTbEryBzvrJjk6yYZ+tPWdJO/sV/1t//uefsT3i339/5jkpiR3J7kqyVMG9lVJXpPkFuCWgbIj+tv7JflAkm1Jvpnkj3cGaz918oUk5yfZDpw1RR8frKqNwEnANrrwIMkLk2wZaNObktzejxZvTnJsklXAm4GT+r59pa/7uSTnJPkC8H3g8L7sPw0cOkneneTeJF9PcuzAis1Jfm1gefC/iV3uz+HpoiS/lGR9v+/1SX5pYN3nkvxpfx99L8lnkhw42f2T5A/6x3p7/9gf3Jd/Azgc+L99OxZPcT9P57l0VpIrklzWt/H6JD8/ya6PAv65qj5bne9V1ceq6lsTtGHCx1rjYegvAEkOBX4T+HJf9FFgC3AwcALw9oGwehfwrqr6N8DPAJf35b/c/96/nwb4Yrr53zcDLweWAp/v9z3oZcBzgBUTNO3ddKO6w4FfAX4X+P2B9c8BbgUOAs4Zpa9V9RDwl8ALhtcleRpwOvDsfsT4G8Dmqvo08Ha6keS+VTUYTL8DnAY8HvjmBIfc2cYDgbcCH09ywAhN3eX+HGrrAcCVwP8ElgDvBK5MsmSg2ivo7q+DgH2AN0x0oCTHAO8Afht4Ut+PSwGq6meAb9GP5Kvqh7tr9DSfSwCrgb+gG51/BPhkkr0n2PX1wM/2L/K/mhGmj3b3WGvmDP357ZNJ7gH+Dvgbuj/IQ4HnA2+qqvur6gbgvXThBvAgcESSA6vqvqr60m72/2rgHVV1U1XtoAvOowZH+/367VX1g8ENkyyiG6md2Y/qNgP/Y6AdAFur6t1VtWN4+ylspQuZYQ8Bi4EVSfauqs1V9Y0p9vX+qtrYt+HBCdbfAfx5P/q8DLgZOH4abZ3M8cAtVfXB/tgfBb4OvGSgzvuq6p/6++ZyutHyRF4JXFxV1/ehfibwi+nny0c0k+cSwHVVdUV/370TeAzw3OGd9+d5Xggs6/tyZ5L3jxD+kz3WmiFDf357WVXtX1VPqar/0ofDwcD2qvreQL1v0v2xAZwKPBX4ej+l8OLd7P8pwLuS3NMHwnYgA/sCuG2SbQ+kG50Ojp4H27G7baeyrG/LT6iqTcBr6aaK7khy6c5pjt2Yqg23109elfCbdPfxbB3Mrv9ZDN8//zJw+/vAZAH5E/uqqvuAu4b2NZWZPJdg4P6rqof58X8Fu6iqL1XVb1fVUrrR+y8D/22Kdk34WGvmDP2FZytwQJLHD5Q9GbgdoKpuqaqT6aYMzgOuSPI4YKLLrd4GvLoPg50/j62qvx+oM9llWu+k+69i8L+CH7Vjim0n1Z8TeAndVNMuquojVfX8/rhF18fdHWuqNixLkoHlJ9PdxwD/Cvz0wLonTmO/W/nJ+2bnvm+foO5UfmJf/eO5ZIb7Gt7vpM+l3qEDx/0p4BB+fP9MqqrWAx+nP3cwkakea82Mob/AVNVtwN8D70jymCTPpBvdfxggyauSLO1HZff0mz1Ed8LsYbr5953+N3Bmkmf02+6X5MQR2/EQ3b/x5yR5fD8l9DpgyrdNTiTJ3kmeTjfH/ES6qYThOk9Lckx/svJ+4Ad93wC+AyzP9N+hcxDwh/3xTwSeDqzr190ArOnXraSb895povtz0DrgqUlekWSvJCfRnRf5q2m2D7q59N9PclTf97cD1/ZTajM21XOp96wkL0/37qvXAj8EdpkyTPL8/mTzQf3yzwIvnaTulI+1Zs7QX5hOBpbTjbg+Aby1qq7u160CNia5j+6k7pp+vvb7dCdTv9BP5zy3qj5BN1K+NMl3ga8Bx02jHf+VbjR8K91c8UeAi6fZl5P6tt4DrKWbtnhWVU00mlwMnEv3X8a/0AX2m/t1f9H/vivJ9dM4/rXAkf0+zwFOqKq7+nV/Qncy/G7gbXT9A2Ci+3Nwp/0+Xkz3zpS7gDcCL66qO6fRtp37+mzflo8B3+7btGa6+5nE7p5L0J1oPYnuPvgd4OWTnBu5hy7kb+wfz0/3+/uzgTrTeaw1Q/FLVCTNRJKzgCOq6lVz3RaNzpG+JDXE0Jekhji9I0kNcaQvSQ2Z1UWyHgkHHnhgLV++fK6bIUnzynXXXXdn/8G33XrUhf7y5cvZsGHDXDdDkuaVJBNdO2oXTu9IUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDHnWfyH2kLD/jygnLN587ju+4lqT5wZG+JDVkwY30JxvRS5Ic6UtSUwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGLLi3bE6XH9qS1BJH+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaMlLoJ1mV5OYkm5KcMcH6X05yfZIdSU4YWndKklv6n1PG1XBJ0vRNGfpJFgEXAMcBK4CTk6wYqvYt4PeAjwxtewDwVuA5wNHAW5M8YfbNliTNxCgj/aOBTVV1a1U9AFwKrB6sUFWbq+qrwMND2/4GcHVVba+qu4GrgVVjaLckaQZGCf1lwG0Dy1v6slGMtG2S05JsSLJh27ZtI+5akjRdo4R+JiirEfc/0rZVdWFVrayqlUuXLh1x15Kk6Rol9LcAhw4sHwJsHXH/s9lWkjRmo4T+euDIJIcl2QdYA6wdcf9XAS9K8oT+BO6L+jJJ0hyYMvSragdwOl1Y3wRcXlUbk5yd5KUASZ6dZAtwIvB/kmzst90O/CndC8d64Oy+TJI0B0b6usSqWgesGyp7y8Dt9XRTNxNtezFw8SzaKEkaEz+RK0kNMfQlqSGGviQ1ZKQ5/RYtP+PKCcs3n3v8Hm6JJI2PI31JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ/xw1jT5oS1J85kjfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUkJFCP8mqJDcn2ZTkjAnWL05yWb/+2iTL+/K9k1yS5MYkNyU5c7zNlyRNx5Shn2QRcAFwHLACODnJiqFqpwJ3V9URwPnAeX35icDiqvo54FnAq3e+IEiS9rxRRvpHA5uq6taqegC4FFg9VGc1cEl/+wrg2CQBCnhckr2AxwIPAN8dS8slSdM2SugvA24bWN7Sl01Yp6p2APcCS+heAP4V+DbwLeC/V9X2WbZZkjRDo4R+JiirEescDTwEHAwcBrw+yeG7HCA5LcmGJBu2bds2QpMkSTMxSuhvAQ4dWD4E2DpZnX4qZz9gO/AK4NNV9WBV3QF8AVg5fICqurCqVlbVyqVLl06/F5KkkYwS+uuBI5MclmQfYA2wdqjOWuCU/vYJwDVVVXRTOsek8zjgucDXx9N0SdJ0TRn6/Rz96cBVwE3A5VW1McnZSV7aV7sIWJJkE/A6YOfbOi8A9gW+Rvfi8b6q+uqY+yBJGtFeo1SqqnXAuqGytwzcvp/u7ZnD2903UbkkaW74iVxJaoihL0kNMfQlqSGGviQ1ZKQTuZra8jOunLB887nH7+GWSNLkHOlLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhfjH6I8wvTJf0aOJIX5IaYuhLUkNGCv0kq5LcnGRTkjMmWL84yWX9+muTLB9Y98wkX0yyMcmNSR4zvuZLkqZjytBPsgi4ADgOWAGcnGTFULVTgbur6gjgfOC8ftu9gA8B/7mqngG8EHhwbK2XJE3LKCP9o4FNVXVrVT0AXAqsHqqzGrikv30FcGySAC8CvlpVXwGoqruq6qHxNF2SNF2jhP4y4LaB5S192YR1qmoHcC+wBHgqUEmuSnJ9kjdOdIAkpyXZkGTDtm3bptsHSdKIRgn9TFBWI9bZC3g+8Mr+928lOXaXilUXVtXKqlq5dOnSEZokSZqJUUJ/C3DowPIhwNbJ6vTz+PsB2/vyv6mqO6vq+8A64Bdm22hJ0syMEvrrgSOTHJZkH2ANsHaozlrglP72CcA1VVXAVcAzk/x0/2LwK8A/jqfpkqTpmvITuVW1I8npdAG+CLi4qjYmORvYUFVrgYuADybZRDfCX9Nve3eSd9K9cBSwrqom/oiqJOkRN9JlGKpqHd3UzGDZWwZu3w+cOMm2H6J726YkaY75iVxJaoihL0kNMfQlqSFeWnmOeMllSXPBkb4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGuIF1x5lvBCbpEeSI31JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhI4V+klVJbk6yKckZE6xfnOSyfv21SZYPrX9ykvuSvGE8zZYkzcSUoZ9kEXABcBywAjg5yYqhaqcCd1fVEcD5wHlD688HPjX75kqSZmOUkf7RwKaqurWqHgAuBVYP1VkNXNLfvgI4NkkAkrwMuBXYOJ4mS5JmapTQXwbcNrC8pS+bsE5V7QDuBZYkeRzwJuBtuztAktOSbEiyYdu2baO2XZI0TaOEfiYoqxHrvA04v6ru290BqurCqlpZVSuXLl06QpMkSTMxypeobAEOHVg+BNg6SZ0tSfYC9gO2A88BTkjyZ8D+wMNJ7q+q98y65ZKkaRsl9NcDRyY5DLgdWAO8YqjOWuAU4IvACcA1VVXAC3ZWSHIWcJ+BL0lzZ8rQr6odSU4HrgIWARdX1cYkZwMbqmotcBHwwSSb6Eb4ax7JRkuSZmak78itqnXAuqGytwzcvh84cYp9nDWD9kmSxshP5EpSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1JCR3qevubf8jCsnLN987vF7uCWS5jNH+pLUEENfkhpi6EtSQwx9SWqIJ3LnuclO8IIneSXtypG+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhriBdcWML9tS9IwR/qS1BBDX5IaMlLoJ1mV5OYkm5KcMcH6xUku69dfm2R5X/7rSa5LcmP/+5jxNl+SNB1Thn6SRcAFwHHACuDkJCuGqp0K3F1VRwDnA+f15XcCL6mqnwNOAT44roZLkqZvlJH+0cCmqrq1qh4ALgVWD9VZDVzS374CODZJqurLVbW1L98IPCbJ4nE0XJI0faOE/jLgtoHlLX3ZhHWqagdwL7BkqM5/AL5cVT8cPkCS05JsSLJh27Zto7ZdkjRNo4R+Jiir6dRJ8gy6KZ9XT3SAqrqwqlZW1cqlS5eO0CRJ0kyMEvpbgEMHlg8Btk5WJ8lewH7A9n75EOATwO9W1Tdm22BJ0syNEvrrgSOTHJZkH2ANsHaozlq6E7UAJwDXVFUl2R+4Ejizqr4wrkZLkmZmytDv5+hPB64CbgIur6qNSc5O8tK+2kXAkiSbgNcBO9/WeTpwBPAnSW7ofw4aey8kSSNJ1fD0/NxauXJlbdiwYcbbT3bpAU3NyzNI81eS66pq5VT1/ESuJDXE0Jekhhj6ktQQQ1+SGmLoS1JD/BIV/YhfuiItfI70Jakhhr4kNcTQl6SGOKevKTnXLy0cjvQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqI79PXjE33W8p8X7809xzpS1JDDH1JaoihL0kNMfQlqSGeyNWc84Ju0p5j6OtRyxcDafwMfe0x032Lp6Txc05fkhriSF/zjtM+0swZ+lowZjJ95AuFWjNS6CdZBbwLWAS8t6rOHVq/GPgA8CzgLuCkqtrcrzsTOBV4CPjDqrpqbK2XZslLSag1U4Z+kkXABcCvA1uA9UnWVtU/DlQ7Fbi7qo5IsgY4DzgpyQpgDfAM4GDg/yV5alU9NO6OSHvCdKeWxjUV5ZSWxmWUkf7RwKaquhUgyaXAamAw9FcDZ/W3rwDekyR9+aVV9UPgn5Ns6vf3xfE0X3p0mO5/DON6J9O4XoR2t40WllFCfxlw28DyFuA5k9Wpqh1J7gWW9OVfGtp22fABkpwGnNYv3pfk5inadCBw5whtn8/s4/w3Z/3LeXtmGxb+Ywjzp49PGaXSKKGfCcpqxDqjbEtVXQhcOEJbuoMlG6pq5aj15yP7OP8t9P6BfZyPRnmf/hbg0IHlQ4Ctk9VJshewH7B9xG0lSXvIKKG/HjgyyWFJ9qE7Mbt2qM5a4JT+9gnANVVVffmaJIuTHAYcCfzDeJouSZquKad3+jn604Gr6N6yeXFVbUxyNrChqtYCFwEf7E/Ubqd7YaCvdzndSd8dwGvG9M6dkaeC5jH7OP8t9P6BfZx30g3IJUkt8No7ktQQQ1+SGjLvQj/JqiQ3J9mU5Iy5bs84JLk4yR1JvjZQdkCSq5Pc0v9+wly2cTaSHJrkr5PclGRjkj/qyxdSHx+T5B+SfKXv49v68sOSXNv38bL+zRDzVpJFSb6c5K/65QXVP4Akm5PcmOSGJBv6sgXzXJ1XoT9wSYjjgBXAyf2lHua79wOrhsrOAD5bVUcCn+2X56sdwOur6unAc4HX9I/bQurjD4FjqurngaOAVUmeS3dJkvP7Pt5Nd8mS+eyPgJsGlhda/3b61ao6auD9+QvmuTqvQp+BS0JU1QPAzktCzGtV9bd073oatBq4pL99CfCyPdqoMaqqb1fV9f3t79GFxjIWVh+rqu7rF/fufwo4hu7SJDDP+5jkEOB44L39clhA/ZvCgnmuzrfQn+iSELtc1mGB+LdV9W3oQhM4aI7bMxZJlgP/HriWBdbHfurjBuAO4GrgG8A9VbWjrzLfn69/DrwReLhfXsLC6t9OBXwmyXX9JWJgAT1X59v19Ee6rIMenZLsC3wMeG1VfbcbKC4c/WdQjkqyP/AJ4OkTVduzrRqPJC8G7qiq65K8cGfxBFXnZf+GPK+qtiY5CLg6ydfnukHjNN9G+i1d1uE7SZ4E0P++Y47bMytJ9qYL/A9X1cf74gXVx52q6h7gc3TnL/bvL00C8/v5+jzgpUk2002rHkM38l8o/fuRqtra/76D7sX7aBbQc3W+hf4ol4RYKAYvbXEK8Jdz2JZZ6ed+LwJuqqp3DqxaSH1c2o/wSfJY4Nfozl38Nd2lSWAe97GqzqyqQ6pqOd3f3TVV9UoWSP92SvK4JI/feRt4EfA1FtJzdb59IjfJb9KNMHZeEuKcOW7SrCX5KPBCuku4fgd4K/BJ4HLgycC3gBOravhk77yQ5PnA54Eb+fF88Jvp5vUXSh+fSXeCbxHdYOryqjo7yeF0I+MDgC8Dr+q/X2Le6qd33lBVL15o/ev784l+cS/gI1V1TpIlLJTn6nwLfUnSzM236R1J0iwY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakh/x90grO/jaMYMQAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGC1JREFUeJzt3XuUHGWdxvHvswkEBQyQhBWSyAQTkbAiSgyoiEoUgyhBN0gianRZcc/Cqkc9mngBxBt4dkFW2UuOoIiXgPGWhWhkQV1FjRkuAmOIDDGQIQoTEsCIXBJ++0e9o0WnZ6Z6pjPTM+/zOWfOdL31VvWvq2uefuftnhpFBGZmloe/Ge4CzMxs6Dj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tC3ppK0TdLBQ3h/L5O0ron7+76kRen22yX9rIn7Pk3SD5u1vwbu96WS7kzPzclDff/WWhz6I5SkDZL+nH6Q75P0JUl7DWJ/bZJC0tjB1BURe0XE+sHso1TTuZKekPTH9PVbSV+QdEDp/n4aEYdU3NdX++sXESdExOVNqH2n4xkRX4uI4we77wE4D/hCem6+W7uy2edSbyTNk3SLpIclbZZ0naS2tK7f59qaw6E/sr0+IvYCXgi8CPjocBUy2BeLPra/MiL2BvYD3gA8E7ix2WGgwmj9eTgI6Oinzy49lyRNB74CvB8YD0wD/gN4stRtSJ7r3I3WkzwrEXEv8H3g7wAkHShphaQtkjolvbOnr6TZktrTaOs+SRemVf+Xvj+YRnwvTv3/QdJaSVslrZJ0UGlfIelMSXcCd5bapqfb4yV9RVK3pLslfbQnWNPUyQ2SLpK0BTi3n8f4RER0AKcC3RThgaRXSOoq1fQhSfem0eI6SXMkzQU+DJyaHtuvU98fS/qUpBuAR4CDU9s/lu5akj4v6SFJd0iaU1qxQdKrSsvl3yZ2Op6100WSXiJpTdr3GkkvKa37saRPpGP0R0k/lDSxt+Mj6Z3pud6SnvsDU/tdwMHA/6Q6xvVznBs5l86VtFzSlanGmyQ9v5ddHwH8LiKui8IfI+JbEXFPnRrqPtfWHA79UUDSVOC1wM2p6RtAF3AgMB/4dCmsLgYujohnAM8Grkrtx6bv+6RpgF+omP/9MPBGYBLw07TvspOBo4CZdUr7PMWo7mDg5cDbgHeU1h8FrAf2Bz5V5bFGxA7ge8DLatdJOgQ4C3hRGjG+BtgQET8APk0xktwrIsrB9FbgDGBv4O46d9lT40TgHODbkvarUOpOx7Om1v2Aa4B/ByYAFwLXSJpQ6vZmiuO1P7A78IF6dyTpOOAzwJuAA9LjWAYQEc8G7iGN5CPisb6KbvBcApgHfJNidP514LuSdquz65uA56YX+VeqwvRRX8+1DZxDf2T7rqQHgZ8BP6H4gZwKHAN8KCIejYhbgC9ShBvAE8B0SRMjYltE/LKP/b8L+ExErI2I7RTBeUR5tJ/Wb4mIP5c3lDSGYqS2JI3qNgD/VqoDYFNEfD4ittdu349NFCFTawcwDpgpabeI2BARd/Wzry9HREeq4Yk66+8HPpdGn1cC64ATG6i1NycCd0bEFem+vwHcAby+1OdLEfHbdGyuohgt13MacFlE3JRCfQnwYqX58ooGci4B3BgRy9OxuxDYAzi6dufpfZ5XAJPTY9ks6csVwr+359oGyKE/sp0cEftExEER8c8pHA4EtkTEH0v97qb4YQM4HXgOcEeaUnhdH/s/CLhY0oMpELYAKu0LYGMv206kGJ2WR8/lOvratj+TUy1PERGdwHspporul7SsZ5qjD/3VcG889aqEd1Mc48E6kJ1/s6g9Pn8o3X4E6C0gn7KviNgGPFCzr/4M5FyC0vGLiCf5628FO4mIX0bEmyJiEsXo/VjgI/3UVfe5toFz6I8+m4D9JO1dansWcC9ARNwZEQsppgwuAJZL2hOod7nVjcC7Uhj0fD0tIn5e6tPbZVo3U/xWUf6t4C919LNtr9J7Aq+nmGraSUR8PSKOSfcbFI+xr/vqr4bJklRafhbFMQb4E/D00rpnNrDfTTz12PTs+946ffvzlH2l53PCAPdVu99ez6Vkaul+/waYwl+PT68iYg3wbdJ7B/X091zbwDj0R5mI2Aj8HPiMpD0kHU4xuv8agKS3SJqURmUPps12ULxh9iTF/HuP/wKWSDosbTte0ikV69hB8Wv8pyTtnaaE3gf0+7HJeiTtJulQijnmZ1JMJdT2OUTScenNykeBP6fHBnAf0KbGP6GzP/DudP+nAIcCK9O6W4AFad0sijnvHvWOZ9lK4DmS3ixprKRTKd4XubrB+qCYS3+HpCPSY/80sDpNqQ1Yf+dScqSkN6r49NV7gceAnaYMJR2T3mzePy0/Fzipl779Ptc2cA790Wkh0EYx4voOcE5EXJvWzQU6JG2jeFN3QZqvfYTizdQb0nTO0RHxHYqR8jJJDwO3Ayc0UMe/UIyG11PMFX8duKzBx3JqqvVBYAXFtMWREVFvNDkOOJ/it4w/UAT2h9O6b6bvD0i6qYH7Xw3MSPv8FDA/Ih5I6z5G8Wb4VuDjFI8PgHrHs7zTtI/XUXwy5QHgg8DrImJzA7X17Ou6VMu3gN+nmhY0up9e9HUuQfFG66kUx+CtwBt7eW/kQYqQvy09nz9I+/tsqU8jz7UNkPxPVMxsICSdC0yPiLcMdy1WnUf6ZmYZceibmWXE0ztmZhnxSN/MLCODukjWrjBx4sRoa2sb7jLMzEaUG2+8cXP6w7c+tVzot7W10d7ePtxlmJmNKJLqXTtqJ57eMTPLiEPfzCwjDn0zs4w49M3MMuLQNzPLiEPfzCwjDn0zs4w49M3MMuLQNzPLSMv9Re5Qa1t8Td32Dec3439fm5m1Fo/0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8tIpdCXNFfSOkmdkhbXWX+spJskbZc0v2bdIkl3pq9FzSrczMwa12/oSxoDXAKcAMwEFkqaWdPtHuDtwNdrtt0POAc4CpgNnCNp38GXbWZmA1FlpD8b6IyI9RHxOLAMmFfuEBEbIuJW4MmabV8DXBsRWyJiK3AtMLcJdZuZ2QBUCf3JwMbScldqq6LStpLOkNQuqb27u7virs3MrFFV/l2i6rRFxf1X2jYilgJLAWbNmlV137tUb/9GsTf+94pmNhJUGel3AVNLy1OATRX3P5htzcysyaqE/hpghqRpknYHFgArKu5/FXC8pH3TG7jHpzYzMxsG/YZ+RGwHzqII67XAVRHRIek8SScBSHqRpC7gFOC/JXWkbbcAn6B44VgDnJfazMxsGFSZ0yciVgIra9rOLt1eQzF1U2/by4DLBlGjmZk1if8i18wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy0il0Jc0V9I6SZ2SFtdZP07SlWn9akltqX03SZdLuk3SWklLmlu+mZk1ot/QlzQGuAQ4AZgJLJQ0s6bb6cDWiJgOXARckNpPAcZFxPOAI4F39bwgmJnZ0Ksy0p8NdEbE+oh4HFgGzKvpMw+4PN1eDsyRJCCAPSWNBZ4GPA483JTKzcysYVVCfzKwsbTcldrq9omI7cBDwASKF4A/Ab8H7gH+NSK2DLJmMzMboCqhrzptUbHPbGAHcCAwDXi/pIN3ugPpDEntktq7u7srlGRmZgNRJfS7gKml5SnApt76pKmc8cAW4M3ADyLiiYi4H7gBmFV7BxGxNCJmRcSsSZMmNf4ozMyskiqhvwaYIWmapN2BBcCKmj4rgEXp9nzg+ogIiimd41TYEzgauKM5pZuZWaP6Df00R38WsApYC1wVER2SzpN0Uup2KTBBUifwPqDnY52XAHsBt1O8eHwpIm5t8mMwM7OKxlbpFBErgZU1bWeXbj9K8fHM2u221Ws3M7Ph4b/INTPLiEPfzCwjDn0zs4w49M3MMlLpjdzRoG3xNcNdgpnZsPNI38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMZPPvEne13v4d44bzTxziSszMeueRvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRiqFvqS5ktZJ6pS0uM76cZKuTOtXS2orrTtc0i8kdUi6TdIezSvfzMwa0W/oSxoDXAKcAMwEFkqaWdPtdGBrREwHLgIuSNuOBb4K/FNEHAa8AniiadWbmVlDqoz0ZwOdEbE+Ih4HlgHzavrMAy5Pt5cDcyQJOB64NSJ+DRARD0TEjuaUbmZmjaoS+pOBjaXlrtRWt09EbAceAiYAzwFC0ipJN0n6YL07kHSGpHZJ7d3d3Y0+BjMzq6hK6KtOW1TsMxY4BjgtfX+DpDk7dYxYGhGzImLWpEmTKpRkZmYDUSX0u4CppeUpwKbe+qR5/PHAltT+k4jYHBGPACuBFw62aDMzG5gqob8GmCFpmqTdgQXAipo+K4BF6fZ84PqICGAVcLikp6cXg5cDv2lO6WZm1qh+L7gWEdslnUUR4GOAyyKiQ9J5QHtErAAuBa6Q1Ekxwl+Qtt0q6UKKF44AVkZE/SuTmZnZLlfpKpsRsZJiaqbcdnbp9qPAKb1s+1WKj22amdkw81/kmpllxKFvZpYRh76ZWUYc+mZmGfG/S9zF/G8UzayVeKRvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUbGDncBuWpbfE3d9g3nnzjElZhZTjzSNzPLSKXQlzRX0jpJnZIW11k/TtKVaf1qSW01658laZukDzSnbDMzG4h+Q1/SGOAS4ARgJrBQ0syabqcDWyNiOnARcEHN+ouA7w++XDMzG4wqI/3ZQGdErI+Ix4FlwLyaPvOAy9Pt5cAcSQKQdDKwHuhoTslmZjZQVUJ/MrCxtNyV2ur2iYjtwEPABEl7Ah8CPt7XHUg6Q1K7pPbu7u6qtZuZWYOqhL7qtEXFPh8HLoqIbX3dQUQsjYhZETFr0qRJFUoyM7OBqPKRzS5gaml5CrCplz5dksYC44EtwFHAfEmfBfYBnpT0aER8YdCVm5lZw6qE/hpghqRpwL3AAuDNNX1WAIuAXwDzgesjIoCX9XSQdC6wzYFvZjZ8+g39iNgu6SxgFTAGuCwiOiSdB7RHxArgUuAKSZ0UI/wFu7JoMzMbmEp/kRsRK4GVNW1nl24/CpzSzz7OHUB9ZmbWRP6LXDOzjIy6a+/0dk0bMzPzSN/MLCsOfTOzjDj0zcwyMurm9Ec6X2ffzHYlj/TNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsI5VCX9JcSeskdUpaXGf9OElXpvWrJbWl9ldLulHSben7cc0t38zMGtFv6EsaA1wCnADMBBZKmlnT7XRga0RMBy4CLkjtm4HXR8TzgEXAFc0q3MzMGldlpD8b6IyI9RHxOLAMmFfTZx5webq9HJgjSRFxc0RsSu0dwB6SxjWjcDMza1yV0J8MbCwtd6W2un0iYjvwEDChps/fAzdHxGO1dyDpDEntktq7u7ur1m5mZg0aW6GP6rRFI30kHUYx5XN8vTuIiKXAUoBZs2bV7tuAtsXX1G3fcP6JQ1yJmY1kVUb6XcDU0vIUYFNvfSSNBcYDW9LyFOA7wNsi4q7BFmxmZgNXJfTXADMkTZO0O7AAWFHTZwXFG7UA84HrIyIk7QNcAyyJiBuaVbSZmQ1Mv6Gf5ujPAlYBa4GrIqJD0nmSTkrdLgUmSOoE3gf0fKzzLGA68DFJt6Sv/Zv+KMzMrJIqc/pExEpgZU3b2aXbjwKn1Nnuk8AnB1mjmZk1if8i18wsIw59M7OMVJresdbV20c5wR/nNLOdeaRvZpYRh76ZWUYc+mZmGXHom5llxKFvZpYRh76ZWUb8kU37C1/J02z0c+hnqK/P9pvZ6ObQH8Uc7mZWy3P6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWXEH9m0AfMfc5mNPB7pm5llxKFvZpYRT+/YkPF0kNnw80jfzCwjDn0zs4x4esf65Qu3mY0eHumbmWXEoW9mlhGHvplZRjynb03X6HsA/iin2dCpFPqS5gIXA2OAL0bE+TXrxwFfAY4EHgBOjYgNad0S4HRgB/DuiFjVtOrNKujrRWikv7D4BdMa1W/oSxoDXAK8GugC1khaERG/KXU7HdgaEdMlLQAuAE6VNBNYABwGHAj8r6TnRMSOZj8Qy0czP03UrNB0+A5cjsduOB9zlTn92UBnRKyPiMeBZcC8mj7zgMvT7eXAHElK7csi4rGI+B3QmfZnZmbDoMr0zmRgY2m5Cziqtz4RsV3SQ8CE1P7Lmm0n196BpDOAM9LiNknrKlW/s4nA5gFuO9RGUq0wDPXqggFvOuhaB3HfA9lP049ts+rvxZCcC016DCPq50wXDKreg6p0qhL6qtMWFftU2ZaIWAosrVBLnyS1R8Sswe5nKIykWmFk1TuSagXXuyuNpFphaOqtMr3TBUwtLU8BNvXWR9JYYDywpeK2ZmY2RKqE/hpghqRpknaneGN2RU2fFcCidHs+cH1ERGpfIGmcpGnADOBXzSndzMwa1e/0TpqjPwtYRfGRzcsiokPSeUB7RKwALgWukNRJMcJfkLbtkHQV8BtgO3DmLv7kzqCniIbQSKoVRla9I6lWcL270kiqFYagXhUDcjMzy4Evw2BmlhGHvplZRkZF6EuaK2mdpE5Ji4e7nlqSLpN0v6TbS237SbpW0p3p+77DWWMPSVMl/UjSWkkdkt6T2lu13j0k/UrSr1O9H0/t0yStTvVemT6E0BIkjZF0s6Sr03Ir17pB0m2SbpHUntpa8lwAkLSPpOWS7kjn8ItbsV5Jh6Rj2vP1sKT3DkWtIz70S5eJOAGYCSxMl39oJV8G5ta0LQaui4gZwHVpuRVsB94fEYcCRwNnpuPZqvU+BhwXEc8HjgDmSjqa4lIgF6V6t1JcKqRVvAdYW1pu5VoBXhkRR5Q+P96q5wIU1wj7QUQ8F3g+xXFuuXojYl06pkdQXLPsEeA7DEWtETGiv4AXA6tKy0uAJcNdV50624DbS8vrgAPS7QOAdcNdYy91f4/iukstXy/wdOAmir8Y3wyMrXeODHONU9IP83HA1RR/wNiStaZ6NgATa9pa8lwAngH8jvQBlVavt1Tf8cANQ1XriB/pU/8yETtd6qEF/W1E/B4gfd9/mOvZiaQ24AXAalq43jRdcgtwP3AtcBfwYERsT11a6Zz4HPBB4Mm0PIHWrRWKv6D/oaQb0+VSoHXPhYOBbuBLafrsi5L2pHXr7bEA+Ea6vctrHQ2hX+lSD9YYSXsB3wLeGxEPD3c9fYmIHVH8mjyF4oJ+h9brNrRV7UzS64D7I+LGcnOdrsNea8lLI+KFFNOnZ0o6drgL6sNY4IXAf0bEC4A/0QJTOX1J79+cBHxzqO5zNIT+SL3Uw32SDgBI3+8f5nr+QtJuFIH/tYj4dmpu2Xp7RMSDwI8p3ovYJ10SBFrnnHgpcJKkDRRXqz2OYuTfirUCEBGb0vf7KeacZ9O650IX0BURq9PycooXgVatF4oX05si4r60vMtrHQ2hX+UyEa2ofOmKRRRz58MuXRL7UmBtRFxYWtWq9U6StE+6/TTgVRRv3v2I4pIg0CL1RsSSiJgSEW0U5+n1EXEaLVgrgKQ9Je3dc5ti7vl2WvRciIg/ABslHZKa5lBcDaAl600W8tepHRiKWof7TYwmvRHyWuC3FHO5HxnueurU9w3g98ATFKOR0ynmcq8D7kzf9xvuOlOtx1BML9wK3JK+XtvC9R4O3JzqvR04O7UfTHGdp06KX53HDXetNXW/Ari6lWtNdf06fXX0/Gy16rmQajsCaE/nw3eBfVu1XooPHjwAjC+17fJafRkGM7OMjIbpHTMzq8ihb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlG/h9vN194yrGc8gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -194,7 +204,7 @@ }, { "cell_type": "code", - "execution_count": 267, + "execution_count": 353, "metadata": {}, "outputs": [], "source": [ @@ -207,7 +217,7 @@ }, { "cell_type": "code", - "execution_count": 333, + "execution_count": 354, "metadata": {}, "outputs": [ { @@ -217,21 +227,23 @@ "Chance is empirical probability that given school is best\n", " effect se Chance\n", "school \n", - "A 28 15 0.2209\n", - "B 8 10 0.1119\n", - "C -3 16 0.0972\n", - "D 7 11 0.1116\n", - "E -1 9 0.0629\n", - "F 1 11 0.0829\n", - "G 18 10 0.1750\n", - "H 12 18 0.1376\n" + "A 28 15 0.191\n", + "B 8 10 0.109\n", + "C -3 16 0.105\n", + "D 7 11 0.124\n", + "E -1 9 0.077\n", + "F 1 11 0.088\n", + "G 18 10 0.186\n", + "H 12 18 0.120\n", + "Only thing that worries me is that Gelman has A best with prob=10%\n" ] } ], "source": [ "print('Chance is empirical probability that given school is best')\n", "p55['Chance']=[sum([best_school(x,i) for x in predictions])/len(predictions) for i in range(8)]\n", - "print(p55)" + "print(p55)\n", + "print(\"Only thing that worries me is that Gelman has A best with prob=10%\")" ] }, {