Skip to content

Commit 42c7f51

Browse files
martin-smiramartin-smira
committed
add CaseStudies models
1 parent a85ab98 commit 42c7f51

38 files changed

+4250
-7
lines changed

‎Bayesian_Cognitive_Modeling/!underDevelopment/!Ability_Stan_old.R‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
# clears workspace:
2+
rm(list=ls(all=TRUE))
3+
4+
library(rstan)
5+
6+
model <- "
7+
# Ability Correlation for ESP Replication
8+
data {
9+
int<lower=1> nsubjs;
10+
int<lower=1> ntrials;
11+
int<lower=0> k[nsubjs, 2];
12+
}
13+
parameters {
14+
vector[2] mu;
15+
vector<lower=0>[2] lambda;
16+
vector[2] thetap[nsubjs];
17+
real<lower=0, upper=1> r;
18+
}
19+
transformed parameters {
20+
vector<lower=0>[2] sigma;
21+
vector[2] theta[nsubjs];
22+
cov_matrix[2] T;
23+
24+
# Reparametrization
25+
sigma[1] <- 1 / (sqrt(lambda[1]));
26+
sigma[2] <- 1 / (sqrt(lambda[2]));
27+
28+
T[1, 1] <- 1 / lambda[1];
29+
T[1, 2] <- r * sigma[1] * sigma[2];
30+
T[2, 1] <- r * sigma[1] * sigma[2];
31+
T[2, 2] <- 1 / lambda[2];
32+
33+
for (i in 1:nsubjs) {
34+
for (j in 1:2) {
35+
theta[i, j] <- Phi(thetap[i, j]);
36+
}
37+
}
38+
}
39+
model {
40+
# Priors
41+
mu ~ normal(0, sqrt(1000));
42+
lambda ~ gamma(.001, .001);
43+
44+
# Data
45+
for (i in 1:nsubjs) {
46+
thetap[i] ~ multi_normal(mu, T);
47+
for (j in 1:2) {
48+
k[i, j] ~ binomial(ntrials, theta[i, j]);
49+
}
50+
}
51+
}"
52+
53+
54+
# Proportion correct on erotic pictures, block 1 and block 2:
55+
prc1.ero <- c(0.6000000, 0.5333333, 0.6000000, 0.6000000, 0.4666667,
56+
0.6666667, 0.6666667, 0.4000000, 0.6000000, 0.6000000,
57+
0.4666667, 0.6666667, 0.4666667, 0.6000000, 0.3333333,
58+
0.4000000, 0.4000000, 0.2666667, 0.3333333, 0.5333333,
59+
0.6666667, 0.5333333, 0.6000000, 0.4000000, 0.4666667,
60+
0.7333333, 0.6666667, 0.6000000, 0.6666667, 0.5333333,
61+
0.5333333, 0.6666667, 0.4666667, 0.3333333, 0.4000000,
62+
0.5333333, 0.4000000, 0.4000000, 0.3333333, 0.4666667,
63+
0.4000000, 0.4666667, 0.4666667, 0.5333333, 0.3333333,
64+
0.7333333, 0.2666667, 0.6000000, 0.5333333, 0.4666667,
65+
0.4000000, 0.5333333, 0.6666667, 0.4666667, 0.5333333,
66+
0.5333333, 0.4666667, 0.4000000, 0.4666667, 0.6666667,
67+
0.4666667, 0.3333333, 0.3333333, 0.3333333, 0.4000000,
68+
0.4000000, 0.6000000, 0.4666667, 0.3333333, 0.3333333,
69+
0.6666667, 0.5333333, 0.3333333, 0.6000000, 0.4666667,
70+
0.4666667, 0.4000000, 0.3333333, 0.4666667, 0.5333333,
71+
0.8000000, 0.4000000, 0.5333333, 0.5333333, 0.6666667,
72+
0.6666667, 0.6666667, 0.6000000, 0.6000000, 0.5333333,
73+
0.3333333, 0.4666667, 0.6666667, 0.5333333, 0.3333333,
74+
0.3333333, 0.2666667, 0.2666667, 0.4666667, 0.6666667)
75+
76+
prc2.ero <- c(0.3333333, 0.6000000, 0.5333333, 0.2666667, 0.6666667,
77+
0.5333333, 0.6666667, 0.4666667, 0.4666667, 0.6666667,
78+
0.4000000, 0.6666667, 0.2666667, 0.4000000, 0.4666667,
79+
0.3333333, 0.5333333, 0.6000000, 0.3333333, 0.4000000,
80+
0.4666667, 0.4666667, 0.6000000, 0.5333333, 0.5333333,
81+
0.6000000, 0.5333333, 0.6666667, 0.6000000, 0.2666667,
82+
0.4666667, 0.4000000, 0.6000000, 0.5333333, 0.4000000,
83+
0.4666667, 0.5333333, 0.3333333, 0.4000000, 0.4666667,
84+
0.8000000, 0.6000000, 0.2000000, 0.6000000, 0.4000000,
85+
0.4000000, 0.2666667, 0.2666667, 0.6000000, 0.4000000,
86+
0.4000000, 0.4000000, 0.4000000, 0.4000000, 0.6666667,
87+
0.7333333, 0.5333333, 0.5333333, 0.3333333, 0.6000000,
88+
0.5333333, 0.5333333, 0.4666667, 0.5333333, 0.4666667,
89+
0.5333333, 0.4000000, 0.4000000, 0.4666667, 0.6000000,
90+
0.6000000, 0.6000000, 0.4666667, 0.6000000, 0.6666667,
91+
0.5333333, 0.4666667, 0.6000000, 0.2000000, 0.5333333,
92+
0.4666667, 0.4000000, 0.5333333, 0.5333333, 0.5333333,
93+
0.5333333, 0.6000000, 0.6666667, 0.4000000, 0.4000000,
94+
0.5333333, 0.8000000, 0.6000000, 0.4000000, 0.2000000,
95+
0.6000000, 0.6666667, 0.4666667, 0.4666667, 0.4666667)
96+
97+
k <- matrix(cbind(as.integer(round(prc1.ero * 60)),
98+
as.integer(round(prc2.ero * 60))), nrow=100)
99+
nsubjs <- nrow(k) # number of participants
100+
ntrials <- 60
101+
102+
103+
data <- list(k=k, nsubjs=nsubjs, ntrials=ntrials) # To be passed on to Stan
104+
105+
# myinits <- list(
106+
# list(r=0, mu=c(0, 0), lambda=c(1, 1),
107+
# thetap=matrix(rnorm(200), 100, 2)))
108+
109+
parameters <- c("r", "mu", "sigma") # Parameters to be monitored
110+
111+
# The following command calls Stan with specific options.
112+
# For a detailed description type "?rstan".
113+
samples <- stan(model_code=model,
114+
data=data,
115+
init=myinits, # If not specified, gives random inits
116+
pars=parameters,
117+
iter=1000,
118+
chains=1,
119+
thin=1,
120+
# warmup = 100, # Stands for burn-in; Don't set to 0 or
121+
# ... low values, it can malfunction; Default = iter/2
122+
# seed = 123 # Setting seed; Default is random seed
123+
)
124+
# Now the values for the monitored parameters are in the "samples" object,
125+
# ready for inspection.
126+
127+
print(samples, digits_summary=2)

