Skip to content

Commit 5eb7b7e

Browse files
martin-smiramartin-smira
committed
Finished models for first 2 chapters
1 parent 40a1a89 commit 5eb7b7e

29 files changed

+3750
-0
lines changed

‎Bayesian_Cognitive_Modeling/!underDevelopment/!Exams_1_Stan (with 2 models).R‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
###############################################################################
7+
# Model 1 - is very simple, but is returning only parameter Phi
8+
# Model 2 - model is more complicated and is able to generate parameters z
9+
###############################################################################
10+
# Your choice: ################################################################
11+
choice <- 2 ################################################################
12+
###############################################################################
13+
14+
if (choice == 1) {
15+
model <- "
16+
# Exam Scores
17+
data {
18+
int<lower=0> p;
19+
int<lower=0> k[p];
20+
int<lower=0> n;
21+
}
22+
transformed data {
23+
real psi;
24+
25+
psi <- .5;
26+
}
27+
parameters {
28+
real<lower=.5, upper=1> phi;
29+
vector<lower=0, upper=1>[p] mix;
30+
}
31+
model {
32+
phi ~ beta(1, 1)T[.5, 1];
33+
34+
for (i in 1:p) {
35+
increment_log_prob(log_sum_exp(log(mix[i]) + binomial_log(k[i], n, phi),
36+
log1m(mix[i]) + binomial_log(k[i], n, psi)));
37+
}
38+
}"
39+
} else if (choice == 2) {
40+
model <- "
41+
# Exam Scores
42+
data {
43+
int<lower=0> p;
44+
int<lower=0> k[p];
45+
int<lower=0> n;
46+
}
47+
transformed data {
48+
real psi;
49+
50+
psi <- .5;
51+
}
52+
parameters {
53+
real<lower=.5, upper=1> phi;
54+
vector<lower=0, upper=1>[p] mix;
55+
}
56+
transformed parameters {
57+
matrix[p, 2] lp_parts;
58+
59+
for (i in 1:p) {
60+
lp_parts[i, 1] <- log(mix[i]) + binomial_log(k[i], n, phi);
61+
lp_parts[i, 2] <- log1m(mix[i]) + binomial_log(k[i], n, psi);
62+
}
63+
}
64+
model {
65+
phi ~ beta(1, 1)T[.5, 1];
66+
67+
for (i in 1:p) {
68+
increment_log_prob(log_sum_exp(lp_parts[i]));
69+
}
70+
}
71+
generated quantities {
72+
vector<lower=0, upper=1>[p] prob;
73+
int<lower=0, upper=1> z[p];
74+
75+
for (i in 1:p) {
76+
prob[i] <- exp(lp_parts[i, 1]) / sum(exp(lp_parts[i]));
77+
z[i] <- bernoulli_rng(prob[i]);
78+
}
79+
}"
80+
}
81+
82+
k <- c(21, 17, 21, 18, 22, 31, 31, 34, 34, 35, 35, 36, 39, 36, 35)
83+
p <- length(k) # number of people
84+
n <- 40 # number of questions
85+
86+
data <- list(p=p, k=k, n=n) # to be passed on to Stan
87+
88+
if (choice == 1) {
89+
myinits <- list(
90+
list(phi=.75, mix=rep(.5, p))) # Initial group assignment (for model 1)
91+
parameters <- c("phi") # parameters to be monitored:
92+
} else if (choice == 2) {
93+
myinits <- list(
94+
list(phi=.75, mix=rep(.5, p))) # Initial group assignment (for model 2)
95+
parameters <- c("phi", "z") # parameters to be monitored:
96+
}
97+
98+
# The following command calls Stan with specific options.
99+
# For a detailed description type "?rstan".
100+
samples <- stan(model_code=model,
101+
data=data,
102+
init=myinits, # If not specified, gives random inits
103+
pars=parameters,
104+
iter=20000,
105+
chains=1,
106+
thin=1,
107+
# warmup = 100, # Stands for burn-in; Don't set to 0 or low
108+
# ... values, it can malfunction; Default = iter/2
109+
# seed = 123 # Setting seed; Default is random seed
110+
)
111+
# Now the values for the monitored parameters are in the "samples" object,
112+
# ready for inspection.
113+
114+
print(samples, digits_summary=3)

