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