Multivariate Normal code does not work since version 3.2

Development
1

Hello,

A few days ago - before the release of PyMC 3.2 - I was testing the support of multivariate Normal distributions with the code below (from here)

import pymc3 as pm
import numpy as np
from matplotlib.patches import Ellipse
import scipy.stats as stats
import theano.tensor as T
import matplotlib.pyplot as plt

def make_ellipse(mu, cov, ci=0.95):
    e, v = np.linalg.eig(cov)
    angle = 180/np.pi * np.arccos(v[0, 0])
    q = stats.chi2(2).ppf(ci)
    e = Ellipse(mu, 2*np.sqrt(q*e[0]), 2*np.sqrt(q*e[1]), angle=angle)
    return e

#%%
np.random.seed(123)
n, p = 200, 2
mu_true = np.random.uniform(-5, 5, 2)
cov_sqrt = np.random.uniform(0,2, (2,2))
cov_true = cov_sqrt.T @ cov_sqrt
y = np.random.multivariate_normal(mu_true, cov_true, n)

#%%
niter = 1000
with pm.Model() as mvn:
    # priors on standard deviation
    sigma = pm.Lognormal('sigma', mu=np.zeros(2), tau=np.ones(2), shape=2)
    # prior on LKJ shape
    nu = pm.Uniform('nu', 0, 5)
    # LKJ prior for correlation matrix as upper triangular vector
    C_triu = pm.LKJCorr('C_triu', n=nu, p=2)
    # convert to matrix form
    C = T.fill_diagonal(C_triu[np.zeros((2,2), 'int')], 1)
    sigma_diag = T.nlinalg.diag(sigma)
    # indduced covariance matrix
    cov = pm.Deterministic('cov', T.nlinalg.matrix_dot(sigma_diag, C, sigma_diag))
    tau = pm.Deterministic('tau', T.nlinalg.matrix_inverse(cov))
    mu = pm.MvNormal('mu', mu=0, tau=tau, shape=2)
    Y_obs_ = pm.MvNormal('Y_obs', mu=mu, tau=tau, observed=y)

    start = pm.find_MAP()

    trace = pm.sample(niter, start=start)

It was running well in PyMC 3.1, but now I get this error:

ValueError: Input dimension mis-match. (input[0].shape[0] = 2, input[1].shape[0] = 1)
Apply node that caused the error: Elemwise{Composite{((-log(i0)) + i1 + (i2 * i0))}}[(0, 0)](Elemwise{exp}.0, TensorConstant{(1,) of -1.0}, (d__logp_nojac/dsigma))
Toposort index: 144
Inputs types: [TensorType(float64, vector), TensorType(float64, vector), TensorType(float64, vector)]
Inputs shapes: [(2,), (1,), (2,)]
Inputs strides: [(8,), (8,), (8,)]
Inputs values: [array([ 1.,  1.]), array([-1.]), array([ 1058.3009407 ,  1131.60014677])]
Outputs clients: [[Join(TensorConstant{0}, (d__logp_nojac/dmu), Rebroadcast{1}.0, Rebroadcast{1}.0, (d__logp_nojac/dsigma_log__))]]

I don’t really get where it comes from, and from the changelog of PyMC, it didn’t seem to me that the support of multivariate distribution changed that much…

I was experimenting with that example before moving to my real problem with many variables and many observables, so I would like to know which version of PyMC to use before starting if possible.

2

I think this might be a bug. Reported to Github issue.

3

FYI, in your case change
sigma = pm.Lognormal('sigma', mu=np.zeros(2), tau=np.ones(2), shape=2)
into
sigma = pm.Lognormal('sigma', mu=np.zeros(2), tau=1, shape=2)
would work.

4

Also, this is something unrelated but you should use the pm.LKJCholeskyCov, and avoid using pm.find_MAP() to get starting value (using the default is better in most cases).

5

Both these made the code work, thanks!