Generating a random matrix with MvNormal

Is it possible to generate a matrix comprised of N draws from a multivariate normal from within the pm.Model() context?

Here is what I’ve tried

import pymc3 as pm

N = 10
Sigma = np.eye(2)

with pm.Model() as model:
    
    X = pm.MvNormal('X', mu=np.zeros(2), cov = Sigma, shape = N)
    betas = pm.Normal('betas', 0, 1, shape = 2)
    y = pm.Deterministic('y', pm.math.dot(X,betas))
    
    prior_pred = pm.sample_prior_predictive(1)

However, I receive the following error:

ValueError: operands could not be broadcast together with shapes (10,) (2,)

Clearly some sort of shape error, but when I try X=pm.MvNormal('X',mu=np.zeros(2), cov=Sigma, shape=(N,2)) I get

ValueError: 1-dimensional argument does not have enough dimensions for all core dimensions ('i_0_0', 'i_0_1')

What’s the problem here? How can I generate what I want from within the model context?

There are some shape weirdness happening here, but you can try:

with pm.Model() as model: 
  X = pm.MvNormal('X', mu=np.zeros(2), cov=Sigma, shape=(N, 2)) 
  betas = pm.Normal('betas', 0, 1, shape=2) 
  y = pm.Deterministic('y', pm.math.dot(X,betas)) 
  prior_pred = pm.sample_prior_predictive(10)  # <== using something > 1

spoke too fast, still getting wrong answer - definitely a bug

Could you raise an issue?

Yep, will raise now.

1 Like

Issue here.

1 Like