# Initialzing NUTS with covariance matrix as mass matrix

**URL:** <https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972>\
**Category:** Questions\
**Created:** [October 15, 2019, 7:56pm UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972 "2019-10-15T19:56:05Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![kappie](https://avatars.discourse-cdn.com/v4/letter/k/dc4da7/32.png) [@kappie](https://discourse.pymc.io/u/kappie)\
**Post date:** [October 15, 2019, 7:56pm UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/1 "2019-10-15T19:56:05Z")

</div>

Hello,

I want to initialize the mass matrix of the NUTS sampler with a covariance matrix calculated from the posterior samples of an MCMC chain that has converged. My small example is:

```
import pandas as pd
import numpy as np
import pymc3 as pm

#### Generate data
np.random.seed(123)
alpha, sigma = 1, 1 # True parameter values
size = 100 # Size of dataset

# Simulate outcome variable
Y = alpha + np.random.randn(size)*sigma

basic_model = pm.Model()
with basic_model:

   # Priors for unknown model parameters
   alpha = pm.Normal('alpha', mu=0, sd=10)
   sigma = pm.HalfNormal('sigma', sd=1)

   # Expected value of outcome
   mu = alpha

   # Likelihood (sampling distribution) of observations
   Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
   step1 = pm.NUTS(vars=[alpha,sigma], is_cov=True)
   trace = pm.sample(500, step=[step1], tune=1000, discard_tuned_samples=False)

```

##### creating the dataframe of posterior samples in the predefined order

```
df = pd.DataFrame()
for myname in trace.varnames:
  if myname in ["alpha", "sigma"]:
    k = pd.DataFrame(trace[500:][myname])
    print(myname)
    df = pd.merge(df, k, left_index=True, right_index=True, how='outer')

# var/cov matrix
covmatrix = np.cov(df.values)

```

then if I run my model again only with

```
step1 = pm.NUTS(vars=[alpha, sigma], is_cov=True, scaling=covmatrix)

```

Here in the sampling step I always get

##### “PositiveDefiniteError: Scaling is not positive definite: Simple check failed. Diagonal contains negatives”

Of course in this small and simple example tuning the NUTS would not be necessary. Nevertheless, it should still work, shouldn’t it?  
It is possible to e.g. invert the calculated covariance matrix with numpy, so it is not singular. Thus, I don’t understand the error message. Can anyone help me to get this code functioning?

---

<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 16, 2019, 6:12am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/2 "2019-10-16T06:12:31Z")

</div>

You can follow the code here:

> <https://github.com/pymc-devs/pymc3/blob/75ec55f08b9dbe35a8d94b97d157138c856c6bf3/pymc3/sampling.py#L1568-L1579>

---

<div class="post-metadata">

**Author:** ![kappie](https://avatars.discourse-cdn.com/v4/letter/k/dc4da7/32.png) [@kappie](https://discourse.pymc.io/u/kappie)\
**Post date:** [October 22, 2019, 6:22am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/3 "2019-10-22T06:22:26Z")

</div>

The problem with pm.trace\_cov() is that it creates a covariance matrix of all parameters. In case I am sampling say one parameter with a different sampler, I won’t be able to use it.  
Thus, I am creating my own covariance matrix using np.cov() of only the relevant parameters.  
Which order do the parameters have to be in?  
The order of the trace? (the order of the columns when saving it as text backend)  
or the order that is displyed as an output to the console, e.g. NUTS: [kappa, s, alpha, gamma, sigmasq] ?

---

<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 22, 2019, 7:55am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/4 "2019-10-22T07:55:02Z")

</div>

Oh right I see, yeah that’s a tricky one, as sigma is a bounded random variable, and PyMC3 samples everything in continuous space, so you need to get the sigma after **tranformed**. Maybe try something like:

```python
df = pd.DataFrame()
for myname in trace.varnames:
  if myname in ["alpha", "sigma_log__"]:
    k = pd.DataFrame(trace[500:][myname])
    print(myname)
    df = pd.merge(df, k, left_index=True, right_index=True, how='outer')

# var/cov matrix
covmatrix = np.cov(df.values)

```

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [October 23, 2019, 9:14am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/5 "2019-10-23T09:14:36Z")

</div>

@junpenglao I don’t think that is correct. You still don’t know the correct order.

You can use the dict -\> array mapping from `model.logp_dlogp_function([vars]).dict_to_array`. That is the order NUTS uses internally.

---

<div class="post-metadata">

**Author:** ![kappie](https://avatars.discourse-cdn.com/v4/letter/k/dc4da7/32.png) [@kappie](https://discourse.pymc.io/u/kappie)\
**Post date:** [October 23, 2019, 9:32am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/6 "2019-10-23T09:32:29Z")

</div>

I assume this corresponds to the order in the list trace.varnames?

---

<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 23, 2019, 10:20am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/7 "2019-10-23T10:20:21Z")

</div>

In your case you need to do `model.logp_dlogp_function([alpha, sigma]).dict_to_array` to get the correct mapping.

---

<div class="post-metadata">

**Author:** ![kappie](https://avatars.discourse-cdn.com/v4/letter/k/dc4da7/32.png) [@kappie](https://discourse.pymc.io/u/kappie)\
**Post date:** [October 23, 2019, 10:37am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/8 "2019-10-23T10:37:25Z")

</div>

@junpenglao when I execute this code I get the error:

TypeError  
Traceback (most recent call last)  
 in   
----\> 1 basic\_model.logp\_dlogp\_function([alpha, sigma]).dict\_to\_array()

TypeError: dict\_to\_array() missing 1 required positional argument: ‘point’

What is the ‘point’ argument I have to pass?

---

<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 23, 2019, 11:42am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/9 "2019-10-23T11:42:06Z")

</div>

it’s the fancy way of saying `dict`, it should have the same structure as `point = basic_model.test_point`

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [October 24, 2019, 10:18am UTC](https://discourse.pymc.io/t/initialzing-nuts-with-covariance-matrix-as-mass-matrix/3972/10 "2019-10-24T10:18:23Z")

</div>

And by the way, you could also try [https://github.com/aseyboldt/covadapt](https://github.com/aseyboldt/covadapt)  
I’m working on further improving and testing this right now, but even as it is it should work better than initializing with the empirical covariance.