Skip to content

Instantly share code, notes, and snippets.

@thorstenwagner

thorstenwagner/catreg.ipynb Secret

Created July 15, 2021 14:39
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save thorstenwagner/4b69a797b8c1508eb2ee07dc9763f9a4 to your computer and use it in GitHub Desktop.
Save thorstenwagner/4b69a797b8c1508eb2ee07dc9763f9a4 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
catreg.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [],
"source": [
"library(rethinking)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Modelling Problem"
]
},
{
"cell_type": "code",
"execution_count": 192,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] 0.25 0.25 0.25 0.25\n"
]
}
],
"source": [
"# simulate observation of 500 category pairs\n",
"N <- 1000\n",
"prob_cat_b <- c(1,1,1,1)\n",
"pa <- softmax(prob_cat_b[1], prob_cat_b[2], prob_cat_b[3], prob_cat_b[4])\n",
"print(pa)"
]
},
{
"cell_type": "code",
"execution_count": 193,
"metadata": {},
"outputs": [],
"source": [
"cat_pos_b <- rep(NA,N)\n",
"set.seed(34302)\n",
"for ( i in 1:N ) cat_pos_b[i] <- sample( 1:4 , size = 1, prob=pa )"
]
},
{
"cell_type": "code",
"execution_count": 194,
"metadata": {},
"outputs": [],
"source": [
"cat_pos_a <- rep(NA,N)\n",
"for ( i in 1:N ) {\n",
" if (cat_pos_b[i] == 1 ) {\n",
" prob_cat_a <- c(1,1,1,5)\n",
" }\n",
" else if ( cat_pos_b[i] == 4 ) {\n",
" prob_cat_a <- c(5,1,1,1)\n",
" } \n",
" else {\n",
" prob_cat_a <- c(1,1,1,1)\n",
" }\n",
" \n",
" pb <- softmax(prob_cat_a[1], prob_cat_a[2], prob_cat_a[3], prob_cat_a[4])\n",
" cat_pos_a[i] <- sample( 1:4 , size = 1, prob=pb ) # I should rename pb/pa\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 195,
"metadata": {},
"outputs": [],
"source": [
"code <- \"\n",
"data{\n",
" int N; // number of amino acid rows\n",
" int K; // number of possible categories\n",
" int cat_pos_a[N]; // outcome\n",
" int posb_1[N]; // predictor indicator for class 1\n",
" int posb_2[N]; // predictor indicator for class 2\n",
" int posb_3[N]; // predictor indicator for class 3\n",
" int posb_4[N]; // predictor indicator for class 4\n",
"}\n",
"\n",
"parameters {\n",
" vector[K-1] a; // intercepts\n",
" vector[K-1] b1;\n",
" vector[K-1] b2;\n",
" vector[K-1] b3;\n",
" vector[K-1] b4;\n",
"\n",
"}\n",
"\n",
"model {\n",
" vector[K] p[N]; // Soft max probabilties for each amino acid row using the regressed values\n",
" vector[K] s[N]; // Regressed values\n",
" a ~ normal(0, 1);\n",
" b1 ~ normal(0, 5);\n",
" b2 ~ normal(0, 5);\n",
" b3 ~ normal(0, 5);\n",
" b4 ~ normal(0, 5);\n",
" \n",
" for(n in 1:N) {\n",
" s[n][1] = a[1] + b1[1]*posb_1[n] + b2[1]*posb_2[n] + b3[1]*posb_3[n] + b4[1]*posb_4[n];\n",
" s[n][2] = a[2] + b1[2]*posb_1[n] + b2[2]*posb_2[n] + b3[2]*posb_3[n] + b4[2]*posb_4[n];\n",
" s[n][3] = a[3] + b1[3]*posb_1[n] + b2[3]*posb_2[n] + b3[3]*posb_3[n] + b4[3]*posb_4[n];\n",
" s[n][4] = 0; // Pivot\n",
" p[n] = softmax( s[n] );\n",
" cat_pos_a[n] ~ categorical( p[n] );\n",
" }\n",
"}\n",
"\""
]
},
{
"cell_type": "code",
"execution_count": 196,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 1).\n",
"Chain 1: \n",
"Chain 1: Gradient evaluation took 0.001035 seconds\n",
"Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.35 seconds.\n",
"Chain 1: Adjust your expectations accordingly!\n",
"Chain 1: \n",
"Chain 1: \n",
"Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)\n",
"Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)\n",
"Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)\n",
"Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)\n",
"Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)\n",
"Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)\n",
"Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)\n",
"Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)\n",
"Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)\n",
"Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)\n",
"Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)\n",
"Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)\n",
"Chain 1: \n",
"Chain 1: Elapsed Time: 47.9183 seconds (Warm-up)\n",
"Chain 1: 49.8576 seconds (Sampling)\n",
"Chain 1: 97.7759 seconds (Total)\n",
"Chain 1: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 2).\n",
"Chain 2: \n",
"Chain 2: Gradient evaluation took 0.000654 seconds\n",
"Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 6.54 seconds.\n",
"Chain 2: Adjust your expectations accordingly!\n",
"Chain 2: \n",
"Chain 2: \n",
"Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)\n",
"Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)\n",
"Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)\n",
"Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)\n",
"Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)\n",
"Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)\n",
"Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)\n",
"Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)\n",
"Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)\n",
"Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)\n",
"Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)\n",
"Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)\n",
"Chain 2: \n",
"Chain 2: Elapsed Time: 50.8036 seconds (Warm-up)\n",
"Chain 2: 55.7507 seconds (Sampling)\n",
"Chain 2: 106.554 seconds (Total)\n",
"Chain 2: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 3).\n",
"Chain 3: \n",
"Chain 3: Gradient evaluation took 0.000678 seconds\n",
"Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.78 seconds.\n",
"Chain 3: Adjust your expectations accordingly!\n",
"Chain 3: \n",
"Chain 3: \n",
"Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)\n",
"Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)\n",
"Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)\n",
"Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)\n",
"Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)\n",
"Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)\n",
"Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)\n",
"Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)\n",
"Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)\n",
"Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)\n",
"Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)\n",
"Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)\n",
"Chain 3: \n",
"Chain 3: Elapsed Time: 48.8396 seconds (Warm-up)\n",
"Chain 3: 53.1494 seconds (Sampling)\n",
"Chain 3: 101.989 seconds (Total)\n",
"Chain 3: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 4).\n",
"Chain 4: \n",
"Chain 4: Gradient evaluation took 0.001245 seconds\n",
"Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 12.45 seconds.\n",
"Chain 4: Adjust your expectations accordingly!\n",
"Chain 4: \n",
"Chain 4: \n",
"Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)\n",
"Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)\n",
"Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)\n",
"Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)\n",
"Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)\n",
"Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)\n",
"Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)\n",
"Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)\n",
"Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)\n",
"Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)\n",
"Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)\n",
"Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)\n",
"Chain 4: \n",
"Chain 4: Elapsed Time: 49.3037 seconds (Warm-up)\n",
"Chain 4: 50.5946 seconds (Sampling)\n",
"Chain 4: 99.8983 seconds (Total)\n",
"Chain 4: \n"
]
}
],
"source": [
"dat_list <- list(N=N, \n",
" K=4, \n",
" cat_pos_a=cat_pos_a, \n",
" posb_1=ifelse ( cat_pos_b==1, 1, 0), \n",
" posb_2=ifelse ( cat_pos_b==2, 1, 0),\n",
" posb_3=ifelse ( cat_pos_b==3, 1, 0),\n",
" posb_4=ifelse ( cat_pos_b==4, 1, 0)\n",
" )\n",
"erics_model <- stan(model_code=code, data=dat_list, chains=4)"
]
},
{
"cell_type": "code",
"execution_count": 197,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<thead><tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th><th scope=col>n_eff</th><th scope=col>Rhat4</th></tr></thead>\n",
"<tbody>\n",
"\t<tr><th scope=row>a[1]</th><td>-0.12322060</td><td>0.9287056 </td><td> -1.590157 </td><td> 1.348636 </td><td>1824.087 </td><td>1.0009527 </td></tr>\n",
"\t<tr><th scope=row>a[2]</th><td>-0.13923145</td><td>0.9480458 </td><td> -1.666414 </td><td> 1.395261 </td><td>2047.836 </td><td>0.9993716 </td></tr>\n",
"\t<tr><th scope=row>a[3]</th><td>-0.16074493</td><td>0.9560072 </td><td> -1.702290 </td><td> 1.378757 </td><td>1830.529 </td><td>1.0004890 </td></tr>\n",
"\t<tr><th scope=row>b1[1]</th><td>-7.99609449</td><td>2.4393442 </td><td>-12.341984 </td><td>-4.661093 </td><td>1985.239 </td><td>1.0001452 </td></tr>\n",
"\t<tr><th scope=row>b1[2]</th><td>-3.29114388</td><td>1.0151762 </td><td> -4.901370 </td><td>-1.652014 </td><td>2037.876 </td><td>0.9993963 </td></tr>\n",
"\t<tr><th scope=row>b1[3]</th><td>-4.33610984</td><td>1.1283548 </td><td> -6.185982 </td><td>-2.560401 </td><td>1980.587 </td><td>1.0008439 </td></tr>\n",
"\t<tr><th scope=row>b2[1]</th><td> 0.30409253</td><td>0.9427030 </td><td> -1.209730 </td><td> 1.823396 </td><td>1909.280 </td><td>1.0007188 </td></tr>\n",
"\t<tr><th scope=row>b2[2]</th><td>-0.12380315</td><td>0.9673072 </td><td> -1.690192 </td><td> 1.426593 </td><td>2066.851 </td><td>0.9993875 </td></tr>\n",
"\t<tr><th scope=row>b2[3]</th><td>-0.02505068</td><td>0.9716536 </td><td> -1.554056 </td><td> 1.545895 </td><td>1888.990 </td><td>1.0003950 </td></tr>\n",
"\t<tr><th scope=row>b3[1]</th><td>-0.04704746</td><td>0.9457446 </td><td> -1.553789 </td><td> 1.436726 </td><td>1918.053 </td><td>1.0009478 </td></tr>\n",
"\t<tr><th scope=row>b3[2]</th><td>-0.06128586</td><td>0.9628079 </td><td> -1.637718 </td><td> 1.469471 </td><td>2089.271 </td><td>0.9996764 </td></tr>\n",
"\t<tr><th scope=row>b3[3]</th><td> 0.01051474</td><td>0.9678124 </td><td> -1.533536 </td><td> 1.576381 </td><td>1863.273 </td><td>1.0004232 </td></tr>\n",
"\t<tr><th scope=row>b4[1]</th><td> 4.10625345</td><td>1.0283142 </td><td> 2.441152 </td><td> 5.723109 </td><td>2026.396 </td><td>1.0003478 </td></tr>\n",
"\t<tr><th scope=row>b4[2]</th><td>-0.16850855</td><td>1.1553958 </td><td> -2.014304 </td><td> 1.695935 </td><td>2368.305 </td><td>0.9995455 </td></tr>\n",
"\t<tr><th scope=row>b4[3]</th><td>-0.13295372</td><td>1.1783590 </td><td> -2.000113 </td><td> 1.739428 </td><td>2151.816 </td><td>1.0008620 </td></tr>\n",
"</tbody>\n",
"</table>\n"
],
"text/latex": [
"\\begin{tabular}{r|llllll}\n",
" & mean & sd & 5.5\\% & 94.5\\% & n\\_eff & Rhat4\\\\\n",
"\\hline\n",
"\ta{[}1{]} & -0.12322060 & 0.9287056 & -1.590157 & 1.348636 & 1824.087 & 1.0009527 \\\\\n",
"\ta{[}2{]} & -0.13923145 & 0.9480458 & -1.666414 & 1.395261 & 2047.836 & 0.9993716 \\\\\n",
"\ta{[}3{]} & -0.16074493 & 0.9560072 & -1.702290 & 1.378757 & 1830.529 & 1.0004890 \\\\\n",
"\tb1{[}1{]} & -7.99609449 & 2.4393442 & -12.341984 & -4.661093 & 1985.239 & 1.0001452 \\\\\n",
"\tb1{[}2{]} & -3.29114388 & 1.0151762 & -4.901370 & -1.652014 & 2037.876 & 0.9993963 \\\\\n",
"\tb1{[}3{]} & -4.33610984 & 1.1283548 & -6.185982 & -2.560401 & 1980.587 & 1.0008439 \\\\\n",
"\tb2{[}1{]} & 0.30409253 & 0.9427030 & -1.209730 & 1.823396 & 1909.280 & 1.0007188 \\\\\n",
"\tb2{[}2{]} & -0.12380315 & 0.9673072 & -1.690192 & 1.426593 & 2066.851 & 0.9993875 \\\\\n",
"\tb2{[}3{]} & -0.02505068 & 0.9716536 & -1.554056 & 1.545895 & 1888.990 & 1.0003950 \\\\\n",
"\tb3{[}1{]} & -0.04704746 & 0.9457446 & -1.553789 & 1.436726 & 1918.053 & 1.0009478 \\\\\n",
"\tb3{[}2{]} & -0.06128586 & 0.9628079 & -1.637718 & 1.469471 & 2089.271 & 0.9996764 \\\\\n",
"\tb3{[}3{]} & 0.01051474 & 0.9678124 & -1.533536 & 1.576381 & 1863.273 & 1.0004232 \\\\\n",
"\tb4{[}1{]} & 4.10625345 & 1.0283142 & 2.441152 & 5.723109 & 2026.396 & 1.0003478 \\\\\n",
"\tb4{[}2{]} & -0.16850855 & 1.1553958 & -2.014304 & 1.695935 & 2368.305 & 0.9995455 \\\\\n",
"\tb4{[}3{]} & -0.13295372 & 1.1783590 & -2.000113 & 1.739428 & 2151.816 & 1.0008620 \\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"| <!--/--> | mean | sd | 5.5% | 94.5% | n_eff | Rhat4 |\n",
"|---|---|---|---|---|---|---|\n",
"| a[1] | -0.12322060 | 0.9287056 | -1.590157 | 1.348636 | 1824.087 | 1.0009527 |\n",
"| a[2] | -0.13923145 | 0.9480458 | -1.666414 | 1.395261 | 2047.836 | 0.9993716 |\n",
"| a[3] | -0.16074493 | 0.9560072 | -1.702290 | 1.378757 | 1830.529 | 1.0004890 |\n",
"| b1[1] | -7.99609449 | 2.4393442 | -12.341984 | -4.661093 | 1985.239 | 1.0001452 |\n",
"| b1[2] | -3.29114388 | 1.0151762 | -4.901370 | -1.652014 | 2037.876 | 0.9993963 |\n",
"| b1[3] | -4.33610984 | 1.1283548 | -6.185982 | -2.560401 | 1980.587 | 1.0008439 |\n",
"| b2[1] | 0.30409253 | 0.9427030 | -1.209730 | 1.823396 | 1909.280 | 1.0007188 |\n",
"| b2[2] | -0.12380315 | 0.9673072 | -1.690192 | 1.426593 | 2066.851 | 0.9993875 |\n",
"| b2[3] | -0.02505068 | 0.9716536 | -1.554056 | 1.545895 | 1888.990 | 1.0003950 |\n",
"| b3[1] | -0.04704746 | 0.9457446 | -1.553789 | 1.436726 | 1918.053 | 1.0009478 |\n",
"| b3[2] | -0.06128586 | 0.9628079 | -1.637718 | 1.469471 | 2089.271 | 0.9996764 |\n",
"| b3[3] | 0.01051474 | 0.9678124 | -1.533536 | 1.576381 | 1863.273 | 1.0004232 |\n",
"| b4[1] | 4.10625345 | 1.0283142 | 2.441152 | 5.723109 | 2026.396 | 1.0003478 |\n",
"| b4[2] | -0.16850855 | 1.1553958 | -2.014304 | 1.695935 | 2368.305 | 0.9995455 |\n",
"| b4[3] | -0.13295372 | 1.1783590 | -2.000113 | 1.739428 | 2151.816 | 1.0008620 |\n",
"\n"
],
"text/plain": [
" mean sd 5.5% 94.5% n_eff Rhat4 \n",
"a[1] -0.12322060 0.9287056 -1.590157 1.348636 1824.087 1.0009527\n",
"a[2] -0.13923145 0.9480458 -1.666414 1.395261 2047.836 0.9993716\n",
"a[3] -0.16074493 0.9560072 -1.702290 1.378757 1830.529 1.0004890\n",
"b1[1] -7.99609449 2.4393442 -12.341984 -4.661093 1985.239 1.0001452\n",
"b1[2] -3.29114388 1.0151762 -4.901370 -1.652014 2037.876 0.9993963\n",
"b1[3] -4.33610984 1.1283548 -6.185982 -2.560401 1980.587 1.0008439\n",
"b2[1] 0.30409253 0.9427030 -1.209730 1.823396 1909.280 1.0007188\n",
"b2[2] -0.12380315 0.9673072 -1.690192 1.426593 2066.851 0.9993875\n",
"b2[3] -0.02505068 0.9716536 -1.554056 1.545895 1888.990 1.0003950\n",
"b3[1] -0.04704746 0.9457446 -1.553789 1.436726 1918.053 1.0009478\n",
"b3[2] -0.06128586 0.9628079 -1.637718 1.469471 2089.271 0.9996764\n",
"b3[3] 0.01051474 0.9678124 -1.533536 1.576381 1863.273 1.0004232\n",
"b4[1] 4.10625345 1.0283142 2.441152 5.723109 2026.396 1.0003478\n",
"b4[2] -0.16850855 1.1553958 -2.014304 1.695935 2368.305 0.9995455\n",
"b4[3] -0.13295372 1.1783590 -2.000113 1.739428 2151.816 1.0008620"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"precis(erics_model,depth=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Making predictions"
]
},
{
"cell_type": "code",
"execution_count": 198,
"metadata": {},
"outputs": [],
"source": [
"post <- extract.samples(erics_model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case A: I want to predict the probabilties for position a, given I observed a 2 for position b. I expect a uniform distribution."
]
},
{
"cell_type": "code",
"execution_count": 199,
"metadata": {},
"outputs": [],
"source": [
"s1 = post$a[,1]+post$b2[,1]\n",
"s2 = post$a[,2]+post$b2[,2]\n",
"s3 = post$a[,3]+post$b2[,3]\n",
"s4= 0"
]
},
{
"cell_type": "code",
"execution_count": 200,
"metadata": {},
"outputs": [],
"source": [
"simulations <- matrix(nrow = 0, ncol = 4) \n",
"for (i in 1:length(s1)){\n",
" s = c(s1[i],s2[i],s3[i],s4);\n",
" u = softmax(s);\n",
" simulations = rbind(simulations, u)\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 201,
"metadata": {},
"outputs": [],
"source": [
"require(datasets)"
]
},
{
"cell_type": "code",
"execution_count": 202,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] 0.3150507 0.2028150 0.2189107 0.2632235\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAA/1BMVEUAAAABAQECAgIHBwcK\nCgoNDQ0PDw8REREWFhYZGRkaGhonJycwMDAzMzM1NTU5OTk+Pj4/Pz9BQUFCQkJERERFRUVI\nSEhKSkpMTExNTU1PT09QUFBTU1NUVFRVVVVYWFhcXFxeXl5fX19gYGBiYmJjY2NmZmZpaWlr\na2tubm5ycnJ3d3eBgYGDg4OEhISIiIiMjIyPj4+ZmZmenp6lpaWqqqqrq6uvr6+ysrK4uLi+\nvr7FxcXMzMzOzs7S0tLb29vc3Nzd3d3f39/g4ODh4eHl5eXn5+fs7Ozt7e3u7u7x8fHz8/P0\n9PT19fX39/f4+Pj5+fn6+vr7+/v8/Pz///9R8t0KAAAACXBIWXMAABJ0AAASdAHeZh94AAAU\nF0lEQVR4nO3ca5MkaVmA4UJFPIGg7gKCrIKgqHgEj6jgAV1FgWX//2+RzZmpTqO6aioz7zc7\nt/u6PsxmxORTkd0T93b2U4fTh8Bmp6e+AHgOhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFB\nQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEh\nQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQB\nIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQE\nASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSE\nBAEhQUBIEBASBIQEASFBQEgQeMqQ3v/ucfzgCb8PPANPGdI3TsfxzhN+H3gGnjKkr73790fx\n1c894feBZ0BIEyGxjZAmQmIbIU2ExDZCmgiJbYQ0ERLbCGkiJLYR0kRIbCOkiZDYRkgTIbGN\nkCZCYhshTYTENptCev2Cz7XjQuLZ2BDS7LXT6x5ASDwb60N66GdtSULi2dgS0mOHSwiJZ0NI\nEyGxjVu7iZDYxrJhIiS2sf6eCIltPCE7ERLbCGkiJLbZ9jvSqz/XPoaQeDY2be1+lpBlA3y4\n9XmkVxFZf/PibXxC9lVCQuKlE9JESGzjJUITIbHN5q3do78i/csvfersF3/+p1ceQEg8G1ue\nRzqHdPE3P/3b75x96/TjK/NC4tkY/4TsPwiJ509IEyGxjZAmQmKbJKSbWzsh8QIIaSIktnFr\nNxES2whpIiS2Gf8OWSHxAoz/zAYh8QKM/xQhIfECjH/RqpB4AYQ0ERLbuLWbCIltLBsmQmIb\n6++JkNjGE7ITIbGNkCZCYhshTYTENkKaCIlthDQREtsIaSIkthHSREhsI6SJkNhGSBMhsY2Q\nJkJiGyFNhMQ2QpoIiW2ENBES2whpIiS2EdJESGwjpImQ2EZIEyGxjZAmQmIbIU2ExDZCmgiJ\nbYQ0ERLbCGkiJLYR0kRIbCOkiZDYRkgTIbGNkCZCYhshTYTENkKaCIlthDQREtsIaSIkthHS\nREhsI6SJkNhGSBMhsY2QJkJiGyFNhMQ2QpoIiW2ENBES2whpIiS2EdJESGwjpImQ2EZIEyGx\njZAmQmIbIU2ExDZCmgiJbYQ0ERLbCGkiJLYR0kRIbCOkiZDYRkgTIbGNkCZCYhshTYTENkKa\nCIlthDQREtsIaSIkthHSREh7+ce/PI5/Dr8uIU2EtJff+OQvH8Uv/H74dQlpIqS9fO6rT/1v\nffbu18KvS0gTIe1FSGsJiRkhrSUkZoS0lpCYEdJaQmJGSGsJiRkhrSUkZoR0Zf4jN88QEjNC\nupg8ffi6o9slCYkZIV1MnqaOPnzz5zVCYkZIF5On1z+Vzv95nJCYEdLFpJBYTkgXk0JiOSFd\nTL7ZM7xl2yAkZoR0MTkL6dZ5QmJGSGsJiRkhrSUkZoT0yOybezq/I3EvIV2OnrcMQuJeQrqY\nfNgzCIl7Celi8uHVQULiXkK6mHz9w+gkJO4npIvJN78fnYTE3YR0OXou6SKk//y9r5z9ppB4\nIKTL0fPy+yKk9//wG2e/LSQeCGktt3bMCGktITEjpLWExIyQbj2IrR13EtKtBxESdxLSWkJi\nRkhrCYkZIT02/PZPtRMS/4+QLkcf3DpNSMwI6WLy3I8PP+FuQrqYPD12eElIzAjpYlJILCek\ni0m3diwnpMtRywYWE9Jjw9bfLCSktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpL\nSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTE\njJDWEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwI\naS0hMSOktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwIaS0hMSOktYTEjJDW\nEhIzQlpLSMwIaS0hMSOktYTEjJDWEhIzQlpLSMwIaS0h7eO/v38c/3P9MoW0lpD28Tun4/jq\n9csU0lpC2sfnv/Dto3jnK9cvU0hrCWkfn3/vqb+JZ18S0gBC2oeQlhJST0glIY0gpH0IaSkh\n9YRUEtIIQtqHkJYSUk9IJSGNIKR9CGkpIfWEVBLSCELah5CWElJPSCUhjSCkfQhpKSH1hFQS\n0ghC2oeQlhJST0glIY0gpH0IaSkh9YRUEtIIQtqHkJYSUk9IJSEtHX7l9klC2oeQljpKSLNP\njbl1mpD2IaSlDhLSQz+3SxLSPoS01GFCeuzwkpD2IaSlhNQTUklIiybd2h2KkJY6SEiWDcci\npKWOEpL196EIaanjhHQXIe1DSEsJqSekkpAWzr65q3vey4Y///Rx/On1yxTSUkcJ6eH3o+cd\n0ld+7Y+O4jOfv36ZQlrqICF9VM/rkp55SF966qs7e09IoQOF9LokIe1ESKUjhfSqJCHtREil\nQ4U0lXQR0gd/852zbwkpI6TSQUKalXQR0r/+yqfOPnn60ZUHENJSQiodLKQP3/LiBrd2HSGV\njhLSnYTUEVJJSD0hlYQ0gpA6QiodLyS/I+1ESCUh9YRUEtIIQuoIqSSknpBKQlo6/CLeISuk\npYS0bPSFfGaDkJYS0qLJl/IpQkJaSkiLJl/K59oJaSkhLZoU0u6EVDpMSG7t9iak0kFCsmzY\nn5BKRwnJ+nt3QiodJ6S7CKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFS\nSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgj\nCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIqCakn\npJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWE\nVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLS\nCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLq\nCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgd\nIZWE1BNSSUgr5j9y8wwhdYRUOkhIUz+n01tLElJHSKUDhXR6U9ON84TUEVLpUCGdD68SUkdI\nJSH1hFQS0qJJIe1OSKXDhPR6z/CWbYOQOkIqHTCkW+cJqSOk0kFCupeQOkIqCaknpJKQRhBS\nR0glIfWEVBLSosm5G+cJqSOkkpB6QioJadnofbNC6gip9HEI6Yfvn/2dkDJCKn0MQvrez83v\n+3505SwhLSWk0lFCuuWfvnv2V34iZYRU+jiENON3pI6QSkLqCakkpBGE1BFS6XgheR5pJ0Iq\nCaknpJKQRhBSR0glIfWEVBLS0uG3f6qdkEpCKh0lpPtesyqkkJBKBwnpoR+v/t6LkEqHCemx\nw0tC6gipJKSekEpCWjTp1m53QiodJCTLhv0JqXSUkKy/dyek0nFCuouQOkIqCaknpJKQRhBS\nR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0gl\nIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkk\npJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE\n1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNS\nSUgjCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIq\nCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkE\nIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWEtHT4ldsnCakjpNJRQjo9uHWakDpC\nKh0kpId+bpckpI6QSocJ6bHDS0LqCKkkpJ6QSkJaNOnWbndCKh0kJMuG/QmpdJSQrL93J6TS\ncUK6i5A6QioJqSekkpBGEFJHSCUh9YRUEtKiydN9azshdYRUElJPSCUhLRu9b1ZIHSGVhNQT\nUklIy0avz/7g3d86+4yQMkIqHSWkG374J3989rtCygip9DEIac6tXUdIJSH1hFQS0ghC6gip\ndLyQPI+0EyGVhNQTUklIIwipI6SSkHpCKglp6bB3yO5LSKWjhOQzG3YnpNJBQvIpQvsTUukw\nIT12eElIHSGVhNQTUklIiybd2u1OSKWDhGTZsD8hlY4SkvX37oRUOk5IdxFSR0glIfWEVBLS\nCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLq\nCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgd\nIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWE\nNIKQOkIqCaknpJKQRhBSR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQ\nekIqCWkEIXWEVBJST0glIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBS\nR0glIfWEVBLSCELqCKkkpJ6QSkIaQUgdIZWE1BNSSUgjCKkjpJKQekIqCWkEIXWEVBJST0gl\nIY0gpI6QSkLqCakkpBGE1BFSSUg9IZWENIKQOkIqCaknpJKQRhBSR0ilQ4V0+sjNM4TUEVLp\nICFN/ZxOby1JSB0hlQ4U0ulNTTfOE1JHSKVDhXQ+vEpIHSGVhNQTUklIiyaFtDshlQ4T0us9\nw1u2DULqCKl0wJBunSekjpBKBwnpXkLqCKkkpJ6QSkJaOvz2p2OFVBJS6SghnR7cOk1IHSGV\nDhLSQz+2dnsRUukwIT12eElIHSGVhNQTUklIiyZv3Nr971/82dkfCCkjpNJBQrq1bPj3X//0\n2a+efnTlAYS0lJBKRwnJ+nt3QiodJ6S7CKkjpJKQekIqCWkEIXWEVDpeSNbfOxFSSUg9IZWE\nNIKQOkIqCaknpJKQlg57HmlfQiodJSRvo9idkEoHCcnbKPYnpNJhQnrs8JKQOkIqCaknpJKQ\nFk26tdudkEoHCcmyYX9CKh0lJOvv3QmpdJyQ7iKkjpBKQuoJqSSkEYTUEVJJSD0hlYQ0gpA6\nQioJqSekkpBGEFJHSCUh9YRUEtIIQuoIqSSknpBKQhpBSB0hlYTUE1JJSCMIqSOkkpB6QioJ\naQQhdYRUElJPSCUhjSCkjpBKQuoJqSSkEYTUEVJJSD0hlYQ0gpA6QioJqSekkpBGEFJHSCUh\n9YRUEtIIQuoIqSSknpBKQhpBSB0hlYTUE1JJSCMIqSOkkpB6QioJaQQhdYRUElJPSCUhjSCk\njpBKQuoJqSSkEYTUEVJJSD0hlYQ0gpA6QioJqSekkpBGEFJHSCUh9YRUEtIIQuoIqSSknpBK\nQhpBSB0hlYTUE1JJSCMIqSOkkpB6QioJaQQhdYRUElJPSCUhjSCkjpBKQuoJqSSkEYTUEVJJ\nSD0hlYQ0gpA6QioJqSekkpBGEFJHSCUh9YRUEtIIQuoIqSSknpBKQhpBSB0hlYTUE1JJSCMI\nqSOkkpB6QioJaQQhdYRUElJPSCUhLR1+5fZJQuoIqXSUkE4Pbp0mpI6QSgcJ6aGf2yUJqSOk\n0mFCeuzwkpA6QioJqSekkpAWTbq1252QSgcJybJhf0IqHSUk6+/dCal0nJDuIqSOkEpC6gmp\nJKSlw27t9iWk0lFCsmzYnZBKBwnJ+nt/QiodJqTHDi8JqSOkkpB6QioJadHkrVu7f/v+2V9f\nD+mz3z6KL98K6Z2nvrqzL9wK6QtPfXVn79wK6ctPfXVnnz1GSLeWDd/7xOwvP/HBlQf45uk4\nvnj96/z6U1/bzHvXL/O9p762ma9fv8wvPvW1zXzz+mUuNmj9/V/vP/iPa+MfvH8c135q/sxP\nnvraZn7ysb/MHz/1tc1c+z/8GuOfkIUXQEgQEBIEkpDe9ioheO6EBAEJQEBIEBj/Ngp4Aca/\njQJegPFvo4AXYPyrv+EFEBIE3NpBwLIBAtbfEFABBIQEASFBQEgQEBIEhAQBIUFASBAQEgSE\nBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIE\nhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBAS\nBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQ\nEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQ\nEBIEhAQBIUFASBAQEgSEBAEhQUBIEBASBIQEASFBQEgQEBIEhAQBIUFASBAQEgSEBAEhQUBI\nEBASBIQEASFBQEgQ+D9iY+wojjvqDQAAAABJRU5ErkJggg==",
"text/plain": [
"plot without title"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"print(colMeans(simulations))\n",
"barplot(colMeans(simulations))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case B: I want to predict the probabilties for position a, given I observed a 1 for position b. I expect a distribution with higher probabilties for category 4. The probabilties for cateogire 1,2,3 should be uniform."
]
},
{
"cell_type": "code",
"execution_count": 203,
"metadata": {},
"outputs": [],
"source": [
"s1 = post$a[,1]+post$b1[,1]\n",
"s2 = post$a[,2]+post$b1[,2]\n",
"s3 = post$a[,3]+post$b1[,3]\n",
"s4= 0"
]
},
{
"cell_type": "code",
"execution_count": 204,
"metadata": {},
"outputs": [],
"source": [
"simulations <- matrix(nrow = 0, ncol = 4) \n",
"for (i in 1:length(s1)){\n",
" s = c(s1[i],s2[i],s3[i],s4);\n",
" u = softmax(s);\n",
" simulations = rbind(simulations, u)\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 205,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] 0.00131950 0.03263794 0.01248308 0.95355948\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAw1BMVEUAAAABAQENDQ0REREY\nGBgZGRk1NTU+Pj4/Pz9ERERFRUVISEhJSUlKSkpMTExNTU1QUFBSUlJTU1NUVFRVVVVXV1dY\nWFhcXFxdXV1eXl5fX19gYGBjY2NmZmZubm53d3d4eHiAgICCgoKDg4OIiIiLi4ulpaWqqqqr\nq6uvr6+wsLCxsbG+vr7FxcXMzMzc3Nzd3d3f39/g4ODh4eHn5+fs7Ozu7u7w8PDx8fH09PT1\n9fX4+Pj5+fn6+vr7+/v8/Pz///8ERfSQAAAACXBIWXMAABJ0AAASdAHeZh94AAARp0lEQVR4\nnO3ca4Odd1XG4QUoClgOilUR8YCKRzzVauVgvv+nUmk7uemzZ2dm5W/WZPW6XiT7xaxn7gR+\nZLozpV4Bb62mB8AGNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkB\nsEFND4ANanoAbFDTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb\n1PQA2KCmB8AGNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFN\nD4ANanoAbFDTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA\n2KCmB8AGNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4AN\nanoAbFDTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCm\nB8AGNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4ANanoA\nbFDTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCmB/Dl\n8tFPX46PD/666uCz4I2+WS/Hnxz8ddXBZ8Eb/fYP/uGl+NYfH/x11cFnwRsJCQ4QEhwgJDhA\nSHCAkOAAIcEBQoIDhAQHCAkOEBIcICQ4QEhwgJDgACHBAUKCA4QEBwgJDhASHCAkOEBIcICQ\n4AAhwQFCggOEBAcICQ4QEhwgJDhASHCAkOAAIcEBQoIDhAQHCAkOEBIcICQ4QEhwgJDgACHB\nAUKCA4QEBwgJDhASHCAkOEBIcICQ4AAhwQFCggOEBAcICQ4QEhwgJDhASHCAkOAAIcEBQoID\nhAQHCAkOEBIcICQ4QEhwgJDgACHBAUKCA4QEBwgJDhASHCAkOEBIcICQ4AAhwQFCggOEBAcI\nCQ4QEhwgJDhASHCAkOAAIcEBQoIDhAQHCAkOEBIcICQ4QEhwgJDgACHduK361Q9v8wy+ZIR0\nPX3t1BrWE9Ll8rM/jz77AZ5CSJfLergXEk8lpMtlPdwLiacS0uXSl3Y8n5Cup95s4NmEdOPW\n2988l5DgACHBAUJ6jk/+6i8f/MUf/b98Ct5PQrr3kC8+5aNvfePB1+sXJz4HOwjp3kPuPeUf\n6+cnPgc7CKlLSAQhdQmJIKRbx0/561ghEYR0PX3adzYIiSCky+VDP/dLEhJBSJfLuvXySkgE\nIV0u69bLKyERhHS59KUdzyek66k3G3g2Id069vY3zySkLiERhNQlJIKQuoREEFKXkAhC6hIS\nQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQu\nIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC\n6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJ\nIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKX\nkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQh\ndQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoRE\nEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlL\nSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQ\nuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYR0Pa1Pf6z7zxASQUjX03r1aUf3SxIS\nQUiXy3r1+k+lOx8nJIKQLpf1Skg8l5Aul/VKSDyXkC6Xn5b08OoxQiII6Xr62r0PExJBSLeO\nn9CRkEhC6hISQUhdQiIIqUtIBCHde8gXn/Kvv/m1B1+tn534HOwgpHsP+eJTfvm3P3nwY38i\n8ZqQunxpRxBSl5AIQrp1/IS/RhISSUjXU9/ZwLMJ6XL50I/vtePJhHS5rFsvr4REENLlsm69\nvBISQUiXS1/a8XxCup56s4FnE9KtY29/80xC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSC\nkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1C\nIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTU\nJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJB\nSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4h\nEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELq\nEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkg\npC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQ\nCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1\nCYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCOnRJ7zhEUIiCOlyme58nJAIQrpcConn\nE9L1tPKnRwmJIKQbt/X6x8cJiSCkW8f/F5GQeAYh3bz+34qExDMI6ZH7++80vBISv0ZIjz1A\nSDyDkLqERBDSc/z33/3kwY+FxGtCuveQLz7lX37jaw++Wj878TnYQUj3HnLvKb60IwipS0gE\nIXUJiSCkW8dv/I7VV0Li1wjpevqkb/4WEklIl8uHfvxrFDyZkC6XdevllZAIQrpc1q2XV0Ii\nCOly6Us7nk9I11NvNvBsQrp17O1vnklIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQ\nCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1\nCYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQ\nUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtI\nBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6\nhEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiII\nqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUk\ngpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhd\nQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE\n1CUkgpC6hEQQUpeQCELqEhJBSF1CIgjpxm3V5y/ufJSQCEK6ntbnJQmJpxLS5bIe/kwSEk8l\npMtlvfq8JCHxVEK6XNanP5aQeDohXS7rs59KSDyZkK6n9dlPJSSeSkjX0/r8ZyHxVELqEhJB\nSF1CIgipS0gEId17yBef8tHvfOPB1+tnJz4HOwjp3kO++JRP/vxPH/yBP5F4TUhdvrQjCKlL\nSAQh3TquN/1t7Csh8WuEdD197d6HCYkgpMvlQz++s4EnE9Llsm69vBISQUiXy7r18kpIBCFd\nLn1px/MJ6XrqzQaeTUi3jr39zTMJqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoSEkFI\nXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCkLiER\nhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AIQuoS\nEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJiSCk\nLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBSl5AI\nQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gEIXUJ\niSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhdQiIIqUtIBCF1CYkgpC4hEYTUJSSCkLqERBBS\nl5AIQuoSEkFIXUIiCKlLSAQhdQmJIKQuIRGE1CUkgpC6hEQQUpeQCELqEhJBSF1CIgipS0gE\nIXUJiSCkLiERhNQlJIKQuoREEFKXkAhC6hISQUhd731I//XPL8cn078Zb01IXe99SD+ol+P3\npn8z3pqQut77kL7/zb9+Kb79wfRvxlsTUtf7H9J3p/8Tf/ChkA4S0rslpJOE1CWkc4R0kpDe\nLSGdJKQuIZ0jpJOE9G4J6SQhdQnpHCGdJKR3S0gnCenW8afuf5CQzhHSSS8lpPjOlXsfJqRz\nhHTSCwnpdT/3SxLSOfdC+vinL8d/PD5TSJfLuvXySkjn3Avpe+/ym2ff4A8fnymky2Xdenkl\npHPuhfTB9/7+pfjg+4/PFNLl0pd279zdkD6cXvfgu0J61ulr9z5MSOcI6aSXEpK3v985IZ30\nckJ6EiGdI6SThPRuCem5hHTMv73+/+v4m0dD+rN3+Ibsm/z+47+W9+RfNf/g29PrHnzzXki/\nO73uwW+9uJAu/5j0T1+J/45+5ZePnP1oqpobvvP4r+6H09vCh4/P/HB6W/jh4zO/M70t/Ojx\nmc9WRx5yecp/fvzavz929suPX447X37+Ynpb+MV7P/Pn09vCY/8L31EHnwVfWjU9ADaotzr+\n1Jkl8B6rtzh97dQaeE9V//KJ32sHXwLVv6xbL+FLqfqXdeslfClV/9KXdvC5eotTbzbAZ+qt\njmUEv1LTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCm\nB8AGNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4ANanoA\nbFDTA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCmB8AG\nNT0ANqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4ANanoAbFDT\nA2CDmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCmB8AGNT0A\nNqjpAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4ANanoAbFDTA2CD\nmh4AG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCmB8AGNT0ANqjp\nAbBBTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2qOkBsEFND4ANanoAbFDTA2CDmh4A\nG9T0ANigpgfABjU9ADao6QGwQU0PgA1qegBsUNMDYIOaHgAb1PQA2KCmB8AGNT0ANqjpAbBB\nTQ+ADWp6AGxQ0wNgg5oeABvU9ADYoKYHwAY1PQA2+B+QlAY1cOXu3gAAAABJRU5ErkJggg==",
"text/plain": [
"plot without title"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"print(colMeans(simulations))\n",
"barplot(colMeans(simulations))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "R",
"language": "R",
"name": "ir"
},
"language_info": {
"codemirror_mode": "r",
"file_extension": ".r",
"mimetype": "text/x-r-source",
"name": "R",
"pygments_lexer": "r",
"version": "3.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment