「線形代数の半歩先」をPythonで写経 ～ ２章 線形回帰、リッジ回帰、ラッソ回帰

第３部「ならべた数に応用を」

書籍の著者　大久保 潤 先生

この記事は、書籍「線形代数の半歩先」の 第３部「ならべた数に応用を」に掲載の「データサイエンスと機械学習」に関する Python写経活動 のドキュメンタリーです。

第３部は機械学習でおなじみの「回帰」と「次元削減」を線形代数と接続しています。
機械学習に慣れた方は書籍の内容がすっと腹落ちするのではないでしょうか？

この記事では、第15話、第16話の 線形回帰、多項式回帰、リッジ回帰、ラッソ回帰 の可視化に取り組みます！
数学素人なので、どうぞお手柔らかにお願いいたします。

では書籍を開いて線形代数の旅に出発です🚀

はじめに

このブログシリーズは、書籍「線形代数の半歩先 データサイエンス・機械学習に挑む前の30話」（講談社サイエンティフィク、「テキスト」と呼びます）の Python 写経の実践を通じて得た個人的な知見を書きます。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

第３部 ならべた数に応用を

テキスト第13話「世界の一部をモデルに写し取る」から始まる第３部では、現実の世界→データの世界→モデルの世界 の流れを押さえつつ、モデルを関数とベクトルに接続するイメージを理解します。
第15話「データに合う関数を探す」では、コスト関数（損失関数）に二乗誤差を使う最小二乗法を行列とベクトルで書き、線形回帰・多項式回帰を実践します。
つづく第16話「学び過ぎはよくない？」では、過学習対策の一つとして正則化を押さえて、リッジ回帰とラッソ回帰を実践します。

では、この記事のコードで利用するライブラリのインポートからスタートしましょう。

### インポート

# 数値計算・表計算
import numpy as np
import pandas as pd

# 線形回帰
from sklearn.linear_model import Ridge, Lasso         # リッジ回帰, Lasso回帰
from sklearn.preprocessing import PolynomialFeatures  # 多項式特徴量の生成

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

１．線形回帰・多項式回帰

① 線形回帰の概要
テキストより、コスト関数に二乗誤差を用いて回帰係数を推定する最小二乗法の推定式（15.11）をお借りします。

$${\widehat{\bm w}}$$ は パラメータである回帰係数の推定量、$${X}$$ は説明変数の行列、$${\bm y}$$ は目的変数のベクトルです。
$${^\top}$$ の記号は行列の転置、$${(X^{\top} X)^{-1}}$$ の $${^{-1}}$$ はかっこで求めた行列の逆行列です。

この式の詳しい導出過程はぜひテキストをお読みください。
ベクトルで関数を微分する方法と応用を深く学べると思います！

◆ ◆ ◆

② 多項式回帰の概要
説明変数に $${x^2, x^3}$$ などの多項式を適用して、非線形のデータに対応できるモデルを検討します。
テキストより「４次の多項式」に関する数式を引用して、主役級の役者の顔ぶれを紹介します！

説明変数 $${\bm x}$$ に $${1, x, x^2, x^3, x^4}$$ の変換を施す関数 $${\bm \phi(\bm x)}$$ を定義します。
小文字のファイです。

最初の $${1}$$ は定数項（切片）に該当します。
４次（$${M=4}$$）の場合、５個（$${M+1}$$）の変数が現れます。

つづいて多項式変換を一般化した式です。
テキストは行列 $${\Phi}$$ をデータ行列と呼んでいます。
$${\Phi}$$ は大文字のファイです。
列数は多項式 $${M+1}$$ 個、行数はデータ点 $${N}$$ 個です。

データ行列 $${\Phi}$$ を用いた多項式回帰のパラメータ推定量 $${\widehat{\bm w}}$$ は次の数式で表されます。

線形回帰のパラメータ推定量 $${\widehat{\bm w}}$$ の式（15.11）とよく似ています。
$${X}$$ と $${\Phi}$$ が入れ替わっています！

◆ ◆ ◆

③ 線形回帰・多項式回帰の実装
多項式回帰のパラメータ推定の数式を用いて、テキストの図 15.4「曲がっていても「線形」？」を描きます。
なお、データセットの値をなるべくテキストに寄せたつもりですが、完全には一致していないため、テキストの図と異なっていることをご了承ください。

### p.129 図15.4 曲がっていても「線形」？

