Problem 5.9.3 (in part)

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<lower=0> 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": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "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
+}