Here’s a minimal working example of an error I’m getting with sample_prior_predictive():

```
import numpy as np
import pymc3 as pm
K = 3
D = 15
mu_0 = np.zeros((D,K))
lambd = 1.0
with pm.Model() as model:
sd_dist = pm.HalfCauchy.dist(beta=2.5)
packedL = pm.LKJCholeskyCov(f"packedL",eta=2, n=D, sd_dist=sd_dist)
L = pm.expand_packed_triangular(D, packedL, lower=True)
Sigma = pm.Deterministic(f"Sigma", L.dot(L.T)) # D x D covariance
mu = pm.MatrixNormal(
f"mu",
mu=mu_0,
rowcov=(1 / lambd) * Sigma,
colcov = np.eye(K), shape=(D,K))
prior = pm.sample_prior_predictive(2)
```

The beginning of the traceback:

```
ValueError Traceback (most recent call last)
<ipython-input-8-37de3b873b9c> in <module>
10 rowcov=(1 / lambd) * Sigma,
11 colcov = np.eye(K), shape=(D,K))
---> 12 prior = pm.sample_prior_predictive(2)
```

And the end:

```
/opt/conda/lib/python3.6/site-packages/pymc3/distributions/multivariate.py in random(self, point, size)
1617 samples.append(mu[j] +
1618 np.matmul(rowchol[j], np.matmul(standard_normal, colchol[j].T)))
-> 1619 samples = np.array(samples).reshape(size + tuple(self.shape))
1620 return samples
1621
ValueError: cannot reshape array of size 180 into shape (2,15,3)
```

I suppose I could instead just put each (15,) mu vector into a length-3 array. But I’d rather not, since my actual model is considerably more complex and I’d rather not have nested arrays.