sklearnで多項式回帰（Polynomial Regression）入門から実践まで【PolynomialFeatures】

こんにちは。現役エンジニアの”はやぶさ”@Cpp_Learningです。

多項式回帰（Polynomial Regression）を実践するときに、sklearnのPolynomialFeaturesがとても便利だったので、使い方を紹介します。

本記事のストーリ

線形回帰（Linear Regression）では上手くフィットしなかった
多項式回帰（Polynomial Regression）にトライ
モデリング自動化とオーバーフィットについて

Contents

1 実践！sklearnで線形回帰
2 実践！sklearnで多項式回帰（2次）
- 2.1 scikit-learnのPolynomialFeatures
3 実践！sklearnで多項式回帰（3次）
- 3.1 推定パラメータ取得
4 モデリング自動化とオーバーフィット
5 まとめ

実践！sklearnで線形回帰

回帰とは？については下記で説明済みのため、本記事では割愛します。

【統計モデリング入門】一般化線形モデル（GLM）を基礎から学ぶ -Pythonによる実践あり-統計モデリング（statistical modelling）の入門記事を書きました。線形モデル（Linear Model）と一般化線形モデル（Generalized Linear Model）の理論から実践まで学べます。Pythonライブラリ statsmodels によるソースコードも公開中です。...

以降からソースコードを交えて説明します。

Import

まずはimportから

import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(0)

import numpy as np

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression

from sklearn.metrics import mean_squared_error, r2_score

from sklearn.preprocessing import PolynomialFeatures

np.random.seed(0)

データセット作成

適当なデータセットを作成します。

x = np.linspace(-5, 5, 20)
y = 2 * x - 2 * (x ** 2) + 0.5 * (x ** 3) + np.random.normal(loc=0, scale=10, size=len(x))

# [a,b,c,..] to [[a], [b], [c]]
x = x[:, np.newaxis]
y = y[:, np.newaxis]

x = np.linspace(-5, 5, 20)

y = 2 * x - 2 * (x ** 2) + 0.5 * (x ** 3) + np.random.normal(loc=0, scale=10, size=len(x))

# [a,b,c,..] to [[a], [b], [c]]

x = x[:, np.newaxis]

y = y[:, np.newaxis]

今回は3次式にノイズを加えたデータにしました。

実際の問題では、データがどんな数式に従うのか不明なので、データにフィットする数式モデルを検討することが重要です。

線形回帰（Linear regression）

初手では線形回帰モデルを適用してみます。つまり 一次式：y = b0 + b1x がデータにフィットすると考え、機械学習でパラメータ b0, b1 を推定します。

# モデル：y = b0 + b1x の b0 と b1 を算出  
model = LinearRegression()
model.fit(x, y)

# 推論
y_pred = model.predict(x)

# 評価
rmse = np.sqrt(mean_squared_error(y, y_pred))
r2 = r2_score(y, y_pred)
print(f'rmse : {rmse}')
print(f'R2 : {r2}')

# 可視化
plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

# モデル：y = b0 + b1x の b0 と b1 を算出

model = LinearRegression()

model.fit(x, y)

# 推論

y_pred = model.predict(x)

# 評価

rmse = np.sqrt(mean_squared_error(y, y_pred))

r2 = r2_score(y, y_pred)

print(f'rmse : {rmse}')

print(f'R2 : {r2}')

# 可視化

plt.scatter(x, y)

plt.plot(x, y_pred, color='r')

plt.show()

rmse : 19.313677119446293
R2 : 0.672751565398112

「RMSE と R2 の評価」と「データと線形回帰モデルの関係」を確認したものが上記です。

はやぶさ

イマイチですね

RMSE や R2 などの評価指標については、以下の記事が参考になります。

精度評価指標と回帰モデルの評価｜Instruction of chemoinformatics

実践！sklearnで多項式回帰（2次）

y = b0 + b1x では上手くいかなかったので、 y = b0 + b1x + b1x^2 の多項式を採用します。

