からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】3.1.4.a:リッジ回帰の実装【PRMLのノート】

はじめに

 『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでのスクラッチ実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。

 この記事は、3.1.4項「正則化最小二乗法」の内容です。リッジ回帰(Ridge回帰)の最尤推定をR言語で実装します。

【数理読解編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

・リッジ回帰の実装

 リッジ回帰(L2正則化)を実装します。線形基底関数モデルについては「【R】3.1.1:線形基底関数モデル【PRMLのノート】 - からっぽのしょこ」を参照してください。
 この記事では、リッジ回帰により観測データから真のモデルを推定(近似)します。

 利用するパッケージを読み込みます。

# 3.1.4項で利用するパッケージ
library(tidyverse)


・モデルの設定とデータの生成

 データを生成する関数を設定して、利用するデータを人工的に生成します。処理の流れは3.1.1項と同じなので、詳しくはそちらを参照してください。

 データ生成関数(推定対象の関数)と基底関数を作成します。

# 真の関数を作成
y_true <- function(x) {
  # 計算式を指定
  y <- sin(2 * pi * x)
  return(y)
}

# 多項式基底関数を作成
phi_poly <- function(x, m) {
  # m乗を計算
  y <- x^m
  return(y)
}

# 計画行列を作成
Phi <- function(x_n, M) {
  # 変数を初期化
  x_nm <- matrix(0, nrow = length(x_n), ncol = M)
  
  # m列目を計算
  for(m in 0:(M-1)) {
    x_nm[, m+1] <- phi_poly(x_n, m)
  }
  return(x_nm)
}

 この例では、多項式基底関数を使います。ガウス基底関数やシグモイド基底関数を使う場合は3.1.1項を参照してください。

 真のモデルの精度パラメータを指定して出力を計算します。

# データの範囲を指定
x_min <- 0
x_max <- 1

# 作図用のxの値を作成
x_vals <- seq(x_min, x_max, by = 0.01)

# 真の精度パラメータを指定
beta_true <- 3
sigma_true <- sqrt(1 / beta_true) # 真の標準偏差

# 真のモデルを計算
model_df <- tidyr::tibble(
  x = x_vals, # x軸の値
  y = y_true(x_vals), # y軸の値
  minus_sigma = y - 2 * sigma_true, # μ - 2σ
  plus_sigma = y + 2 * sigma_true # μ + 2σ
)


 設定したモデルに従って、データを生成します。

# データ数を指定
N <- 100

# (観測)データを生成
x_n <- runif(n = N, min = x_min, max = x_max) # 入力
t_n <- y_true(x_n) + rnorm(n = N, mean = 0, sd = sigma_true) # 出力

# 観測データをデータフレームに格納
data_df <- tidyr::tibble(
  x_n = x_n, # 入力
  t_n = t_n # 出力
)

 $N$個の観測データ(入力値$\mathbf{x} = \{x_1, \cdots, x_N\}$と目標値$\mathbf{t} = \{t_1, \cdots, t_N\}$)を作成できました。

・最尤推定

 正則化しない場合と正則化する場合の最尤推定を行いパラメータ$\mathbf{w} = (w_0, w_1, \cdots, w_{M-1})^{\top}$を求めて、結果を比較します。

 パラメータの数(数基底関数の数)$M$を指定して、入力値を基底関数により変換します。

# 重みパラメータの次元数(基底関数の数)を指定
M <- 6

# 基底関数により入力を変換
phi_x_nm <- Phi(x_n, M)

 計画行列

