# A Sequential Monte Carlo model was used to fit the subjects transitive inference (TI) process

The task was: there were a total of 7 picture stimuli with different potential ratings (L=7), and 2 stimuli were presented in each pair for subjects to choose from (choice[nSubjects,nTrials]), with a feedback of 1 if the subject chose the higher-ranked one, and 0 if they did the opposite (reward[nSubjects,nTrials]). Nsmc(10000) represents the number of states I want to use the SMC process to build a hierarchical model fitting the process by which subjects learn the potential rank of a stimulus. Here is my stan code

``````data {
int<lower=1> nTrials;
int<lower=1> nSubjects;
int<lower=1> Nsmc;
int<lower=1> L;
int<lower=1,upper=7> nchoice[nSubjects,nTrials];
int<lower=1,upper=7> choice[nSubjects,nTrials];
int<lower=0,upper=1> reward[nSubjects,nTrials];
}

parameters {
// sigma
real  sigma_mu;
real <lower=0> sigma_sd;

// tau
real <lower=0,upper=3> tau_mu;
real <lower=0> tau_sd;

vector <lower=0,upper=1>[nSubjects] sigma_raw;
vector <lower=0,upper=3>[nSubjects] tau_raw;

}

transformed parameters {
vector <lower=0,upper=1>[nSubjects] sigma;
vector <lower=0,upper=3>[nSubjects] tau;

sigma = Phi_approx(sigma_mu+sigma_sd*sigma_raw);
tau = Phi_approx(tau_mu+tau_sd*tau_raw)*3;

}
model {

sigma_sd ~ cauchy(0,1);
tau_sd ~ cauchy(0,3);
sigma_raw ~ normal(0,1);
tau_raw ~ normal(0,1);

for (s in 1:nSubjects){
vector[2] stim;
matrix[Nsmc,L] x;
for(i in 1:Nsmc){
for(j in 1:L){
x[i,j] ~ normal(0,1);
}
}
vector[Nsmc] w_old;
w_old = rep_vector(1.0/Nsmc, Nsmc);
vector[Nsmc] w;
w = rep_vector(0, Nsmc);

matrix[Nsmc,2] p;
real Neff;
vector[Nsmc] ix;

for (t in 1:nTrials){
matrix[Nsmc,L] delta;
for(i in 1:Nsmc){
for(j in 1:L){
delta[i,j] ~ normal(0,sigma);
}
}

if(reward[s,t] > 0){

for(n in 1:Nsmc){
stim[1]=x[n,choice[s,t]];
stim[2]=x[n,nchoice[s,t]];
p[n,]=softmax(tau[s]*stim);
1 ~ categorical(p);
}
x = x+delta;
w = p[,1]*w_old;
w = w/sum(w);
Neff = 1.0/sum(w_old*w_old);
if (Neff < 0.25*Nsmc) {
ix[Nsmc] ~ categorical_logit(w);
x = x[ix,];
w_old = rep_vector(1.0/Nsmc,Nsmc);
} else {
w_old = w;
}

}else{
for(n in 1:Nsmc){
stim[1]=x[n,nchoice[s,t]];
stim[2]=x[n,choice[s,t]];
p[n,]=softmax(tau[s]*stim);
1 ~ categorical(p);
}
x = x+delta;
w = p[,1]*w_old;
w = w/sum(w);
Neff = 1.0/sum(w_old*w_old);
if (Neff < 0.25*Nsmc) {
ix[Nsmc] ~ categorical_logit(w);
x = x[ix,];
w_old = rep_vector(1.0/Nsmc,Nsmc);
} else {
w_old = w;
}

}

}
}
}

``````

The reported error message is:Ill-typed arguments supplied to assignment operator =: lhs has type row_vector and rhs has type vector

This code I adapted from the r-code provided in the paper. Here is the r-code SMC function from the paper. I found that stan doesn’t seem to be able to sample random numbers. In the code, N represents the number of states.the state variable is initialized a normal distribution with fixed
initial variance , and zero mean(sig0). The state process is described as a Gaussian random walk
with evolution variance sig (free parameter).resp is feedback.stim represents the paired two stimuli.

``````smc <- function(stims,L,N,Sig,Sig0,Luce,RespBool=FALSE,Resp=NULL,Counterfactual=TRUE) {
g_smc <- function(x,a,b,Luce) {
xa <- x[,a]
xb <- x[,b]
y <- inv_logit(Luce*(xa-xb))
return(y)
}

inv_logit <- function( x ) {
p <- 1/(1+exp(-x))
p <- ifelse(x==Inf,1,p)
return(p)
}

npairs <- L*(L-1)/2
pairs <- matrix(0,npairs,2)
dex <- 0
for (i in 1:(L-1)) {
for (j in 1:(L-i)) {
dex <- dex+1
pairs[dex,1] <- j
pairs[dex,2] <- j+i
}
}

ntrials <- length(stims[,1])
xbar <- matrix(0,L,ntrials)
pairp <- matrix(0,npairs,ntrials)
correct_prob <- rep(0,ntrials)
if (RespBool==FALSE) {
Resp <- rep(0,ntrials)
}

w_old <- rep(1/N,N)
w <- rep(0,N)

x <- rnorm(N*L,0,Sig0)
dim(x) <- c(N,L)

for (i in 1:ntrials) {
xbar[,i] <- apply(x,2,mean)
delta <- rnorm(N*L,0,Sig)
dim(delta) <- c(N,L)
x <- x+delta

a <- stims[i,1]
b <- stims[i,2]
p <- max(c(0,min(c(1,sum((x[,a]>x[,b])*w)))))
if (Counterfactual==TRUE) {
for (j in 1:npairs) {
pairp[j,i] <- max(c(0,min(c(1,sum((x[,pairs[j,1]]>x[,pairs[j,2]])*w)))))
}
}

correct_prob[i] <- p
if (RespBool==FALSE) {
Resp[i] <- sample(c(1,0),1,prob=c(p,1-p))
}

w <- g_smc(x,a,b,Luce)*w_old
w <- w/sum(w)

Neff <- 1/sum(w_old^2)
if (Neff < 0.25*N) {
ix <- sample(1:N,N,replace = TRUE,prob=w)
x <- x[ix,]
w_old <- rep(1/N,N)
} else {
w_old <- w
}
}
output <- list("xbar" = xbar, "choices" = correct_prob, "resps" = Resp, "pairp" = pairp)
return(output)
}
``````

Sorry, @haohao, but Stan doesn’t implement SMC.

I mean using stan to model SMC cognitive processes, not using SMC as a proxy for HMC to estimate model parameters.