From 54ba9a1a66e2acb7c8f7b21d526056e04f1e9457 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Sat, 21 Apr 2018 16:40:37 -0400 Subject: [PATCH] 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,