‎Bayesian_Cognitive_Modeling/!underDevelopment/!TwoCountryQuiz_Stan.R‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
8+
data {
9+
int nx;
10+
int nz;
11+
int<lower=0, upper=1> k[nx, nz];
12+
}
13+
parameters {
14+
real<lower=0, upper=1> alpha; # Match
15+
real<lower=0, upper=alpha> beta; # Mismatch
16+
real<lower=0, upper=1> x_prob[nx];
17+
real<lower=0, upper=1> z_prob[nz];
18+
}
19+
transformed parameters {
20+
real lp_parts[4, nx, nz];
21+
22+
for (i in 1:nx) {
23+
for (j in 1:nz) {
24+
lp_parts[1, i, j] <- exp(bernoulli_log(0, x_prob[i]) +
25+
bernoulli_log(0, z_prob[j]) +
26+
bernoulli_log(k[i, j], alpha));
27+
lp_parts[2, i, j] <- exp(bernoulli_log(1, x_prob[i]) +
28+
bernoulli_log(1, z_prob[j]) +
29+
bernoulli_log(k[i, j], alpha));
30+
lp_parts[3, i, j] <- exp(bernoulli_log(0, x_prob[i]) +
31+
bernoulli_log(1, z_prob[j]) +
32+
bernoulli_log(k[i, j], beta));
33+
lp_parts[4, i, j] <- exp(bernoulli_log(1, x_prob[i]) +
34+
bernoulli_log(0, z_prob[j]) +
35+
bernoulli_log(k[i, j], beta));
36+
}
37+
}
38+
}
39+
model {
40+
for (i in 1:nx) {
41+
for (j in 1:nz) {
42+
increment_log_prob(log(lp_parts[1, i, j] + lp_parts[2, i, j] +
43+
lp_parts[3, i, j] + lp_parts[4, i, j]));
44+
}
45+
}
46+
}
47+
generated quantities {
48+
int x[nx];
49+
int z[nz];
50+
real prob_x[nx];
51+
real prob_z[nz];
52+
53+
for (i in 1:nx) {
54+
for (j in 1:nz) {
55+
prob_x[i] <- (sum(lp_parts[2, i]) + sum(lp_parts[4, i])) /
56+
(sum(lp_parts[1, i]) + sum(lp_parts[2, i]) +
57+
sum(lp_parts[3, i]) + sum(lp_parts[4, i]));
58+
}
59+
x[i] <- bernoulli_rng(prob_x[i]);
60+
}
61+
}"
62+
63+
64+
dset <- 1
65+
66+
if (dset==1) {
67+
k <- c(1,0,0,1,1,0,0,1,
68+
1,0,0,1,1,0,0,1,
69+
0,1,1,0,0,1,0,0,
70+
0,1,1,0,0,1,1,0,
71+
1,0,0,1,1,0,0,1,
72+
0,0,0,1,1,0,0,1,
73+
0,1,0,0,0,1,1,0,
74+
0,1,1,1,0,1,1,0)
75+
k <- matrix(k, nrow=8, byrow=T)
76+
}
77+
78+
79+
nx <- nrow(k)
80+
nz <- ncol(k)
81+
82+
data <- list(nx=nx, nz=nz, k=k) # To be passed on to Stan
83+
84+
myinits <-list(
85+
list(z=runif(nz), x=runif(nx), alpha=0.75, beta=0.25))
86+
87+
parameters <- c("alpha", "beta", "x") # Parameters to be monitored
88+
# parameters <- c("z", "x", "alpha", "beta") # Parameters to be monitored
89+
90+
# The following command calls Stan with specific options.
91+
# For a detailed description type "?rstan".
92+
samples <- stan(model_code=model,
93+
data=data,
94+
# init=myinits, # If not specified, gives random inits
95+
# pars=parameters,
96+
iter=10000,
97+
chains=1,
98+
thin=1,
99+
# warmup = 100, # Stands for burn-in; Don't set to 0 or
100+
# ... low values, it can malfunction; Default = iter/2
101+
# seed = 123 # Setting seed; Default is random seed
102+
)
103+
# Now the values for the monitored parameters are in the "samples" object,
104+
# ready for inspection.
105+
106+
print(samples, digits=3)

