Régression linéaire en Python (statmodels, scikit-learn, PyMC3)

introduction

Il existe plusieurs bibliothèques en Python que vous pouvez utiliser pour effectuer une régression linéaire. J'avais personnellement besoin de faire une régression linéaire, et j'ai recherché une méthode pour cela, alors j'aimerais la partager sous forme de mémo. La bibliothèque utilisée est la suivante:

• statmodels
• scikit-learn
• PyMC3

Préparation des données

Tout d'abord, préparez les données appropriées. Cette fois, nous utiliserons des données artificielles avec du bruit $ \ epsilon $ ajouté à la formule suivante.

Les valeurs estimées ici sont $ \ 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()

Lors de l'utilisation de Statsmodels

import statsmodels.api as sm

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

Lors de l'utilisation de 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]

Lors de l'utilisation de PyMC3

Dans les deux cas ci-dessus, les paramètres ont été estimés en utilisant l'estimation la plus probable, mais dans PyMC3, la distribution de probabilité des paramètres est calculée par la méthode de Monte Carlo en chaîne de Markov (MCMC) basée sur le théorème bayésien.

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:])

Essayez-le dans un cas un peu plus compliqué

Essayons un cas un peu plus compliqué pour voir tout ce que PyMC3 peut faire. La formule originale est:

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]

en conclusion

Je ne l'ai pas expliqué en particulier car il peut être écrit en short code, mais vous pouvez voir que tous sont capables d'estimer bien les paramètres. Si vous avez des questions, n'hésitez pas à commenter.

** (Addition) ** Corrigé car l'image de sortie était incorrecte.