# Multivariate Normal code does not work since version 3.2

**URL:** <https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456>\
**Category:** Development\
**Created:** [October 20, 2017, 2:44pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456 "2017-10-20T14:44:31Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![viiv](https://avatars.discourse-cdn.com/v4/letter/v/f9ae1b/32.png) [@viiv](https://discourse.pymc.io/u/viiv)\
**Post date:** [October 20, 2017, 2:44pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456/1 "2017-10-20T14:44:31Z")

</div>

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](https://github.com/cliburn/sta-663-2016/blob/master/lectures/16C_PyMC3.ipynb))

```
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.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [October 20, 2017, 4:49pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456/2 "2017-10-20T16:49:17Z")

</div>

I think this might be a bug. Reported to [Github issue](https://github.com/pymc-devs/pymc3/issues/2644).

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [October 20, 2017, 4:53pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456/3 "2017-10-20T16:53:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [October 20, 2017, 5:48pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456/4 "2017-10-20T17:48:40Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![viiv](https://avatars.discourse-cdn.com/v4/letter/v/f9ae1b/32.png) [@viiv](https://discourse.pymc.io/u/viiv)\
**Post date:** [October 30, 2017, 3:54pm UTC](https://discourse.pymc.io/t/multivariate-normal-code-does-not-work-since-version-3-2/456/5 "2017-10-30T15:54:29Z")

</div>

Both these made the code work, thanks!