Reshape error inside sample_prior_predictive

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.

Hm, that shape error might be a bug, right?
In the mean time you could just transform a iid standard normal array using this: https://en.wikipedia.org/wiki/Matrix_normal_distribution#Drawing_values_from_the_distribution
Since you already have the cholesky of the rowcol, that should be faster. It is a differerent parametrization though, so the sampler might like it more or less. (often more :slight_smile: )

The shape bug is now here: