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.