Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

BDA/BDA 5.9.3.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
500 lines (500 sloc) 24.5 KB
Raw Blame
Open in GitHub Desktop
{
"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": "markdown",
"metadata": {},
"source": [
"## First part"
]
},
{
"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",
" 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": 355,
"metadata": {},
"outputs": [],
"source": [
"answers=sm.sampling(data=dict({'means':p55['effect'],'se':p55['se']}),iter=5000)"
]
},
{
"cell_type": "code",
"execution_count": 356,
"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] 11.66 0.17 8.25 -2.16 6.27 10.51 16.07 30.71 2223 1.0\n",
"theta[1] 8.02 0.1 6.27 -4.44 4.18 7.82 11.76 21.02 3674 1.0\n",
"theta[2] 6.36 0.12 7.68 -10.59 2.17 6.78 11.15 20.84 3806 1.0\n",
"theta[3] 7.79 0.11 6.43 -5.29 4.02 7.62 11.79 20.79 3689 1.0\n",
"theta[4] 5.19 0.13 6.38 -8.66 1.46 5.72 9.44 16.84 2357 1.0\n",
"theta[5] 6.3 0.12 6.62 -7.68 2.38 6.54 10.57 18.98 3277 1.0\n",
"theta[6] 10.99 0.15 6.78 -0.97 6.41 10.4 14.99 25.86 2181 1.0\n",
"theta[7] 8.61 0.13 7.97 -7.29 4.13 8.21 12.91 25.66 3946 1.0\n",
"mu 8.2 0.11 5.3 -1.85 5.01 7.96 11.32 18.88 2354 1.0\n",
"tau 6.84 0.2 5.58 0.93 2.99 5.46 9.11 20.69 751 1.01\n",
"results[0] 11.6 0.18 12.19 -8.79 4.97 9.96 16.83 40.26 4692 1.0\n",
"results[1] 7.97 0.13 10.69 -14.26 2.73 7.74 13.36 29.84 6547 1.0\n",
"results[2] 6.4 0.15 11.84 -18.96 1.18 6.95 12.45 28.63 5970 1.0\n",
"results[3] 7.79 0.14 11.1 -14.73 2.66 7.66 13.29 30.0 6438 1.0\n",
"results[4] 5.24 0.16 11.09 -19.44 0.35 6.13 11.03 24.67 4830 1.0\n",
"results[5] 6.36 0.15 10.98 -17.1 1.36 6.66 11.93 28.06 5695 1.0\n",
"results[6] 10.93 0.16 11.19 -8.59 4.94 9.77 16.05 35.59 4800 1.0\n",
"results[7] 8.64 0.14 11.63 -13.75 2.97 8.25 14.03 33.19 6700 1.0\n",
"lp__ -17.68 0.41 5.52 -27.62 -21.48 -18.16 -14.18 -5.46 179 1.04\n",
"\n",
"Samples were drawn using NUTS at Sat Apr 21 16:45:42 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": 357,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFzBJREFUeJzt3X+0XWV95/H3p4lEBQsSwigBvKHgj9gfViOoVetAa6GosQ6UoLbUocVZI2MddWlwWkVaFbo6UkeZdliCxZ9A8UczJUodaTvWKpMAVk2REmOEEH8EAigqQuA7f+yderw5N/fcm0su9z7v11p33bOf/ey9n7PPvp/znGfvs2+qCklSG35qthsgSdp7DH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+ppRSe5OcsRe3N5zktw4g+v7ZJLT+se/k+QfZ3DdL0vytzO1vils95eS3NS/Ni/e29vXQ4uhP0cl2Zzkh/0f8reTvC/JfnuwvrEklWThnrSrqvarqk17so6BNp2d5L4k3+t//jXJe5I8dmB7n62qJ4y4rg9OVq+qTqiqS2ag7bvsz6r6UFU9f0/XPQ3nAO/pX5tPjJ8508fSRJKsTPLFJN9NcluSzyQZ6+dN+lprZhj6c9sLq2o/4KnA04E/mK2G7OmbxW6Wv6yqHgUcCPwG8Bjg2pkOg3Tm69/D44ANk9R5UI+lJEcC7wdeB+wPLAP+J/DAQLW98lq3br4e5E2pqluBTwI/C5DkkCRrkmxPsjHJ7+2sm+ToJOv73ta3k7yzn/V/+9939j2+Z/b1/2OSG5LckeSqJI8bWFcleVWSm4CbBsqO7B/vn+T9SbYl+UaSP9gZrP3QyeeSnJ9kO3D2JM/xvqraAJwCbKMLD5I8L8mWgTa9McmtfW/xxiTHJTkeeBNwSv/c/rmv+/dJ3pbkc8APgCP6st8d2HSSvDvJXUm+muS4gRmbk/zKwPTgp4ld9uf44aIkz0qyrl/3uiTPGpj390n+qN9H30vyt0kOmmj/JPm9/rXe3r/2h/TlXwOOAP53345Fk+znqRxLZye5IsllfRuvS/ILE6z6KcDXq+oz1fleVX20qm4e0oahr7VmhqE/DyQ5DPh14Pq+6CPAFuAQ4CTg7QNh9S7gXVX108DPAJf35c/tfx/QDwN8Pt3475uAlwBLgM/26x70YuAYYPmQpr2brld3BPDLwG8DrxiYfwywCTgYeNsoz7Wq7gf+GnjO+HlJngCcCTy97zH+GrC5qj4FvJ2uJ7lfVQ0G028BZwCPAr4xZJM723gQ8BbgY0kOHKGpu+zPcW09ELgS+B/AYuCdwJVJFg9Ueynd/joY2Ad4/bANJTkWeAfwm8Bj++dxKUBV/QxwM31Pvqp+tLtGT/FYAlgJ/BVd7/zDwCeSPGzIqq8Dnti/yf/7jDB8tLvXWtNn6M9tn0hyJ/CPwD/Q/UEeBjwbeGNV3VNVXwTeSxduAPcBRyY5qKrurqov7Gb9rwTeUVU3VNUOuuB8ymBvv5+/vap+OLhgkgV0PbWz+l7dZuC/D7QDYGtVvbuqdoxffhJb6UJmvPuBRcDyJA+rqs1V9bVJ1vWXVbWhb8N9Q+Z/B/izvvd5GXAjcOIU2jqRE4GbquoD/bY/AnwVeOFAnfdV1b/2++Zyut7yMC8DLq6q6/pQPwt4Zvrx8hFN51gCuLaqruj33TuBhwPPGL/y/jzP84Cl/XO5LclfjhD+E73WmiZDf257cVUdUFWPq6r/3IfDIcD2qvreQL1v0P2xAZwOPB74aj+k8ILdrP9xwLuS3NkHwnYgA+sCuGWCZQ+i650O9p4H27G7ZSeztG/LT6iqjcBr6IaKvpPk0p3DHLsxWRturZ+8K+E36PbxnjqEXT9ZjN8/3xp4/ANgooD8iXVV1d3A7ePWNZnpHEswsP+q6gF+/KlgF1X1har6zapaQtd7fy7w3yZp19DXWtNn6M8/W4EDkzxqoOxw4FaAqrqpqk6lGzI4D7giyb7AsNut3gK8sg+DnT+PqKp/Gqgz0W1ab6P7VDH4qeDf2jHJshPqzwm8kG6oaRdV9eGqena/3aJ7jrvb1mRtWJokA9OH0+1jgO8DjxyY95gprHcrP7lvdq771iF1J/MT6+pfz8XTXNf49U54LPUOG9juTwGH8uP9M6GqWgd8jP7cwTCTvdaaHkN/nqmqW4B/At6R5OFJfp6ud/8hgCQvT7Kk75Xd2S92P90Jswfoxt93+gvgrCRP7pfdP8nJI7bjfrqP8W9L8qh+SOi1wKSXTQ6T5GFJnkQ3xvwYuqGE8XWekOTY/mTlPcAP++cG8G1gLFO/Qudg4NX99k8GngSs7ed9EVjVz1tBN+a907D9OWgt8PgkL02yMMkpdOdF/maK7YNuLP0VSZ7SP/e3A9f0Q2rTNtmx1Htakpeku/rqNcCPgF2GDJM8uz/ZfHA//UTgRRPUnfS11vQZ+vPTqcAYXY/r48BbqurT/bzjgQ1J7qY7qbuqH6/9Ad3J1M/1wznPqKqP0/WUL03yXeArwAlTaMd/oesNb6IbK/4wcPEUn8spfVvvBNbQDVs8raqG9SYXAefSfcr4Fl1gv6mf91f979uTXDeF7V8DHNWv823ASVV1ez/vD+lOht8BvJXu+QEwbH8OrrRfxwvorky5HXgD8IKqum0Kbdu5rs/0bfko8M2+Taumup4J7O5Ygu5E6yl0++C3gJdMcG7kTrqQ/3L/en6qX9+fDNSZymutaYr/REXSdCQ5Gziyql4+223R6OzpS1JDDH1JaojDO5LUEHv6ktSQPbpJ1oPhoIMOqrGxsdluhiTNKddee+1t/RffdushF/pjY2OsX79+tpshSXNKkmH3jtqFwzuS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQh9w3ch8sY6uvHFq++dyZ+B/XkjQ3zLvQnyjcJUkO70hSUwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyEihn+T4JDcm2Zhk9ZD5z01yXZIdSU4aN++0JDf1P6fNVMMlSVM3aegnWQBcAJwALAdOTbJ8XLWbgd8BPjxu2QOBtwDHAEcDb0ny6D1vtiRpOkbp6R8NbKyqTVV1L3ApsHKwQlVtrqovAQ+MW/bXgE9X1faqugP4NHD8DLRbkjQNo4T+UuCWgektfdkoRlo2yRlJ1idZv23bthFXLUmaqlFCP0PKasT1j7RsVV1YVSuqasWSJUtGXLUkaapGCf0twGED04cCW0dc/54sK0maYaOE/jrgqCTLkuwDrALWjLj+q4DnJ3l0fwL3+X2ZJGkWTBr6VbUDOJMurG8ALq+qDUnOSfIigCRPT7IFOBn4X0k29MtuB/6I7o1jHXBOXyZJmgUj/bvEqloLrB1X9uaBx+vohm6GLXsxcPEetFGSNEP8Rq4kNcTQl6SGGPqS1BBDX5IaYuhLUkNGunpnPhtbfeXQ8s3nnriXWyJJDz57+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1ZKTQT3J8khuTbEyyesj8RUku6+dfk2SsL39YkkuSfDnJDUnOmtnmS5KmYtLQT7IAuAA4AVgOnJpk+bhqpwN3VNWRwPnAeX35ycCiqvo54GnAK3e+IUiS9r5RevpHAxuralNV3QtcCqwcV2clcEn/+ArguCQBCtg3yULgEcC9wHdnpOWSpCkbJfSXArcMTG/py4bWqaodwF3AYro3gO8D3wRuBv60qrbvYZslSdM0SuhnSFmNWOdo4H7gEGAZ8LokR+yygeSMJOuTrN+2bdsITZIkTccoob8FOGxg+lBg60R1+qGc/YHtwEuBT1XVfVX1HeBzwIrxG6iqC6tqRVWtWLJkydSfhSRpJKOE/jrgqCTLkuwDrALWjKuzBjitf3wScHVVFd2QzrHp7As8A/jqzDRdkjRVk4Z+P0Z/JnAVcANweVVtSHJOkhf11S4CFifZCLwW2HlZ5wXAfsBX6N483ldVX5rh5yBJGtHCUSpV1Vpg7biyNw88vofu8szxy909rFySNDv8Rq4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDRrrLZovGVl85tHzzuSfu5ZZI0syxpy9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JashIoZ/k+CQ3JtmYZPWQ+YuSXNbPvybJ2MC8n0/y+SQbknw5ycNnrvmSpKmYNPSTLAAuAE4AlgOnJlk+rtrpwB1VdSRwPnBev+xC4IPAf6qqJwPPA+6bsdZLkqZklJ7+0cDGqtpUVfcClwIrx9VZCVzSP74COC5JgOcDX6qqfwaoqtur6v6ZabokaapGCf2lwC0D01v6sqF1qmoHcBewGHg8UEmuSnJdkjcM20CSM5KsT7J+27ZtU30OkqQRjRL6GVJWI9ZZCDwbeFn/+zeSHLdLxaoLq2pFVa1YsmTJCE2SJE3HKKG/BThsYPpQYOtEdfpx/P2B7X35P1TVbVX1A2At8NQ9bbQkaXoWjlBnHXBUkmXArcAq4KXj6qwBTgM+D5wEXF1VleQq4A1JHgncC/wy3YneOWts9ZVDyzefe+JebokkTd2koV9VO5KcCVwFLAAurqoNSc4B1lfVGuAi4ANJNtL18Ff1y96R5J10bxwFrK2q4akpSXrQjdLTp6rW0g3NDJa9eeDxPcDJEyz7QbrLNiVJs8xv5EpSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNWSk0E9yfJIbk2xMsnrI/EVJLuvnX5NkbNz8w5PcneT1M9NsSdJ0TBr6SRYAFwAnAMuBU5MsH1ftdOCOqjoSOB84b9z884FP7nlzJUl7YuEIdY4GNlbVJoAklwIrgX8ZqLMSOLt/fAXwniSpqkryYmAT8P0Za/VD0NjqK4eWbz73xL3cEkma2CjDO0uBWwamt/RlQ+tU1Q7gLmBxkn2BNwJv3d0GkpyRZH2S9du2bRu17ZKkKRol9DOkrEas81bg/Kq6e3cbqKoLq2pFVa1YsmTJCE2SJE3HKMM7W4DDBqYPBbZOUGdLkoXA/sB24BjgpCR/AhwAPJDknqp6zx63XJI0ZaOE/jrgqCTLgFuBVcBLx9VZA5wGfB44Cbi6qgp4zs4KSc4G7jbwJWn2TBr6VbUjyZnAVcAC4OKq2pDkHGB9Va0BLgI+kGQjXQ9/1YPZaEnS9IzS06eq1gJrx5W9eeDxPcDJk6zj7Gm0T5I0g/xGriQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkNGusumps//nSvpocSeviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0ZKfSTHJ/kxiQbk6weMn9Rksv6+dckGevLfzXJtUm+3P8+dmabL0maiklDP8kC4ALgBGA5cGqS5eOqnQ7cUVVHAucD5/XltwEvrKqfA04DPjBTDZckTd0oPf2jgY1Vtamq7gUuBVaOq7MSuKR/fAVwXJJU1fVVtbUv3wA8PMmimWi4JGnqRvl3iUuBWwamtwDHTFSnqnYkuQtYTNfT3+k/ANdX1Y+m39z5w3+jKGk2jBL6GVJWU6mT5Ml0Qz7PH7qB5AzgDIDDDz98hCZJkqZjlOGdLcBhA9OHAlsnqpNkIbA/sL2fPhT4OPDbVfW1YRuoqgurakVVrViyZMnUnoEkaWSjhP464Kgky5LsA6wC1oyrs4buRC3AScDVVVVJDgCuBM6qqs/NVKMlSdMzaehX1Q7gTOAq4Abg8qrakOScJC/qq10ELE6yEXgtsPOyzjOBI4E/TPLF/ufgGX8WkqSRjDKmT1WtBdaOK3vzwON7gJOHLPfHwB/vYRslSTPEb+RKUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQka7T197jjdgkPZjs6UtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkP8Ru4c4Td1Jc0Ee/qS1BBDX5IaYuhLUkMMfUlqiCdy57iJTvCCJ3kl7cqeviQ1xNCXpIYY+pLUEMf057HdjfcP4zkAaf6zpy9JDTH0JakhDu/o33h/H2n+Gyn0kxwPvAtYALy3qs4dN38R8H7gacDtwClVtbmfdxZwOnA/8OqqumrGWq+9wnMD0vwxaegnWQBcAPwqsAVYl2RNVf3LQLXTgTuq6sgkq4DzgFOSLAdWAU8GDgH+T5LHV9X9M/1E9NDhJwbpoWuUnv7RwMaq2gSQ5FJgJTAY+iuBs/vHVwDvSZK+/NKq+hHw9SQb+/V9fmaar7lkqp8YZsru3mxm6g3KNzrNFaOE/lLgloHpLcAxE9Wpqh1J7gIW9+VfGLfs0vEbSHIGcEY/eXeSG0dq/Y8dBNw2xWVa4H4Bct4uRZPulyHLzNS2H6o8VoabS/vlcaNUGiX0M6SsRqwzyrJU1YXAhSO0Zagk66tqxXSXn6/cL8O5X3blPhluPu6XUS7Z3AIcNjB9KLB1ojpJFgL7A9tHXFaStJeMEvrrgKOSLEuyD92J2TXj6qwBTusfnwRcXVXVl69KsijJMuAo4P/NTNMlSVM16fBOP0Z/JnAV3SWbF1fVhiTnAOurag1wEfCB/kTtdro3Bvp6l9Od9N0BvOpBunJn2kND85z7ZTj3y67cJ8PNu/2SrkMuSWqBt2GQpIYY+pLUkDkf+kmOT3Jjko1JVs92e2ZDksOS/F2SG5JsSPL7ffmBST6d5Kb+96Nnu62zIcmCJNcn+Zt+elmSa/r9cll/gUJTkhyQ5IokX+2Pm2e2frwk+a/9389XknwkycPn47Eyp0N/4BYRJwDLgVP7Wz+0Zgfwuqp6EvAM4FX9flgNfKaqjgI+00+36PeBGwamzwPO7/fLHXS3EWnNu4BPVdUTgV+g2z/NHi9JlgKvBlZU1c/SXbSy85Yy8+pYmdOhz8AtIqrqXmDnLSKaUlXfrKrr+sffo/sDXkq3Ly7pq10CvHh2Wjh7khwKnAi8t58OcCzd7UKgwf2S5KeB59JddUdV3VtVd+LxshB4RP9do0cC32QeHitzPfSH3SJil9s8tCTJGPCLwDXAv6uqb0L3xgAcPHstmzV/BrwBeKCfXgzcWVU7+ukWj5kjgG3A+/phr/cm2ZeGj5equhX4U+BmurC/C7iWeXiszPXQH+k2D61Ish/wUeA1VfXd2W7PbEvyAuA7VXXtYPGQqq0dMwuBpwJ/XlW/CHyfhoZyhunPX6wEltHdEXhfumHj8eb8sTLXQ9/bPPSSPIwu8D9UVR/ri7+d5LH9/McC35mt9s2SXwJelGQz3dDfsXQ9/wP6j/DQ5jGzBdhSVdf001fQvQm0fLz8CvD1qtpWVfcBHwOexTw8VuZ66I9yi4h5rx+nvgi4oareOTBr8PYYpwF/vbfbNpuq6qyqOrSqxuiOjaur6mXA39HdLgTa3C/fAm5J8oS+6Di6b823fLzcDDwjySP7v6ed+2TeHStz/hu5SX6drve28xYRb5vlJu11SZ4NfBb4Mj8eu34T3bj+5cDhdAf1yVW1fVYaOcuSPA94fVW9IMkRdD3/A4HrgZf3//OhGUmeQndyex9gE/AKuk5gs8dLkrcCp9BdDXc98Lt0Y/jz6liZ86EvSRrdXB/ekSRNgaEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGvL/ARuZcyj6Jn9hAAAAAElFTkSuQmCC\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": 358,
"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": 359,
"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.2047\n",
"B 8 10 0.1160\n",
"C -3 16 0.1010\n",
"D 7 11 0.1075\n",
"E -1 9 0.0766\n",
"F 1 11 0.0854\n",
"G 18 10 0.1761\n",
"H 12 18 0.1327\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)\n",
"print(\"Only thing that worries me is that Gelman has A best with prob=10%\")"
]
},
{
"cell_type": "code",
"execution_count": 360,
"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.59 0.62 0.60 0.66 0.64 0.52 0.58 \n",
"B 0.41 1.00 0.54 0.51 0.57 0.56 0.42 0.49 \n",
"C 0.38 0.46 1.00 0.47 0.54 0.51 0.39 0.45 \n",
"D 0.40 0.49 0.53 1.00 0.56 0.54 0.41 0.48 \n",
"E 0.34 0.43 0.46 0.44 1.00 0.48 0.35 0.41 \n",
"F 0.36 0.44 0.49 0.46 0.52 1.00 0.38 0.43 \n",
"G 0.48 0.58 0.61 0.59 0.65 0.62 1.00 0.57 \n",
"H 0.42 0.51 0.55 0.52 0.59 0.57 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": "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<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],se[i]);\n",
" }\n",
"}\n",
"'''\n",
"sm2=pystan.StanModel(model_code=stan_code_2)"
]
},
{
"cell_type": "code",
"execution_count": 341,
"metadata": {},
"outputs": [],
"source": [
"answers=sm2.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,
"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.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}