準備

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

pacman::p_load(tidyverse,
               broom)

次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。LinuxにはIPAexフォント がインストールされていることを想定している（IPAex はインストールすれば maxOSやWindows でも使える）。

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))


最小二乗法

単回帰

単回帰で最小二乗法を実行する。

まず、データ \(\{(x_n, y_n)\}_{n=1}^N\) をランダムに生成する。

set.seed(2021-10-08) # この資料と異なるデータを得るためにはseedの値を変える
beta0 <- rnorm(1, mean = 0, sd = 2)  # 切片
beta1 <- rnorm(1, mean = 0, sd = 1)  # 傾き
N <- 50  # サンプルサイズ
x <- rnorm(N, mean = 0, sd = 1)  # 説明変数の観測値
mu <- beta0 + beta1 * x          # 線形関数 f(x)
y <- rnorm(N, mean = mu, sd = 1) # ノイズを加えて応答変数の観測値を得る 
myd <- tibble(x = x, y = y)  # データを tibble に

得られたデータの散布図を作る。

scat1 <- ggplot(myd, aes(x = x, y = y)) +
  geom_point()
plot(scat1)

最小二乗法。

b1 <- with(myd, sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2))
b0 <- with(myd, mean(y) - b1 * mean(x))
cbind(b0, b1)  # データから学習した係数
##            b0         b1
## [1,] 2.835935 -0.9265309

得られた直線を散布図に描き足す。

scat1a <- scat1 +
  geom_abline(intercept = b0, 
              slope = b1, 
              color = "dodgerblue")
plot(scat1a)

この直線は、予測値 \(\hat{y} = b_0 + b_1 x\) の集合のはずだ。 それを確認しよう。

myd <- myd %>% 
  mutate(y_hat = b0 + b1 * x)

scat1b <- scat1a +
  geom_line(data = myd, 
            aes(y = y_hat), 
            color = "tomato")
plot(scat1b)

予測値をつないだ線（新たに加えた赤い直線）が、さきほど求めた回帰直線（青い線）に重なっていることがわかる（当たり前だが）。

この学習における「損失」すなわち残差平方和を確認しよう。

e <- with(myd, y - y_hat)
sum(e^2)
## [1] 43.01426

二乗せずに残差の和を求めると、

sum(e)
## [1] -9.992007e-15

となり、（実質的に）ゼロになっていることがわかる。

残差の和をゼロにする直線は、他にも（無数に）ある。 たとえば、y_hat_b を次のように作ってみる。

d1 <- 2
d0 <- with(myd, mean(y) - d1 * mean(x))
myd <- myd %>% 
  mutate(y_hat_b = d0 + d1 * x)
scat1c <- scat1a +
  geom_line(data = myd,
            aes(y = y_hat_b),
            color = "orange")
plot(scat1c)

このように、y_hat2 はデータにうまくフィットしていない。しかし、その残差の和は、

e2 <- with(myd, y_hat_b - y)
sum(e2)
## [1] 1.154632e-14

となり、ほぼゼロである。

それに対し、残差平方和を求めると、

sum(e2^2)
## [1] 385.4724

となって、上で求めた回帰直線の残差平方和 43.0142596 よりもかなり大きいことがわかる。

計量経済学で学習した（する）とおり、実際に回帰分析を行う場合には lm()（または、estimatr::lm_robust()）を使う。

fit1 <- lm(y ~ x, data = myd)

結果は、broom::tidy() を使うと見やすい。

tidy(fit1) %>% 
  select(term, estimate)
## # A tibble: 2 × 2
##   term        estimate
##   <chr>          <dbl>
## 1 (Intercept)    2.84 
## 2 x             -0.927

y_hatは、fitted.values として結果に含まれている。

fit1$fitted.values
##         1         2         3         4         5         6         7         8 
## 3.0830425 1.8412688 2.1660640 4.0433629 2.1299525 1.6168763 2.6861920 2.3373218 
##         9        10        11        12        13        14        15        16 
## 3.8783028 3.6927392 3.1203250 2.0230169 2.1703641 3.8693058 4.1877619 2.5059036 
##        17        18        19        20        21        22        23        24 
## 2.6588048 3.1444056 2.5771941 3.6079043 0.5923387 3.4755991 3.2670257 2.8450502 
##        25        26        27        28        29        30        31        32 
## 2.6665862 1.9714502 2.5899234 3.8245133 2.5329407 1.9041414 2.0691229 2.1630957 
##        33        34        35        36        37        38        39        40 
## 1.8314990 4.0860341 2.0653530 1.4540032 3.1686931 2.6899617 2.8190824 3.9344093 
##        41        42        43        44        45        46        47        48 
## 4.0346036 2.2647769 2.0727798 2.7629147 1.7519134 2.8017257 1.8069491 4.1439777 
##        49        50 
## 2.1176556 2.1113411

同様に、残差は residuals である。

(fit1$residuals)^2 %>% 
  sum()
## [1] 43.01426

実習課題

  1. 残差の和をゼロにする直線は他にもあることを確かめなさい。
  2. 残和の和をゼロにする直線をどのように選んでも、最小二乗法による回帰直線よりも残差平方和を小さくできる直線がみつけられなそうなことを確かめなさい（回帰直線の傾きを少しずつ変えると、残差平方和が大きくなってしまうことを確認しなさい）。