$$ \boldsymbol{\Phi} = \begin{pmatrix} \phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_{M-1}(x_1) \\ \phi_0(x_2) & \phi_1(x_2) & \cdots & \phi_{M-1}(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(x_N) & \phi_1(x_N) & \cdots & \phi_{M-1}(x_N) \end{pmatrix} \tag{3.16} $$

phi_x_nmとします。また、$\boldsymbol{\phi}(x_n)$は$\boldsymbol{\Phi}$の$n$行目です。

・正則化しない場合

 まずは、正則化を行わない場合の重みパラメータと分散パラメータの最尤解を計算します。

# 重みパラメータの最尤解を計算
w_ml_m <- solve(t(phi_x_nm) %*% phi_x_nm) %*% t(phi_x_nm) %*% t_n %>% 
  as.vector()

# 分散パラメータの最尤解を計算
sigma2_ml <- sum((t_n - phi_x_nm %*% w_ml_m)^2) / N

# 精度パラメータの最尤解を計算
beta_ml <- 1 / sigma2_ml

# 確認
w_ml_m
beta_ml
## [1]    0.07435094    2.49476598   23.49882292 -113.08672718  137.94215416
## [6]  -50.62912633
## [1] 3.153926

 重みパラメータ、分散パラメータ、精度パラメータをそれぞれ次の式で計算します。

$$ \begin{align} \mathbf{w}_{\mathrm{ML}} &= (\boldsymbol{\Phi}^{\top} \boldsymbol{\Phi})^{-1} \boldsymbol{\Phi}^{\top} \mathsf{t} \tag{3.15}\\ \sigma^2_{\mathrm{ML}} &= \frac{1}{N} \sum_{n=1}^N \Bigl( t_n - \mathbf{w}_{\mathrm{ML}}^{\top} \phi(\mathbf{x}_n) \Bigr)^2 \tag{3.21}\\ \beta_{\mathrm{ML}} &= \frac{1}{\sigma^2_{\mathrm{ML}}} \end{align} $$

 ただし、$t_n - \mathbf{w}^{\top} \phi(\mathbf{x}_n)$の計算を$N$個同時に行うため、$\mathbf{t} - \boldsymbol{\Phi} \mathbf{w}$で計算しています。

 推定した線形回帰のパラメータを使って回帰曲線を計算します。

# 推定したパラメータによるモデルを計算
ml_df <- tidyr::tibble(
  x = x_vals, # x軸の値
  t = Phi(x_vals, M) %*% w_ml_m %>% 
    as.vector(), # y軸の値
  minus_sigma = t - 2 * sqrt(sigma2_ml), # μ - 2σ
  plus_sigma = t + 2 * sqrt(sigma2_ml) # μ + 2σ
)

# 確認
head(ml_df)
## # A tibble: 6 x 4
##       x      t minus_sigma plus_sigma
##   <dbl>  <dbl>       <dbl>      <dbl>
## 1  0    0.0744      -1.05        1.20
## 2  0.01 0.102       -1.02        1.23
## 3  0.02 0.133       -0.993       1.26
## 4  0.03 0.167       -0.959       1.29
## 5  0.04 0.205       -0.921       1.33
## 6  0.05 0.245       -0.882       1.37

 未知の入力値(x軸の値)$x_{*}$に対する予測値(y軸の値)$t_{*}$を次の式で計算します。

$$ \begin{aligned} \mathbf{t} &= \boldsymbol{\Phi} \mathbf{w}_{\mathrm{ML}} \\ t_{*} &= \mathbf{w}_{\mathrm{ML}}^{\top} \boldsymbol{\phi}(x_{*}) \end{aligned} $$

 さらに、分散パラメータsigma2_mlを使って、推定した精度を求めます。

 ここまでは、3.1.1項と同様の処理です。

・正則化する場合

 続いて、正則化を行う場合の重みパラメータと分散パラメータの最尤解を計算します。

# 正則化係数を指定
lambda <- 0.5

# リッジ回帰の重みパラメータの最尤解を計算
w_ridge_m <- solve(lambda * diag(M) + t(phi_x_nm) %*% phi_x_nm) %*% t(phi_x_nm) %*% t_n %>% 
  as.vector()

# リッジ回帰の分散パラメータの最尤解を計算
sigma2_ridge <- sum((t_n - phi_x_nm %*% w_ridge_m)^2) / N

# リッジ回帰の精度パラメータの最尤解を計算
beta_ridge <- 1 / sigma2_ridge

# 確認
w_ridge_m
beta_ridge
## [1]  0.7343945 -0.5998903 -1.7423368 -0.7757288  0.4985178  1.5943355
## [1] 2.235324

 重みパラメータ、分散パラメータ、精度パラメータをそれぞれ次の式で計算します。

$$ \begin{align} \mathbf{w}_{\mathrm{Ridge}} &= \Bigl( \lambda \mathbf{I} + \boldsymbol{\Phi}^{\top} \boldsymbol{\Phi} \Bigr)^{-1} \boldsymbol{\Phi}^{\top} \mathsf{t} \tag{3.28}\\ \sigma^2_{\mathrm{Ridge}} &= \frac{1}{N} \sum_{n=1}^N \Bigl( t_n - \mathbf{w}_{\mathrm{Ridge}}^{\top} \phi(\mathbf{x}_n) \Bigr)^2 \tag{3.21}\\ \beta_{\mathrm{Ridge}} &= \frac{1}{\sigma^2_{\mathrm{Ridge}}} \end{align} $$

 こちらも、$N$個同時に計算(処理)します。
 diag()で単位行列を作成できます。

 推定したリッジ回帰のパラメータを使って回帰曲線を計算します。

# 推定したパラメータによるモデルを計算
ridge_df <- tidyr::tibble(
  x = x_vals, # x軸の値
  t = Phi(x_vals, M) %*% w_ridge_m %>% 
    as.vector(), # y軸の値
  minus_sigma = t - 2 * sqrt(sigma2_ridge), # μ - 2σ
  plus_sigma = t + 2 * sqrt(sigma2_ridge) # μ + 2σ
)

# 確認
head(ridge_df)
## # A tibble: 6 x 4
##       x     t minus_sigma plus_sigma
##   <dbl> <dbl>       <dbl>      <dbl>
## 1  0    0.734      -0.603       2.07
## 2  0.01 0.728      -0.609       2.07
## 3  0.02 0.722      -0.616       2.06
## 4  0.03 0.715      -0.623       2.05
## 5  0.04 0.708      -0.630       2.05
## 6  0.05 0.700      -0.638       2.04

 先ほどと同様に、y軸の値を計算します。

$$ \begin{aligned} \mathbf{t} &= \boldsymbol{\Phi} \mathbf{w}_{\mathrm{Ridge}} \\ t_{*} &= \mathbf{w}_{\mathrm{Ridge}}^{\top} \boldsymbol{\phi}(x_{*}) \end{aligned} $$

 以上がリッジ回帰の推定処理です。推定したパラメータの各要素の値が小さくなっているのを確認できます。

・推定結果の確認

 リッジ回帰により推定した回帰曲線(近似曲線)を作図します。

# 推定したパラメータによるモデルを作図
ggplot() + 
  geom_point(data = data_df, aes(x = x_n, y = t_n)) + # 観測データ
  geom_line(data = model_df, aes(x = x, y = y), color = "cyan4") + # 真のモデル
  geom_ribbon(data = model_df, aes(x = x, ymin = minus_sigma, ymax = plus_sigma), 
              fill = "cyan4", alpha = 0.1, color = "cyan4", linetype = "dotted") + # 真のノイズ範囲
  geom_line(data = ml_df, aes(x = x, y = t), color = "blue", linetype ="dashed") + # 推定したモデル:(正則化なし)
  geom_line(data = ridge_df, aes(x = x, y = t), color = "blue") + # 推定したモデル:(L2正則化)
  geom_ribbon(data = ridge_df, aes(x = x, ymin = minus_sigma, ymax = plus_sigma), 
              fill = "blue", alpha = 0.1, color = "blue", linetype = "dotted") + # 推定したノイズ範囲
  #ylim(c(-5, 5)) + # y軸の表示範囲
  labs(title = "Ridge Regression", 
       subtitle = paste0("N=", N, ", M=", M, ", q=", 2, ", lambda=", lambda, 
                         ", w=(", paste0(round(w_ridge_m, 2), collapse = ", "), ")", 
                         ", beta=", round(beta_ridge, 2)), 
       x = "x", y = "t")

f:id:anemptyarchive:20211221172237p:plain
リッジ回帰

 緑の実線が真のモデル、青の実線がリッジ回帰(L2正則化)による近似曲線、青の破線が正則化なしの近似曲線です。
 この例だと、正則化しない方がうまく近似できていますね。
 $\lambda$の値が大きいほど、正則化の効果が強まります。

 正則化の効果を見るために、観測データ数を減らしパラメータ数を増やして、過学習が起きやすい設定にしてみます。

f:id:anemptyarchive:20211221172252p:plain
過学習とリッジ回帰

 正則化していない方(青の破線)はデータに過剰に適合しているのが分かります。正則化している方(青の実線)は過剰適合を回避できています。

 以上でリッジ回帰を実装できました。

・おまけ:正則化係数と回帰曲線の関係

 最後に、正則化係数と回帰曲線(近似曲線)の関係をアニメーションで確認します。

・作図コード(クリックで展開)

 アニメーション(gif画像)の作成にgganimateパッケージを利用します。

# 追加パッケージ
library(gganimate)


 for()ループで正則化係数$\lambda$の値を変更しながら繰り返しパラメータと回帰曲線を計算します。

# 使用するlambdaの値を作成
lambda_vec <- seq(0, 0.05, by = 0.001)

# lambdaごとに最尤解を計算
anime_ridge_df <- dplyr::tibble()
for(lambda in lambda_vec) {
  # 重みパラメータの最尤解を計算
  w_ridge_m <- solve(lambda * diag(M) + t(phi_x_nm) %*% phi_x_nm) %*% t(phi_x_nm) %*% t_n %>% 
    as.vector()
  
  # 分散パラメータの最尤解を計算
  sigma2_ridge <- sum((t_n - t(w_ridge_m) %*% t(phi_x_nm))^2) / N
  
  # 精度パラメータの最尤解を計算
  beta_ridge <- 1 / sigma2_ridge
  
  # 推定したパラメータによるモデルを計算
  tmp_ridge_df <- tidyr::tibble(
    x = x_vals, # x軸の値
    t = Phi(x_vals, M) %*% w_ridge_m %>% 
      as.vector(), # y軸の値
    minus_sigma = t - 2 * sqrt(sigma2_ridge), # μ - 2σ
    plus_sigma = t + 2 * sqrt(sigma2_ridge), # μ + 2σ
    label = paste0(
      "lambda=", lambda, 
      ", w=(", paste0(round(w_ridge_m, 2), collapse = ", "), ")", 
      ", beta=", round(beta_ridge, 2)
    ) %>% 
      as.factor() # フレーム切替用のラベル
  )
  
  # 結果を結合
  anime_ridge_df <- rbind(anime_ridge_df, tmp_ridge_df)
}

# 確認
head(anime_ridge_df)
## # A tibble: 6 x 5
##       x       t minus_sigma plus_sigma label                                    
##   <dbl>   <dbl>       <dbl>      <dbl> <fct>                                    
## 1  0    -12.2       -13.1      -11.2   lambda=0, w=(-12.18, 589.6, -9177.01, 67~
## 2  0.01  -7.14       -8.10      -6.18  lambda=0, w=(-12.18, 589.6, -9177.01, 67~
## 3  0.02  -3.56       -4.52      -2.60  lambda=0, w=(-12.18, 589.6, -9177.01, 67~
## 4  0.03  -1.13       -2.09      -0.172 lambda=0, w=(-12.18, 589.6, -9177.01, 67~
## 5  0.04   0.411      -0.549      1.37  lambda=0, w=(-12.18, 589.6, -9177.01, 67~
## 6  0.05   1.29        0.329      2.25  lambda=0, w=(-12.18, 589.6, -9177.01, 67~

 $\lambda$として使う値を作成してlambda_valsとします。
 lambda_valsから値を順番に取り出して、パラメータの最尤解を求めて回帰曲線を計算します。計算結果はanime_ridge_dfに追加していきます。

 作図してgif画像として出力します。

# 推定したパラメータによるモデルを作図
anime_graph <- ggplot() + 
  geom_point(data = data_df, aes(x = x_n, y = t_n)) + # 観測データ
  geom_line(data = model_df, aes(x = x, y = y), color = "cyan4") + # 真のモデル
  geom_ribbon(data = model_df, aes(x = x, ymin = minus_sigma, ymax = plus_sigma), 
              fill = "cyan4", alpha = 0.1, color = "cyan4", linetype = "dotted") + # 真のノイズ範囲
  geom_line(data = ml_df, aes(x = x, y = t), color = "blue", linetype ="dashed") + # 推定したモデル:(正則化なし)
  geom_line(data = anime_ridge_df, aes(x = x, y = t), color = "blue") + # 推定したモデル:(L2正則化)
  geom_ribbon(data = anime_ridge_df, aes(x = x, ymin = minus_sigma, ymax = plus_sigma), 
              fill = "blue", alpha = 0.1, color = "blue", linetype = "dotted") + # 推定したノイズ範囲
  gganimate::transition_manual(label) + # フレーム
  ylim(c(-4, 4)) + # y軸の表示範囲
  labs(title = "Ridge Regression", 
       subtitle = "{current_frame}", 
       x = "x", y = "t")

# gif画像に変換
gganimate::animate(anime_graph, nframes = length(lambda_vec), fps = 10)


f:id:anemptyarchive:20211221172355g:plain
リッジ回帰の正則化の様子

 正則化係数が大きくなるほど、パラメータの値が小さくなっていき、回帰曲線がシンプルになっていくことから、正則化の効果が強まるのを確認できます。
 $\lambda = 0$のとき正則化しない最尤解と一致しているのが分かります。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.

おわりに

 こっちは問題なくできた問題はラッソ回帰...

 先ほど配信されたモーニング娘。'21の新曲MVをどうぞ!

 してほしーい🔫

  • 2021.12.21:追記して記事を分割しました。


【次節の内容】

www.anarchive-beta.com

【関連する内容】

www.anarchive-beta.com