Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
finished off (sort of) 5.9.3
  • Loading branch information
jet08013 committed Apr 21, 2018
1 parent 93bb64e commit 54ba9a1
Showing 1 changed file with 190 additions and 11 deletions.
201 changes: 190 additions & 11 deletions BDA 5.9.3.ipynb
Expand Up @@ -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",
Expand Down Expand Up @@ -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<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",
"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,
Expand Down

0 comments on commit 54ba9a1

Please sign in to comment.