前回記事

前回まで

線形回帰モデルの定義と最小二乗法、勾配降下法という2つのアルゴリズムでの解き方を学びました。
今回は実際にコードを書いて線形回帰を実行していきます。

データ準備

今回使うデータはペンギンデータセットです。

ペンギンデータセットは3種類のペンギンの生体情報がデータになったもので、irisデータセットにかわる機械学習の入門的データセットと呼べるものです。
以下にデータに含まれる項目の説明を記載します。

species：ペンギンの種類を表す特徴量です。Adelie（アデリーペンギン）、Chinstrap（ヒゲペンギン）、Gentoo（ジェンツーペンギン）の3種類が含まれています。
island：ペンギンの生息地の島を表す特徴量です。Torgersen、Biscoe、Dreamの3種類が含まれています。
bill_length_mm：ペンギンのくちばしの長さ（単位：mm）を表す数値特徴量です。
bill_depth_mm：ペンギンのくちばしの奥行き（単位：mm）を表す数値特徴量です。
flipper_length_mm：ペンギンのひれの長さ（単位：mm）を表す数値特徴量です。
body_mass_g：ペンギンの体重（単位：g）を表す数値特徴量です。
sex：ペンギンの性別を表す特徴量です。Male、Femaleの2種類が含まれています。
year：データが収集された年を表す数値特徴量です。

実際にデータをロードして散布図プロットしてみます。

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# ペンギンデータセットをロード
penguins = sns.load_dataset('penguins')

