からっぽのしょこ

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

【R】3.5:線形回帰の例【緑ベイズ入門のノート】

はじめに

 『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

 この記事は3.5節の内容です。ベイズ推論を用いて線形回帰モデルのパラメータの学習と未観測データの予測をRで実装します。

 省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。

【数式読解編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

Rでやってみよう

 人工的に生成したデータを用いて、パラメータを推定してみましょう。

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

# 利用パッケージ
library(tidyverse)
library(mvtnorm)

 mvtnormパッケージは、多次元ガウス分布に関するパッケージです。多次元ガウス分布に従う乱数生成関数rmvnorm()を使います。

・観測モデルの設定

 この例では、入力値$x_n$に対して$\mathbf{x}_n = \{1, x_n, x_n^2, x_n^3\}$としたベクトルを用います。このとき次元数$M = 4$であり、観測モデル$p(\mathbf{y} | \mathbf{X} ,\mathbf{w})$のパラメータを$\mathbf{w} = (w_1, w_2, \cdots, w_4)$とします。また平均パラメータ$0$、分散パラメータ$\sigma^2 = \lambda^{-1}$の1次元のガウス分布に従うノイズ成分$\epsilon_n$を加えた

$$ \begin{align} y_n &= \mathbf{w}^{\top} \mathbf{x}_n + \epsilon_n \\ &= w_1 + w_2 x_n + w_3 x_n^2 + w_4 x_n^3 + \epsilon_n \tag{3.141} \end{align} $$

を出力値$y_n$とします。ここで$\lambda$を精度パラメータと呼びます。

 まずは、$x_n$から$\mathbf{x}_n$を作成する関数を定義しておきます。

# (1, x, ..., x^(n-1))のベクトル作成関数を定義
x_vector <- function(x_n, M) {
  x_mn <- matrix(NA, nrow = M, ncol = length(x_n))
  for(m in 1:M) {
    x_mn[m, ] <- x_n^(m-1)
  }
  x_mn
}

# 確認
x_vector(2, 5)
##      [,1]
## [1,]    1
## [2,]    2
## [3,]    4
## [4,]    8
## [5,]   16

 $\mathbf{x}_n$は縦ベクトルであるため、$M$行1列のマトリクスを返します。また$x^0 = 1$であることを利用しています。

 次に、観測モデルのパラメータを設定します。

# 観測モデルのパラメータを指定
M_truth <- 4
sigma <- 1
lambda <- 1 / sigma^2
w_m <- sample(x = seq(-1, 1, by = 0.1), size = M_truth) %>% 
  matrix()

# 確認
w_m
##      [,1]
## [1,]  0.3
## [2,] -0.8
## [3,]  0.5
## [4,]  0.2

 w_mも$M$行1列のマトリクスで縦ベクトルを表現します。

 観測モデルを設定できたので、グラフを作成してモデルを確認しましょう。

# 作図用のx軸の値
x_line <- seq(-3, 3, by = 0.01)
y_line <- t(w_m) %*% x_vector(x_line, M_truth) %>% 
  as.vector()

# ノイズを含まない観測モデルを作図
model_df <- tibble(
  x = x_line, 
  y = y_line
)
ggplot(model_df, aes(x, y)) + 
  geom_line() + 
  labs(title = "Observation Model", 
       subtitle = paste0("w=(", paste0(round(w_m, 1), collapse = ', '), ")"))

f:id:anemptyarchive:20201113174131p:plain
観測モデル


 設定した観測モデルに従って乱数(入力データ)を生成します。

# データをサンプリング
N <- 50
smp_x_n <- sample(seq(min(x_line), max(x_line), by = 0.01), size = N, replace = TRUE)
y_1n <- t(w_m) %*% x_vector(smp_x_n, M_truth) + rnorm(n = N, mean = 0, sd = sqrt(1 / lambda))