scikit-learnのPolynomialFeatures

scikit-learnのPolynomialFeaturesを使えば、簡単に多項式や交互作用の特徴量を生成できます。

例えば degree=2 をセットすると、特徴量 X を以下のように変換してくれます。

【origin】

X = [[x1], [x2], [x3]]

【transform】

X = [[1, x1, x1^2], [1, x2, x2^2], [1, x3, x3^2]]

あとは先ほど同様、多項式：y = b0 + b1x + b2x^2 のパラメータ b0～b2 を機械学習で推定します。

# 2次元の特徴量に変換
polynomial_features= PolynomialFeatures(degree=2)
x_poly = polynomial_features.fit_transform(x)

# y = b0 + b1x + b2x^2 の b0～b2 を算出
model = LinearRegression()
model.fit(x_poly, y)
y_pred = model.predict(x_poly)

# 評価
rmse = np.sqrt(mean_squared_error(y, y_pred))
r2 = r2_score(y, y_pred)
print(f'rmse : {rmse}')
print(f'R2 : {r2}')

# 可視化
plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

# 2次元の特徴量に変換

polynomial_features= PolynomialFeatures(degree=2)

x_poly = polynomial_features.fit_transform(x)

# y = b0 + b1x + b2x^2 の b0～b2 を算出

model = LinearRegression()

model.fit(x_poly, y)

y_pred = model.predict(x_poly)

# 評価

rmse = np.sqrt(mean_squared_error(y, y_pred))

r2 = r2_score(y, y_pred)

print(f'rmse : {rmse}')

print(f'R2 : {r2}')

# 可視化

plt.scatter(x, y)

plt.plot(x, y_pred, color='r')

plt.show()

rmse : 11.007765829915302
R2 : 0.8936968806347743

はやぶさ

イイ感じですね

実践！sklearnで多項式回帰（3次）

最後に degree=3 を試してみます。

# 3次元の特徴量に変換
polynomial_features= PolynomialFeatures(degree=3)
x_poly = polynomial_features.fit_transform(x)

# y = b0 + b1x + b2x^2 + b3x^3  の b0～b3 を算出
model = LinearRegression()
model.fit(x_poly, y)
y_pred = model.predict(x_poly)

# 評価
rmse = np.sqrt(mean_squared_error(y, y_pred))
r2 = r2_score(y, y_pred)
print(f'rmse : {rmse}')
print(f'R2 : {r2}')

# 可視化
plt.scatter(x, y)
plt.plot(x, y_pred, color='r')
plt.show()

# 3次元の特徴量に変換

polynomial_features= PolynomialFeatures(degree=3)

x_poly = polynomial_features.fit_transform(x)

# y = b0 + b1x + b2x^2 + b3x^3 の b0～b3 を算出

model = LinearRegression()

model.fit(x_poly, y)

y_pred = model.predict(x_poly)

# 評価

rmse = np.sqrt(mean_squared_error(y, y_pred))

r2 = r2_score(y, y_pred)

print(f'rmse : {rmse}')

print(f'R2 : {r2}')

# 可視化

plt.scatter(x, y)

plt.plot(x, y_pred, color='r')

plt.show()

rmse : 7.304125160357986
R2 : 0.9531958828956855

データセットを3次式から作成したので、degree=3 が高評価なのは明らかですね。

推定パラメータ取得

機械学習で推定したパラメータ b0～bn は以下のコードで取得できます。

# y = b0 + b1 * x + b2 * x ** 2 + b3 * x ** 3
b0 = model.coef_[0, 0]
b1 = model.coef_[0, 1]
b2 = model.coef_[0, 2]
b3 = model.coef_[0, 3]
intercept = model.intercept_[0]

print("data:  y = 0.0 + 2x - 2x^2 + 0.5x^3 + random.normal")
print(f'model: y = {b0} + {b1}x + {b2}x^2 + {b3}x^3 + {intercept}')