# NaN値を除去
penguins = penguins.dropna(subset=['bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'body_mass_g'])

# カラーマップを作成
color_map = {'Adelie': 'red', 'Chinstrap': 'green', 'Gentoo': 'blue'}
colors = penguins['species'].map(color_map)

# 散布図行列を作成
scatter_matrix = pd.plotting.scatter_matrix(penguins,c=colors, figsize=(10, 10))

# グラフのタイトルを追加
plt.suptitle("Penguins Dataset - Scatter Matrix", fontsize=16)

# グラフのラベルが重ならないように調整
plt.tight_layout()

# 図を保存
plt.savefig('scatter_matrix.png')

# 散布図行列を表示
plt.show()

まずseabornを使ってペンギンデータセットをロードします。ロードしたデータはpandas.DataFrame形式で返却されます。
データにはNaN値が含まれているのでdropnaでNaN値が含まれるレコードを削除します。
pandas.plotting.scatter_matrixにデータを渡すと数値項目を抽出し散布図マトリックスを作成してくれます。（体格成分は各項目のヒストグラムになります。）

最小二乗法

ペンギンデータセットを使って最小二乗法により線形回帰を実行するコードを作成します。
データはジェンツーペンギンのくちばしの長さと体重を使います。

# ジェンツーペンギンのデータを抽出
gentoo_penguins = penguins[penguins['species'] == 'Gentoo']

# くちばしの長さと体重のデータを抽出
bill_length = gentoo_penguins['bill_length_mm']
body_mass = gentoo_penguins['body_mass_g']

# 散布図をプロット
plt.scatter(bill_length, body_mass)
plt.xlabel('Bill Length (mm)')
plt.ylabel('Body Mass (g)')
plt.title('Gentoo Penguins - Bill Length vs Body Mass')

# 図を保存
plt.savefig('scatter_gentoo.png')

# グラフを表示
plt.show()

penguinsデータフレームの項目'species'で'Gentoo'を指定しジェンツーペンギンのみのデータセットにします。
'bill_length_mm'と'body_mass_gを抽出し2次元プロットすると、くちばしが長い個体ほど体重も重い関係が確認できます。

続けて最小二乗法による線形回帰を実行します。scikit-learnではsklearn.linear_model.LinearRegressionを使って最小二乗法による学習を行う事ができます。

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

# 線形回帰モデルの訓練
reg = LinearRegression()
reg.fit(bill_length.values.reshape(-1, 1), body_mass)

# 線形回帰の結果（パラメータ）を出力
slope = reg.coef_[0]
intercept = reg.intercept_
print("Slope: ", slope)
print("Intercept: ", intercept)

# 平均二乗誤差を計算
predictions = reg.predict(bill_length.values.reshape(-1, 1))
mse = mean_squared_error(body_mass, predictions)
print("Mean Squared Error: ", mse)

# 散布図をプロット
plt.scatter(bill_length, body_mass, label='Data')
plt.plot(bill_length, reg.predict(bill_length.values.reshape(-1, 1)), color='red', label='Linear Regression')
plt.xlabel('Bill Length (mm)')
plt.ylabel('Body Mass (g)')
plt.title('Gentoo Penguins - Bill Length vs Body Mass')
plt.legend()

# グラフを表示
plt.show()

結果

Slope:  109.45916306693874
Intercept:  -123.82793265390046
Mean Squared Error:  139195.62409167975

まずLinearRegression()によりインスタンスを作成します。コード中は記載していませんがパラメタを指定することもできます。

fit_intercept :モデルの切片（定数項$\theta_0$）を計算するかどうかを指定します。デフォルトはTrueで、Falseに設定すると、切片は計算されません。
copy_X:入力データをコピーして使うかどうかを指定します。デフォルトはTrueで入力データは変更されませんが、Falseにすると入力データが標準化されます。（元のデータが更新される。）
n_jobs:使用するジョブ数を指定します。デフォルトは'None'で、一つのCPUコアを使用しますが、数値を設定するとその数だけCPUコアを使用します。
positive:パラメータの符号を指定します。デフォルトはFalseで符号は指定しませんが、Trueにするとパラメータの符号が正になります。

インスタンス作成後、LinearRegression.fitというメソッドで最小二乗法の計算を行います。引数は以下の通りです。

第一引数X : 線形回帰の予測変数（または独立変数、特徴量）
第二引数y : 線形回帰の応答変数（または従属変数、目的変数）
第三引数sample_weight（オプション） : データを均一ではなく重みをつけて学習させたい場合に使用します。一部のデータが重要な不均衡なデータセットの場合等に使用しますが基本的には指定は不要です。

学習後はインスタンス変数coef_に切片以外のパラメタ（$\theta_1\cdots$\theta_n）が、intercept_に切片が入ります。これらのパラメタから各入力に対する予測値を計算することもできますがLinearRegression.predictメソッドを使っても計算ができます。（引数は入力データXです）

平均二乗誤差については、sklearn.metrics.mean_squared_errorメソッドを用いて計算します。
なおLinearRegressionクラスでは正則化は使えないので、後述するRidge等のクラスを使用する必要があります。

正則化

sklearnでL2回帰を実行するにはsklearn.linear_model.Ridgeクラスを使用します。

from sklearn.linear_model import Ridge

# 線形回帰モデルの訓練
rg = Ridge()
rg.fit(bill_length.values.reshape(-1, 1), body_mass)

# 線形回帰の結果（パラメータ）を出力
slope = rg.coef_[0]
intercept = rg.intercept_
print("Slope: ", slope)
print("Intercept: ", intercept)

# 平均二乗誤差を計算
predictions = rg.predict(bill_length.values.reshape(-1, 1))
mse = mean_squared_error(body_mass, predictions)
print("Mean Squared Error: ", mse)

結果

Slope:  109.36478032657529
Intercept:  -119.3442920830239
Mean Squared Error:  139195.70801157953

先程と同様にRidge()でインスタンスを生成します。主なパラメータは下記の通りです。

alpha：正則化の強度を制御するパラメータで、デフォルトは1.0です。alphaが大きいほど、正則化の強度が高くなり、係数の値は0に近づきます。alphaが小さい場合、正則化の効果は弱まり、リッジ回帰は標準の最小二乗回帰に近くなります。

その後Ridge.fitでL2正則化の学習を行います。coef_とintercept_で学習したパラメータを確認すると最小二乗法の場合とほぼ同じ結果が得られました。

外れ値がある場合の比較

今回のペンギンデータセットの場合、最小二乗法とL2正則化での結果に大きな差がありませんでした。これは正則化で大きなペナルティとなるような外れ値が少なかったためと考えられます。そこで、人工的に外れ値データを追加し最小二乗法とL2正則化での学習結果の違いを確認していきます。

import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_squared_error

# データセットのロード
penguins = sns.load_dataset('penguins')

# NaN値を除去
penguins = penguins.dropna(subset=['bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'body_mass_g'])

# Gentooペンギンのデータのみを抽出
gentoo = penguins[penguins.species == 'Gentoo']

# くちばしの長さと体重データを取得
X = gentoo[['bill_length_mm']]
y = gentoo['body_mass_g']

# 外れ値の数
num_outliers = 5

# 外れ値を追加 (ここではくちばしの長さ 50 mm ±5、体重 9000 g ±500 付近のランダムデータを追加)
np.random.seed(0)
outlier_X = pd.DataFrame({'bill_length_mm': 50 + np.random.randn(num_outliers) * 5})
outlier_y = pd.Series(9000 + np.random.randn(num_outliers) * 500)
X = pd.concat([X, outlier_X], ignore_index=True)
y = pd.concat([y, outlier_y], ignore_index=True)

# Linear Regressionモデルの作成と訓練
linear_model = LinearRegression()
linear_model.fit(X, y)

# Ridge Regressionモデルの作成と訓練
ridge_model = Ridge(alpha=1000)
ridge_model.fit(X, y)

# 元のデータと外れ値の予測
y_pred_linear = linear_model.predict(X)
y_pred_ridge = ridge_model.predict(X)

# 平均二乗誤差を計算
mse_lm = mean_squared_error(y, y_pred_linear)
mse_rm = mean_squared_error(y, y_pred_ridge)

# プロット作成
plt.figure(figsize=(10, 6))
plt.scatter(X[:-num_outliers], y[:-num_outliers], color='blue', label='Original data')  # 元のデータ
plt.scatter(X[-num_outliers:], y[-num_outliers:], color='red', label='Outliers')  # 外れ値
plt.plot(X, y_pred_linear, color='green', label='Linear Regression')  # Linear Regressionの結果
plt.plot(X, y_pred_ridge, color='purple', label='Ridge Regression')  # Ridge Regressionの結果
plt.xlabel('Bill Length (mm)')
plt.ylabel('Body Mass (g)')
plt.title('Effect of Outliers on Linear and Ridge Regression')
plt.legend()

# 図を保存
plt.savefig('scatter_lr_rr.png')

# グラフを表示
plt.show()

# パラメータの比較
slope_lm = linear_model.coef_[0]
intercept_lm = linear_model.intercept_
print("Slope in Linear Regression: ", slope_lm)
print("Intercept in Linear Regression: ", intercept_lm)
print("Mean Squared Error in Linear Regression: ", mse_lm)

slope_rm = ridge_model.coef_[0]
intercept_rm = ridge_model.intercept_
print("Slope in Ridge Regression: ", slope_rm)
print("Intercept in Ridge Regression: ", intercept_rm)
print("Mean Squared Error in Ridge Regression: ", mse_rm)

結果

Slope in Linear Regression:  184.55200078405588
Intercept in Linear Regression:  -3607.5834391099033
Mean Squared Error in Linear Regression:  383452.26191625826
Slope in Ridge Regression:  115.45017982770396
Intercept in Ridge Regression:  -298.60117698488193
Mean Squared Error in Ridge Regression:  445778.9623524508

先程と違う点として、くちばしの長さ 50 mm ±5、体重 9000 g ±500 付近にランダムな人工データを5点追加しています。このデータで学習すると、回帰式は先程より外れ値側に引っ張られたものになりますが、L2回帰のほうがより影響が少ないものになっています。（紫の直線のほうが緑の直線よりもとのデータに当てはまりが良い）

また平均二乗誤差を比較すると最小二乗法のほうが小さくなっている事がわかります。L2回帰では誤差は大きくなっていますが外れ値の影響は小さくなっていることが確認できます。