### データセットの作成
## 設定と準備
# サンプルサイズの設定
n = 20
# 乱数生成器の初期化
rng = np.random.default_rng(seed=56)
## データセットの作成
# xの作成： xn ～ Uniform[-0.4, 1.9]
x = rng.uniform(low=-0.4, high=1.9, size=n)
# yの作成： yn = -sin(4x) + 1.5x + 2 + 𝜖, 𝜖 ～ Normal(0, 0.2²)
eps = rng.normal(loc=0, scale=0.2, size=n)
y = -np.sin(4*x) + 1.5*x + 2 + eps


### 最小二乗法によるパラメータ推定
## 関数の定義
# 特徴量行列 Φ 作成関数の定義 ※式（15.23）
Phi = lambda x, deg: np.array([x**i for i in range(deg + 1)]).T
# 最小二乗法による線形回帰のパラメータwの推定関数の定義  ※式（15.28）
lsm = lambda Phi, y: np.linalg.inv(Phi.T @ Phi) @ Phi.T @ y
# 線形回帰の予測値算出関数の定義 ※式（15.24）
predict = lambda phi, w: phi @ w

## 予測値の算出
# x軸の値の作成
x_val = np.linspace(-1, 2, 101)
# 線形回帰（次数1の多項式）の予測値算出
degree = 1
f1_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))
# 次数4の多項式回帰の予測値算出
degree = 4
f2_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))


### 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)

## [1] 線形回帰（次数1の多項式）の最尤推定
# データセットの散布図の描画
ax1.plot(x, y, 'o', color='royalblue', alpha=0.7)
# 多項式の曲線の描画
ax1.plot(x_val, f1_pred, color='tab:red')
# 修飾
ax1.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', title='[1] $f(x)=w_0 + w_1x$')

## [2] 次数4の多項式の最尤推定
# データセットの散布図の描画 
ax2.plot(x, y, 'o', color='royalblue', alpha=0.7)
# 多項式の曲線の描画
ax2.plot(x_val, f2_pred, color='tab:green')
# 修飾
ax2.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', 
        title=f'[2] $f(x)=w_0 + w_1x + w_2 x^2 + w_3 x^3 + w_4 x^4$')
fig.suptitle('図 15.4 曲がっていても「線形」？', fontsize=14)

plt.show()

【実行結果】
左側の単回帰では非線形のデータ点をうまく捉えきれていません。
一方で右側の４次多項式回帰では非線形のデータ点をいい感じに捉えています！

■ コードの補足説明
コードの「最小二乗法によるパラメータ推定」の箇所で $${\widehat{w}}$$ を推定しています。
関数は lambda 関数で書いてみました。

• データ行列（特徴量行列）$${\Phi}$$ を作成する関数定義

  • 式（15.23）に基づく式（15.27）の実装です。

  • 変数 $${x}$$ に対して、$${1, x, x^2,  \ldots, x^{\text{deg}}}$$ の値を作成します。

# 特徴量行列 Φ 作成関数の定義 ※式（15.23）
Phi = lambda x, deg: np.array([x**i for i in range(deg + 1)]).T

• 最小二乗法でパラメータ推定する関数定義

  • 式（15.28）の実装です。

# 最小二乗法による線形回帰のパラメータwの推定関数の定義  ※式（15.28）
lsm = lambda Phi, y: np.linalg.inv(Phi.T @ Phi) @ Phi.T @ y

• 線形回帰（多項式回帰）で予測値を算出する関数定義

  • 式（15.24）の実装です。

# 線形回帰の予測値算出関数の定義 ※式（15.24）
predict = lambda phi, w: phi @ w

• 単回帰（１次の多項式回帰）による予測値の算出

# 線形回帰（次数1の多項式）の予測値算出
degree = 1
f1_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))

• ４次の多項式回帰による予測値の算出

# 次数4の多項式回帰の予測値算出
degree = 4
f2_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))

２．リッジ回帰・ラッソ回帰

① 次数を大きくすると・・・
テキストの図 16.1「柔軟なモデルがよい、とは限らない⋯」を実装します。
データセットがテキストと異なるため、次数 $${12}$$ で実装しています。

### p.133 図16.1 柔軟なモデルがよい、とは限らない⋯ ★テキストは13次の多項式

### データセットの作成
## 設定と準備
# サンプルサイズの設定
n = 20
# 乱数生成器の初期化
rng = np.random.default_rng(seed=56)  # 0, 6, 31, 33, 35, 43, 46, 50, 52
## データセットの作成
# xの作成： xn ～ Uniform[-0.4, 1.9]
x = rng.uniform(low=-0.4, high=1.9, size=n)
# yの作成： yn = -sin(4x) + 1.5x + 2 + 𝜖, 𝜖 ～ Normal(0, 0.2²)
eps = rng.normal(loc=0, scale=0.2, size=n)
y = -np.sin(4*x) + 1.5*x + 2 + eps


