Skip to content

Instantly share code, notes, and snippets.

@thorstenwagner

thorstenwagner/catreg.ipynb Secret

Created July 15, 2021 14:13
Show Gist options
  • Save thorstenwagner/3d95f0bdf73e5da516d44079f2a5d676 to your computer and use it in GitHub Desktop.
Save thorstenwagner/3d95f0bdf73e5da516d44079f2a5d676 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": 156,
"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_a <- c(1,1,1,1)\n",
"pa <- softmax(prob_cat_a[1], prob_cat_a[2], prob_cat_a[3], prob_cat_a[4])\n",
"print(pa)"
]
},
{
"cell_type": "code",
"execution_count": 157,
"metadata": {},
"outputs": [],
"source": [
"cat_pos_a <- rep(NA,N)\n",
"set.seed(34302)\n",
"for ( i in 1:N ) cat_pos_a[i] <- sample( 1:4 , size = 1, prob=pa )"
]
},
{
"cell_type": "code",
"execution_count": 158,
"metadata": {},
"outputs": [],
"source": [
"cat_pos_b <- rep(NA,N)\n",
"for ( i in 1:N ) {\n",
" if (cat_pos_a[i] == 1 ) {\n",
" prob_cat_b <- c(1,1,1,5)\n",
" }\n",
" else if ( cat_pos_a[i] == 4 ) {\n",
" prob_cat_b <- c(5,1,1,1)\n",
" } \n",
" else {\n",
" prob_cat_b <- c(1,1,1,1)\n",
" }\n",
" \n",
" pb <- softmax(prob_cat_b[1], prob_cat_b[2], prob_cat_b[3], prob_cat_b[4])\n",
" cat_pos_b[i] <- sample( 1:4 , size = 1, prob=pb )\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 159,
"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": 160,
"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.000817 seconds\n",
"Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.17 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: 61.1961 seconds (Warm-up)\n",
"Chain 1: 61.0941 seconds (Sampling)\n",
"Chain 1: 122.29 seconds (Total)\n",
"Chain 1: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 2).\n",
"Chain 2: \n",
"Chain 2: Gradient evaluation took 0.000595 seconds\n",
"Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.95 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: 62.5527 seconds (Warm-up)\n",
"Chain 2: 49.3767 seconds (Sampling)\n",
"Chain 2: 111.929 seconds (Total)\n",
"Chain 2: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 3).\n",
"Chain 3: \n",
"Chain 3: Gradient evaluation took 0.00066 seconds\n",
"Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.6 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: 62.9365 seconds (Warm-up)\n",
"Chain 3: 56.2621 seconds (Sampling)\n",
"Chain 3: 119.199 seconds (Total)\n",
"Chain 3: \n",
"\n",
"SAMPLING FOR MODEL '4ff463eaa934e574e79b49a9bb230404' NOW (CHAIN 4).\n",
"Chain 4: \n",
"Chain 4: Gradient evaluation took 0.000696 seconds\n",
"Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 6.96 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: 62.0976 seconds (Warm-up)\n",
"Chain 4: 73.5381 seconds (Sampling)\n",
"Chain 4: 135.636 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": 161,
"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.1333182 </td><td>0.9395081 </td><td> -1.6509627</td><td> 1.3561493 </td><td>2029.099 </td><td>1.003792 </td></tr>\n",
"\t<tr><th scope=row>a[2]</th><td> 0.2132939 </td><td>0.9619450 </td><td> -1.3014405</td><td> 1.7491328 </td><td>1584.769 </td><td>1.004521 </td></tr>\n",
"\t<tr><th scope=row>a[3]</th><td> 0.2124571 </td><td>0.9298039 </td><td> -1.3103088</td><td> 1.6882524 </td><td>1485.001 </td><td>1.001244 </td></tr>\n",
"\t<tr><th scope=row>b1[1]</th><td>-8.1399452 </td><td>2.5878861 </td><td>-12.9382444</td><td>-4.7239700 </td><td>1933.191 </td><td>1.001000 </td></tr>\n",
"\t<tr><th scope=row>b1[2]</th><td>-1.3788134 </td><td>0.9681511 </td><td> -2.9259794</td><td> 0.1955229 </td><td>1607.110 </td><td>1.004258 </td></tr>\n",
"\t<tr><th scope=row>b1[3]</th><td>-1.7251122 </td><td>0.9437645 </td><td> -3.2079879</td><td>-0.2025233 </td><td>1519.283 </td><td>1.001214 </td></tr>\n",
"\t<tr><th scope=row>b2[1]</th><td> 0.8292764 </td><td>1.1237247 </td><td> -0.9673302</td><td> 2.6289144 </td><td>2064.360 </td><td>1.002502 </td></tr>\n",
"\t<tr><th scope=row>b2[2]</th><td> 2.3905512 </td><td>1.0966038 </td><td> 0.6575043</td><td> 4.1698508 </td><td>1745.774 </td><td>1.003277 </td></tr>\n",
"\t<tr><th scope=row>b2[3]</th><td> 2.4528399 </td><td>1.0577231 </td><td> 0.7875985</td><td> 4.1739539 </td><td>1661.334 </td><td>1.000736 </td></tr>\n",
"\t<tr><th scope=row>b3[1]</th><td>-0.2562160 </td><td>1.2478642 </td><td> -2.2755439</td><td> 1.7250921 </td><td>2338.915 </td><td>1.002619 </td></tr>\n",
"\t<tr><th scope=row>b3[2]</th><td> 2.4678523 </td><td>1.1031657 </td><td> 0.7468782</td><td> 4.2462882 </td><td>1704.054 </td><td>1.003890 </td></tr>\n",
"\t<tr><th scope=row>b3[3]</th><td> 2.5050899 </td><td>1.0729435 </td><td> 0.8478881</td><td> 4.2248268 </td><td>1606.760 </td><td>1.000925 </td></tr>\n",
"\t<tr><th scope=row>b4[1]</th><td> 4.0276696 </td><td>1.0333838 </td><td> 2.3604726</td><td> 5.7026892 </td><td>2088.128 </td><td>1.002616 </td></tr>\n",
"\t<tr><th scope=row>b4[2]</th><td> 2.3888018 </td><td>1.0614906 </td><td> 0.6866429</td><td> 4.0829525 </td><td>1634.832 </td><td>1.003456 </td></tr>\n",
"\t<tr><th scope=row>b4[3]</th><td> 2.3871083 </td><td>1.0303443 </td><td> 0.7443001</td><td> 4.0304591 </td><td>1579.773 </td><td>1.000682 </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.1333182 & 0.9395081 & -1.6509627 & 1.3561493 & 2029.099 & 1.003792 \\\\\n",
"\ta{[}2{]} & 0.2132939 & 0.9619450 & -1.3014405 & 1.7491328 & 1584.769 & 1.004521 \\\\\n",
"\ta{[}3{]} & 0.2124571 & 0.9298039 & -1.3103088 & 1.6882524 & 1485.001 & 1.001244 \\\\\n",
"\tb1{[}1{]} & -8.1399452 & 2.5878861 & -12.9382444 & -4.7239700 & 1933.191 & 1.001000 \\\\\n",
"\tb1{[}2{]} & -1.3788134 & 0.9681511 & -2.9259794 & 0.1955229 & 1607.110 & 1.004258 \\\\\n",
"\tb1{[}3{]} & -1.7251122 & 0.9437645 & -3.2079879 & -0.2025233 & 1519.283 & 1.001214 \\\\\n",
"\tb2{[}1{]} & 0.8292764 & 1.1237247 & -0.9673302 & 2.6289144 & 2064.360 & 1.002502 \\\\\n",
"\tb2{[}2{]} & 2.3905512 & 1.0966038 & 0.6575043 & 4.1698508 & 1745.774 & 1.003277 \\\\\n",
"\tb2{[}3{]} & 2.4528399 & 1.0577231 & 0.7875985 & 4.1739539 & 1661.334 & 1.000736 \\\\\n",
"\tb3{[}1{]} & -0.2562160 & 1.2478642 & -2.2755439 & 1.7250921 & 2338.915 & 1.002619 \\\\\n",
"\tb3{[}2{]} & 2.4678523 & 1.1031657 & 0.7468782 & 4.2462882 & 1704.054 & 1.003890 \\\\\n",
"\tb3{[}3{]} & 2.5050899 & 1.0729435 & 0.8478881 & 4.2248268 & 1606.760 & 1.000925 \\\\\n",
"\tb4{[}1{]} & 4.0276696 & 1.0333838 & 2.3604726 & 5.7026892 & 2088.128 & 1.002616 \\\\\n",
"\tb4{[}2{]} & 2.3888018 & 1.0614906 & 0.6866429 & 4.0829525 & 1634.832 & 1.003456 \\\\\n",
"\tb4{[}3{]} & 2.3871083 & 1.0303443 & 0.7443001 & 4.0304591 & 1579.773 & 1.000682 \\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"| <!--/--> | mean | sd | 5.5% | 94.5% | n_eff | Rhat4 |\n",
"|---|---|---|---|---|---|---|\n",
"| a[1] | -0.1333182 | 0.9395081 | -1.6509627 | 1.3561493 | 2029.099 | 1.003792 |\n",
"| a[2] | 0.2132939 | 0.9619450 | -1.3014405 | 1.7491328 | 1584.769 | 1.004521 |\n",
"| a[3] | 0.2124571 | 0.9298039 | -1.3103088 | 1.6882524 | 1485.001 | 1.001244 |\n",
"| b1[1] | -8.1399452 | 2.5878861 | -12.9382444 | -4.7239700 | 1933.191 | 1.001000 |\n",
"| b1[2] | -1.3788134 | 0.9681511 | -2.9259794 | 0.1955229 | 1607.110 | 1.004258 |\n",
"| b1[3] | -1.7251122 | 0.9437645 | -3.2079879 | -0.2025233 | 1519.283 | 1.001214 |\n",
"| b2[1] | 0.8292764 | 1.1237247 | -0.9673302 | 2.6289144 | 2064.360 | 1.002502 |\n",
"| b2[2] | 2.3905512 | 1.0966038 | 0.6575043 | 4.1698508 | 1745.774 | 1.003277 |\n",
"| b2[3] | 2.4528399 | 1.0577231 | 0.7875985 | 4.1739539 | 1661.334 | 1.000736 |\n",
"| b3[1] | -0.2562160 | 1.2478642 | -2.2755439 | 1.7250921 | 2338.915 | 1.002619 |\n",
"| b3[2] | 2.4678523 | 1.1031657 | 0.7468782 | 4.2462882 | 1704.054 | 1.003890 |\n",
"| b3[3] | 2.5050899 | 1.0729435 | 0.8478881 | 4.2248268 | 1606.760 | 1.000925 |\n",
"| b4[1] | 4.0276696 | 1.0333838 | 2.3604726 | 5.7026892 | 2088.128 | 1.002616 |\n",
"| b4[2] | 2.3888018 | 1.0614906 | 0.6866429 | 4.0829525 | 1634.832 | 1.003456 |\n",
"| b4[3] | 2.3871083 | 1.0303443 | 0.7443001 | 4.0304591 | 1579.773 | 1.000682 |\n",
"\n"
],
"text/plain": [
" mean sd 5.5% 94.5% n_eff Rhat4 \n",
"a[1] -0.1333182 0.9395081 -1.6509627 1.3561493 2029.099 1.003792\n",
"a[2] 0.2132939 0.9619450 -1.3014405 1.7491328 1584.769 1.004521\n",
"a[3] 0.2124571 0.9298039 -1.3103088 1.6882524 1485.001 1.001244\n",
"b1[1] -8.1399452 2.5878861 -12.9382444 -4.7239700 1933.191 1.001000\n",
"b1[2] -1.3788134 0.9681511 -2.9259794 0.1955229 1607.110 1.004258\n",
"b1[3] -1.7251122 0.9437645 -3.2079879 -0.2025233 1519.283 1.001214\n",
"b2[1] 0.8292764 1.1237247 -0.9673302 2.6289144 2064.360 1.002502\n",
"b2[2] 2.3905512 1.0966038 0.6575043 4.1698508 1745.774 1.003277\n",
"b2[3] 2.4528399 1.0577231 0.7875985 4.1739539 1661.334 1.000736\n",
"b3[1] -0.2562160 1.2478642 -2.2755439 1.7250921 2338.915 1.002619\n",
"b3[2] 2.4678523 1.1031657 0.7468782 4.2462882 1704.054 1.003890\n",
"b3[3] 2.5050899 1.0729435 0.8478881 4.2248268 1606.760 1.000925\n",
"b4[1] 4.0276696 1.0333838 2.3604726 5.7026892 2088.128 1.002616\n",
"b4[2] 2.3888018 1.0614906 0.6866429 4.0829525 1634.832 1.003456\n",
"b4[3] 2.3871083 1.0303443 0.7443001 4.0304591 1579.773 1.000682"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"precis(erics_model,depth=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Making predictions"
]
},
{
"cell_type": "code",
"execution_count": 184,
"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 roughly equal distribution category 2 and 3."
]
},
{
"cell_type": "code",
"execution_count": 185,
"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": 186,
"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": 187,
"metadata": {},
"outputs": [],
"source": [
"require(datasets)"
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] 0.06798879 0.43437440 0.46161577 0.03602103\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAA6lBMVEUAAAABAQENDQ0ODg4R\nERETExMUFBQbGxsiIiImJiYnJycpKSkzMzM0NDQ8PDw+Pj4/Pz9AQEBCQkJFRUVGRkZISEhK\nSkpMTExNTU1PT09QUFBTU1NVVVVYWFhbW1tcXFxfX19gYGBhYWFjY2NmZmZnZ2dpaWltbW1u\nbm55eXl8fHyIiIiLi4uNjY2QkJCVlZWYmJiZmZmfn5+lpaWrq6u7u7u+vr6/v7/FxcXMzMzO\nzs7c3Nzd3d3e3t7h4eHi4uLm5ubq6urr6+vs7Ozu7u7w8PDx8fH19fX4+Pj5+fn6+vr7+/v9\n/f3///+E1kStAAAACXBIWXMAABJ0AAASdAHeZh94AAATzklEQVR4nO3dWXcj13WGYTjzbCUe\nI6sdO87gxJknZ56sJJat/v9/JxEkUb1EkLs38W3UDvQ8F+y6wCkckHybhVOHi6fXwNVOR08A\n7oGQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkC\nhAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFB\ngJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAk\nCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKE\nBAFCggAh3Yv/+Ic9fnT0J+P2hHQv3j3t8e2jPxm3J6R78fVXf7fFu+8d/cm4PSHdCyEdSkj3\nQkiHEtK9ENKhhHQvhHQoId0LIR1KSPdCSIcS0r0Q0qGEdC+EdCgh3QshHUpI90JIhxLSvRDS\noYR0L4R0KCHdCyEdSkj3QkiHEtK9ENKhhHQvhHQoId0LIR1KSPdCSIcS0r0Q0qGEdC+EdCgh\n3QshHUpI90JIhxLSvRDSoYR0L4R0KCHdCyEdSkj3QkiHEtK9ENKhhHQvhHQoId0LIR1KSPdC\nSIcS0r0Q0qGEdC+EdCgh3QshHUpIpQ/f3+PDp6cppEMJqfT7R/9F1jf8ztPTFNKhhFR6752/\n2OKrX396mkI6lJBK77179Dfmg1dC2kpIJSF1CYkLhNQlJC4QUpeQuEBIXULiAiF1CYkLhNQl\nJC4QUpeQuEBIXULiAiF1CYkLhNQlJC4QUpeQuEBIXULiAiF1CYkLhNQlJC4QUpeQuEBIXULi\nAiF1CYkLhNQlJC4QUpeQuEBIXULiAiF1CYkLhNQlJC4QUpeQuEBIXULiAiF1CYkLhNQlJC4Q\nUpeQuEBIXULiAiF1CYkLhNQlJC4QUpeQuEBIXUJ60RnuvUUhdQmpNfJNwRmtI6QuIbVGCunm\nhLTWFQl8ks99V/RaSH1Cao49ffbxjgmpS0jdwR9FJKTbEdJa11XwUUlCuh0hrXVtBfe+0vBa\nSH1CeskJhHQ7Qlrr3isIEFKXkFI++Ku/fPDnfzzyFLcjpC4hvfQknz/Lv/3qLz34xdMHiec4\njpC6hPTSkzx3lr8XUoyQ1pp/jySkHCGtJaSSkLqE1B38NjtWhZQjpLWu2rT6Vpu/hZQjpLWu\n+jWKR0eXCClHSGtdE9Klw8eElCOktYRUElKXkFojXdrdnJDWsthQElKXkLqDLX/flpDWckO2\nJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlCWktIJSF1CWmCkHKEtJaQSkLqEtIE\nIeUIaS0hlYTUJaQJQsoR0lpCKgmpS0gThJQjpLWEVBJSl5AmCClHSGsJqSSkLiFNEFKOkNYS\nUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpByhLSWkEpC6hLSBCHlCGktIZWE1CWk\nCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakkpC4hTRBSjpDWElJJSF1CmiCkHCGt\nJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh5QhpLSGVhNQlpAlCyhHSWkIqCalL\nSBOElCOktYRUElKXkCYIKUdIawmpJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlC\nWktIJSF1CWmCkHKEtJaQSkLqEtIEIeUIaS0hlYTUJaQJQsoR0lpCKgmpS0gThJQjpLWEVBJS\nl5AmCClHSGsJqSSkLiFNEFKOkNYSUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpBy\nhLSWkEpC6hLSBCHlCGktIZWE1CWkCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakk\npC4hTRBSjpDWElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh\n5QhpLSGVhNQlpAlCyhHSWkIqCalLSBOElCOktYRUElKXkCYIKUdIawmpJKQuIU0QUo6Q1hJS\nSUhdQpogpBwhrSWkkpC6hDRBSDlCWktIJSF1CWmCkHKEtJaQSkLqEtIEIeUIaS0hlYTUJaQJ\nQsoR0lpCKgmpS0gThJQjpLWuCOl0On388flzCClHSGu9PKTTOaHTqSpJSDlCWuuakD79cfR8\nSULKEdJa14X08YfXQroRIa0lpJKQuoTUGnm6dPiYkHKEtNbVq3bVWyQhBQlprWvuIz2E9Oyj\nhJQjpLXckC0JqUtIE4SUI6S1ZkL6rz/8vQe/LaQYIa0VCenRu6T//PZ7D37r9OPEcxxHSF1C\neulJLH/fhpDW8h6pJKQuIU0QUo6Q1roqpFO991tISUJa67qdDae3SElIOUJa66rfR3p0dImQ\ncoS0lk2rJSF1Cak1Ukg3J6S1XNqVhNQlpN5Qiw23JqS1LH+XhNQlpAlCyhHSWkIqCalLSBOE\nlCOktYRUElKXkCYIKUdIawmpJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlCWktI\nJSF1CWmCkHKEtJaQSkLqEtIEIeUIaS0hlYTUJaQJQsoR0lpCKgmpS0gThJQjpLWEVBJSl5Am\nCClHSGsJqSSkLiFNEFKOkNYSUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpByhLSW\nkEpC6hLSBCHlCGktIZWE1CWkCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakkpC4h\nTRBSjpDWElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh5Qhp\nLSGVhNQlpAlCyhHSWkIqCalLSBOElCOktYRUElKXkCYIKUdIawmpJKQuIU0QUo6Q1hJSSUhd\nQpogpBwhrSWkkpC6hDRBSDlCWktIJSF1CWmCkHKEtJaQSkLqEtIEIeUIaS0hlYTUJaQJQsoR\n0lpCKgmpS0gThJQjpLWEVBJSl5AmCClHSGsJqSSkLiFNEFKOkNYSUklIXUKaIKQcIa0lpJKQ\nuoQ0QUg5QlpLSCUhdQlpgpByhLSWkEpC6hLSBCHlCGktIZWE1CWkCULKEdJaQioJqUtIE4SU\nI6S1hFQSUpeQJggpR0hrCakkpC4hTRBSjpDWElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0gl\nIXUJaYKQcoS0lpBKQuoS0gQh5QhpLSGVhNQlpAlCyhHSWkIqCalLSBOElCOktYRUElKXkCYI\nKUdIawmpJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlCWktIJSF1CWmCkHKEtNY1\nIZ1Op08PnnmUkHKEtNYVIZ1On5YkpBsR0lovD+mjej4pSUg3IqS1rgvpk5KEdCNCWuvKkD4u\nSUg3IqS1rg3pXJKQbkRIa1212PDJPych3YiQ1ro+pNdCuhkhreWGbElIXUKaIKQcIa0lpJKQ\nuoT00pN8/iz//hu/9uBXTj9OPMdxhNQlpJee5PNn+Z8//cGDP/ATKUZIa7m0KwmpS0gThJQj\npLWuCul0qu7GvhZSkpDWuvbXKOqUhJQjpLWu2Wt3enR0iZByhLTW1ZtWP3f4mJByhLSWkEpC\n6hJSa6RLu5sT0loWG0pC6hJSd7Dl79sS0lpuyJaE1CWkCULKEdJaQioJqUtIE4SUI6S1hFQS\nUpeQJggpR0hrCakkpC4hTRBSjpDWElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0glIXUJaYKQ\ncoS0lpBKQuoS0gQh5QhpLSGVhNQlpAlCyhHSWkIqCalLSBOElCOktYRUElKXkCYIKUdIawmp\nJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlCWktIJSF1CWmCkHKEtJaQSkLqEtIE\nIeUIaS0hlYTUJaQJQsoR0lpCKgmpS0gThJQjpLWEVBJSl5AmCClHSGsJqSSkLiFNEFKOkNYS\nUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpByhLSWkEpC6hLSBCHlCGktIZWE1CWk\nCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakkpC4hTRBSjpDWElJJSF1CmiCkHCGt\nJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh5QhpLSGVhNQlpAlCyhHSWkIqCalL\nSBOElCOktYRUElKXkCYIKUdIawmpJKQuIU0QUo6Q1hJSSUhdQpogpBwhrSWkkpC6hDRBSDlC\nWktIJSF1CWmCkHKEtJaQSkLqEtIEIeUIaS0hlYTUJaQJQsoR0lpCKgmpS0gThJQjpLWEVBJS\nl5AmCClHSGsJqSSkLiFNEFKOkNYSUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpBy\nhLSWkEpC6hLSBCHlCGktIZWE1CWkCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakk\npC4hTRBSjpDWElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh\n5QhprZeHdHrTM48TUo6Q1hJSSUhdQuoNfbuxQsoR0lpCKgmpS0i9oUK6NSGtZdWuJKQuIU0Q\nUo6Q1hJSSUhdQpogpBwhrRUJ6dGyw7/8ws89+BkhxQhprZmQPvzh3z74MyHFCGktl3YlIXUJ\naYKQcoS01lUh1RvtXgspSUhrXbWz4W32rAopSEhrXbX7+9HRJULKEdJa14R06fAxIeUIaS0h\nlYTUJaTWSJd2NyektSw2lITUJaTuYMvftyWktdyQLQmpS0gThJQjpLWEVBJSl5AmCClHSGsJ\nqSSkLiFNEFKOkNYSUklIXUKaIKQcIa0lpJKQuoQ0QUg5QlpLSCUhdQlpgpByhLSWkEpC6hLS\nBCHlCGktIZWE1CWkCULKEdJaQioJqUtIE4SUI6S1hFQSUpeQJggpR0hrCakkpC4hTRBSjpDW\nElJJSF1CmiCkHCGtJaSSkLqENEFIOUJaS0glIXUJaYKQcoS0lpBKQuoS0gQh5QhprSND+scf\n7PHDp1+AkLqENOHpkH73Z395i5//zadfgJC6hDTh6ZC++7Wjv+IPvvPO0y9ASF1CmiCkHCGt\nJaQzISUJaYKQcoS0lpDOhJQkpAlCyhHSWkI6E1KSkCYIKUdIawnpTEhJQpogpBwhrSWkMyEl\nCWmCkHKEtJaQzoSUJKQJQsoR0lpCOhNSkpAmCClHSGsJ6UxISUKaIKQcIa0lpDMhJQlpgpBy\nhLSWkM6ElCSkCULKEdJaQjoTUpKQJggpR0hrCelMSElCmiCkHCGtJaQzISUJaYKQcoS0lpDO\nhJQkpAlCyhHSWkI6E1KSkCYIKUdIawnpTEhJQpogpBwhrSWkMyElCWmCkHKEtJaQzoSUJKQJ\nQsoR0lpCOhNSkpAmCClHSGsJ6UxISUKaIKScOwjp1WmPP3qrb+C3I6QzISU9F9I7X/2TLX79\nu2/1Dfx2hHQmpKRnQ/rO0bN78DUhxQkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQj\npCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5Qkp\nSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQj\npCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5QkpSUgThJQjpCQh5Qkp\nSUgThJQjpKQ9IX3yN22ff5CQcoSUtCWkN/489HMPE1KOkJKWhPRZP8+XJKQcISWtCenS4WNC\nyhFSkpDyhJQkpNZIl3Y3J6SkJSFZbLg9ISVtCcny980JKWlPSG9FSDlCShJSnpCShBTzr//8\n4G+eDunLf73Ft54L6StHz+7BN54L6RtHz+7BV54L6VtHz+7Bl9eF9Oht0j996Y2ViC/99Ilh\n3z/t8c2nX933jp7bG149Pc1XR8/tDd97eprfPHpub/j+09Nsmwnp9X+//5kfPTXsp+/v8dRP\nzf/zk6Pn9oaf/L+f5gdHz+0NT/0P/xLz75HgC0BIEDB/Hwm+AOZ3NsAXwPxeO/gCmN/9DV8A\nQoIAl3YQYLEBAix/Q4AKIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQB\nQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFC\nggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBA\nSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIE\nCAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQL+F8xpTonZXF/eAAAAAElFTkSuQmCC",
"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"
]
},
{
"cell_type": "code",
"execution_count": 189,
"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": 190,
"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": 191,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] 0.0008275309 0.2038204867 0.1444212091 0.6509307733\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAA21BMVEUAAAABAQEHBwcNDQ0R\nEREZGRkaGhoiIiImJiY1NTU+Pj4/Pz9ERERFRUVISEhKSkpMTExNTU1QUFBSUlJTU1NVVVVW\nVlZYWFhcXFxdXV1fX19gYGBjY2NmZmZubm53d3d8fHx+fn6CgoKDg4Oenp6lpaWqqqqrq6ut\nra2urq6vr6+xsbG4uLi6urq7u7u+vr7CwsLFxcXMzMzU1NTc3Nzd3d3g4ODh4eHn5+fo6Ojp\n6ens7Ozu7u7v7+/x8fHy8vLz8/P19fX39/f4+Pj5+fn6+vr7+/v+/v7////AsNXjAAAACXBI\nWXMAABJ0AAASdAHeZh94AAASJElEQVR4nO3ba59kV1nG4ULxgBITQAwKgorg+YQgHiKoIXG+\n/ydSSKb7zuzuml5P7r26ZnJdL7rrxX5W7er5/ad2rd59eQF8bpfnPgF4GwgJCoQEBUKCAiFB\ngZCgQEhQICQoEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAk\nKBASFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqE\nBAVCggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIUCAkKhAQFQoICIUGB\nkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCgQEhQICQo\nEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBASW/3rj27H\nB8XXJSS2eudyO/64+LqExFa/+0f/cCt+79vF1yUkthISFAgJCoQEBUKCAiFBgZCgQEhQICQo\nEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBASFAgJCoQE\nBUKCAiFBgZCgQEhQICQoEBIUCAkKhHQcvXzy9SJGnkxIx9HLi086UhJPJqTD5OXF/btS63R4\n2wnpMCkk1gnpMCkk1gnpMPlJSXeP4AmEdBy91zsf3nJCemhYRywSEhQICQqEtOJ///Hv7/zd\nX57yFLyZhHRtkVdX+cmvxU7E5aPGc/B2ENK1Ra6t8k+XXzSeg7eDkKaERBDSlJAIQnpo+Cm/\nRhISQUjH0afd2SAkgpAOk3f9XC9JSAQhHSYvDz08EhJBSIdJIbFOSIdJl3asE9Jx1GYDy4T0\n0LDtbxYJaUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQp\nIRGENCUkgpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhC\nmhISQUhTQiIIaUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJ\nIKQpIRGENCUkgpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKU\nkAhCmhISQUhTQiIIaUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQh\nTQmJIKQpIRGENCUkgpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoRE\nENKUkAhCmhISQUiPzP/S1SOERBDSYfLy4tOOrpckJIKQDpOXX3X04uXXxwiJIKTD5OXTd6W7\nbw8TEkFIh0khsU5Ih0khsU5Ih8mX+wyv2W0QEkFIh8kI6dpxQiIIaUpIBCFNCYkgpCU//+DO\n3wqJe0K6tsirq/z4S5fwYeM5eDsI6doih1X++Ud3/tw7EveENOUzEkFIU0IiCOmh4dff+y0k\nPkNIx9F71w4TEkFIh8m7ftwixJMJ6TB5eejhkZAIQjpMCol1QjpMurRjnZCOozYbWCakh4Zt\nf7NISFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpIBCFNCYkg\npCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCmhEQQ0pSQ\nCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpIBCFN\nCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCmhEQQ\n0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpI\nBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCm\nhEQQ0pSQCEJ6dIXXLCEkgpAOk+nKcUIiCOkwKSTWCek4eslvjxISQUgPzF7uvz5OSAQhPTT8\ny4iExAIhPTj9/xUJiQVCemT++k7DCyHxGUJ6bAEhsUBIU0IiCGlKSAQhXVvk1VV++tu/dec3\nLx82noO3g5CuLfLqKh//1Q/vfN87EveENOXSjiCkKSERhPTQ8GvvWH0hJD5DSMfRJ938LSSS\nkA6Td/34MwqeTEiHyctDD4+ERBDSYVJIrBPSYdKlHeuEdBy12cAyIT00bPubRUKaEhJBSFNC\nIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpIBCFNCYkgpCkhEYQ0\nJSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCmhEQQ0pSQCEKaEhJB\nSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpIBCFNCYkgpCkh\nEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCmhEQQ0pSQCEKa\nEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhTQiIIaUpIBCFNCYkg\npCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUkgpCmhEQQ0nH0\ncvnk6/U1hEQQ0mHyVwldLq8rSUgEIR0mLy/fjq6XJCSCkA6Tl5dfXgiJpxLSYVJIrBPSYfLy\n0MMjIRGEdBz9tB+bDTydkB6YveS3xwiJIKQpIRGENCUkgpCmhEQQ0rVFXl3l37/21TtfuXzY\neA7eDkK6tsirq/z8T75351vekbgnpCmXdgQhTQmJIKSHhl9/77eQ+AwhHUfvXTtMSAQhHSbv\n+vFnFDyZkA6TblplnZAOk0JinZAOky7tWCek46jNBpYJ6aFh298sEtKUkAhCmhISQUhTQiII\naUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGENCUk\ngpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhISQUhT\nQiIIaUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQpIRGE\nNCUkgpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENKUkAhCmhIS\nQUhTQiIIaUpIBCFNCYkgpCkhEYQ0JSSCkKaERBDSlJAIQpoSEkFIU0IiCGlKSAQhTQmJIKQp\nIRGENCUkgpCmhEQQ0pSQCEKaEhJBSFNCIghpSkgEIU0JiSCkKSERhDQlJIKQpoREENIDs5fL\nywdXjhISQUjH0cvLkoTEUwnpMHm5e08SEk8lpMPk5cXLkoTEUwnpMPnpm9FFSDydkA6TLz8f\nXYTEkwnpOHpX0tsd0v/85Hb893P/MD43IR1H7za/3+6Qvn25HX/w3D+Mz01IU298SN985we3\n4r13n/uH8bkJaerND+nrz/0vfud9IRW9ASH9159+7863hFQjpKbbC+nwGenf3v3qna9cPmw8\nx/MRUpOQri3yVm82CKlJSFNC6hFSk5D2ElKTkB4a/sT1g4TUI6SmWwkpfk947TAh9Qip6UZC\nuu/neklC6hFS082E9NDDIyH1CKlJSHsJqUlIh0mXdtsJqelGQrLZsJ+Qmm4lJNvf2wmp6XZC\nehIh9QipSUh7CalJSFNC6hFSk5D2ElKTkKaE1COkJiHtJaQmIU0JqUdITULaS0hNQpoSUo+Q\nmoS0l5CahDQlpB4hNQlpLyE1CWlKSD1CahLSXkJqEtKUkHqE1CSkvYTUJKQpIfUIqUlIewmp\nSUhTQuoRUpOQ9hJSk5CmhNQjpCYh7SWkJiFNCalHSE1C2ktITUKaElKPkJqEtJeQmoQ0JaQe\nITUJaS8hNQlpSkg9QmoS0l5CahLSlJB6hNQkpL2E1CSkKSH1CKlJSHsJqUlIU0LqEVKTkPYS\nUpOQpoTUI6QmIe0lpCYhTQmpR0hNQtrrDQnpP390O372+GkKaUpIPddC+sbldvzh46cppCkh\n9VwL6d1v/M2tePebj5+mkKaE1HM1pPef++zufF1IJxBSj5CahLSXkFYJ6QxC6hFSk5D2EtIq\nIZ1BSD1CahLSXkJaJaQzCKlHSE1C2ktIq4R0BiH1CKlJSHsJaZWQziCkHiE1CWkvIa0S0hmE\n1COkJiHtJaRVQjqDkHqE1CSkvYS0SkhnEFKPkJqEtJeQVgnpDELqEVKTkPYS0iohnUFIPUJq\nEtJeQlolpDMIqUdITULaS0irhHQGIfUIqUlIewlplZDOIKQeITUJaS8hrRLSGYTUI6QmIe0l\npFVCWppMV44TUo+QmoS0l5BWCWlt9GmzQuoRUpOQ9hLSKiGtjQppNyE13UpITySkHiE1CWkv\nIa0S0hmE1COkpjcgpI//+od3vi+kGiE13V5Ih22Hf/mNL9/59cuHjed4PkJaJaTpIn4hu4eQ\nmm4vpKuE1COkJiHtJaRVQlodfu2Ndi+E1CSkplsJ6Wn3rAqpSEhNNxLSfT/u/t5FSE03E9JD\nD4+E1COkJiHtJaRVQlqadGm3nZCabiQkmw37CanpVkKy/b2dkJpuJ6QnEVKPkJqEtJeQVgnp\nDELqEVKTkPYS0iohnUFIPUJqEtJeQlolpDMIqUdITULaS0irhHQGIfUIqUlIewlplZDOIKQe\nITUJaS8hrRLSGYTUI6QmIe0lpFVCOoOQeoTUJKS9hLRKSGcQUo+QmoS0l5BWCekMQuoRUpOQ\n9hLSKiGdQUg9QmoS0l5CWiWkMwipR0hNQtpLSKuEdAYh9QipSUh7CWmVkM4gpB4hNQlpLyGt\nEtIZhNQjpCYh7SWkVUI6g5B6hNQkpL2EtEpIZxBSj5CahLSXkFYJ6QxC6hFSk5D2EtIqIZ1B\nSD1CahLSXkJaJaQzCKlHSE1C2ktIq4R0BiH1CKlJSHsJaZWQziCkHiE1CWkvIa0S0hmE1COk\nJiHtJaRVQjqDkHqE1CSkvYS0SkhnEFKPkJqEtJeQVgnpDELqEVKTkPYS0iohnUFIPUJqEtJe\nQlolpDMIqUdITULaS0irhHQGIfUIqUlIewlplZDOIKQeITUJaS8hrRLSGYTUI6QmIe0lpFVC\nOoOQeoTUJKS9hLRKSGcQUo+QmoS0l5BWCekMQuoRUpOQ9hLSKiGdQUg9QmoS0l5CWiWkMwip\nR0hNtxPS5RPXDxJSj5CabiWky71rhwmpR0hNNxLSfT/XSxJSj5Cabiakhx4eCalHSE1C2ktI\nq4S0NOnSbjshNd1ISDYb9hNS062EZPt7OyE13U5ITyKkHiE1CWkvIa0SUs1Pf3LnLx4N6c++\nfDu+8/hr+eY7P7gV710L6b3nPrs771wL6fef++zu/M7NhXT4mPTjL8VOxJc+fmTsu5fb8bXH\nX913nvvcwvuPn+b7z31u4cp/S1977nML3338NJedE9KLn31w7z8eG/v4g9tx5fLzo+c+t/DR\nG3+av3jucwuP/Q8/cf5nJPgCEBIUnP97JPgCOP/OBvgCOP9eO/gCOP/ub/gCEBIUuLSDApsN\nUGD7GwpUAAVCggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIUCAkKhAQF\nQoICIUGBkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCg\nQEhQICQoEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBAS\nFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIUCAkKhAQFQoICIUGBkKBASFAgJCgQEhQICQqEBAVC\nggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIUCAkKhAQFQoICIUGBkKBA\nSFAgJCgQEhQICQqEBAVCggIhQYGQoEBIUCAkKBASFAgJCoQEBUKCAiFBgZCgQEhQICQoEBIU\nCAkKhAQFQoICIUGBkKBASFDwf8/p7Tl0yJqTAAAAAElFTkSuQmCC",
"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": []
}
],
"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