Skip to content

Commit 3b7ee0a

Browse files
idoba2030idoba2030
committed
confirmatory_two_alphas
1 parent 9495114 commit 3b7ee0a

File tree

9 files changed

+814
-0
lines changed

9 files changed

+814
-0
lines changed

‎functions/working_model.rdata‎

13 Bytes
Binary file not shown.

‎stan_modeling/models/confirmatory_two_alphas/README.html‎

Lines changed: 419 additions & 0 deletions
Large diffs are not rendered by default.

‎stan_modeling/models/confirmatory_two_alphas/README.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
## Model two alphas
2+
This is a model which has two separate alpha parameters depending on whether the prediction error was positive or negative.
3+
4+
Read ["The computational roots of positivity and
5+
confirmation biases in reinforcement learning"](https://www.sciencedirect.com/science/article/pii/S1364661322000894)
6+
for more information.
7+
8+
### Parameter recovery
9+
![Blue dashed lines indicate the true population parameter, black dot the posterior median, and black lines the 89% and 97% CI](parameter_recovery.png)

‎stan_modeling/models/confirmatory_two_alphas/confirmatory_two_alphas.R‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
#### simulate Rescorla-Wagner block for participant ----
2+
sim.block = function(subject,parameters,cfg){
3+
print(paste('subject',subject))
4+
5+
#pre-allocation
6+
7+
#set parameters
8+
alpha_confirmatory = parameters['alpha_confirmatory']
9+
alpha_disconfirmatory = parameters['alpha_disconfirmatory']
10+
beta = parameters['beta']
11+
12+
13+
#set initial var
14+
Narms = cfg$Narms
15+
Nraffle = cfg$Nraffle
16+
Nblocks = cfg$Nblocks
17+
Ntrials_perblock = cfg$Ntrials_perblock
18+
expvalues = cfg$rndwlk
19+
rownames(expvalues)=c('ev1','ev2','ev3','ev4')
20+
Qval = as.matrix(t(rep(0.5,Narms)))
21+
colnames(Qval) =sapply(1:Narms, function(n) {paste('Qbandit',n,sep="")})
22+
df =data.frame()
23+
24+
for (block in 1:Nblocks){
25+
26+
Qval = as.matrix(t(rep(0.5,Narms)))
27+
28+
for (trial in 1:Ntrials_perblock){
29+
30+
#computer offer
31+
raffle = sample(1:Narms,Nraffle,prob=rep(1/Narms,Narms))
32+
#raffle = sort(raffle)
33+
34+
#players choice
35+
p = exp(beta*Qval[raffle]) / sum(exp(beta*Qval[raffle]))
36+
choice = sample(raffle,1,prob=p)
37+
unchosen = raffle[choice!=raffle]
38+
39+
#outcome
40+
reward = sample(0:1,1,prob=c(1-expvalues[choice,trial],expvalues[choice,trial]))
41+
42+
#save trial's data
43+
44+
#create data for current trials
45+
dfnew=data.frame(
46+
subject = subject,
47+
block = block,
48+
trial = trial,
49+
first_trial_in_block = (trial==1)*1,
50+
choice = choice,
51+
selected_offer = (choice==raffle[2])*1,
52+
unchosen = unchosen,
53+
offer1 = raffle[1],
54+
offer2 = raffle[2],
55+
expval_ch = expvalues[choice,trial],
56+
expval_unch = expvalues[raffle[choice!=raffle],trial],
57+
reward = reward
58+
)
59+
60+
dfnew=cbind(dfnew,Qval)
61+
dfnew=cbind(dfnew,t(t(expvalues)[trial,]))
62+
63+
#bind to the overall df
64+
df=rbind(df,dfnew)
65+
66+
67+
68+
#updating Qvalues
69+
if(reward==1){
70+
Qval[choice] = Qval[choice] + alpha_confirmatory*(reward - Qval[choice])
71+
Qval[unchosen] = Qval[unchosen] + alpha_confirmatory*((1-reward) - Qval[unchosen])
72+
}
73+
else{
74+
Qval[choice] = Qval[choice] + alpha_disconfirmatory*(reward - Qval[choice])
75+
Qval[unchosen] = Qval[unchosen] + alpha_disconfirmatory*((1-reward) - Qval[unchosen])
76+
}
77+
}
78+
}
79+
return (df)
80+
}

‎stan_modeling/models/confirmatory_two_alphas/confirmatory_two_alphas.exe‎

5.14 MB
Binary file not shown.

‎stan_modeling/models/confirmatory_two_alphas/confirmatory_two_alphas.stan‎

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
data {
2+
3+
//General fixed parameters for the experiment/models
4+
int<lower = 1> Nsubjects;
5+
int<lower = 1> Nblocks;
6+
int<lower = 1> Ntrials;
7+
int<lower = 1> Ntrials_per_subject[Nsubjects];
8+
int<lower = 4> Narms;
9+
int<lower = 2> Nraffle;
10+
11+
12+
//Behavioral data:
13+
int<lower = 0> choice[Nsubjects,Ntrials];
14+
int<lower = 0> unchosen[Nsubjects,Ntrials];
15+
int<lower = 0> reward[Nsubjects,Ntrials];
16+
int<lower = 0> offer1[Nsubjects,Ntrials];
17+
int<lower = 0> offer2[Nsubjects,Ntrials];
18+
int<lower = 0> selected_offer[Nsubjects,Ntrials];
19+
int<lower = 0> first_trial_in_block[Nsubjects,Ntrials];
20+
21+
}
22+
23+
24+
transformed data{
25+
int<lower = 1> Nparameters=3;
26+
vector[Narms] Qvalue_initial;
27+
Qvalue_initial = rep_vector(0.5, Narms);
28+
}
29+
30+
31+
32+
33+
parameters {
34+
//population level parameters
35+
vector [Nparameters] population_locations;
36+
vector<lower=0>[Nparameters] population_scales;
37+
38+
//individuals level
39+
vector[Nsubjects] alpha_confirmatory_random_effect;
40+
vector[Nsubjects] alpha_disconfirmatory_random_effect;
41+
vector[Nsubjects] beta_random_effect;
42+
}
43+
44+
45+
transformed parameters {
46+
47+
vector<lower=0, upper=1>[Nsubjects] alpha_confirmatory;
48+
vector<lower=0, upper=1>[Nsubjects] alpha_disconfirmatory;
49+
vector [Nsubjects] beta;
50+
matrix [Ntrials,Nsubjects] p_ch_action;
51+
matrix [Ntrials,Nsubjects] Qdiff_external;
52+
matrix [Ntrials,Nsubjects] Qval1_external;
53+
matrix [Ntrials,Nsubjects] Qval2_external;
54+
matrix [Ntrials,Nsubjects] Qval3_external;
55+
matrix [Ntrials,Nsubjects] Qval4_external;
56+
matrix [Ntrials,Nsubjects] PE_external;
57+
58+
59+
60+
//RL
61+
for (subject in 1:Nsubjects) {
62+
//internal variabels
63+
real Qdiff;
64+
real PE;
65+
real Qval[Narms];
66+
67+
//set indvidual parameters
68+
alpha_confirmatory[subject] = inv_logit(population_locations[1] + population_scales[1] * alpha_confirmatory_random_effect[subject]);
69+
alpha_disconfirmatory[subject] = inv_logit(population_locations[2] + population_scales[2] * alpha_disconfirmatory_random_effect[subject]);
70+
beta[subject] = (population_locations[3] + population_scales[3] * beta_random_effect [subject]);
71+
72+
//likelihood estimation
73+
for (trial in 1:Ntrials_per_subject[subject]){
74+
75+
//reset Qvalues (first trial only)
76+
if (first_trial_in_block[subject,trial] == 1) {
77+
Qval = rep_array(0.5, Narms);
78+
}
79+
80+
//calculate probability for each action
81+
Qdiff = Qval[offer2[subject,trial]]- Qval[offer1[subject,trial]];
82+
83+
p_ch_action[trial,subject] = inv_logit(beta[subject]*Qdiff);
84+
85+
//update Qvalues
86+
87+
if(reward[subject,trial]==1){
88+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_confirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
89+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_confirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
90+
}
91+
else{
92+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_disconfirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
93+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_disconfirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
94+
}
95+
#appened to external variabels
96+
Qdiff_external[trial,subject] = Qdiff;
97+
Qval1_external[trial,subject] = Qval[1];
98+
Qval2_external[trial,subject] = Qval[2];
99+
Qval3_external[trial,subject] = Qval[3];
100+
Qval4_external[trial,subject] = Qval[4];
101+
PE_external[trial,subject] = PE;
102+
103+
104+
}
105+
106+
}
107+
108+
}
109+
110+
111+
model {
112+
113+
// population level
114+
population_locations ~ normal(0,2);
115+
population_scales ~ cauchy(0,2);
116+
117+
// indvidual level
118+
alpha_confirmatory_random_effect ~ std_normal();
119+
alpha_disconfirmatory_random_effect ~ std_normal();
120+
beta_random_effect ~ std_normal();
121+
122+
123+
for (subject in 1:Nsubjects){
124+
for (trial in 1:Ntrials_per_subject[subject]){
125+
target+= bernoulli_logit_lpmf(selected_offer[subject,trial] | beta[subject] * Qdiff_external[trial,subject]);
126+
}
127+
}
128+
}

‎stan_modeling/models/confirmatory_two_alphas/confirmatory_two_alphas_loo.stan‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
data {
2+
3+
//General fixed parameters for the experiment/models
4+
int<lower = 1> Nsubjects;
5+
int<lower = 1> Nblocks;
6+
int<lower = 1> Ntrials;
7+
int<lower = 1> Ntrials_per_subject[Nsubjects];
8+
int<lower = 4> Narms;
9+
int<lower = 2> Nraffle;
10+
int<lower = 0> fold[Nsubjects,Ntrials];
11+
real testfold;
12+
13+
//Behavioral data:
14+
int<lower = 0> choice[Nsubjects,Ntrials];
15+
int<lower = 0> unchosen[Nsubjects,Ntrials];
16+
int<lower = 0> reward[Nsubjects,Ntrials];
17+
int<lower = 0> offer1[Nsubjects,Ntrials];
18+
int<lower = 0> offer2[Nsubjects,Ntrials];
19+
int<lower = 0> selected_offer[Nsubjects,Ntrials];
20+
int<lower = 0> first_trial_in_block[Nsubjects,Ntrials];
21+
22+
23+
}
24+
25+
26+
transformed data{
27+
int<lower = 1> Nparameters=3;
28+
vector[Narms] Qvalue_initial;
29+
Qvalue_initial = rep_vector(0.5, Narms);
30+
}
31+
32+
33+
34+
35+
parameters {
36+
//population level parameters
37+
vector [Nparameters] population_locations;
38+
vector<lower=0>[Nparameters] population_scales;
39+
40+
//individuals level
41+
vector[Nsubjects] alpha_confirmatory_random_effect;
42+
vector[Nsubjects] alpha_disconfirmatory_random_effect;
43+
vector[Nsubjects] beta_random_effect;
44+
}
45+
46+
47+
transformed parameters {
48+
49+
vector<lower=0, upper=1>[Nsubjects] alpha_confirmatory;
50+
vector<lower=0, upper=1>[Nsubjects] alpha_disconfirmatory;
51+
vector [Nsubjects] beta;
52+
matrix [Ntrials,Nsubjects] p_ch_action;
53+
matrix [Ntrials,Nsubjects] Qdiff_external;
54+
matrix [Ntrials,Nsubjects] Qval1_external;
55+
matrix [Ntrials,Nsubjects] Qval2_external;
56+
matrix [Ntrials,Nsubjects] Qval3_external;
57+
matrix [Ntrials,Nsubjects] Qval4_external;
58+
matrix [Ntrials,Nsubjects] PE_external;
59+
60+
61+
62+
//RL
63+
for (subject in 1:Nsubjects) {
64+
//internal variabels
65+
real Qdiff;
66+
real PE;
67+
real Qval[Narms];
68+
69+
//set indvidual parameters
70+
alpha_confirmatory[subject] = inv_logit(population_locations[1] + population_scales[1] * alpha_confirmatory_random_effect[subject]);
71+
alpha_disconfirmatory[subject] = inv_logit(population_locations[2] + population_scales[2] * alpha_disconfirmatory_random_effect[subject]);
72+
beta[subject] = (population_locations[3] + population_scales[3] * beta_random_effect [subject]);
73+
74+
//likelihood estimation
75+
for (trial in 1:Ntrials_per_subject[subject]){
76+
if(fold[subject,trial]!=testfold){
77+
78+
//reset Qvalues (first trial only)
79+
if (first_trial_in_block[subject,trial] == 1) {
80+
Qval = rep_array(0.5, Narms);
81+
}
82+
83+
//calculate probability for each action
84+
Qdiff = Qval[offer2[subject,trial]]- Qval[offer1[subject,trial]];
85+
86+
p_ch_action[trial,subject] = inv_logit(beta[subject]*Qdiff);
87+
88+
//update Qvalues
89+
if(reward[subject,trial]==1){
90+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_confirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
91+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_confirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
92+
}
93+
else{
94+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_disconfirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
95+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_disconfirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
96+
}
97+
//appened to external variabels
98+
Qdiff_external[trial,subject] = Qdiff;
99+
Qval1_external[trial,subject] = Qval[1];
100+
Qval2_external[trial,subject] = Qval[2];
101+
Qval3_external[trial,subject] = Qval[3];
102+
Qval4_external[trial,subject] = Qval[4];
103+
PE_external[trial,subject] = PE;
104+
}
105+
}
106+
107+
}
108+
109+
}
110+
111+
112+
model {
113+
114+
// population level
115+
population_locations ~ normal(0,2);
116+
population_scales ~ cauchy(0,2);
117+
118+
// indvidual level
119+
alpha_confirmatory_random_effect ~ std_normal();
120+
alpha_disconfirmatory_random_effect ~ std_normal();
121+
beta_random_effect ~ std_normal();
122+
123+
124+
for (subject in 1:Nsubjects){
125+
for (trial in 1:Ntrials_per_subject[subject]){
126+
if(fold[subject,trial]!=testfold){ //fit parameters only for training set
127+
target+= bernoulli_logit_lpmf(selected_offer[subject,trial] | beta[subject] * Qdiff_external[trial,subject]);
128+
}
129+
}
130+
}
131+
}
132+
133+
generated quantities {
134+
135+
matrix[Ntrials,Nsubjects] log_lik;
136+
vector[Narms] Qval;
137+
vector[Nraffle] Qoffer;
138+
real PE;
139+
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
140+
//Likelihood function per subject per trial
141+
log_lik=rep_matrix(0,Ntrials,Nsubjects);
142+
for (subject in 1:Nsubjects) {
143+
144+
for (trial in 1:Ntrials_per_subject[subject]){
145+
146+
if(fold[subject,trial] == testfold) {
147+
148+
//reset Qvalues (first trial only)
149+
if (first_trial_in_block[subject,trial] == 1) {
150+
Qval =rep_vector(0.5, Narms);
151+
}
152+
//offer values
153+
Qoffer[1]=Qval[offer1[subject,trial]];
154+
Qoffer[2]=Qval[offer2[subject,trial]];
155+
log_lik[trial,subject] = bernoulli_logit_lpmf(selected_offer[subject, trial] | beta[subject] * Qoffer);
156+
157+
158+
159+
//update Qvalues
160+
if(reward[subject,trial]==1){
161+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_confirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
162+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_confirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
163+
}
164+
else{
165+
Qval[choice[subject,trial]] = Qval[choice[subject,trial]]+alpha_disconfirmatory[subject]*(reward[subject,trial] - Qval[choice[subject,trial]]);
166+
Qval[unchosen[subject,trial]] = Qval[unchosen[subject,trial]]+alpha_disconfirmatory[subject]*((1-reward[subject,trial]) - Qval[unchosen[subject,trial]]);
167+
}
168+
169+
170+
}
171+
}
172+
}
173+
}

0 commit comments

Comments
 (0)