Skip to content

Commit 7017e58

Browse files
martin-smiramartin-smira
committed
duplicate all in one R-scripts to separate folder
1 parent 712427b commit 7017e58

File tree

74 files changed

+8283
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

74 files changed

+8283
-0
lines changed

‎Bayesian_Cognitive_Modeling/! rstan - all in one/CaseStudies/BART/BART_1_Stan.R‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
// BART Model of Risky Decision-Making
8+
data {
9+
int<lower=1> ntrials;
10+
real<lower=0,upper=1> p;
11+
int<lower=1> options[ntrials];
12+
int d[ntrials,30];
13+
}
14+
parameters {
15+
// Priors
16+
real<lower=0,upper=10> gplus;
17+
real<lower=0,upper=10> beta;
18+
}
19+
transformed parameters {
20+
real<lower=0> omega;
21+
22+
// Optimal Number of Pumps
23+
omega <- -gplus / log1m(p); // log1m(p) equals log(1 - p), but faster
24+
}
25+
model {
26+
// Choice Data
27+
for (j in 1:ntrials) {
28+
for (k in 1:options[j]) {
29+
real theta;
30+
theta <- 1 - inv_logit(-beta * (k - omega));
31+
d[j,k] ~ bernoulli(theta);
32+
}
33+
}
34+
}"
35+
36+
p <- .15 # (Belief of) bursting probability
37+
ntrials <- 90 # Number of trials for the BART
38+
39+
Data <- matrix(data=as.numeric(as.matrix(read.table("GeorgeSober.txt"))[-1, ]),
40+
ntrials, 8)
41+
d <- matrix(-99, ntrials, 30) # Data in binary format
42+
43+
cash <-(Data[, 7] != 0) * 1 # Cash or burst?
44+
npumps <- Data[, 6] # Nr. of pumps
45+
options <- cash + npumps # Nr. of decision possibilities
46+
47+
for (j in 1:ntrials) {
48+
if (npumps[j] > 0) {d[j, 1:npumps[j]] <- rep(0, npumps[j])}
49+
if (cash[j] == 1) {d[j, (npumps[j] + 1)] <- 1}
50+
}
51+
52+
# To be passed on to Stan:
53+
data <- list(ntrials=ntrials, p=p, options=options, d=d)
54+
55+
myinits <- list(
56+
list(gplus=1.2, beta=.5))
57+
58+
parameters <- c("gplus", "beta") # Parameters to be monitored
59+
60+
# The following command calls Stan with specific options.
61+
# For a detailed description type "?rstan".
62+
samples <- stan(model_code=model,
63+
data=data,
64+
init=myinits, # If not specified, gives random inits
65+
pars=parameters,
66+
iter=5000,
67+
chains=1,
68+
thin=1,
69+
warmup=2000, # Stands for burn-in; Default = iter/2
70+
# seed=123 # Setting seed; Default is random seed
71+
)
72+
# Now the values for the monitored parameters are in the "samples" object,
73+
# ready for inspection.
74+
75+
gplus <- extract(samples)$gplus
76+
beta <- extract(samples)$beta
77+
78+
#################### PLOT RESULTS
79+
80+
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
81+
layout (matrix (1:3, 1, 3))
82+
83+
hist(npumps, main = " ", xlab="Number of Pumps", ylab="Frequency", breaks=c(0:7),
84+
col="lightblue", axes=F)
85+
axis(2)
86+
axis(1, at=seq(.5,6.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7"))
87+
88+
plot(density(gplus), xlab = expression (gamma^'+'),
89+
main = " ", bty = 'n', lwd=2, ylab="Posterior Density")
90+
plot(density(beta), xlab = expression (beta),
91+
main = " ", bty = 'n', lwd=2, lab = c(5,3,5), ylab="Posterior Density")
92+

‎Bayesian_Cognitive_Modeling/! rstan - all in one/CaseStudies/BART/BART_2_Stan.R‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
// BART Model of Risky Decision-Making
8+
data {
9+
int<lower=1> nconds;
10+
int<lower=1> ntrials;
11+
real<lower=0,upper=1> p;
12+
int<lower=1> options[nconds,ntrials];
13+
int d[nconds,ntrials,30];
14+
}
15+
parameters {
16+
vector<lower=0>[nconds] gplus;
17+
vector<lower=0>[nconds] beta;
18+
19+
// Priors
20+
real<lower=0,upper=10> mug;
21+
real<lower=0,upper=10> sigmag;
22+
real<lower=0,upper=10> mub;
23+
real<lower=0,upper=10> sigmab;
24+
}
25+
transformed parameters {
26+
vector<lower=0>[nconds] omega;
27+
28+
// Optimal Number of Pumps
29+
omega <- -gplus / log1m(p);
30+
}
31+
model {
32+
for (i in 1:nconds) {
33+
gplus[i] ~ normal(mug, sigmag)T[0,];
34+
beta[i] ~ normal(mub, sigmab)T[0,];
35+
}
36+
// Choice Data
37+
for (i in 1:nconds) {
38+
for (j in 1:ntrials) {
39+
for (k in 1:options[i,j]) {
40+
real theta;
41+
theta <- 1 - inv_logit(-beta[i] * (k - omega[i]));
42+
d[i,j,k] ~ bernoulli(theta);
43+
}
44+
}
45+
}
46+
}"
47+
48+
p <- .15 # (Belief of) bursting probability
49+
ntrials <- 90 # Number of trials for the BART
50+
51+
################### READ IN THE DATA
52+
53+
Data <- list(
54+
matrix(data=as.numeric(as.matrix(read.table("GeorgeSober.txt"))[-1, ]),
55+
ntrials, 8),
56+
matrix(data=as.numeric(as.matrix(read.table("GeorgeTipsy.txt"))[-1, ]),
57+
ntrials, 8),
58+
matrix(data=as.numeric(as.matrix(read.table("GeorgeDrunk.txt"))[-1, ]),
59+
ntrials, 8)
60+
)
61+
62+
nconds <- length(Data)
63+
cash <- npumps <- matrix (, nconds, ntrials) # Cashes and nr. of pumps
64+
d <- array (-99, c(nconds, ntrials, 30)) # Data in binary format
65+
66+
for (i in 1:nconds) {
67+
cash[i, ] <- (Data[[i]][, 7] != 0) * 1 # Cash or burst?
68+
npumps[i, ] <- Data[[i]][, 6] # Nr. of pumps
69+
70+
for (j in 1:ntrials) {
71+
if (npumps[i, j] > 0) {d[i, j, 1:npumps[i, j]] <- rep(0, npumps[i, j])}
72+
if (cash[i, j] == 1) {d[i, j, (npumps[i, j] + 1)] <- 1}
73+
}
74+
}
75+
options <- cash + npumps # Nr. of decision possibilities
76+
77+
# To be passed on to Stan:
78+
data <- list(nconds=nconds, ntrials=ntrials, p=p, options=options, d=d)
79+
80+
myinits <- list(
81+
list(gplus=rep(1.2, 3), beta=rep(.5, 3), mug=1.2, sigmag=.1, mub=.8, sigmab=.8),
82+
list(gplus=rep(1.2, 3), beta=rep(.5, 3), mug=1.5, sigmag=.2, mub=1, sigmab=1.2))
83+
84+
# Parameters to be monitored:
85+
parameters <- c("beta", "gplus", "mub", "mug", "sigmab", "sigmag")
86+
87+
# The following command calls Stan with specific options.
88+
# For a detailed description type "?rstan".
89+
samples <- stan(model_code=model,
90+
data=data,
91+
init=myinits, # If not specified, gives random inits
92+
pars=parameters,
93+
iter=1500,
94+
chains=2,
95+
thin=1,
96+
warmup=500, # Stands for burn-in; Default = iter/2
97+
seed=1234 # Setting seed; Default is random seed
98+
)
99+
# Now the values for the monitored parameters are in the "samples" object,
100+
# ready for inspection.
101+
102+
gplus.sober <- extract(samples)$gplus[,1]
103+
gplus.tipsy <- extract(samples)$gplus[,2]
104+
gplus.drunk <- extract(samples)$gplus[,3]
105+
106+
beta.sober <- extract(samples)$beta[,1]
107+
beta.tipsy <- extract(samples)$beta[,2]
108+
beta.drunk <- extract(samples)$beta[,3]
109+
110+
#################### PLOT SOME RESULTS
111+
windows()
112+
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
113+
layout (matrix (1:9, 3, 3, byrow = T))
114+
115+
hist(npumps[1,], xlab=" ", main = "Sober: # pumps", breaks=c(0:max(npumps[1,])),
116+
xlim = range(c(0:9)), col="lightblue", axes=F)
117+
axis(2)
118+
axis(1, at=seq(.5,8.5,by=1),
119+
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
120+
121+
plot(density(gplus.sober), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
122+
main = expression (paste ("Sober: Posterior ", gamma^'+')), bty = 'n')
123+
plot (density(beta.sober), xlab = expression (beta), bty = 'n',
124+
main = expression (paste ("Sober: Posterior ", beta)), xlim = c(0.2,1.4))
125+
126+
hist(npumps[2,], xlab=" ", main = "Tipsy: # pumps", breaks=c(0:max(npumps[2,])),
127+
xlim = range(c(0:9)), col="lightblue", axes=F)
128+
axis(2)
129+
axis(1, at=seq(.5,8.5,by=1),
130+
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
131+
plot(density (gplus.tipsy), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
132+
main = expression (paste ("Tipsy: Posterior ", gamma^'+')), bty = 'n')
133+
plot(density (beta.tipsy), xlab = expression (beta), xlim = c(0.2,1.4),
134+
main = expression (paste ("Tipsy: Posterior ", beta)), bty = 'n')
135+
136+
hist(npumps[3,], xlab="Number of Pumps", main = "Drunk: # pumps",
137+
breaks=c(0:max(npumps[3,])), xlim = range(c(0:9)), col="lightblue", axes=F)
138+
axis(2)
139+
axis(1, at=seq(.5,8.5,by=1),
140+
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
141+
plot(density(gplus.drunk), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
142+
main = expression(paste ("Drunk: Posterior ", gamma^'+')), bty = 'n')
143+
plot(density(beta.drunk), xlab = expression (beta), xlim = c(0.2,1.4),
144+
main = expression(paste ("Drunk: Posterior ", beta)), bty = 'n')
145+

‎Bayesian_Cognitive_Modeling/! rstan - all in one/CaseStudies/BART/BillDrunk.txt‎

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
pres.bl block gr.fact prob. trial pumps cash total
2+
1 3 0.25 20 1 4 0.99 0.99
3+
1 3 0.25 20 2 1 0 0.99
4+
1 3 0.25 20 3 1 0 0.99
5+
1 3 0.25 20 4 4 0.99 1.98
6+
1 3 0.25 20 5 5 1.24 3.22
7+
1 3 0.25 20 6 1 0 3.22
8+
1 3 0.25 20 7 5 1.24 4.46
9+
1 3 0.25 20 8 6 1.55 6.01
10+
1 3 0.25 20 9 5 1.24 7.25
11+
1 3 0.25 20 10 6 0 7.25
12+
1 3 0.25 20 11 5 0 7.25
13+
1 3 0.25 20 12 5 1.24 8.49
14+
1 3 0.25 20 13 4 0 8.49
15+
1 3 0.25 20 14 4 0 8.49
16+
1 3 0.25 20 15 6 0 8.49
17+
1 3 0.25 20 16 2 0 8.49
18+
1 3 0.25 20 17 5 0 8.49
19+
1 3 0.25 20 18 1 0 8.49
20+
1 3 0.25 20 19 1 0 8.49
21+
1 3 0.25 20 20 6 1.55 10.04
22+
1 3 0.25 20 21 5 1.24 11.28
23+
1 3 0.25 20 22 1 0 11.28
24+
1 3 0.25 20 23 1 0.5 11.78
25+
1 3 0.25 20 24 2 0 11.78
26+
1 3 0.25 20 25 3 0 11.78
27+
1 3 0.25 20 26 3 0 11.78
28+
1 3 0.25 20 27 4 0 11.78
29+
1 3 0.25 20 28 8 2.42 14.2
30+
1 3 0.25 20 29 3 0 14.2
31+
1 3 0.25 20 30 11 0 14.2
32+
2 2 0.1765 15 1 5 0 0
33+
2 2 0.1765 15 2 6 1.12 1.12
34+
2 2 0.1765 15 3 1 0 1.12
35+
2 2 0.1765 15 4 5 0 1.12
36+
2 2 0.1765 15 5 5 0.95 2.07
37+
2 2 0.1765 15 6 2 0 2.07
38+
2 2 0.1765 15 7 1 0.5 2.57
39+
2 2 0.1765 15 8 2 0.59 3.16
40+
2 2 0.1765 15 9 3 0.69 3.85
41+
2 2 0.1765 15 10 1 0 3.85
42+
2 2 0.1765 15 11 4 0.81 4.66
43+
2 2 0.1765 15 12 3 0 4.66
44+
2 2 0.1765 15 13 5 0 4.66
45+
2 2 0.1765 15 14 2 0 4.66
46+
2 2 0.1765 15 15 1 0 4.66
47+
2 2 0.1765 15 16 13 3.48 8.14
48+
2 2 0.1765 15 17 15 4.81 12.95
49+
2 2 0.1765 15 18 1 0 12.95
50+
2 2 0.1765 15 19 2 0 12.95
51+
2 2 0.1765 15 20 1 0.5 13.45
52+
2 2 0.1765 15 21 4 0 13.45
53+
2 2 0.1765 15 22 4 0 13.45
54+
2 2 0.1765 15 23 10 0 13.45
55+
2 2 0.1765 15 24 1 0 13.45
56+
2 2 0.1765 15 25 7 0 13.45
57+
2 2 0.1765 15 26 6 1.12 14.57
58+
2 2 0.1765 15 27 1 0 14.57
59+
2 2 0.1765 15 28 1 0.5 15.07
60+
2 2 0.1765 15 29 6 0 15.07
61+
2 2 0.1765 15 30 4 0 15.07
62+
3 1 0.1111 10 1 1 0 0
63+
3 1 0.1111 10 2 4 0 0
64+
3 1 0.1111 10 3 1 0.5 0.5
65+
3 1 0.1111 10 4 1 0.5 1
66+
3 1 0.1111 10 5 6 0 1
67+
3 1 0.1111 10 6 5 0.77 1.77
68+
3 1 0.1111 10 7 6 0.86 2.63
69+
3 1 0.1111 10 8 3 0 2.63
70+
3 1 0.1111 10 9 2 0 2.63
71+
3 1 0.1111 10 10 1 0.5 3.13
72+
3 1 0.1111 10 11 2 0.56 3.69
73+
3 1 0.1111 10 12 12 1.63 5.32
74+
3 1 0.1111 10 13 3 0.62 5.94
75+
3 1 0.1111 10 14 3 0 5.94
76+
3 1 0.1111 10 15 14 2.01 7.95
77+
3 1 0.1111 10 16 5 0.77 8.72
78+
3 1 0.1111 10 17 11 0 8.72
79+
3 1 0.1111 10 18 4 0 8.72
80+
3 1 0.1111 10 19 3 0 8.72
81+
3 1 0.1111 10 20 3 0 8.72
82+
3 1 0.1111 10 21 3 0 8.72
83+
3 1 0.1111 10 22 2 0 8.72
84+
3 1 0.1111 10 23 1 0.5 9.22
85+
3 1 0.1111 10 24 1 0.5 9.72
86+
3 1 0.1111 10 25 1 0.5 10.22
87+
3 1 0.1111 10 26 3 0.62 10.84
88+
3 1 0.1111 10 27 5 0.77 11.61
89+
3 1 0.1111 10 28 7 0 11.61
90+
3 1 0.1111 10 29 1 0 11.61
91+
3 1 0.1111 10 30 6 0.86 12.47

0 commit comments

Comments
 (0)