62 views

Skip to first unread message

May 25, 2017, 1:51:33 AM5/25/17

to stan-...@googlegroups.com

Hi Stan users,

I have a question regarding a multinomial item response theory model. I am trying to transform the following BUGS code to a STAN code.

BUGS code here.

for(i in 1:lengthy){

yy[i] ~ dcat(prob[nidx[i],midx[i],1:K[midx[i]]])

}

for (i in 1:lengthy){

for (k in 1:K[midx[i]]){

p[nidx[i],midx[i],k] <- exp(beta[midx[i],k]*x[nidx[i]] + beta2[midx[i],k]*x2[nidx[i]] - alpha[midx[i],k])

}

prob[nidx[i],midx[i],1:K[midx[i]] ] <- p[nidx[i],midx[i],1:K[midx[i]] ]/sum(p[nidx[i],midx[i],1:K[midx[i]]])

}

Through the second for loop, I will obtain an array (number of individuals by number of questions by number of categories in questions) of probabilities. These probabilities will go through dcat to the response yy.

I tried to code this in Stan because it is faster to implement this IRT model. But I got the following error message:

attempt to assign variable in wrong block. left-hand-side variable origin=parameter

Here is the Stan code that I wrote.

library(rstan)

stan.code <- "data {

int<lower=1> J; // number of Individuals

int<lower=1> M; // number of Questions

int<lower=1> N; // number of observations

int<lower=1> K; // maximum number of categories

int<lower=1, upper=J> j[N]; // Individual for observation n

int<lower=1, upper=M> m[N]; // Question for observation n

int<lower=1, upper=K> y[N]; // Answer of observation n

int<lower=1, upper=K> k[M]; // Maximum categories of question m

int<lower=1, upper=K> q[N]; // Maximum categories of observation n

int<lower=1> D; //no. of dimensions

}

parameters {

matrix[M,K] alpha;

matrix[M,K] beta;

matrix[M,K] beta2;

real theta[J];

real theta2[J];

real<lower=0> p[J, M, K];

real<lower=0> prob[J, M, K];

}

model {

to_vector(alpha) ~ normal(0,1);

for (i in 1:M){

alpha[i,1] ~ normal(0,0.000001);

for (g in (k[i]+1):K)

alpha[i,g] ~ normal(999,0.000001);}

to_vector(beta) ~ normal(0,1);

for (i in 1:M){

beta[i,1] ~ normal(0,0.000001);

for (g in (k[i]+1):K)

beta[i,g] ~ normal(999,0.000001);}

to_vector(beta2) ~ normal(0,1);

for (i in 1:M){

beta2[i,1] ~ normal(0,0.000001);

for (g in (k[i]+1):K)

beta2[i,g] ~ normal(999,0.000001);}

theta ~ normal(0,0.1);

theta2 ~ normal(0,0.1); // priors

for (n in 1:N){

for (h in 1:q[n]){

p[j[n], m[n], h] = exp(theta[j[n]] * beta[m[n],h] + theta2[j[n]] * beta2[m[n],h] - alpha[m[n],h])}

prob[j[n], m[n], 1:q[n]] = p[j[n], m[n], 1:q[n]]/sum(p[j[n], m[n], 1:q[n]])}

for (n in 1:N){

y[n] ~ categorical(prob[j[n], m[n], 1:q[n]])

}

}

"

J <- dim(jdata$y)[1] # number of Individuals

M <- dim(jdata$y)[2] # number of Questions

N <- length(jdata$y) # number of observations

K <- max(jdata$K) # maximum number of categories

j <- rep(1:J,M) # individual for observation n

m <- rep(1:M, each=J) # question for observation n

y <- as.vector(jdata$y) # answer of observation n

k <- jdata$K # maximum categories of question m

q <- rep(jdata$K, each=J) # maximum categories of observation n

D <- 2 # no of dimensions

# deleting missing values

miss <- which(is.na(y))

N <- N - length(miss)

j <- j[-miss]

m <- m[-miss]

y <- y[-miss]

q <- q[-miss]

## data and initial values

stan.data <- list(J=J, M=M, N=N, K=K, j=j, m=m, y=y, k=k, q=q, D=D)

starting=list(theta=startvalues+rnorm(length(startvalues)),

theta2=rnorm(length(startvalues)))

stan.fit <- stan(model_code=stan.code, data=stan.data, iter=500, warmup=200,

chains=1, thin=2, init=starting)

Why is my code wrong? How should I fix the code to implement the same thing as BUGS?

Thank you in advance!

John

May 25, 2017, 3:17:06 PM5/25/17

to Stan users mailing list

That error just means you are trying to overwrite a parameter (as defined in parameters block), which you can't do.

From just a quick glance, that seems to be p[j[n], m[n], h], which you're defining as a parameter, but then assigning to in the model block.

I'm guessing this is due to some confusion about what the parameters block is for. Parameters{} isn't for listing "things I want estimated", it's for things you want "sampled" (really, for defining ficticious particles that are simulated in multidimensional space). If p[] is a COMPUTED parameter, then you should compute it in the transformed parameters block. If it's just a utility variable that you don't need to save or do inference on, but is merely used in computing the posterior, then declare and define it directly in the model block.

From just a quick glance, that seems to be p[j[n], m[n], h], which you're defining as a parameter, but then assigning to in the model block.

I'm guessing this is due to some confusion about what the parameters block is for. Parameters{} isn't for listing "things I want estimated", it's for things you want "sampled" (really, for defining ficticious particles that are simulated in multidimensional space). If p[] is a COMPUTED parameter, then you should compute it in the transformed parameters block. If it's just a utility variable that you don't need to save or do inference on, but is merely used in computing the posterior, then declare and define it directly in the model block.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu