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?

Était-ce utile?

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]));
}
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top