‎Bayesian_Cognitive_Modeling/!underDevelopment/!TwoCountryQuiz_Stan_2.R‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
8+
data {
9+
int nx;
10+
int nz;
11+
int<lower=0, upper=1> k[nx, nz];
12+
}
13+
parameters {
14+
real<lower=0, upper=1> alpha; # Match
15+
real<lower=0, upper=alpha> beta; # Mismatch
16+
real<lower=0, upper=1> x_prob[nx];
17+
real<lower=0, upper=1> z_prob[nz];
18+
}
19+
transformed parameters {
20+
vector[4] lp_parts[nx, nz];
21+
for (i in 1:nx) {
22+
for (j in 1:nz) {
23+
lp_parts[i, j, 1] <- log1m(x_prob[i]) + log1m(z_prob[j]) +
24+
bernoulli_log(k[i, j], alpha);
25+
lp_parts[i, j, 2] <- log(x_prob[i]) + log(z_prob[j]) +
26+
bernoulli_log(k[i, j], alpha);
27+
lp_parts[i, j, 3] <- log1m(x_prob[i]) + log(z_prob[j]) +
28+
bernoulli_log(k[i, j], beta);
29+
lp_parts[i, j, 4] <- log(x_prob[i]) + log1m(z_prob[j]) +
30+
bernoulli_log(k[i, j], beta);
31+
}
32+
}
33+
}
34+
model {
35+
for (i in 1:nx) {
36+
for (j in 1:nz) {
37+
increment_log_prob(log_sum_exp(lp_parts[i, j]));
38+
}
39+
}
40+
}
41+
generated quantities {
42+
int x[nx];
43+
int z[nz];
44+
real prob_x[nx];
45+
real prob_z[nz];
46+
47+
for (i in 1:nx) {
48+
for (j in 1:nz) {
49+
prob_x[i] <- (sum(exp(lp_parts[2, i])) + sum(exp(lp_parts[4, i]))) /
50+
(sum(exp(lp_parts[1, i])) + sum(exp(lp_parts[2, i])) +
51+
sum(exp(lp_parts[3, i])) + sum(exp(lp_parts[4, i])));
52+
}
53+
x[i] <- bernoulli_rng(prob_x[i]);
54+
}
55+
}"
56+
57+
58+
dset <- 1
59+
60+
if (dset==1) {
61+
k <- c(1,0,0,1,1,0,0,1,
62+
1,0,0,1,1,0,0,1,
63+
0,1,1,0,0,1,0,0,
64+
0,1,1,0,0,1,1,0,
65+
1,0,0,1,1,0,0,1,
66+
0,0,0,1,1,0,0,1,
67+
0,1,0,0,0,1,1,0,
68+
0,1,1,1,0,1,1,0)
69+
k <- matrix(k, nrow=8, byrow=T)
70+
}
71+
72+
73+
nx <- nrow(k)
74+
nz <- ncol(k)
75+
76+
data <- list(nx=nx, nz=nz, k=k) # To be passed on to Stan
77+
78+
myinits <-list(
79+
list(z_prob=runif(nz), x_prob=runif(nx), alpha=0.75, beta=0.25))
80+
81+
parameters <- c("alpha", "beta", "x") # Parameters to be monitored
82+
# parameters <- c("z", "x", "alpha", "beta") # Parameters to be monitored
83+
84+
# The following command calls Stan with specific options.
85+
# For a detailed description type "?rstan".
86+
samples <- stan(model_code=model,
87+
data=data,
88+
init=myinits, # If not specified, gives random inits
89+
pars=parameters,
90+
iter=10000,
91+
chains=1,
92+
thin=1,
93+
# warmup = 100, # Stands for burn-in; Don't set to 0 or
94+
# ... low values, it can malfunction; Default = iter/2
95+
# seed = 123 # Setting seed; Default is random seed
96+
)
97+
# Now the values for the monitored parameters are in the "samples" object,
98+
# ready for inspection.
99+
100+
print(samples, digits=3)

‎Bayesian_Cognitive_Modeling/GettingStarted/Rate_1_Stan.R‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
# When you work through the code for the first time,
2+
# execute each command one at a time to better understand
3+
# what it does.
4+
5+
# clears workspace:
6+
rm(list=ls(all=TRUE))
7+
8+
library(rstan)
9+
10+
model <- "
11+
// Inferring a Rate
12+
data {
13+
int<lower=1> n;
14+
int<lower=0> k;
15+
}
16+
parameters {
17+
real<lower=0,upper=1> theta;
18+
}
19+
model {
20+
// Prior Distribution for Rate Theta
21+
theta ~ beta(1, 1);
22+
23+
// Observed Counts
24+
k ~ binomial(n, theta);
25+
}"
26+
27+
k <- 5
28+
n <- 10
29+
30+
data <- list(k=k, n=n) # to be passed on to Stan
31+
32+
myinits <- list(
33+
list(theta=.1), # chain 1 starting value
34+
list(theta=.9)) # chain 2 starting value
35+
36+
# parameters to be monitored:
37+
parameters <- c("theta")
38+
39+
# The following command calls Stan with specific options.
40+
# For a detailed description type "?rstan".
41+
samples <- stan(model_code=model,
42+
data=data,
43+
init=myinits, # If not specified, gives random inits
44+
pars=parameters,
45+
iter=40000,
46+
chains=2,
47+
thin=1,
48+
# warmup = 100, # Stands for burn-in; Default = iter/2
49+
# seed = 123 # Setting seed; Default is random seed
50+
)
51+
# Now the values for the monitored parameters are in the "samples" object,
52+
# ready for inspection.
53+
54+
# The commands below are useful for a quick overview:
55+
print(samples) # a rough summary
56+
print(summary(samples)) # more detailed summary
57+
plot(samples) # a visual representation
58+
traceplot(samples) # traceplot
59+
60+
as.array(samples)[1:15,,2] # array: sample, chain, parameter
61+
62+
# Collect posterior and prior samples across all chains:
63+
theta <- extract(samples)$theta
64+
65+
# Now let's plot a histogram for theta.
66+
# NB. Some the plots will not look good in RStudio.
67+
# First, some options to make the plot look better:
68+
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
69+
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
70+
Nbreaks <- 80
71+
y <- hist(theta, Nbreaks, plot=F)
72+
plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
73+
xlim=c(0,1), xlab="Rate", ylab="Posterior Density")
74+

0 commit comments

Comments
 (0)