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