Multiple draws from same multivariate normal

I would like to create a model like the below, where each row in the response matrix is drawn from the same multivariate normal distribution.

data {
  int<lower=1> N;
  matrix[N, 3] ys;
}

parameters {
  cov_matrix[3] covmat;
  vector[3] means;
}

model {
  ys ~ multi_normal(means, covmat);
}

but I can only get it working by using arrays of vectors instead of matrixes and by broadcasting the mean vector manually as below


data {
  int<lower=1> N;
  array[N] vector[3] ys;
}

parameters {
  cov_matrix[3] covmat;
  vector[3] means;
}

transformed parameters {
  array[N] vector[3] b_means;
  
  b_means[,1] = to_array_1d(rep_vector(means[1], N));
  b_means[,2] = to_array_1d(rep_vector(means[2], N));
  b_means[,3] = to_array_1d(rep_vector(means[3], N));
}

model {
  ys ~ multi_normal(b_means, covmat)
}

Is this really the correct way to specify this model in stan 2.29? The signatures for multi_normal seem extremely limited

See discussion here.

Oops, I mean here.

Thanks for the link! The discussion really helped me understand the function