基底関数の数を増やす

上の単回帰における基底関数ベクトルは、 \[ \phi(x) = (x^0, x^1) = (1, x) \] で、この基底関数とパラメタベクトル\(\beta = (\beta_0, \beta_1)^\mathsf{T}\)を利用して、 \[ y = \phi(x) \beta + \varepsilon \] という仮定のもとで最小二乗学習を行った。

では、基底関数ベクトルを拡張し、 \[ \phi(x) = (x^0, x^1, \dots, x^J) \] の\(J\)を大きくすると、何が起こるだろうか。確かめてみよう。

\(J\) を指定すると、 \[ y_i = \beta_0 + \beta_1 x_i + \cdots + \beta_J x_i^J + \varepsilon_i \] を仮定して最小二乗学習を行い、予測値 \(\hat{y}\) を返す関数を作る。

ls_basis <- function(data, J = 1) {
  X <- matrix(NA, ncol = J, nrow = nrow(data))
  for (j in 1:J) {
    X[, j] <- data$x^j
  }
  fit <- lm(myd$y ~ X)
  fit$fitted.values
}

さらに、これを使って図を描くところまで関数にする。

plot_ls <- function(data, J) {
  y_hat <- ls_basis(data = data, J = J)
  df <- tibble(x = data$x,
               y = data$y,
               y_hat = y_hat)
  p <- ggplot(df, aes(x = x)) +
    geom_point(aes(y = y)) +
    geom_line(aes(y = y_hat), color = "dodgerblue") +
    labs(title = paste("J =", J))
  plot(p)
}

\(J = 2\)の場合

plot_ls(data = myd, J = 2)

\(J = 5\)の場合

plot_ls(data = myd, J = 5)

\(J = 10\)の場合

plot_ls(data = myd, J = 10)

実習課題

  • さまざまな\(J\)の値について、上の関数を実行し、基底関数ベクトルの選び方によって最小二乗学習の結果がどのように変化するか調べなさい。
    • \(J\) はどこまで大きくできるだろうか？
    • \(J\) をあ無限に大きくできないはなぜだろうか？
    • サンプルサイズ\(N\)を変えると、\(J\)と損失の関係はどう変化するだろうか？
  • 上で作った関数を、\(J\) の値ごとに残差平方和が計算（確認）できるように改良し、\(J\)の大きくすると残差平方和がどのように変化するか調べなさい。

複数の説明変数を利用する

まず、データ\(\{(x_n, y_n)\}_{n=1}^N\)をランダムに生成する関数を作る。

gen_data <- function(N, D, seed = 1) {
  # N: サンプルサイズ
  # D: 説明変数の数
  # seed 乱数の種
  set.seed(seed)
  beta <- matrix(rnorm(D + 1, mean = 0, sd = 1),
                 ncol = 1)
  X <- matrix(rnorm(N * D, mean = 0, sd = 1),
              nrow = N)
  X <- cbind(1, X)  # 定数項のための第1列を加える
　y <- X %*% beta
  df <- cbind(y, X)
  colnames(df) <- c("y", paste0("x", 0:D))
  as_tibble(df)
}

N=50、D=5でデータを作る。

D5 <- gen_data(N = 50, D = 5, seed = 123)
glimpse(D5)
## Rows: 50
## Columns: 7
## $ y  <dbl> -4.41568188, -0.54360764, 2.72159017, -0.19759786, 0.04956472, -1.1…
## $ x0 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ x1 <dbl> 0.4609162, -1.2650612, -0.6868529, -0.4456620, 1.2240818, 0.3598138…
## $ x2 <dbl> -1.548752804, 0.584613750, 0.123854244, 0.215941569, 0.379639483, -…
## $ x3 <dbl> -0.78490447, -1.66794194, -0.38022652, 0.91899661, -0.57534696, 0.6…
## $ x4 <dbl> 0.56298953, -0.37243876, 0.97697339, -0.37458086, 1.05271147, -1.04…
## $ x5 <dbl> -0.78860284, -0.59461727, 1.65090747, -0.05402813, 0.11924524, 0.24…

最小二乗法でパラメタを求める。

y <- D5 %>% 
  pull(y) %>% 
  as.matrix(ncol = 1)
X <- D5[, -1] %>% 
  as.matrix()

beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
row.names(beta_hat) <- paste0("b", 0:5)
beta_hat
##           [,1]
## b0 -0.56047565
## b1 -0.23017749
## b2  1.55870831
## b3  0.07050839
## b4  0.12928774
## b5  1.71506499

lm() で求める。

fit_D5 <- lm(y ~ x1 + x2 + x3 + x4 + x5, 
             data = D5)
tidy(fit_D5) %>% 
  select(term, estimate)
## # A tibble: 6 × 2
##   term        estimate
##   <chr>          <dbl>
## 1 (Intercept)  -0.560 
## 2 x1           -0.230 
## 3 x2            1.56  
## 4 x3            0.0705
## 5 x4            0.129 
## 6 x5            1.72


授業の内容に戻る