Multivariate linear regression in pymc3
<p>I've recently started learning <code>pymc3</code> after exclusively using <code>emcee</code> for ages and I'm running into some conceptual problems. </p>

<p>I'm practising with Chapter 7 of <a href="https://arxiv.org/abs/1008.4686" rel="nofollow">Hogg's Fitting a model to data</a>. This involves a mcmc fit to a straight line with arbitrary 2d uncertainties. I've accomplished this quite easily in <code>emcee</code>, but <code>pymc</code> is giving me some problems. </p>

<p>It essentially boils down to using a multivariate gaussian likelihood.</p>

<p>Here is what I have so far. </p>

<pre><code>from pymc3 import  *

import numpy as np
import matplotlib.pyplot as plt

size = 200
true_intercept = 1
true_slope = 2

true_x = np.linspace(0, 1, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * true_x
# add noise

# here the errors are all the same but the real world they are usually not!
std_y, std_x = 0.1, 0.1 
y = true_regression_line + np.random.normal(scale=std_y, size=size)
x = true_x + np.random.normal(scale=std_x, size=size)

y_err = np.ones_like(y) * std_y
x_err = np.ones_like(x) * std_x

data = dict(x=x, y=y)

with Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
    # Define priors
    intercept = Normal('Intercept', 0, sd=20)
    gradient = Normal('gradient', 0, sd=20)


    # Define likelihood
    likelihood = MvNormal('y', mu=intercept + gradient * x,
                        tau=1./(np.stack((y_err, x_err))**2.), observed=y)

    # start the mcmc!
    start = find_MAP() # Find starting value by optimization
    step = NUTS(scaling=start) # Instantiate MCMC sampling algorithm
    trace = sample(2000, step, start=start, progressbar=False) # draw 2000 posterior samples using NUTS sampling
</code></pre>

<p>This raises the error: <code>LinAlgError: Last 2 dimensions of the array must be square</code></p>

<p>So I'm trying  to pass <code>MvNormal</code> the measured values for x and y (<code>mu</code>s)  and their associated measurement uncertainties (<code>y_err</code> and <code>x_err</code>). But it appears that it is not liking the 2d <code>tau</code> argument.</p>

<p>Any ideas? This must be possible</p>

<p>Thanks</p>

<p>You may try by adapting the following model. Is a "regular" linear regression. But <code>x</code> and <code>y</code> have been replaced by Gaussian distributions. Here I am assuming not only the measured values of the input and output variables but also a reliable estimation of the their error (for example as provided by a measurement device). If you do not trust those error values you may instead try to estimate them from the data.</p>

<pre><code>with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)
    obs_x = pm.Normal('obs_x', mu=x, sd=x_err, shape=len(x))
    obs_y = pm.Normal('obs_y', mu=y, sd=y_err, shape=len(y))

    likelihood = pm.Normal('y', mu=intercept + gradient * obs_x,
                    sd=epsilon, observed=obs_y)

    trace = pm.sample(2000)
</code></pre>

<p>If you are estimating the error from the data it could be reasonable to assume they could be correlated and hence, instead of using two separate Gaussian you can use a Multivariate Gaussian. In such a case you will end up with a model like the following:</p>

<pre><code>df_data = pd.DataFrame(data)
cov = df_data.cov()

with pm.Model() as model:
    intercept = pm.Normal('intercept', 0, sd=20)
    gradient = pm.Normal('gradient', 0, sd=20)
    epsilon = pm.HalfCauchy('epsilon', 5)

    obs_xy = pm.MvNormal('obs_xy', mu=df_data, tau=pm.matrix_inverse(cov), shape=df_data.shape)

    yl = pm.Normal('yl', mu=intercept + gradient * obs_xy[:,0],
                    sd=epsilon, observed=obs_xy[:,1])

mu, sds, elbo = pm.variational.advi(n=20000)
step =  pm.NUTS(scaling=model.dict_to_array(sds), is_cov=True)
trace = pm.sample(1000, step=step, start=mu)
</code></pre>

<p>Notice that in the previous model the covariance matrix was computed from the data. If you are going to do that then I think is better to go with the first model, but if instead you are going to estimate the covariance matrix then the second model could be a sensible approach.</p>

<p>For the second model I use ADVI to initialize it. ADVI can be a good way to initialize models, often it works much better than find_MAP().</p>

<p>You may also want to check this <a href="https://github.com/davidwhogg/DataAnalysisRecipes" rel="nofollow">repository</a> by David Hogg. And the book <a href="http://xcelab.net/rm/statistical-rethinking/" rel="nofollow">Statistical Rethinking</a> where McElreath discuss the problem of doing linear regression including the errors in the input and output variables.</p>

python statistics mcmc pymc3