# 確認
dim(y_1n)
## [1]  1 50

 データの生成に対して分布の仮定をおいていないため、sample()で一様な確率に従いランダムに値を生成します。ただし$x_n$がとり得る範囲は、グラフの描画範囲にしておきます。サンプリングしたデータsmp_x_nは、計算に用いないためsmple()の返り値のままのベクトル形式です。

 観測モデル(3.141)に従い、出力値$\mathbf{y} = \{y_1, \cdots, y_N\}$を計算します。計算結果は、1行$N$列のマトリクス(横ベクトル)になります。

 観測データの散布図を先ほどの観測モデルのグラフと重ねて確認しておきます。

# 観測データの散布図を作成
sample_df <- tibble(
  x = smp_x_n, 
  y = as.vector(y_1n)
)
ggplot() + 
  geom_point(data = sample_df, aes(x, y)) + # 観測データ
  geom_line(data = model_df, aes(x, y), color = "blue") + # ノイズなし観測モデル
  labs(title = "Observation Model", 
       subtitle = paste0("M=", M_truth, ", N=", N, ", sigma=", round(sqrt(1 / lambda), 1)))

f:id:anemptyarchive:20201113174159p:plain
観測データ

 以上が観測モデルに関する設定です。

 続いて、事前分布$p(\mathbf{w})$のパラメータを設定します。

# 事前分布のパラメータを指定
M <- 4
m_m <- matrix(rep(0, M), nrow = M, ncol = 1)
sigma_mm <- diag(M)
lambda_mm <- solve(sigma_mm^2)

# xのベクトルを作成
x_mn <- x_vector(smp_x_n, M)
x_mline <- x_vector(x_line, M)

# 確認
dim(x_mn)
## [1]  4 50

 事前分布は多次元ガウス分布を仮定するのでした。次元数をMに指定します。この値はパラメータ$\mathbf{w}$の要素数に対応します。

 この値に合うように、平均パラメータ$\mathbf{m} = (m_1, \cdots, m_M)$、精度行列パラメータ$\boldsymbol{\Lambda}^{-1} = \boldsymbol{\Sigma} = (\sigma_{11}^2, \cdots, \sigma_{MM}^2)$を指定する必要があります。
 この例では、平均m_mを全ての要素が0の$M$行1列のマトリクス(縦ベクトル)、標準偏差と相関係数の行列$(\sigma_{11}, \cdots, \sigma_{MM})$を$M$行$M$列の単位行列としてその2乗の逆行列を精度行列lambda_mmとします。$M$行$M$列の単位行列はdiag(M)で作成し、逆行列はsolve()で計算します。この場合lambda_mmも$M$行$M$列の単位行列になるので、直接指定しても問題ありません。

 またサンプリングしたデータ$\{x_1, \cdots, x_N\}$をもとに、入力値$\mathbf{X} = \{\mathbf{x}_1, \cdots, \mathbf{x}_N\}$に変換します。

 設定した事前分布から$\mathbf{w}$をサンプリングしたモデルを確認しましょう。

# 事前分布からのwのサンプリングによるモデルを比較
smp_model_df <- tibble()
for(i in 1:5) {
  tmp_df <- tibble(
    x = x_line, 
    y = as.vector(
      mvtnorm::rmvnorm(n = 1, mean = m_m, sigma = solve(lambda_mm)) %*% x_mline
    ), 
    smp_num = as.factor(i)
  )
  smp_model_df <- rbind(smp_model_df, tmp_df)
}
ggplot() + 
  geom_line(data = smp_model_df, aes(x, y, color = smp_num)) + # 事前分布からサンプリングwによるモデル
  geom_line(data = model_df, aes(x, y), color = "blue", linetype = "dotted") + # 観測モデル
  ylim(min(min(y_1n), min(y_line)), max(max(y_1n), max(y_line))) + 
  geom_line() + 
  labs(title = "Sampling from Prior Distribution")

f:id:anemptyarchive:20201114233523p:plain
事前分布からサンプリングしたパラメータによるモデル

 この時点では当然観測モデルとは程遠いグラフになります。

・事後分布の計算

 観測データy_nx_mnを用いて、事後分布$p(\mathbf{w} | \mathbf{y}, \mathbf{X})$(のパラメータ)を求めます。

