How to make a variable in the function equal to the sum of that variable in previous trials in jags

Question

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) 

Answer 1

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")

---

Edit: Answer to edited question.

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.