### 最小二乗法によるパラメータ推定
## 関数の定義
# 特徴量行列 Φ 作成関数の定義 ※式（15.23）
Phi = lambda x, deg: np.array([x**i for i in range(deg + 1)]).T
# 最小二乗法による線形回帰のパラメータwの推定関数の定義  ※式（15.28）
lsm = lambda Phi, y: np.linalg.inv(Phi.T @ Phi) @ Phi.T @ y
# 線形回帰の予測値算出関数の定義 ※式（15.24）
predict = lambda phi, w: phi @ w

## 予測値の算出
# x軸の値の作成
x_val = np.linspace(-1, 2, 301)
# 次数4の多項式の予測値算出
degree = 4
f1_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))
# 次数12の多項式回帰の予測値算出
degree = 12
f2_pred = predict(Phi(x_val, degree), lsm(Phi(x, degree), y))


### 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)

## [1] 線形回帰（次数1の多項式）の最尤推定
# データセットの散布図の描画
ax1.plot(x, y, 'o', color='royalblue', alpha=0.7)
# 多項式の曲線の描画
ax1.plot(x_val, f1_pred, color='tab:red')
# 修飾
ax1.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', title='[1] 4次までの多項式')

## [2] 次数4の多項式の最尤推定
# データセットの散布図の描画 
ax2.plot(x, y, 'o', color='royalblue', alpha=0.7)
# 多項式の曲線の描画
ax2.plot(x_val, f2_pred, color='tab:green')
# 修飾
ax2.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', 
        title=f'[2] 12次までの多項式')
fig.suptitle('図 16.1 柔軟なモデルがよい、とは限らない⋯', fontsize=14)

plt.show()

【実行結果】

右側の $${12}$$ 次の多項式回帰の結果は・・・

• 回帰の曲線（緑色）がデータ点（青い点）にピタッと張り付いています。

• ところどころ、不自然な線になっています。
  たとえば、右から３点目と２点目の間の上昇、右端の極端な下降、などです。

テキストでは右側のチャートを眺めつつ、とても大切なことを記載しています。引用します。

■ 様々な多項式回帰を動かして確かめよう！
「乱数生成器の初期化」のコード中の「seed=」の数値を変えるとデータセットが変化します。
「予測値」のコード中の「degree=」の数値（次数）を変えるとチャートの曲線が変化します。
ぜひ、コードを変えてデータと次数の関係や過学習を体感してください！

◆ ◆ ◆

② リッジ回帰・ラッソ回帰の概要
過学習＝学び過ぎを防ぐ方法としてテキストはリッジ回帰とラッソ回帰を紹介しています。
数式や過学習を防げる理由などはぜひテキストで確かめてください！

ざっくり二つの回帰は「正則化」という罰則項を適用しています。
正則化の適用により、リッジ回帰では回帰係数（パラメータ）が過度に大きくならないように、ラッソ回帰では特定の説明変数の回帰係数（パラメータ）をゼロにするようになります。
正則化の強さはパラメータ $${\alpha}$$ でコントロールします。
なお$${\alpha}$$は、パラメータ推定の対象にならないので、人為的に値を決めてやる必要があります。

③ リッジ回帰・ラッソ回帰の実装
テキストの図 16.4「リッジ回帰とラッソ回帰の結果」を実装します。
テキストと同様に、scikit-learn ライブラリを使います。
多項式のデータ行列・特徴量行列（$${\Phi}$$）の作成も scikit-learn ライブラリ を使います。PolynomialFeatures です。

（参考）
パラメータ推定の数式に関しては、リッジ回帰はテキストの式（16.8）を用いて計算可能ですが、ラッソ回帰は正則化に「絶対値」（$${L1}$$正則化）を用いていて微分できないので、反復計算（数値計算）のテクニックが必要とされています。

では実装しましょう。
チャートの見栄え上、次数 $${14}$$ で作成しています（テキストは $${13}$$）。

### p.143 図16.4 リッジ回帰とラッソ回帰 ★テキストは次数13

### データセットの作成
## 設定と準備
# サンプルサイズの設定
n = 20
# 乱数生成器の初期化
rng = np.random.default_rng(seed=56)  # 0, 6, 31, 33, 35, 43, 46, 50, 52
## 学習データの作成
# xの作成： xn ～ Uniform[-0.4, 1.9]
x = rng.uniform(low=-0.4, high=1.9, size=n)
# yの作成： yn = -sin(4x) + 1.5x + 2 + 𝜖, 𝜖 ～ Normal(0, 0.2²)
eps = rng.normal(loc=0, scale=0.2, size=n)
y = -np.sin(4*x) + 1.5*x + 2 + eps