# 事後分布のパラメータを計算
lambda_hat_mm <- lambda * x_mn %*% t(x_mn) + lambda_mm
m_hat_m <- solve(lambda_hat_mm) %*% (lambda * x_mn %*% t(y_1n) + lambda_mm %*% m_m)

 事後分布(多次元ガウス分布)の精度行列パラメータ$\hat{\boldsymbol{\Lambda}}^{-1} = \hat{\boldsymbol{\Sigma}} = (\hat{\sigma}_{11}^2, \cdots, \hat{\sigma}_{MM}^2)$、平均パラメータ$\hat{\mathbf{m}} = (\hat{m}_1, \cdots, \hat{m}_M)$を、次の式で計算します。

$$ \begin{align} \hat{\boldsymbol{\Lambda}} &= \lambda \sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top} + \boldsymbol{\Lambda} \\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}^{-1} \left( \lambda \sum_{n=1}^N y_n \mathbf{x}_n + \boldsymbol{\Lambda} \mathbf{m} \right) \tag{3.148} \end{align} $$

 ただし$\sum_{n=1}^N y_n \mathbf{x}_n$の計算を効率よく処理するために、$\mathbf{X} \mathbf{y}^{\top}$で計算しています。計算結果をそれぞれlambda_hat_mmm_hat_mとします。

 事後分布からも$\mathbf{w}$をサンプリングしてモデルを確認しましょう。

# 事後分布による観測モデルのパラメータのサンプリングによるモデルの変化を比較
smp_model_df <- tibble()
for(i in 1:5) {
  tmp_df <- tibble(
    x = x_line, 
    y = as.vector(
      mvtnorm::rmvnorm(n = 1, mean = m_hat_m, sigma = solve(lambda_hat_mm)) %*% x_mline
    ), 
    smp_num = as.factor(i)
  )
  smp_model_df <- rbind(smp_model_df, tmp_df)
}
ggplot() + 
  geom_line(data = smp_model_df, aes(x, y, color = smp_num)) + # 事後分布からサンプリングしたwによるモデル
  geom_line(data = model_df, aes(x, y), color = "blue", linetype = "dotted") + # 観測モデル
  geom_point(data = sample_df, aes(x, y)) + # 観測データ
  ylim(min(min(y_1n), min(y_line)), max(max(y_1n), max(y_line))) + 
  labs(title = "Sampling from Posterior Distribution", 
       subtitle = paste0("M=", M))

f:id:anemptyarchive:20201113174259p:plain
事後分布からサンプリングしたパラメータによるモデル

 観測データから学習(パラメータを更新)することで、どれも観測モデルに近い形になっていることが分かります。

・予測分布の計算

 最後に、事後分布パラメータm_mlambda_hat_mmを用いて、予測分布$(y_{*} | \mathbf{x}_{*}, \mathbf{y}, \mathbf{X})$(のパラメータ)を求めます。

# 予測分布のパラメータを計算
mu_star_hat_line <- t(m_hat_m) %*% x_mline %>% 
  as.vector()
sigma2_star_hat_line <- rep(NA, length(x_line))
for(i in seq_along(x_line)) {
  sigma2_star_hat_line[i] <- 1 / lambda + t(x_mline[, i]) %*% solve(lambda_hat_mm) %*% matrix(x_mline[, i])
}

 予測分布(1次元のガウス分布)の平均パラメータ$\hat{\mu}_{*}$、精度パラメータ$\hat{\lambda}_{*}^{-1} = \hat{\sigma}^2$を、次の式で計算します。

$$ \begin{aligned} \hat{\mu}_{*} &= \hat{\mathbf{m}}^{\top} \mathbf{x}_{*} \\ \hat{\lambda}_{*}^{-1} &= \lambda^{-1} + \mathbf{x}_{*}^{\top} \hat{\boldsymbol{\Lambda}}^{-1} \mathbf{x}_{*} \end{aligned} $$

 グラフの描画範囲全体の未知の入力データ$\mathbf{X}_{*}$に対して計算する必要があるため、x_mlineを使います。計算結果をそれぞれsigma2_star_hat_linemu_star_hat_lineとします。

 作図用に、平均mu_star_hat_lineと平均から標準偏差sqrt(sigma2_star_hat_line)を足し引きしたものをまとめたデータフレームを作成します。

