In my model, I have 10 options (from 1-10) for each subject and each trial to choose (expectancy
). I calculated the value for each option based on the rule in the graph below, so the value for each option updated based on the difference between shock
and v
in every trial (multiply alpha
). Then, I used softmax rule to transform v
for each option to a certain probability with the same function in this threat: JAGS errors: "Resolving undeclared variables" and "Invalid vector argument to exp".
I guess the problem here is I can't make jags update the value for the same choice.
data: expectancy
= number from 1-10 in each trial. shock
=number either 1 or 0 in each trial. (I provided example data below)
The second plot is how this be done in stan with 2 choices/1 subject situation.
RW_model <- function(){
# data
for(i in 1:nsubjects) # for each person
{
# initial value for v
v [i,1,expectancy[i,1]] <- 0
for (j in 2:ntrials) # for each trial
{
# expectancy chosen
expectancy[i,j] ~ dcat(mu[i,j,1:10])
predk[i,j] ~ dcat(mu[i,j,1:10])
# softmax rule to calculate values of each expectancy for each subject
# tau is the value sensitivity parameter
mu[i,j,1:10] <- exp_v[i,j,1:10] / sum(exp_v[i,j,1:10])
exp_v[i,j,expectancy[i,j-1]] <- exp(v[i,j,expectancy[i,j-1]]/tau[i])
# prediction error: difference between feedback and learned values of the chosen expectancy
pe [i,j-1] <- shock [i,j-1] - v [i,j-1,expectancy[i,j-1]]
# value updating process for expectancy
v [i,j,expectancy[i,j-1]] <- v [i,j-1,expectancy[i,j-1]] + alpha [i] * pe [i,j-1]
}
}
# priors
for (i in 1:nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
# example data/ initial value/ parameters
nsubjects <- 42
ntrials <- 14
shock <- matrix(c(0,0,1,1,0,0,1,1,0,0,1,0,1,0),nrow=42,ncol = 14,byrow = T)
expectancy <- matrix(c(1,2,3,4,5,6,7,7,8,8,7,10,10,00),nrow=42,ncol = 14,byrow = T)
data <- list('shock','nsubjects','ntrials','expectancy')
myinits <- list(list(tau = runif (42,0,3),
alpha = runif (42,0,1)))
parameters <- c("tau",'alpha','v','predk')
# jags sampling
samples <- jags(data, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
Since there are no stochastic nodes in the setup, I'm not sure what the MCMC simulation will get you, but here is some code that works. From what I can see, the main problem is that when the table says c(Vmax-Vn)
it is not using c()
as a concatenation operator. c
is a constant defined as 0.3.
dat <- list(
v = c(0, rep(NA, 8)),
vn = rep(NA, 8),
vmax=1,
c=0.3,
Ntrials = 9
)
jmod <- "model{
for(i in 2:Ntrials){
v[i] <- c*(vmax - v[(i-1)]) + v[(i-1)]
}
}"
library(runjags)
out <- run.jags(jmod, data=dat, monitor="v")
You'll have to see if this is doing what you want, but it produces a result. As far as I can tell, this is a direct extension of the Stan code to multiple subjects and multiple choices per trial ported to JAGS. First, the model:
RW_model <- function(){
for(i in 1:Nsubjects){
for(j in 1:Ntrials){
expectancy[i,j] ~ dcat(p[i, 1:Nposs, j])
mu[i,1:Nposs, j] <- tau[i] * v[i,1:Nposs, j]
for(k in 1:Nposs){
q[i,k,j] <- exp(mu[i,k,j])
p[i,k,j] <- q[i,k,j]/sum(q[i,,j])
}
pe[i,j] <- shock[i,j] - v[i,expectancy[i,j], j]
for(k in 1:Nposs){
v[i,k,j+1] <- v[i,k,j] + ifelse(expectancy[i,j] == k,
alpha[i] * pe[i,j],
0)
}
}
}
for (i in 1:Nsubjects){
tau [i] ~ dunif (0,3)
alpha [i] ~ dunif (0,1)
}
}
Next, we could make up some data. The data here are for 20 subjects, 5 possible choices in expectancy
and 14 trials.
set.seed(40120)
# v has to be a Nsubjects x Nchoices x Ntrials+1 matrix
v <- array(NA, dim=c(20, 5, 15))
# The first trial of v is initialized to 0 fo all subjects and choices
v[,,1] <- matrix(0, nrow=20, ncol=5)
# expectancy and shock are both Nsubjects x Ntrials matrices
expectancy <- matrix(sample(1:5, 20*14, replace=TRUE), ncol=14)
shock <- matrix(sample(c(0,1), 20*14, replace=TRUE), ncol=14)
dlist <- list(
Nsubjects = 20,
Ntrials = 14,
Nposs = 5,
expectancy = expectancy,
shock = shock,
v=v
)
Finally, we can specify initial values and run the model.
myinits <- list(list(tau = runif (nrow(expectancy),0,3),
alpha = runif (nrow(expectancy)0,1)))
parameters <- c("tau",'alpha')
library(R2jags)
# jags sampling
samples <- jags(dlist, inits=myinits, parameters,
model.file = RW_model,
n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
You can see from the output that it generated samples and estimates different tau
and alpha
parameters by subject.
samples$BUGSoutput
Inference for Bugs model at "/var/folders/55/q9y1hbcx13b5g50kks_p0mb00000gn/T//Rtmp32zeG4/model1726bdee0c99.txt", fit using jags,
1 chains, each with 1000 iterations (first 500 discarded)
n.sims = 500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5%
alpha[1] 0.6 0.3 0.1 0.4 0.7 0.9 1.0
alpha[2] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[3] 0.3 0.2 0.0 0.1 0.2 0.4 0.9
alpha[4] 0.4 0.3 0.0 0.1 0.3 0.6 0.9
alpha[5] 0.4 0.3 0.0 0.1 0.3 0.6 1.0
alpha[6] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[7] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[8] 0.3 0.2 0.0 0.1 0.2 0.5 0.9
alpha[9] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[10] 0.4 0.2 0.0 0.2 0.4 0.5 0.9
alpha[11] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[12] 0.3 0.3 0.0 0.1 0.2 0.5 0.9
alpha[13] 0.3 0.3 0.0 0.1 0.3 0.5 0.9
alpha[14] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[15] 0.5 0.3 0.0 0.2 0.5 0.7 1.0
alpha[16] 0.4 0.3 0.0 0.1 0.3 0.7 0.9
alpha[17] 0.4 0.3 0.0 0.2 0.4 0.6 0.9
alpha[18] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[19] 0.4 0.3 0.0 0.2 0.4 0.7 1.0
alpha[20] 0.3 0.3 0.0 0.0 0.2 0.5 0.9
deviance 909.3 5.1 901.1 905.8 908.8 912.5 919.9
tau[1] 1.1 0.6 0.2 0.7 1.1 1.5 2.3
tau[2] 0.8 0.7 0.0 0.3 0.6 1.2 2.7
tau[3] 1.0 0.8 0.0 0.3 0.8 1.6 2.8
tau[4] 1.1 0.7 0.0 0.4 0.9 1.6 2.7
tau[5] 0.9 0.7 0.0 0.3 0.7 1.3 2.7
tau[6] 1.0 0.8 0.0 0.3 0.9 1.6 2.8
tau[7] 1.0 0.7 0.1 0.5 0.9 1.5 2.7
tau[8] 1.1 0.8 0.0 0.4 0.9 1.7 2.9
tau[9] 1.0 0.8 0.0 0.3 0.8 1.6 2.7
tau[10] 1.6 0.8 0.1 0.9 1.7 2.3 2.9
tau[11] 1.1 0.9 0.0 0.3 0.9 1.8 2.8
tau[12] 1.0 0.8 0.0 0.4 0.8 1.6 2.9
tau[13] 0.9 0.8 0.0 0.3 0.6 1.4 2.7
tau[14] 1.1 0.8 0.0 0.5 0.9 1.6 2.7
tau[15] 1.2 0.7 0.1 0.7 1.1 1.6 2.7
tau[16] 1.0 0.8 0.1 0.4 0.9 1.5 2.8
tau[17] 1.1 0.8 0.1 0.4 0.8 1.6 2.9
tau[18] 1.2 0.8 0.1 0.5 1.1 1.7 2.7
tau[19] 1.2 0.8 0.1 0.5 1.0 1.8 2.8
tau[20] 1.0 0.9 0.0 0.3 0.8 1.5 2.9
DIC info (using the rule, pD = var(deviance)/2)
pD = 12.9 and DIC = 922.1
DIC is an estimate of expected predictive error (lower deviance is better).
Since I don't know what to expect from this model, you'll have to verify that it is actually doing what you expect, but it seems like this is a good place to start.