1
$\begingroup$

I'm quite new to applying the techniques of Bayesian regression. I was experimenting with a simple example where I was comparing the results of Bayesian with the actual values. In this simple example, I just used univariate linear regression with a normally distributed error term and X that was uniformly sampled

Here is the data:

data = pd.DataFrame(
        {"X": np.random.RandomState(42).choice(map(lambda x: float(x)/10000.0, 
        np.arange(10000)), 10000, replace=False)})

data["Y"] = 5 + 3*data["X"] + np.random.RandomState(42).normal(0.0, 0.5, 10000)

description: X is a uniform sample drawn from [0,1,2,3...9999]/10000 Y = 5 + 3 * X + N(0,0.5)

I used Pymc3 to find the posterior distributions of the intercept and coefficient term. I used a normally distributed prior.

The complete code is as follows:

import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


data = pd.DataFrame(
        {"X": np.random.RandomState(42).choice(map(lambda x: float(x)/10000.0, np.arange(10000)), 10000, replace=False)})

data["Y"] = 5 + 3*data["X"] + np.random.RandomState(42).normal(0.0, 0.5, 10000)


with pm.Model() as normal_model:
    # The prior for the data likelihood is a Normal Distribution
    family = pm.glm.families.Normal()

    # Creating the model requires a formula and data (and optionally a family)
    pm.GLM.from_formula("Y~X", data=data, family=family)

    start = pm.find_MAP()

    # Perform Markov Chain Monte Carlo sampling letting PyMC3 choose the algorithm
    trace = pm.sample(5000, start=start,tune =1000, random_seed=42, progressbar=True)

pm.traceplot(trace[500:])
plt.show()

#Some convergence plots
fig, axes = plt.subplots(2, 5, figsize=(14,6))
axes = axes.ravel()
for i in range(10):
    axes[i].hist(beta_trace[500*i:500*(i+1)])
plt.tight_layout()
plt.show()


z = geweke(trace, intervals=15)
print(z[0]['X'])
plt.scatter(*(z[0]['X']).T)
plt.hlines([-1,1], 0, 3000, linestyles='dotted')
plt.xlim(0, 3000)
plt.show()

In this case the result was quite close to the actual values.

Estimates and GWEKE convergence for the X coefficient:

Bayesian Estimates

GEWEKE convergence

However, now I change the X vector to:

data = pd.DataFrame(
        {"X": np.random.RandomState(42).choice(map(lambda x: float(x)/100.0, np.arange(10000)), 10000, replace=False)})
data["Y"] = 5 + 3*data["X"] + np.random.RandomState(42).normal(0.0, 0.5, 10000)

description: X is a uniform sample drawn from [0,1,2,3...9999]/100 Y = 5 + 3 * X + N(0,0.5)

Notice that I have used float(x)/100.0 instead of float(x)/10000.0.

In this case although the intercept term is converging to the actual value, the mode of the coefficient is around 0.

Estimates and GWEKE convergence for the X coefficient:

Bayesian Estimates

GEWEKE convergence I looked at the convergence stats and it seems that the samples are converging in both cases. Is this because the gaussian prior is incorrect and not suitable in the second case? Even if I increase the observations and samples, the X coefficient is nowhere close to the actual one. If the prior is incorrect what could be a good prior in this case

Can you please help me understand why results of Bayesian regression coming out to be incorrect? As long as X was uniformly varying between 0-1 it worked fine.

Thanks

$\endgroup$
  • $\begingroup$ Could you add a description of your model, data, and results (that can be understood without reading the code) $\endgroup$ – Juho Kokkala Jul 12 '18 at 9:23
  • 1
    $\begingroup$ We especially need to see the priors, and the exact values of the priors' parameters. Note that 'prior for data likelihood' should just be 'data likelihood'. $\endgroup$ – Frank Harrell Jul 12 '18 at 12:19
  • $\begingroup$ Have a closer look at your graphs: there is a tiny "+3" on the bottom right of the trace plot for X. Type: print(pm.summary(trace[500:])). Voilà, it worked. $\endgroup$ – COOLSerdash Jul 12 '18 at 13:03
  • $\begingroup$ As a remedy for the future, set the offset in matplotlib off. This worked for me: plt.rcParams["axes.formatter.useoffset"] = False $\endgroup$ – COOLSerdash Jul 12 '18 at 13:21
  • $\begingroup$ Thanks! I should have just had a closer look at the plots rather than spending so much time!! They did the offset it one case but not in the other. Lol $\endgroup$ – user3427916 Jul 12 '18 at 13:30
2
$\begingroup$

The x-axis in the traceplot for X in the second example has an offset (see the little +3 on the bottom right?). Setting this offset off reveals that the procedure worked as intended.

pm.traceplot(trace[500:])
plt.rcParams["axes.formatter.useoffset"] = False
plt.show()

Traceplot without offset

We can also print a summary of the posterior distribution:

print(pm.summary(trace[500:]))

               mean        sd  mc_error   hpd_2.5  hpd_97.5    n_eff      Rhat
Intercept  4.987825  0.009945  0.000107  4.969082  5.008067   9537.0  0.999896
X          3.000222  0.000172  0.000002  2.999893  3.000566   9310.0  0.999895
sd         0.501767  0.003523  0.000031  0.494985  0.508702  12527.0  1.000146
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.