# 予測分布を計算
predict_df <- tibble(
  x = x_line, 
  E_y = mu_star_hat_line, # 予測分布の期待値
  minus_sigma_y = E_y - sqrt(sigma2_star_hat_line), # 予測分布の期待値-sigma
  plus_sigma_y = E_y + sqrt(sigma2_star_hat_line) # 予測分布の期待値+sigma
)
head(predict_df)
## # A tibble: 6 x 4
##       x   E_y minus_sigma_y plus_sigma_y
##   <dbl> <dbl>         <dbl>        <dbl>
## 1 -3     2.16         0.949         3.37
## 2 -2.99  2.17         0.968         3.38
## 3 -2.98  2.19         0.986         3.39
## 4 -2.97  2.20         1.00          3.40
## 5 -2.96  2.21         1.02          3.41
## 6 -2.95  2.23         1.04          3.41


 予測分布を作図します。

# 作図
ggplot() + 
  geom_line(data = predict_df, aes(x, E_y), color = "orange") + # 予測分布の期待値
  geom_line(data = predict_df, aes(x, plus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値+sigma
  geom_line(data = predict_df, aes(x, minus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値-sigma
  geom_line(data = model_df, aes(x, y), color = "blue", linetype = "dotted") + # 観測モデル
  geom_point(data = sample_df, aes(x, y)) + # 観測データ
  ylim(min(min(y_1n), min(y_line)), max(max(y_1n), max(y_line))) + 
  labs(title = "Predictive Distribution", 
       subtitle = paste0("M=", M, ", N=", N), 
       y = "y")

f:id:anemptyarchive:20201113174336p:plain
予測分布

 以上が推論処理です。

 いくつか値を変えて試してみましょう。

 データ数$N$を減らしてパラメータ数$M$を増やした場合は、次のようになります。

f:id:anemptyarchive:20201113174402p:plain
予測分布:過学習

 モデルへの当てはまりが悪くなり、またx軸の各値に対するy軸の値の標準偏差(緑の破線の縦方向の範囲)も広くなっているのが分かります。

 次に描画範囲を広くして、観測データにない範囲の入力値に対する予測分布を確認してみましょう。

f:id:anemptyarchive:20201113174427p:plain
予測分布:過学習

 $M$の値が真の値よりも大きいと、観測データの範囲を越えたとたんに予測できていません。

f:id:anemptyarchive:20201113174443p:plain
予測分布

 $M$が真の値と一致していると、(ある程度)予測ができています。

・おまけ

 gganimateパッケージを利用して、gif画像を作成するコードです。

 観測データの増加が予測精度へ与える影響を確認します。

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

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

# 観測モデルのパラメータを指定
M_truth <- 4
sigma <- 1
lambda <- 1 / sigma^2
w_m <- sample(x = seq(-1, 1, by = 0.1), size = M_truth) %>% 
  matrix()

# 作図用のx軸の値
x_line <- seq(-3, 3, by = 0.01)
y_line <- t(w_m) %*% x_vector(x_line, M_truth) %>% 
  as.vector()

# 観測モデルのデータフレームを作成
model_df <- tibble(
  x = x_line, 
  y = y_line
)
# サンプルサイズを指定
N <- 100

# 事前分布のパラメータを指定
M <- 10
m_m <- matrix(rep(0, M), nrow = M, ncol = 1)
sigma_mm <- diag(M)
lambda_mm <- solve(sigma_mm^2)

# xのベクトルを作成
x_mline <- x_vector(x_line, M)

# 推論
smp_x_n <- rep(NA, N)
y_1n <- rep(NA, N)
eps_n <- rep(NA, N)
sample_df <- tibble()
predict_df <- tibble()
for(n in 1:N) {
  
  # データをサンプリング
  smp_x_n[n] <- sample(seq(min(x_line), max(x_line), by = 0.01), size = 1, replace = TRUE)
  y_1n[n] <- t(w_m) %*% x_vector(smp_x_n[n], M_truth) + rnorm(n = 1, mean = 0, sd = sqrt(1 / lambda))
  
  # 観測データのデータフレームを作成
  sample_df <- tibble(
    x = smp_x_n[1:n], 
    y = as.vector(y_1n[1:n]), 
    iteration = n
  ) %>% 
    rbind(sample_df, .)
  
  # 観測データからxベクトルを作成
  x_mn <- x_vector(smp_x_n[n], M)
  
  # 事後分布のパラメータを計算
  old_lambda_mm <- lambda_mm
  lambda_mm <- lambda * x_mn %*% t(x_mn) + lambda_mm
  m_m <- solve(lambda_mm) %*% (lambda * x_mn %*% t(y_1n[n]) + old_lambda_mm %*% m_m)
  
  # 予測分布のパラメータを計算
  sigma2_star_line <- rep(NA, length(x_line))
  for(i in seq_along(x_line)) {
    sigma2_star_line[i] <- 1 / lambda + t(x_mline[, i]) %*% solve(lambda_mm) %*% matrix(x_mline[, i])
  }
  mu_star_line <- t(m_m) %*% x_mline %>% 
    as.vector()
  
  # 予測分布を計算
  predict_df <- tibble(
    x = x_line, 
    E_y = mu_star_line, 
    minus_sigma_y = E_y - sqrt(sigma2_star_line), 
    plus_sigma_y = E_y + sqrt(sigma2_star_line), 
    iteration = n
  ) %>% 
    rbind(predict_df, .)
}

# 作図
predict_graph <- ggplot() + 
  geom_line(data = predict_df, aes(x, E_y), color = "orange") + # 予測分布の期待値
  geom_line(data = predict_df, aes(x, plus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値+sigam
  geom_line(data = predict_df, aes(x, minus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値-sigma
  geom_line(data = model_df, aes(x, y), color = "blue", linetype = "dotted") + # 観測モデル
  geom_point(data = sample_df, aes(x, y)) + # 観測データ
  transition_manual(iteration) + # フレーム
  ylim(min(min(y_1n), min(y_line)), max(max(y_1n), max(y_line))) + 
  labs(title = "Predictive Distribution", 
       subtitle = paste("M=", M, ", N={current_frame}", sep = ""), 
       y = "y")

# gif画像を作成
animate(predict_graph, nframes = N, fps = 10)

 異なる点のみを簡単に解説します。

 各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはfor()で1データずつ処理します。よって観測データ数Nがイタレーション数になります。

 一度の処理で事後分布のパラメータを計算するのではなく、事前分布のパラメータに対して繰り返し観測データの情報を与えることで更新(上書き)していきます。$\sum_{n=1}^N$の計算は、N回繰り返し計算することで実行されます。
 ただし事後分布のパラメータの計算(3.148)の計算において、更新前(事前分布)のパラメータ$\boldsymbol{\Lambda}$を使うため、old_lambda_mmとして値を一時的に保存しておきます。

f:id:anemptyarchive:20201113181041g:plain
予測分布の推移:過学習

 (わざわざ学習してるっぽい振る舞いをさせていますが、これが良くないことは先ほど確認した通りです。)

 同じ観測データを用いたとき、パラメータの数(事前分布の次元数)の違いが予測精度へ与える影響を確認します。

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

# サンプルサイズを指定
N <- 10

# データをサンプリング
smp_x_n <- sample(seq(min(x_line), max(x_line), by = 0.01), size = N, replace = TRUE)
y_1n <- t(w_m) %*% x_vector(smp_x_n, M_truth) + rnorm(n = N, mean = 0, sd = sqrt(1 / lambda)) 

# 観測データのデータフレームを作成
sample_df <- tibble(
  x = smp_x_n, 
  y = as.vector(y_1n)
)

# 事前分布のパラメータを指定
max_M <- 15

# 推論
predict_df <- tibble()
for(m in 1:max_M) {
  
  # 事前分布のパラメータを生成
  m_m <- matrix(rep(0, m), nrow = m, ncol = 1)
  sigma_mm <- diag(m)
  lambda_mm <- solve(sigma_mm^2)
  
  # 観測データからxベクトルを作成
  x_mn <- x_vector(smp_x_n, m)
  x_mline <- x_vector(x_line, m)
  
  # 事後分布のパラメータを計算
  lambda_hat_mm <- lambda * x_mn %*% t(x_mn) + lambda_mm
  m_hat_m <- solve(lambda_hat_mm) %*% (lambda * x_mn %*% t(y_1n) + lambda_mm %*% m_m)
  
  # 予測分布のパラメータを計算
  mu_star_hat_line <- t(m_hat_m) %*% x_mline %>% 
     as.vector()
  sigma2_star_hat_line <- rep(NA, length(x_line))
  for(i in seq_along(x_line)) {
     sigma2_star_hat_line[i] <- 1 / lambda + t(x_mline[, i]) %*% solve(lambda_hat_mm) %*% matrix(x_mline[, i])
   }
  
  # 予測分布を計算
  predict_df <- tibble(
    x = x_line, 
    E_y = mu_star_hat_line, 
    minus_sigma_y = E_y - sqrt(sigma2_star_hat_line), 
    plus_sigma_y = E_y + sqrt(sigma2_star_hat_line), 
    iteration = m
  ) %>% 
    rbind(predict_df, .)
}

# 作図
predict_graph <- ggplot() + 
  geom_line(data = predict_df, aes(x, E_y), color = "orange") + # 予測分布の期待値
  geom_line(data = predict_df, aes(x, plus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値+sigam
  geom_line(data = predict_df, aes(x, minus_sigma_y), color = "#00A968", linetype = "dashed") + # 予測分布の期待値-sigma
  geom_line(data = model_df, aes(x, y), color = "blue", linetype = "dotted") + # 観測モデル
  geom_point(data = sample_df, aes(x, y)) + # 観測データ
  transition_manual(iteration) + # フレーム
  ylim(min(min(y_1n), min(y_line)), max(max(y_1n), max(y_line))) + 
  labs(title = "Predictive Distribution", 
       subtitle = paste0("M={current_frame}", ", N=", N), 
       y = "y")

# gif画像を作成
animate(predict_graph, nframes = max_M, fps = 5)

 こちらはMをイタレーション数として、推論処理をM回繰り返します。

f:id:anemptyarchive:20201113181246g:plain
予測分布の変化:過学習


参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

 実装すると数式の間違いや理解が足りていなかったところに気付けてとても良いです。まぁ理解を深める工程の分先へは進めませんがね。
 そういえばこれまでのモデルでは、予測分布のパラメータの計算に入力$\mathbf{x}_{*}$自体を使わなかったので気にもならなかったのですが、1つのデータに対して予測するので$\sum_{n=1}^N$の計算が不要みたいな流れでしたが、じゃあ複数のデータ(x_mline)の場合はどうなるんだ?結果の分布(グラフ)がちゃんとそれっぽくなってるんでこのままでいいんだろうなとして考えると、事後分布では観測した全てのデータの影響を考慮するため$N$個のデータとして扱う必要があるが、予測分布においては未知のデータ間の影響を考慮しようがないため独立に計算すればいいのだろうと理解しました。いかがでしょうか?
 勉強をするほど疑問が出てくるのは、(うまく回っているうちは)良いことなんだと思ってます。

 そんなただでさえ思ってるより進捗が出せない状況なのですが、この内容をPythonでも実装してみることにしました。

 最後に昨日公開されたモーニング娘。'20さんの新曲です!ぜひBGMに♪

 エビエビエビデンス🦐

【次節の内容】

www.anarchive-beta.com