必要なパッケージを読み込む。
::p_load(tidyverse,
pacman broom)
次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。LinuxにはIPAexフォント がインストールされていることを想定している(IPAex はインストールすれば maxOSやWindows でも使える)。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
<- "Yu Gothic"
my_font else {
} <- "Japan1"
my_font
}else if (capabilities("aqua")) {
} <- "HiraginoSans-W3"
my_font else {
} <- "IPAexGothic"
my_font
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))
単回帰で最小二乗法を実行する。
まず、データ \(\{(x_n, y_n)\}_{n=1}^N\) をランダムに生成する。
set.seed(2021-10-08) # この資料と異なるデータを得るためにはseedの値を変える
<- rnorm(1, mean = 0, sd = 2) # 切片
beta0 <- rnorm(1, mean = 0, sd = 1) # 傾き
beta1 <- 50 # サンプルサイズ
N <- rnorm(N, mean = 0, sd = 1) # 説明変数の観測値
x <- beta0 + beta1 * x # 線形関数 f(x)
mu <- rnorm(N, mean = mu, sd = 1) # ノイズを加えて応答変数の観測値を得る
y <- tibble(x = x, y = y) # データを tibble に myd
得られたデータの散布図を作る。
<- ggplot(myd, aes(x = x, y = y)) +
scat1 geom_point()
plot(scat1)
最小二乗法。
<- with(myd, sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2))
b1 <- with(myd, mean(y) - b1 * mean(x))
b0 cbind(b0, b1) # データから学習した係数
## b0 b1
## [1,] 2.835935 -0.9265309
得られた直線を散布図に描き足す。
<- scat1 +
scat1a 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)
<- scat1a +
scat1b geom_line(data = myd,
aes(y = y_hat),
color = "tomato")
plot(scat1b)
予測値をつないだ線(新たに加えた赤い直線)が、さきほど求めた回帰直線(青い線)に重なっていることがわかる(当たり前だが)。
この学習における「損失」すなわち残差平方和を確認しよう。
<- with(myd, y - y_hat)
e sum(e^2)
## [1] 43.01426
二乗せずに残差の和を求めると、
sum(e)
## [1] -9.992007e-15
となり、(実質的に)ゼロになっていることがわかる。
残差の和をゼロにする直線は、他にも(無数に)ある。 たとえば、y_hat_b
を次のように作ってみる。
<- 2
d1 <- with(myd, mean(y) - d1 * mean(x))
d0 <- myd %>%
myd mutate(y_hat_b = d0 + d1 * x)
<- scat1a +
scat1c geom_line(data = myd,
aes(y = y_hat_b),
color = "orange")
plot(scat1c)
このように、y_hat2
はデータにうまくフィットしていない。しかし、その残差の和は、
<- with(myd, y_hat_b - y)
e2 sum(e2)
## [1] 1.154632e-14
となり、ほぼゼロである。
それに対し、残差平方和を求めると、
sum(e2^2)
## [1] 385.4724
となって、上で求めた回帰直線の残差平方和 43.0142596 よりもかなり大きいことがわかる。
計量経済学で学習した(する)とおり、実際に回帰分析を行う場合には lm()
(または、estimatr::lm_robust()
)を使う。
<- lm(y ~ x, data = myd) fit1
結果は、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
として結果に含まれている。
$fitted.values fit1
## 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
である。
$residuals)^2 %>%
(fit1sum()
## [1] 43.01426
実習課題
上の単回帰における基底関数ベクトルは、 \[ \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}\) を返す関数を作る。
<- function(data, J = 1) {
ls_basis <- matrix(NA, ncol = J, nrow = nrow(data))
X for (j in 1:J) {
<- data$x^j
X[, j]
}<- lm(myd$y ~ X)
fit $fitted.values
fit }
さらに、これを使って図を描くところまで関数にする。
<- function(data, J) {
plot_ls <- ls_basis(data = data, J = J)
y_hat <- tibble(x = data$x,
df y = data$y,
y_hat = y_hat)
<- ggplot(df, aes(x = x)) +
p 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)
実習課題
まず、データ\(\{(x_n, y_n)\}_{n=1}^N\)をランダムに生成する関数を作る。
<- function(N, D, seed = 1) {
gen_data # N: サンプルサイズ
# D: 説明変数の数
# seed 乱数の種
set.seed(seed)
<- matrix(rnorm(D + 1, mean = 0, sd = 1),
beta ncol = 1)
<- matrix(rnorm(N * D, mean = 0, sd = 1),
X nrow = N)
<- cbind(1, X) # 定数項のための第1列を加える
X <- X %*% beta
y <- cbind(y, X)
df colnames(df) <- c("y", paste0("x", 0:D))
as_tibble(df)
}
N=50
、D=5
でデータを作る。
<- gen_data(N = 50, D = 5, seed = 123)
D5 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…
最小二乗法でパラメタを求める。
<- D5 %>%
y pull(y) %>%
as.matrix(ncol = 1)
<- D5[, -1] %>%
X as.matrix()
<- solve(t(X) %*% X) %*% t(X) %*% y
beta_hat 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()
で求める。
<- lm(y ~ x1 + x2 + x3 + x4 + x5,
fit_D5 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