‎Bayesian_Cognitive_Modeling/!underDevelopment/!PsychophysicalFunction1_Stan (2 in 1).R‎

Lines changed: 236 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,236 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
###########################################################
7+
# Choose a model: 1 or 2 - the second one is simplified
8+
modelChoice <- 2
9+
###########################################################
10+
11+
if (modelChoice == 1) {
12+
model <- "
13+
# Logistic Psychophysical Function
14+
data {
15+
int nsubjs;
16+
int nstim[nsubjs];
17+
int n[nsubjs, 28];
18+
int r[nsubjs, 28];
19+
int x[nsubjs, 28];
20+
vector[nsubjs] xmean;
21+
}
22+
parameters {
23+
real mua;
24+
real mub;
25+
real<lower=0,upper=1000> sigmaa;
26+
real<lower=0,upper=1000> sigmab;
27+
vector[nsubjs] alpha;
28+
vector[nsubjs] beta;
29+
}
30+
model {
31+
# Priors
32+
mua ~ normal(0, sqrt(1000));
33+
mub ~ normal(0, sqrt(1000));
34+
35+
alpha ~ normal(mua, sigmaa);
36+
beta ~ normal(mub, sigmab);
37+
38+
for (i in 1:nsubjs) {
39+
for (j in 1:nstim[i]) {
40+
41+
real thetalim;
42+
real lthetalim;
43+
real ltheta;
44+
45+
ltheta <- alpha[i] + beta[i] * (x[i, j] - xmean[i]);
46+
47+
if (ltheta < -999)
48+
lthetalim <- -999;
49+
else if (ltheta > 999)
50+
lthetalim <- 999;
51+
else
52+
lthetalim <- ltheta;
53+
54+
thetalim <- inv_logit(lthetalim);
55+
56+
r[i, j] ~ binomial(n[i, j], thetalim);
57+
}
58+
}
59+
}"
60+
}
61+
if (modelChoice == 2) {
62+
model <- "
63+
# Logistic Psychophysical Function
64+
data {
65+
int nsubjs;
66+
int nstim[nsubjs];
67+
int n[nsubjs,28];
68+
int r[nsubjs,28];
69+
int x[nsubjs,28];
70+
vector[nsubjs] xmean;
71+
}
72+
parameters {
73+
real mua;
74+
real mub;
75+
real<lower=0,upper=1000> sigmaa;
76+
real<lower=0,upper=1000> sigmab;
77+
vector[nsubjs] alpha;
78+
vector[nsubjs] beta;
79+
}
80+
model {
81+
real theta;
82+
83+
# Priors
84+
mua ~ normal(0, inv_sqrt(.001));
85+
mub ~ normal(0, inv_sqrt(.001));
86+
87+
alpha ~ normal(mua, sigmaa);
88+
beta ~ normal(mub, sigmab);
89+
90+
for (i in 1:nsubjs) {
91+
for (j in 1:nstim[i]) {
92+
theta <- inv_logit(alpha[i] + beta[i] * (x[i,j] - xmean[i]));
93+
r[i,j] ~ binomial(n[i,j], theta);
94+
}
95+
}
96+
}"
97+
}
98+
99+
x <- as.matrix(read.table("data_x.txt", sep="\t"))
100+
x[is.na(x)] = -5 # transforming because Stan can't handle NAs
101+
102+
n <- as.matrix(read.table("data_n.txt", sep="\t"))
103+
n[is.na(n)] = -5
104+
105+
r <- as.matrix(read.table("data_r.txt", sep="\t"))
106+
r[is.na(r)] = -5
107+
108+
rprop <- as.matrix(read.table("data_rprop.txt", sep="\t"))
109+
110+
xmean <- c(318.888, 311.0417, 284.4444, 301.5909,
111+
296.2000, 305.7692, 294.6429, 280.3571)
112+
nstim <- c(27, 24, 27, 22, 25, 26, 28, 28)
113+
nsubjs <- 8
114+
115+
# to be passed on to Stan
116+
data <- list(x=x, xmean=xmean, n=n, r=r, nsubjs=nsubjs, nstim=nstim)
117+
118+
# myinits <- list( # Doesn't work with this initial list
119+
# list(alpha=runif(nsubjs, -2, 2), beta=runif(nsubjs, 0, .5),
120+
# mua=0, mub=0, sigmaa=1, sigmab=1),
121+
# list(alpha=runif(nsubjs, -2, 2), beta=runif(nsubjs, 0, .5),
122+
# mua=0, mub=0, sigmaa=1, sigmab=1),
123+
# list(alpha=runif(nsubjs, -2, 2), beta=runif(nsubjs, 0, .5),
124+
# mua=0, mub=0, sigmaa=1, sigmab=1))
125+
126+
parameters <- c("alpha", "beta") # Parameters to be monitored
127+
128+
# The following command calls Stan with specific options.
129+
# For a detailed description type "?rstan".
130+
samples <- stan(model_code=model,
131+
data=data,
132+
init=0, # !!! has to be set on 0 for this example,
133+
pars=parameters,
134+
iter=16000,
135+
chains=1,
136+
thin=10,
137+
# warmup = 100, # Stands for burn-in; Default = iter/2
138+
# seed = 123 # Setting seed; Default is random seed
139+
)
140+
# Now the values for the monitored parameters are in the "samples" object,
141+
# ready for inspection.
142+
143+
x[x == -5] = NA # transforming back to NAs
144+
n[n == -5] = NA
145+
r[r == -5] = NA
146+
147+
# Extracting the parameters
148+
alpha = extract(samples)$alpha
149+
beta = extract(samples)$beta
150+
alphaMAP = c(rep(0,nsubjs))
151+
betaMAP = c(rep(0,nsubjs))
152+
alpha_sel = matrix(NA,20,8)
153+
beta_sel = matrix(NA,20,8)
154+
155+
# Constructing MAP-estimates and alpha/beta range
156+
for (i in 1:nsubjs)
157+
{
158+
alphaMAP[i] <- density(alpha[,i])$x[which(density(alpha[,i])$y ==
159+
max(density(alpha[,i])$y))]
160+
betaMAP[i] <- density(beta[,i])$x[which(density(beta[,i])$y ==
161+
max(density(beta[,i])$y))]
162+
alpha_sel[,i] <- sample(alpha[,i],20)
163+
beta_sel[,i] <- sample(beta[,i],20)
164+
}
165+
166+
############################## PSYCHOMETRIC FUNCTIONS ##########################
167+
168+
# only the MAP estimate; use this to plot psychometric functions
169+
F1 <- function(X,s)
170+
{
171+
exp(alphaMAP[s] + betaMAP[s]*(X - xmean[s]))/
172+
(1+exp(alphaMAP[s] + betaMAP[s]*(X - xmean[s])))
173+
}
174+
175+
F1inv <- function(Y,s)
176+
{
177+
(log(-Y/(Y-1))-alphaMAP[s])/betaMAP[s]
178+
}
179+
180+
# function for all the posterior alpha/beta values; use this to calculate JND
181+
# posterior
182+
F2 <- function(X,s)
183+
{
184+
exp(alpha[,s] + beta[,s]*(X - xmean[s]))/
185+
(1+exp(alpha[,s] + beta[,s]*(X - xmean[s])))
186+
}
187+
F2inv <- function(Y,s)
188+
{
189+
(log(-Y/(Y-1))-alpha[,s])/beta[,s]
190+
}
191+
192+
# function for 20 grabbed posterior alpha/beta values; use this to plot
193+
# overlapping sigmoids to visualize variance
194+
F3 <- function(X,s,g)
195+
{
196+
exp(alpha_sel[g,s] + beta_sel[g,s]*(X - xmean[s]))/
197+
(1+exp(alpha_sel[g,s] + beta_sel[g,s]*(X - xmean[s])))
198+
}
199+
200+
##################################### JND/PSE calculation #####################
201+
JND <- F2inv(0.84,c(1:nsubjs))-F2inv(0.5,c(1:nsubjs))
202+
JNDmap <- F1inv(0.84,c(1:nsubjs))-F1inv(0.5,c(1:nsubjs))
203+
204+
PSE <- F2inv(0.5,c(1:nsubjs))+xmean
205+
PSEmap <- F1inv(0.5,c(1:nsubjs))+xmean
206+
207+
################## PLOTS ####################
208+
209+
### Figure 12.2
210+
211+
dev.new(width=10,height=5)
212+
layout(matrix(1:nsubjs,2,4,byrow=T))
213+
par(mar=c(1,2,2,0),oma=c(5,5,1,1))
214+
for (i in 1:nsubjs)
215+
{
216+
scale <- seq(x[i,1],x[i,nstim[i]], by=.1)
217+
plot(x[i,],rprop[i,],main=paste("Subject",as.character(i)),xlab="",ylab="",
218+
pch=15,col="dark grey",ylim=c(0,1),yaxt="n",xaxt="n")
219+
lines(scale,F1(scale,i),type="l")
220+
segments(x0=x[i,1],x1=PSEmap[i]+JNDmap[i],y0=0.84,lty=2)
221+
segments(x0=x[i,1],x1=PSEmap[i],y0=0.5,lty=2)
222+
segments(y0=0,y1=0.84,x0=PSEmap[i]+JNDmap[i],lty=2)
223+
segments(y0=0,y1=0.5,x0=PSEmap[i],lty=2)
224+
if (i==1 | i==5)
225+
{
226+
axis(2,las=1,yaxp=c(0,1,2))
227+
axis(2,at=0.84,las=1)
228+
}
229+
if (i>4) axis(1)
230+
}
231+
mtext("Proportion 'Long' Response",side=2,line=2,outer=T,cex=1.4)
232+
mtext("Test Interval (ms)",side=1,outer=T,line=3,cex=1.4)
233+
234+
### WARNING: Do not close R window.
235+
236+
# NOTE: Answers to the exercises can be found in PsychometricFunction1_Answers.R

0 commit comments

Comments
 (0)