Linear regression in Python (statmodels, scikit-learn, PyMC3)

Introduction

There are several libraries in Python that you can use to do linear regression. I personally needed to do linear regression, and I researched a method for that, so I would like to share it as a memo. The library used is as follows:

• statmodels
• scikit-learn
• PyMC3

Data preparation

First, prepare appropriate data. This time, we will use artificial data with noise $ \ epsilon $ added to the following formula.

The values estimated here are $ \ beta_0, \ beta_1, \ beta_2 $.

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate random data
beta = [1.2, 0.5]
# prep data
x1 = np.random.random(size=1000)*5 
x2 = np.random.random(size=1000)*10 
x = np.transpose([x1, x2])
y = np.dot(x, beta)  + np.random.normal(scale=0.5, size=1000)
#data = dict(x=x, y=y)
data = dict(x1=x1, x2=x2, y=y)
df = pd.DataFrame(data)
# visualisation
plt.scatter(x1, y, color='b')
plt.scatter(x2, y, color='orange')
# 3D
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x1, x2, y)
plt.show()

When using Statsmodels

import statsmodels.api as sm

x = sm.add_constant(x)
results = sm.OLS(endog=y, exog=x).fit()
results.summary()

When using scikit-learn

from sklearn import linear_model
# compare different regressions
lreg = linear_model.LinearRegression()
lreg.fit(x, y)
print("Linear regression: \t", lreg.coef_)
breg = linear_model.BayesianRidge()
breg.fit(x, y)
print("Bayesian regression: \t", breg.coef_)
ereg = linear_model.ElasticNetCV()
ereg.fit(x, y)
print("Elastic net:  \t\t", ereg.coef_)
print("true parameter values: \t", beta)

Linear regression: 	 [ 1.18232244  0.49832431]
Bayesian regression: 	 [ 1.18214701  0.49830501]
Elastic net:  		 [ 1.17795855  0.49756084]
true parameter values: 	 [1.2, 0.5]

When using PyMC3

In the above two cases, the parameters were estimated using maximum likelihood estimation, but in PyMC3, the probability distribution of the parameters is calculated by Markov chain Monte Carlo method (MCMC) based on Bayes' theorem.

import pymc3 as pm

with pm.Model() as model_robust:
    family = pm.glm.families.Normal()
    #pm.glm.glm('y ~ x', data, family=family)
    pm.glm.glm('y ~ x1+x2', df, family=family)
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    #step = pm.Metropolis()
    trace_robust = pm.sample(25000, step)

Optimization terminated successfully.
         Current function value: 764.008811
         Iterations: 19
         Function evaluations: 30
         Gradient evaluations: 30


 100%|██████████| 25000/25000 [00:32<00:00, 761.52it/s]

# show results
pm.traceplot(trace_robust[2000:])
#pm.summary(trace_robust[2000:])
plt.show()
pm.df_summary(trace_robust[2000:])

Try it in a slightly more complicated case

Let's try a little more complicated case to see how much PyMC3 can do. The original formula is:

y = \beta_0 + \beta_1 x_1 + \beta_2 x_2+ \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5 + \epsilon

# Generate random data
beta = [1.2, 0.5, 9.8, 0.2, 1.08]
# prep data
x1 = np.random.random(size=1000)*5 
x2 = np.random.random(size=1000)*10 
x3 = np.random.random(size=1000)
x4 = np.random.normal(size=1000)*2 
x5 = np.random.random(size=1000)*60 
x = np.transpose([x1, x2, x3, x4, x5])
y = np.dot(x, beta)  + np.random.normal(scale=0.5, size=1000)
#data = dict(x=x, y=y)
data = dict(x1=x1, x2=x2, x3=x3, x4=x4, x5=x5, y=y)
df = pd.DataFrame(data)

with pm.Model() as model_robust:
    family = pm.glm.families.Normal()
    pm.glm.glm('y ~ x1+x2+x3+x4+x5', df, family=family)
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    #step = pm.Metropolis()
    trace_robust = pm.sample(25000, step)

# show results
pm.traceplot(trace_robust[2000:])
plt.show()
print("true parameter values are:", beta)
pm.df_summary(trace_robust[2000:])

true parameter values are: [1.2, 0.5, 9.8, 0.2, 1.08]

in conclusion

I didn't explain it because it can be written in short code, but you can see that all of them are able to estimate the parameters well. If you have any questions please feel free to comment.

** (Addition) ** Corrected because the output image was incorrect.