### 多項式特徴量の作成
# 次数14の多項式特徴量（定数項を含む）を生成するためのインスタンスの作成
poly = PolynomialFeatures(degree=14, include_bias=True)
# 学習データを次数14の多項式特徴量に変換
x_poly = poly.fit_transform(x.reshape(-1, 1))
# 予測に用いるデータを作成して、次数14の多項式特徴量に変換
x_val = np.linspace(-1, 2, 301)
x_val_poly = poly.transform(x_val.reshape(-1, 1))


### リッジ回帰とラッソ回帰の実行
## リッジ回帰の実行
ridge = Ridge(alpha=0.01, fit_intercept=False)
ridge.fit(x_poly, y)
f1_pred = ridge.predict(x_val_poly)
## ラッソ回帰の実行
lasso = Lasso(alpha=0.01, fit_intercept=False, tol=1e-2)
lasso.fit(x_poly, y)
f2_pred = lasso.predict(x_val_poly)


### 描画
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)

## [1] リッジ回帰
# データセットの散布図の描画
ax1.plot(x, y, 'o', color='royalblue', alpha=0.7)
# リッジ回帰の予測値の曲線の描画
ax1.plot(x_val, f1_pred, color='tab:green')
# 修飾
ax1.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', title='[1] リッジ回帰')

## [2] ラッソ回帰
# データセットの散布図の描画 
ax2.plot(x, y, 'o', color='royalblue', alpha=0.7)
# ラッソ回帰の予測値の曲線の描画
ax2.plot(x_val, f2_pred, color='tab:red')
# 修飾
ax2.set(xlim=(-0.5, 2), ylim=(-1, 6), yticks=range(0, 7, 2),
        xlabel='$x$', ylabel='$y$', 
        title=f'[2] ラッソ回帰')
fig.suptitle('図 16.2 リッジ回帰とラッソ回帰の結果（次数14）', fontsize=14)

plt.show()

【実行結果】

大きな次数を使っていますが、$${12}$$ 次の多項式回帰と比べて、不自然な動きは少ない感じです。

推定した回帰係数を見てみましょう。
列名の数値は次数です。０は定数項（切片）です。
いずれも小数点第３位で丸めています。

■ リッジ回帰

# リッジ回帰の係数の表示
print('リッジ回帰の係数:')
display(pd.DataFrame(ridge.coef_, columns=['w']).T.round(3))

【実行結果】

■ ラッソ回帰

# ラッソ回帰の係数の表示
print('ラッソ回帰の係数:')
display(pd.DataFrame(lasso.coef_, columns=['w']).T.round(3))

【実行結果】

（参考）$${12}$$ 次の多項式回帰

$${12}$$ 次の多項式回帰と比べて、リッジ回帰もラッソ回帰も回帰係数の絶対値が小さくなっています。
また、ラッソ回帰は $${9}$$ 次以降、回帰係数がゼロになっています。

■ コードの補足説明
リッジ回帰、ラッソ回帰はコードの「リッジ回帰とラッソ回帰の実行」の箇所で実行しています。
３ステップでモデルの学習と予測ができます。
ライブラリの威力ですね！

リッジ回帰を例にしますと・・・
・リッジ回帰の設定（インスタンスの生成）

ridge = Ridge(alpha=0.01, fit_intercept=False)

「ridge」という名のモデルを作っています。
引数 alpha で正則化パラメータ $${\alpha}$$ を設定します。
テキストと同じ $${0.01}$$ を使いました。

・リッジ回帰の学習

ridge.fit(x_poly, y)

引数には学習データの説明変数（ここでは x_poly）と目的変数（ここでは y）を与えます。

・学習済みモデルで予測

f1_pred = ridge.predict(x_val_poly)

モデル名「ridge」に predict をすることで予測ができます。
引数には 予測値を求めたいデータの「説明変数」（ここでは x_val_poly）を与えます。

おまけ：リッジ回帰の深堀り

以前、テキストの式（16.8）と同じ数式でリッジ回帰のパラメータ推定を実践した記事を書いています。

この記事では「数式に基づくリッジ回帰」と「Python の scikit-learn ライブラリと statsmodels ライブラリを用いたリッジ回帰」の留意点をまとめています。
留意点は「切片（定数項）に正則化を適用するかどうか」です。
よかったらお読みください！

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。