Lineare Regression in Python (Statmodelle, Scikit-Learn, PyMC3)

Einführung

In Python gibt es mehrere Bibliotheken, mit denen Sie eine lineare Regression durchführen können. Ich persönlich musste eine lineare Regression durchführen und habe eine Methode dafür recherchiert, daher möchte ich sie als Memo teilen. Die verwendete Bibliothek lautet wie folgt:

• statmodels
• scikit-learn
• PyMC3

Datenaufbereitung

Bereiten Sie zunächst die entsprechenden Daten vor. Dieses Mal verwenden wir künstliche Daten mit Rauschen $ \ epsilon $, das der folgenden Formel hinzugefügt wird.

Die hier geschätzten Werte sind $ \ 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()

Bei Verwendung von Statistikmodellen

import statsmodels.api as sm

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

Bei Verwendung von 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]

Bei Verwendung von PyMC3

In den beiden obigen Fällen wurden die Parameter unter Verwendung der wahrscheinlichsten Schätzung geschätzt, aber in PyMC3 wird die Wahrscheinlichkeitsverteilung der Parameter durch die Markov-Ketten-Monte-Carlo-Methode (MCMC) basierend auf dem Bayes'schen Theorem berechnet.

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

Versuchen Sie es in einem etwas komplizierteren Fall

Versuchen wir einen etwas komplizierteren Fall, um zu sehen, wie viel PyMC3 leisten kann. Die ursprüngliche Formel lautet:

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]

abschließend

Ich habe es nicht besonders erklärt, weil es in Kurzcode geschrieben werden kann, aber Sie können sehen, dass alle in der Lage sind, die Parameter gut abzuschätzen. Wenn Sie Fragen haben, können Sie diese gerne kommentieren.

** (Ergänzung) ** Korrigiert, weil das Ausgabebild falsch war.