準備

必要なパッケージを読み込む。

pacman::p_load(tidyverse,
               knitr)

次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。

if (.Platform$OS.type == "windows") { 
  if (require(fontregisterer)) {
    my_font <- "Yu Gothic"
  } else {
    my_font <- "Japan1"
  }
} else if (capabilities("aqua")) {
  my_font <- "HiraginoSans-W3"
} else {
  my_font <- "IPAexGothic"
}

theme_set(theme_gray(base_size = 9,
                     base_family = my_font))


1つの説明変数

母回帰の設定

以下の母回帰関数を仮定する。 \[ \begin{aligned} Y_i &= \beta_0 + \beta_1 X_i + \beta_2 X_i^2 + \varepsilon_i \\ \varepsilon_i & \sim \mbox{Normal}(0, \sigma) \end{aligned} \]

誤差項 \(\varepsilon\) の標準偏差は\(\sigma = 1\) として、\(\boldsymbol{\beta}\) はランダムに決定する。

set.seed(2021-10-11)
sigma <- 1
beta0 <- runif(1, min = -1, max = 1)
beta1 <- runif(1, min = -1, max = 1)
beta2 <- runif(1, min = -1, max = 1)

標本抽出

上の設定で、データ（観測値、訓練データ）を抽出する関数を作る。

sample_data <- function(N = 50, beta0, beta1, beta2, sigma = 1) {
  tibble(x = rnorm(N, mean = 0, sd = 1)) %>% 
    mutate(y = beta0 + beta1 * x + beta2 * x^2 + rnorm(N, mean = 0, sd = sigma))
}

\(N = 50\) で標本抽出を行う。

training <- sample_data(N = 50, 
                        beta0 = beta0, 
                        beta1 = beta1,
                        beta2 = beta2)

回帰モデル

このデータから \(\hat{f}(x)\) を学習するためのモデルために、\(x\)の基底関数を使った以下の4つのモデルを考える。

models_s <- alist(m1 = y ~ x,
                  m2 = y ~ x + I(x^2),
                  m3 = y ~ x + I(x^2) + I(x^3),
                  m4 = y ~ x + I(x^2) + I(x^3) + I(x^4)) %>% 
  map(formula)

それぞれのモデルを使って最小二乗学習を行う。

fits_s <- models_s %>% 
  map(lm, data = training)

予測値 \(\hat{y}\) をデータフレームに加える

trained <- fits_s %>% 
  map(predict) %>% 
  as_tibble() %>% 
  bind_cols(training)

学習結果を使い、\(x \in [-4, 4]\) の範囲で \(\hat{f}(x)\)のグラフを観測値の散布図に重ねて描く。

res_df <- tibble(x = seq(from = -4, to = 4, length.out = 100))
yhats <- fits_s %>% 
  map(predict, newdata = res_df) %>% 
  as_tibble() %>% 
  bind_cols(res_df) %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat")
p <- ggplot() +
  geom_point(data = training, aes(x = x, y = y)) +
  geom_line(data = yhats, aes(x = x, y = y_hat, color = model))
plot(p)

訓練データに対するMSE を計算する。

trained %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat") %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2)) %>% 
  kable(digits = 3)
model MSE
m1 1.559
m2 1.373
m3 1.373
m4 1.364

m4のMSEが最も小さく、データには最もよくフィットしていることがわかる。

念のため、決定係数を確認する。

lapply(1:4, function(x) summary(fits_s[[x]])$r.squared) %>% 
  kable(digits = 3)
x
0.007
x
0.126
x
0.126
x
0.131

やはり、m4 の決定係数が最も大きい。

予測

新たなデータを抽出し、抽出した \(x\)に学習モデルを当てはめたとき、抽出した\(y\)を正しく予測できるか確かめよう。 まず、データを抽出する。

newdf <- sample_data(N = 50,
                     beta0 = beta0,
                     beta1 = beta1,
                     beta2 = beta2)

学習モデルを使って \(\hat{y}\)を計算する。

pred <- fits_s %>% 
  map(predict, newdata = newdf) %>% 
  as_tibble() %>% 
  mutate(y = newdf$y) %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat")

（検証データに対する）予測のMSEを計算する。

pred %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2)) %>% 
  kable(digits = 4)
model MSE
m1 1.1511
m2 0.8699
m3 0.8697
m4 0.9842

m3 のMSEが最も小さい。すなわち、この特定のデータに関する予測は、m3が最も良い。

別の標本を抽出して、もう1度確かめよう。

newdf <- sample_data(N = 50,
                     beta0 = beta0,
                     beta1 = beta1,
                     beta2 = beta2)

