Question
I am trying to migrate some code from JAGS to Stan. Say I have the following dataset:
N <- 10
nchoices <- 3
ncontrols <- 3
toydata <- list("y" = rbinom(N, nchoices - 1, .5),
"controls" = matrix(runif(N*ncontrols), N, ncontrols),
"N" = N,
"nchoices" = nchoices,
"ncontrols" = ncontrols)
and that I want to run a multinomial logit with the following code (taken from section 9.5 of the documentation):
data {
int N;
int nchoices;
int y[N];
int ncontrols;
vector[ncontrols] controls[N];
}
parameters {
matrix[nchoices, ncontrols] beta;
}
model {
for (k in 1:nchoices)
for (d in 1:ncontrols)
beta[k,d] ~ normal(0,100);
for (n in 1:N)
y[n] ~ categorical(softmax(beta * controls[n]));
}
I now want to fix the first row of beta
to zero. In JAGS I would simply declare in the model block that
for (i in 1:ncontrols) {
beta[1,i] <- 0
}
but I am not sure about how to do this in Stan. I have tried many combinations along the lines of section 6.2 of the documentation (Partially Known Parameters) like, for instance,
parameters {
matrix[nchoices, ncontrols] betaNonObs;
}
transformed parameters {
matrix[nchoices, ncontrols] beta;
for (i in 1:ncontrols) beta[1][i] <- 0
for (k in 2:nchoices) beta[k] <- betaNonObs[k - 1]
}
but none of them work. Any suggestions?
La solution
It would be helpful to mention the error message. In this case, if beta is declared to be a matrix, then the syntax you want is the R-like syntax
beta[1,i] <- 0.0; // you also omitted the semicolon
To answer your broader question, I believe you were on the right track with your last approach. I would create a matrix of parameters in the parameters block called free_beta and copy those elements to another matrix declared in the model block called beta that has one extra row at the top for the fixed zeros. Like
data {
int N;
int nchoices;
int y[N];
int ncontrols;
vector[ncontrols] controls[N];
}
parameters {
matrix[nchoices-1, ncontrols] free_beta;
}
model {
// copy free beta into beta
matrix[nchoices,ncontrols] beta;
for (d in 1:ncontrols)
beta[1,d] <- 0.0;
for (k in 2:nchoices)
for (d in 1:ncontrols)
beta[k,d] <- free_beta[k-1,d];
// priors on free_beta, which execute faster this way
for (k in 1:(nchoices-1))
row(free_beta,k) ~ normal(0.0, 100.0);
// likelihood
for (n in 1:N)
y[n] ~ categorical(softmax(beta * controls[n]));
}