シリーズPython⑩ リッジ回帰を少々Pythonで～書籍「ガウス過程と機械学習」インスパイア
はじめに
この記事は、Pythonの機械学習ライブラリで「リッジ回帰」を実装する際の豆知識を取り扱います。
具体的には「scikit-learn」ライブラリと「statsmodels」ライブラリのリッジ回帰の使い方を記載しています。
記事を書くきっかけは、書籍「ガウス過程と機械学習」の第１章でリッジ回帰を学んだときのある出来事からです・・・
書籍では行列計算を用いてリッジ回帰の推定を行います。
「Pythonライブラリのリッジ回帰は書籍の結果と一致するはず」
と軽いノリで２つのライブラリを動かしたところ、書籍と全然違う結果が出力されて、ドツボにはまりました。。。
Webサイトを右往左往して調べたところ、２つのライブラリの使い方が分かりましたので、メモ代わりに記事にしました。
ガウス過程と機械学習の紹介
この記事は書籍「ガウス過程と機械学習」（講談社、以下「テキスト」と呼びます）の学びをきっかけにして書きました。
テキストは2019年3月に発売され、その名の通り、ガウス過程をとことん学べる書籍です。
数式がたくさん並んでいますがご安心下さい。
丁寧に数式の解き方が説明されているとともに、グラフ・図をふんだんに用いて、ガウス諸々のお気持ちをイメージしやすくなっています。
そして圧巻は「第０章」。
章タイトルは「たった５分でガウス過程法が分かってしまう」です！
ガウス過程を読み解きやすくする仕掛けがふんだんに盛り込まれているのです！
またロジックを補完するPythonのサンプルコードは、書籍のWebサイトで公開されています。
ただ、個人的にはコードの分量がちょっと物足りないかな、と感じています。
そこで、リッジ回帰のPython実装に寄り道しました・・・ところ、ドツボにハマったのです・・・
まえがき
この記事は、出典に記載の書籍に掲載された数式、文章及びデータを引用し、適宜、掲載文章を改変して書いています。
Python のコードは趣味的に作成されたものであり、テキストの著者の意図を反映したものではなく、処理の正確性は担保しておりません。
誤りや改善点がありましたら、ぜひ教えてください。
【出典】
「ガウス過程と機械学習」（第10刷、持橋大地・大羽成征 著、講談社）
第１章 1.4節 リッジ回帰
1. 重回帰とリッジ回帰
① 重回帰
テキストは線形回帰モデルの回帰係数を行列計算で算出するアプローチを採っています。
【公式1.4 重回帰モデルの回帰係数の推定】
$$
\boldsymbol{w} = (\boldsymbol{X}^{\top} \boldsymbol{X})^{-1} \boldsymbol{X}^{\top} \boldsymbol{y}
$$
$${\boldsymbol{X}}$$は説明変数の行列、$${\boldsymbol{y}}$$は目的変数のベクトル、$${\boldsymbol{w}}$$は回帰係数のベクトル、$${\top}$$は行列の転置、$${^{-1}}$$は逆行列です。
じつは、逆行列$${(\boldsymbol{X}^{\top} \boldsymbol{X})^{-1}}$$が存在しない場合、上記公式では計算ができません。
また、逆行列$${(\boldsymbol{X}^{\top} \boldsymbol{X})^{-1}}$$の値次第で、計算結果が不安定になることもあるそうです。
② リッジ回帰
そこでテキストは、計算を可能にするための手段になりうる手法として、リッジ回帰を紹介しています。
リッジ回帰のアイデアは、行列$${(\boldsymbol{X}^{\top} \boldsymbol{X})}$$にパラメータ$${\alpha}$$を加算することで逆行列が存在するように仕向けることです。
$${\alpha}$$はリッジ回帰の正則化パラメータと呼ばれるものです。
テキストの公式をお借りします。
【公式1.6 リッジ回帰の回帰係数の推定】
$$
\boldsymbol{w} = (\boldsymbol{X}^{\top} \boldsymbol{X} + \alpha \boldsymbol{I})^{-1} \boldsymbol{X}^{\top} \boldsymbol{y}
$$
$${\boldsymbol{I}}$$は単位行列です（ numpy の eye() です）。
2. テキストの公式を用いたリッジ回帰の実装
最初に公式1.6 をPythonコードに変換して、リッジ回帰を実践します。
まずは、この記事のコードで使用するライブラリをインポートします。
statsmodels は formula api を利用します。
### インポート
# 数値計算
import numpy as np
import pandas as pd
# リッジ回帰
from sklearn.linear_model import Ridge
import statsmodels.formula.api as smf
# 評価指標
from sklearn.metrics import root_mean_squared_error続いて、テキストの 公式1.6 と 例7 の以下の条件を引用させていただき、リッジ回帰を実装します。
・計画行列$${\boldsymbol{X}}$$
・$${\boldsymbol{y}=(1,\ 2,\ 3)^{\top}}$$
・$${\alpha=0.1}$$
### 例7 リッジ回帰による回帰係数wの推定
## 設定
alpha = 0.1 # 正則化パラメータα
## データの作成
# 計画行列Xの作成 ※第1列は定数項
X = np.array([[1, 2, 4],
[1, 3, 6],
[1, 4, 8]])
# yの作成 shape=(3, 1)
y = np.array([[1, 2, 3]]).T
## リッジ回帰の実行
# 回帰係数の推定 公式1.6
w = np.linalg.inv(X.T @ X + alpha * np.eye(len(X))) @ X.T @ y
print('回帰係数: ', w.T[0])
# 学習データによる予測
y_pred = (X @ w).T[0]
print('予測値 train: ', y_pred)
# RMSEによる評価
print('RMSE train: ', np.sqrt(sum((y.T[0] - y_pred)**2) / len(y)))【実行結果】
予測値は真値$${[1,\ 2,\ 3]}$$に近い感じがいたします。
書籍のリッジ回帰の公式を用いたこの結果をベースにして、Pythonの２つのライブラリの実装を確認していきます。
3. Pythonライブラリのリッジ回帰 その１ テキストに合わせる
scikit-learn と statsmodels でリッジ回帰を実装します。
「その１」では、書籍の公式を用いた結果と一致する方法を探ります。
① scikit-learn
scikit-learn のリッジ回帰は「Ridge クラス」で実装します。
ひとまずコードから。
### scikit-learnでリッジ回帰 ※計画行列に定数項を設定
## データの作成
# 説明変数Xの作成 ※定数項を含む
X = np.array([[1, 2, 4],
[1, 3, 6],
[1, 4, 8]])
# 目的変数の作成
y = np.array([1, 2, 3])
## リッジ回帰モデルの学習
reg = Ridge(alpha=0.1, fit_intercept=False)
reg.fit(X, y)
print('回帰係数: ', reg.coef_)
## 学習データによる予測
# 予測
y_pred = reg.predict(X)
print('予測値 train: ', y_pred)
# RMSEによる評価
print('RMSE train: ', root_mean_squared_error(y, y_pred))【実行結果】
書籍の結果と合っています。
【コードのポイント】
「定数項」の扱いにご注目ください。
このコードでは、まず、「説明変数$${X}$$の中に定数項を含めています」。
そして、リッジ回帰 Ridge の切片パラメータ fit_intercept に False を設定して、ロジックによる定数項考慮をやめるようにしています。
【次なる試練】
上のコードはテキストに合わせる目的で、説明変数$${X}$$の中に定数項を含めました。
しかしscikit-learn を通常利用するケースでは、一般的に「説明変数自体に定数項を含めない」＆「ロジックで定数項を考慮する」方法でモデリングすると思います。
したがいまして、上のコードのような方法を採るのは例外的でしょう。
一般的な方法は「4. その２ 一般的な使用法」節で検討します。
② statsmodels
statsmodelsのリッジ回帰専用機能を探したのですが、見つけることができませんでした。
どうやら、ElasticNet回帰を転用して、リッジ回帰を実装する模様です。
この転用の際、ElasticNet回帰のお気持ち（数式）を理解しつつ、リッジ回帰に変換する必要があるようです。
変換の際には、次のサイトの情報を活用させていただきました。
ありがとうございます！
ではコードに進みます。
ols（最小二乗法）に対して「fit_regularized()」を適用することで、リッジ回帰を実行できます。
### statsmodelsのリッジ回帰 ※定数項にpenaltyを設定
# 説明変数と目的変数をデータフレーム化 ※定数項を含む
df = pd.DataFrame(dict(x1=[2, 3, 4], x2=[4, 6, 8], y=[1, 2, 3]))
# リッジ回帰の実行 ★penaltyを追加
alpha = 0.1
penalty = np.array([alpha, alpha, alpha]) / len(df) # 定数項(第1要素)にalphaを設定
result = (smf.ols(formula='y ~ 1 + x1 + x2', data=df)
.fit_regularized(alpha=penalty, L1_wt=0)) # リッジの場合、L1_wt=0
# 結果の表示
print('回帰係数 : ', result.params)
print('予測値 train: ', result.fittedvalues)
print('RMSE train: ', root_mean_squared_error(df.y, result.fittedvalues))【実行結果】
書籍の結果と合っています。
【コードのポイント】
ポイントは「penalty」（正則化）の工夫です。
ElasticNet 回帰とリッジ回帰は、正則化のスケールが異なるようです。
そこで、先程の参考サイトにならって、正則化パラメータ$${\alpha}$$に次の工夫をします。
・説明変数（定数項を含む）ごとに正則化パラメータ$${\alpha}$$を設定
・正則化の値には$${\alpha/n}$$（$${n}$$: データの数）を設定
コードでは、penalty 変数に上述の正則化パラメータを設定しました。
penalty = np.array([alpha, alpha, alpha]) / len(df)penalty 変数を「fit_regularized()」の引数 alpha に設定します。
result = (smf.ols(formula='y ~ 1 + x1 + x2', data=df)
.fit_regularized(alpha=penalty, L1_wt=0)) # リッジの場合、L1_wt=0なお、ElasticNet 回帰を使用するにあたって、L1正則化（ラッソ回帰のアレです）を off にしたいので、「fit_regularized()」の L1重み引数 L1_wt に 0 を設定します。
4. Pythonライブラリのリッジ回帰 その２ 一般的な使用法
私の小さな経験論で恐縮ですが、機械学習を適用する場面で、説明変数自体に定数項を含める経験はありませんでした。
ですので、前節「3. その１ テキストに合わせる」の実践内容を正直「納得できない」と感じていました。
ここからが本番だ！と意気込んで、進めて行きたいと思います！
① scikit-learn
scikit-learn の一般的な使い方「説明変数自体に定数項を含めない」ケースを実装します。
まずはコードから。
### scikit-learnでリッジ回帰 ※切片をパラメータで設定
## データの作成
# 説明変数Xの作成 ※定数項を含まない
X = np.array([[2, 4],
[3, 6],
[4, 8]])
# 目的変数の作成
y = np.array([1, 2, 3])
## リッジ回帰モデルの学習
reg = Ridge(alpha=0.1, fit_intercept=True) # fit_interceptで定数項を追加
reg.fit(X, y)
print('回帰係数: ', np.hstack([reg.intercept_, reg.coef_]))
## 学習データによる予測
# 予測
y_pred = reg.predict(X)
print('予測値 train: ', y_pred)
# RMSEによる評価
print('RMSE train: ', root_mean_squared_error(y, y_pred))【実行結果】
書籍の結果とは異なります。
こちらの方がRMSEが小さく、予測値が真値に近いです。
【コードのポイント】
その１との違いは次の２点です。
・説明変数$${X}$$自体に定数項を含めない
・Ridge() クラスの引数 fit_intercept で定数項を含めるよう設定
この scikit-learn の結果と一致するように statsmodels を実装してみましょう。
② statsmodels
まずはコードから。
### statsmodelsのリッジ回帰 ※定数項にpenaltyを未設定
# 説明変数と目的変数をデータフレーム化
df = pd.DataFrame(dict(x1=[2, 3, 4], x2=[4, 6, 8], y=[1, 2, 3]))
# リッジ回帰の実行 ★penaltyを追加
alpha = 0.1
penalty = np.array([0, alpha, alpha]) / len(df) # 定数項(第1要素)に0を設定
result = (smf.ols(formula='y ~ 1 + x1 + x2', data=df)
.fit_regularized(alpha=penalty, L1_wt=0)) # リッジの場合、L1_wt=0
# 結果の表示
print('回帰係数 : ', result.params)
print('予測値 train: ', result.fittedvalues)
print('RMSE train: ', root_mean_squared_error(df.y, result.fittedvalues))【実行結果】
書籍の結果とは一致しませんが、statsmodels の結果とは一致しております！
【コードのポイント】
penalty 変数のうち、定数項に対する正則化パラメータを 0 にしました。
penalty = np.array([0, alpha, alpha]) / len(df) # 定数項(第1要素)に0を設定【推察～まとめ】
この statsmodels の設定＆結果からの想像ですが、scikit-learn の一般的な使い方の場合、定数項に対する正則化は行われない、と考えられます。
以上で記事はおしまい、、、
にしようと思っていましたが、上記【推察～まとめ】に関連して、数学の達人から素敵なアドバイスをいただきましたので、ぜひご紹介させてください！
5. 一般的なリッジ回帰
株式会社すうがくぶんか の 内場崇之 先生（@utaka233）より、次のようなメッセージを X 上でいただきました。
先生、ありがとうございます！
書籍で勉強していると、ridge推定量は以下の公式
— 内場 崇之(Takayuki Uchiba) (@utaka233) June 12, 2024
w^ = (X^TX + αI)^{-1}X^Ty
の導出が掲載されています。でも、これはw0も正則化項に加えた場合の式で、実際にはsklearnのRidgeクラスが採用しているように切片w0を正則化項に入れない場合をよく用います。これは、numpyで実装するなら以下の通りです。 https://t.co/KL5qrFnxkq pic.twitter.com/oN9Y7Wy3UX
リッジ回帰の実用上、「切片$${w_0}$$を正則化項に入れない」とのことです。
「4. その２ 一般的な使用法」節の最後の推察は、リッジ回帰の一般的な考え方であることが分かりました。
明快な答えをいただき、誠にありがとうございます！
さっそく先生に教わった numpy によるリッジ回帰の実装を試してみましょう。
### 切片w0に正則化を適用しない
# https://x.com/utaka233/status/1800935766828417142
## 設定
alpha = 0.1 # 正則化パラメータα
## データの作成
# 計画行列Xの作成 ※第1列は定数項
X = np.array([[1, 2, 4],
[1, 3, 6],
[1, 4, 8]])
# yの作成 shape=(3, 1)
y = np.array([[1, 2, 3]]).T
## リッジ回帰の実行
# 行列Iの作成：(0,0)番目の要素に0をセット⇒切片w0に正則化を適用しないようにする
I = np.diag([0, 1, 1])
print('行列I:')
print(I)
# 回帰係数の推定
w = np.linalg.inv(X.T @ X + alpha * I) @ X.T @ y
print('\n回帰係数: ', w.T[0])
# 学習データによる予測
y_pred = (X @ w).T[0]
print('予測値 train: ', y_pred)
# RMSEによる評価
print('RMSE train: ', np.sqrt(sum((y.T[0] - y_pred)**2) / len(y)))【実行結果】
「4. その２ 一般的な使用法」節の scikit-learn や statsmodels の結果と一致しました！
【コードのポイント】
定数項（切片）に正則化を適用しないようにするため、行列$${I}$$（単位行列ではありません）を導入します。
行列$${I}$$の内容は、上の実行結果をご覧ください。
この行列を正則化パラメータ$${\alpha}$$に乗じることで、定数項（切片）に正則化を適用しないように制御します。
回帰係数の推定について、書籍との違いを比較してみましょう。
まずは書籍から。
# 回帰係数の推定 公式1.6
w = np.linalg.inv(X.T @ X + alpha * np.eye(len(X))) @ X.T @ y「np.eye(len(X))」は 次に示す「$${3 \times 3}$$の単位行列」です。
$$
\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{pmatrix}
$$
定数項（切片）に効く$${(0,\ 0)}$$番目の要素は「$${1}$$」になっており、定数項（切片）に正則化を適用しています。
次に内場先生のコードです。
定数項（切片）に効く$${(0,\ 0)}$$番目の要素が「$${0}$$」の 行列$${I}$$になっており、定数項（切片）に正則化を適用しないようにしています。
# 回帰係数の推定
w = np.linalg.inv(X.T @ X + alpha * I) @ X.T @ y
内場先生の回帰係数ベクトルの推定に関する投稿もぜひご参考にしてくださいね。
定数項（切片）$${\boldsymbol{\beta_0}}$$を「０に近づける必要はない」、「本当は$${\boldsymbol{\beta_0}}$$を除く」の言葉が力強いです！
実は、昨日のお話は講座『統計・機械学習のための線形代数と微分積分』でも紹介しています o(^ ^)o 少しでも受講生の皆さんの参考になっていればいいなと思っています！#すうがくぶんか pic.twitter.com/LwT0o5Sjxs
— 内場 崇之(Takayuki Uchiba) (@utaka233) June 13, 2024
先生、ありがとうございました！
まとめ
① 切片の正則化
結局のところ、リッジ回帰の切片に正則化を適用するかどうかが論点でした。
教科書的な書籍の公式・数式では「切片に正則化を適用する」場合があります。
一方で、Pythonの機械学習ライブラリ scikit-learn のリッジ回帰は「切片に正則化を適用しない」です。
リッジ回帰を実践する際には相違点に留意しておくのが良さそうです。
ちなみに、説明変数と目的変数の両方に標準化または中心化を施すと、理論上、切片が$${0}$$になります。
標準化・中心化する場合には、切片の正則化を気にしなくても良さそうです。
もしかすると、教科書的な書籍では、データの標準化・中心化に配慮して、公式・数式を掲載しているのかもしれません。
以下はデータの標準化・中心化と切片の推定値に関するコード例です。
画像ですみません。。。
② statsmodels のリッジ回帰
続いて、Pythonの統計解析ライブラリ statsmodels でリッジ回帰を実装する際の留意点です。
ElasticNet 回帰機能「.fit_regularized()」を代用的に使用するため、次の２点に留意します。
正則化パラメータには、スケーリングした正則化パラメータ値$${\alpha / n}$$を説明変数の個数、設定します。
$${n}$$はデータの個数（サンプルサイズ）です。定数項（切片）に正則化を適用しない場合には、定数項の正則化パラメータに$${0}$$を設定します。
スケーリングした正則化パラメータ penalty 作成のコード例です。
定数項（切片）に正則化を適用しない方法＝ scikit-learn のリッジ回帰と同じ結果を得られると思います。
alpha = 0.1 # リッジ回帰の正則化パラメータの値を設定
n = x.shape[0] # 説明変数xのデータ数（サンプルサイズ）を設定
penalty = np.array([0] + np.tile(alpha, n-1).tolist()) / nそうそう、fit_regularized() で推定した回帰結果は「summary()」で見ることができないので、ご注意下さい。未実装だそうです (^_^;)
③ 参考サイトのご紹介
最後に、定数項を正則化しない理由に言及された記事を見つけましたので、引用させていただきます。
以上で記事はおしまいです（本当です）。
みなさま、ありがとうございました。
それでは本論である「ガウス過程」の学びを頑張ります！
おわり
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？