PyMC3 Bayesian Linear Regression prediction with sklearn.datasets
<p>I've been trying to implement <strong>Bayesian Linear Regression</strong> models using <code>PyMC3</code> with <em>REAL DATA</em> (i.e. not from linear function + gaussian noise) from the datasets in <code>sklearn.datasets</code>.  I chose the regression dataset with the smallest number of attributes (i.e. <code>load_diabetes()</code>) whose shape is <code>(442, 10)</code>; that is, <code>442 samples</code> and <code>10 attributes</code>. </p>

<p>I believe I got the model working, the posteriors look decent enough to try and predict with to figure out how this stuff works but...I realized I have no idea how to predict with these Bayesian Models!  I'm trying to avoid using the <code>glm</code> and <code>patsy</code> notation because it's difficult for me to understand what is actually going on when using that. </p>

<p>I tried following: 
<a href="https://stackoverflow.com/questions/33249122/generating-predictions-from-inferred-parameters-in-pymc3">Generating predictions from inferred parameters in pymc3</a>
and also <a href="http://pymc-devs.github.io/pymc3/posterior_predictive/" rel="noreferrer">http://pymc-devs.github.io/pymc3/posterior_predictive/</a> but my model is either extremely terrible at predicting or I'm doing it wrong. </p>

<p>If I actually am doing the prediction correctly (which I'm probably not) then can anyone help me <em>optimize</em> my model. I don't know if least <code>mean squared error</code>, <code>absolute error</code>, or anything like that works in Bayesian frameworks.  Ideally, I would like to get an array of number_of_rows = the amount of rows in my <code>X_te</code> attribute/data test set, and the number of columns to be samples from the posterior distribution.</p>

<pre><code>import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from scipy import stats, optimize
from sklearn.datasets import load_diabetes
from sklearn.cross_validation import train_test_split
from theano import shared

np.random.seed(9)

%matplotlib inline

#Load the Data
diabetes_data = load_diabetes()
X, y_ = diabetes_data.data, diabetes_data.target

#Split Data
X_tr, X_te, y_tr, y_te = train_test_split(X,y_,test_size=0.25, random_state=0)

#Shapes
X.shape, y_.shape, X_tr.shape, X_te.shape
#((442, 10), (442,), (331, 10), (111, 10))

#Preprocess data for Modeling
shA_X = shared(X_tr)

#Generate Model
linear_model = pm.Model()

with linear_model: 
    # Priors for unknown model parameters    
    alpha = pm.Normal("alpha", mu=0,sd=10)
    betas = pm.Normal("betas", mu=0,#X_tr.mean(), 
                               sd=10, 
                               shape=X.shape[1])
    sigma = pm.HalfNormal("sigma", sd=1)

    # Expected value of outcome
    mu = alpha + np.array([betas[j]*shA_X[:,j] for j in range(X.shape[1])]).sum()

    # Likelihood (sampling distribution of observations)
    likelihood = pm.Normal("likelihood", mu=mu, sd=sigma, observed=y_tr)

    # Obtain starting values via Maximum A Posteriori Estimate
    map_estimate = pm.find_MAP(model=linear_model, fmin=optimize.fmin_powell)

    # Instantiate Sampler
    step = pm.NUTS(scaling=map_estimate)

    # MCMC
    trace = pm.sample(1000, step, start=map_estimate, progressbar=True, njobs=1)


#Traceplot
pm.traceplot(trace)
</code></pre>

<p><a href="https://i.stack.imgur.com/VcMXs.png" rel="noreferrer"><img src="https://i.stack.imgur.com/VcMXs.png" alt="enter image description here"></a></p>

<pre><code># Prediction
shA_X.set_value(X_te)
ppc = pm.sample_ppc(trace, model=linear_model, samples=1000)

#What's the shape of this? 
list(ppc.items())[0][1].shape #(1000, 111) it looks like 1000 posterior samples for the 111 test samples (X_te) I gave it

#Looks like I need to transpose it to get `X_te` samples on rows and posterior distribution samples on cols

for idx in [0,1,2,3,4,5]:
    predicted_yi = list(ppc.items())[0][1].T[idx].mean()
    actual_yi = y_te[idx]
    print(predicted_yi, actual_yi)
# 158.646772735 321.0
# 160.054730647 215.0
# 149.457889418 127.0
# 139.875149489 64.0
# 146.75090354 175.0
# 156.124314452 275.0 
</code></pre>