Let’s say we have a matrix:

and a matrix:

where each row of a matrix is supposed to be a mean vector of some m-dimensional Multivariate Normal Distribution with Covariance matrix .

And we need to define Multivariate Normal Distributions with different mean vectors but with the same covariance matrix inside of a Model. How to achieve this without explicit "for loop" which is a quite performance overhead?

A rough code:

```
with pm.Model():
Sigma = pm.Deterministic(...)
M = pm.Deterministic(...)
# unoptimal solution
for i in xrange(n):
_Y = pm.MvNormal('_Y_' + str(i), mu=mu_[:, i], cov=sd_, shape=(m, ))
# Is it possible to vectorize last lines using pm.MvNormal or pm.MatrixNormal ??
```

Thanks,

Best Regards

Alex