学習モデルを使って \(\hat{y}\)を計算する。

pred <- fits_s %>% 
  map(predict, newdata = newdf) %>% 
  as_tibble() %>% 
  mutate(y = newdf$y) %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat")

（検証データに対する）予測のMSEを計算する。

pred %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2)) %>% 
  kable(digits = 4)
model MSE
m1 1.4201
m2 1.2605
m3 1.2614
m4 1.3300

今回は、m2の予測誤差が最も小さい。すなわち、この特定のデータに関する予測は、m2が最も良い。

バイアスと分散のトレードオフ

バイアス\(\mathbb{E}[\hat{f}(x)] - f(x)\) と分散 \(\mathbb{V}[\hat{f}(x)] = \mathbb{E}\left[( \hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2 \right]\) を計算するために、標本抽出と学習を繰り返し、それぞれのモデルについて1,000個の \(\hat{f}\)を手に入れよう。そして、それぞれについて、次の\(x\) (x_test)から\(\hat{f}(x)\)を求める。

x_test <- c(-3, -1.5, 0, 1.5, 3)
f_x <- beta0 + beta1 * x_test + beta2 * x_test^2
names(f_x) <- paste("x =", x_test)

4つのモデルそれぞれについて、1,000個の\(\hat{f}(x)\)をシミュレートする。

n_trials <- 1000
m1_res <- m2_res <- m3_res <- m4_res <-
  matrix(NA, ncol = length(x_test), nrow = n_trials,
         dimnames = list(NULL, names(f_x)))
for (i in 1:n_trials) {
  df <- sample_data(N = 20, beta0, beta1, beta2)
  fits <- models_s %>% 
    map(lm, data = df)
  preds <- fits %>% 
    map(predict, newdata = tibble(x = x_test))
  m1_res[i, ] <- preds[[1]]
  m2_res[i, ] <- preds[[2]]
  m3_res[i, ] <- preds[[3]]
  m4_res[i, ] <- preds[[4]]
}

バイアスを計算する。

bias_m1 <- colMeans(m1_res) - f_x
bias_m2 <- colMeans(m2_res) - f_x
bias_m3 <- colMeans(m3_res) - f_x
bias_m4 <- colMeans(m4_res) - f_x
rbind(bias_m1, bias_m2, bias_m3, bias_m4)
##              x = -3   x = -1.5        x = 0     x = 1.5      x = 3
## bias_m1  2.19747712 0.35192828 -0.252391655  0.38451731 2.26265518
## bias_m2  0.05829404 0.01443367 -0.006526547 -0.00458661 0.02025348
## bias_m3 -0.02720893 0.02115474 -0.004223326  0.02161864 0.22364241
## bias_m4  0.17380162 0.01977264 -0.008354480  0.07458547 0.91192772

ある程度複雑なモデルのほうが、バイアスが小さくなる、すなわちデータへの当てはまりがよくなることがわかる。ただし、観測データの範囲の外側については、m4のバイアスが大きくなることが見て取れる。

次に、分散を計算する。

var_m1 <- apply(m1_res, 2, var)
var_m2 <- apply(m2_res, 2, var)
var_m3 <- apply(m3_res, 2, var)
var_m4 <- apply(m4_res, 2, var)
rbind(var_m1, var_m2, var_m3, var_m4)
##              x = -3  x = -1.5      x = 0   x = 1.5        x = 3
## var_m1    0.9341418 0.2756142 0.05892877 0.2840856    0.9510848
## var_m2    4.8659515 0.3771458 0.08286757 0.3787034    4.8921944
## var_m3   51.6297154 0.8047427 0.09438242 1.1259269   66.0720183
## var_m4 1632.0171083 6.6139021 0.12992508 7.7465269 1724.3047314

複雑なモデルになるほど、分散が大きいことがわかる。つまり、複雑なモデルが学習する\(\hat{f}\)は、データによって大きく変動する。したがって、新しいデータに対する予測はうまくいかない可能性が高い。 観測データの範囲の外側については、その傾向が特に顕著である。


複数の説明変数

母回帰の設定

以下の母回帰関数を仮定する。 \[ \begin{aligned} Y_i &= \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \varepsilon_i \\ \varepsilon_i & \sim \mbox{Normal}(0, \sigma) \end{aligned} \]

誤差項 \(\varepsilon\) の標準偏差は\(\sigma = 1\) として、\(\boldsymbol{\beta}\) を適当に決める

beta0 <- 0.5
beta1 <- 0.8
beta2 <- -0.6

標本抽出

上の設定で、データ（観測値、訓練データ）を抽出する関数を作る。

sample_data2 <- function(N = 50, 
                         beta0 = 0.5, 
                         beta1 = 0.8, 
                         beta2 = -0.6, 
                         sigma = 1) {
  tibble(x1 = rnorm(N, mean = 0, sd = 1),
         x2 = rnorm(N, mean = 0, sd = 1),
         x3 = rnorm(N, mean = 0, sd = 1),
         x4 = rnorm(N, mean = 0, sd = 1),
         x5 = rnorm(N, mean = 0, sd = 1),
         x6 = rnorm(N, mean = 0, sd = 1),
         x7 = rnorm(N, mean = 0, sd = 1),
         x8 = rnorm(N, mean = 0, sd = 1),
         x9 = rnorm(N, mean = 0, sd = 1),
         x10 = rnorm(N, mean = 0, sd = 1)) %>% 
    mutate(y = beta0 + beta1 * x1 + beta2 * x2 + rnorm(N, mean = 0, sd = sigma))
}

\(N = 50\) で標本抽出を行う。

training <- sample_data2(N = 50)

回帰モデル

このデータから \(\hat{f}(x)\) を学習するためのモデルために、\(x\)の基底関数を使った以下の4つのモデルを考える。

models_m <- alist(m1 = y ~ x1,
                  m2 = y ~ x1 + x2,
                  m3 = y ~ x1 + x2 + x3,
                  m4 = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10) %>% 
  map(formula)

それぞれのモデルを使って最小二乗学習を行う。

fits_m <- models_m %>% 
  map(lm, data = training)

予測値 \(\hat{y}\) をデータフレームに加える

trained <- fits_m %>% 
  map(predict) %>% 
  as_tibble() %>% 
  bind_cols(training)

訓練データに対するMSE を計算する。

trained %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat") %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2))
## # A tibble: 4 × 2
##   model   MSE
##   <chr> <dbl>
## 1 m1    1.10 
## 2 m2    0.782
## 3 m3    0.756
## 4 m4    0.654

m4のMSEが最も小さく、データには最もよくフィットしていることがわかる。

念のため、決定係数を確認する。

lapply(1:4, function(x) summary(fits_m[[x]])$r.squared)
## [[1]]
## [1] 0.3277038
## 
## [[2]]
## [1] 0.5208046
## 
## [[3]]
## [1] 0.5366069
## 
## [[4]]
## [1] 0.5992128

やはり、m4 の決定係数が最も大きい。 つまり、yとは無関係のx3 〜 x10を使ったモデルのほうが、データへの当てはまりはよい。

予測

新たなデータを抽出し、抽出した \(x\)に学習モデルを当てはめたとき、抽出した\(y\)を正しく予測できるか確かめよう。 まず、データを抽出する。

newdf <- sample_data2(N = 50)

学習モデルを使って \(\hat{y}\)を計算する。

pred <- fits_m %>% 
  map(predict, newdata = newdf) %>% 
  as_tibble() %>% 
  mutate(y = newdf$y) %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat")

（検証データに対する）予測のMSEを計算する。

pred %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2))
## # A tibble: 4 × 2
##   model   MSE
##   <chr> <dbl>
## 1 m1    0.925
## 2 m2    0.587
## 3 m3    0.618
## 4 m4    0.904

m2 のMSEが最も小さい。すなわち、この特定のデータに関する予測は、m2が最も良い。

別の標本を抽出して、もう1度確かめよう。

newdf <- sample_data2(N = 50)

学習モデルを使って \(\hat{y}\)を計算する。

pred <- fits_m %>% 
  map(predict, newdata = newdf) %>% 
  as_tibble() %>% 
  mutate(y = newdf$y) %>% 
  pivot_longer(cols = m1:m4,
               names_to = "model",
               values_to = "y_hat")

（検証データに対する）予測のMSEを計算する。

pred %>% 
  group_by(model) %>% 
  summarize(MSE = mean((y - y_hat)^2))
## # A tibble: 4 × 2
##   model   MSE
##   <chr> <dbl>
## 1 m1    0.874
## 2 m2    0.751
## 3 m3    0.851
## 4 m4    1.07

今回も、m2の予測誤差が最も小さい。すなわち、この特定のデータに関する予測は、m2が最も良い。

複雑なモデルにすればデータへの当てはまりは良くなるが、それは予測に貢献しないばかりか、かえって予測性能が低くなるようだ。つまり、モデルを複雑（柔軟）にすると、過剰適合が起こる。

実習課題



授業の内容に戻る