How to generate a matrix from multionmial distribution

Hi, I am following Jim Savage rank ordered logit approach generally. However, instead of his beta_individual definition, I want to generate beta_individual (a matrix) from a multinomial distribution. Then I get the error because stan can only generate a vector or row_vector rather than a matrix under the distribution. Any suggestions? (The sigma equation is not correct while I just list an example to show sigma is a ixp matrix.)

Stan code for the ranked choice model

// saved as ranked_rcl.stan
functions {
  real rank_logit_lpmf(int[] rank_order, vector delta) {
    // We reorder the raw utilities so that the first rank is first, second rank second... 
    vector[rows(delta)] tmp = delta[rank_order];
    real out;
    // ... and sequentially take the log of the first element of the softmax applied to the remaining
    // unranked elements.
    for(i in 1:(rows(tmp) - 1)) {
      if(i == 1) {
        out = tmp[1] - log_sum_exp(tmp);
      } else {
        out += tmp[i] - log_sum_exp(tmp[i:]);
    // And return the log likelihood of observing that ranking
data {
  int N; // number of rows
  int T; // number of inidvidual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  int rank_order[N]; // The vector describing the index (within each task) of the first, second, third, ... choices. 
  // In R, this is order(-utility) within each task
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  matrix[I, P] beta_individual ;
transformed parameters {
  // here we use the reparameterization discussed on slide 30
  matrix[I, P] sigma=X2*Gamma';
model {
  // priors on the parameters
  tau ~ normal(0, .5);
  beta ~ normal(0, 1);
  to_vector(z) ~ normal(0, 1);
  to_vector(Gamma) ~ normal(0, 1);
  **beta_individual ~ multi_normal(beta,sigma);**
  // log probabilities of each choice in the dataset
  for(t in 1:T) {
    vector[K] utilities; // tmp vector holding the utilities for the task/individual combination
    // add utility from product attributes with individual part-worths/marginal utilities
    utilities = X[start[t]:end[t]]*beta_individual[task_individual[t]]';
    rank_order[start[t]:end[t]] ~ rank_logit(utilities);

So a draw from the multinomial distribution is an array of integers.

Is it that you want to reshape this array of N integers to look like a 2d array?

Or do you want to arrange N arrays of M integers beside each other to form an array with N * M elements?

Hi, thanks for responding. I want a matrix I X P from the multinomial distribution like N= I X P then I row P column, while I am confused with how to make it.

A multinomial distribution is defined for a 1d array of integers though, so how does this 1d array map to a matrix or how is this matrix construct from 1d arrays?