情報科学3
制約付き最小二乗学習
矢内 勇生
2021-10-14
準備
必要なパッケージを読み込む。
pacman::p_load(tidyverse,
broom,
rgl,
glmnet,
ISLR,
rsample)次に、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))制約付き最小二乗学習とは？
データの準備
母数を適当に設定し、その母集団からサイズ\(N\)の標本を抽出する。
set.seed(2021-10-14)
N <- 100
beta0 <- 0.8
beta1 <- 2.5
beta2 <- 1.4
x <- runif(N, min = -2, max = 2)
z <- -0.5 * x + rnorm(N, mean = 0, sd = 1) # xとzに相関をもたせる
mu <- beta0 + beta1 * x + beta2 * z
y <- rnorm(N, mean = mu, sd = 2)
sim_d <- tibble(y = y, x = x, z = z)抽出した標本の各変数を中心化する。
sim_d <- sim_d %>%
mutate(y_c = y - mean(y),
x_c = x - mean(x),
z_c = z - mean(z))制約なしの最小二乗学習
中心化した変数を用いて最小二乗学習を行う。切片なしを指定するために、formula に-1を加える
fit1 <- lm(y_c ~ -1 + x_c + z_c,
data = sim_d)
tidy(fit1) %>%
select(term:estimate)## # A tibble: 2 × 2
## term estimate
## <chr> <dbl>
## 1 x_c 2.46
## 2 z_c 1.70ここで、さまざまな \(\hat{\boldsymbol{\beta}} = (\hat{\beta}_1, \hat{\beta}_2)\) について、損失（残差平方和） \[ S(\boldsymbol{\beta}) = \sum_{n=1}^N \left( y_n - (\beta_1 x_n + \beta_2 z_n) \right)^2 \] を計算してみよう。
LOSS <- function(beta1, beta2) {
sum((sim_d$y_c - (beta1 * sim_d$x_c + beta2 * sim_d$z_c))^2)
}
b1 <- seq(-1.5, 5, length.out = 500)
b2 <- seq(-1.5, 5, length.out = 500)
LOSS_D <- expand_grid(beta1 = b1, beta2 = b2)
LOSS_D$loss <- LOSS_D %>%
pmap(.f = LOSS) %>%
unlist()
loss_matrix <- LOSS_D %>%
pull(loss) %>%
matrix(ncol = length(b1),
byrow = TRUE)計算結果を可視化する。3-D プロットはウェブページに表示できないので、各自で確認されたい。
persp3d(x = b1, y = b2, z = loss_matrix,
xlab = "beta1", ylab = "beta2", zlab = "Loss",
color = "dodgerblue")最小二乗法で推定された点で、損失（残差平方和）が最小化されていることが確認できる。
同じ結果を等高線図 (contour plot) にする。
cont1 <- ggplot(LOSS_D) +
geom_contour_filled(aes(x = beta1, y = beta2, z = loss)) +
geom_vline(xintercept = coef(fit1)[1],
linetype = "dashed",
color = "gray") +
geom_hline(yintercept = coef(fit1)[2],
linetype = "dashed",
color = "gray") +
scale_fill_brewer(palette = "RdBu",
name = "loss") +
labs(x = expression(beta[1]),
y = expression(beta[2])) +
coord_fixed()
plot(cont1)
この図の中央ほど損失が小さく、図の外側に近づくほど損失が大きくなることがわかる。
パラメタに制約を課さない場合には、パラメタ空間内の任意の点を選べるので、残差平方和\(S(\boldsymbol{\beta})\) を最小化する点（2本の点線が交わる点）をパラメタの推定値と利用する。
制約付きの最小二乗学習
パラメタ空間に制約を付けた以下の損失関数を最小化する問題を考える。 \[ L(\boldsymbol{\beta}) = S(\boldsymbol{\beta}) + \lambda |R(\boldsymbol{\beta})| \]
- \(\lambda\)：正則化パラメタ
- \(R(\boldsymbol{\beta})\)：正則化項
パラメタ空間に\(L_2\)制約を課す
正則化項に\(L_2\)ノルムの二乗を使う。\(\boldsymbol{\beta}\)の\(L_2\)ノルムは、 \[ \|\boldsymbol{\beta}\|_2 = \sqrt{\boldsymbol{\beta}^\top \boldsymbol{\beta}} = \sqrt{\beta_1^2 + \beta_2^2} \] である。
これを使うと、損失関数は、 \[ L(\boldsymbol{\beta}) = S(\boldsymbol{\beta}) + \lambda \| (\boldsymbol{\beta})\|_2^2 \] と、なる。
\(\lambda\)がある程度大きければ、\(L\)を最小化するために\(\|\boldsymbol{\beta}\|_2\) を小さくすることが必要になる。
\(\|\boldsymbol{\beta}\|_2\) を小さくするという制約は、どのようなことを意味するのだろうか。たとえば、 \(\|\boldsymbol{\beta}\|_2^2 \leq 1\) というのはどのような制約だろうか。 \[ \|\boldsymbol{\beta}\|_2^2 = \beta_1^2 + \beta_2^2 \leq 1 \] なので、パラメタ空間が原点を中心とする半径1の円の内側（境界線を含む）に制限されるということである。
この状況を図に表す。
L2 <- tibble(x = seq(-1, 1, length.out = 500),
y_p = sqrt(1 - x^2),
y_n = -1 * y_p) %>%
pivot_longer(cols = y_p:y_n,
names_to = "sign",
values_to = "y")
L2_1 <- cont1 +
geom_path(data = L2, aes(x = x, y = y, group = sign)) +
geom_vline(xintercept = 0,
linetype = "dotted",
color = "black") +
geom_hline(yintercept = 0,
linetype = "dotted",
color = "black") +
coord_fixed()
plot(L2_1)
この例では、\(\|\boldsymbol{\beta}\|_2 \leq 1\) という制約により、\(\hat{\beta}_1\)も\(\hat{\beta}_2\)も0に近づくことがわかる。
実際の分析では、\(\lambda\)の値によって制約の強さを調整する。 また、パラメタの次元は2よりも大きいことがほとんである。
パラメタ空間に\(L_1\)制約を課す
正則化項に\(L_1\)ノルムを使う。\(\boldsymbol{\beta}\)の\(L_1\)ノルムは、 \[ \|\boldsymbol{\beta}\|_1 = |\beta_1| + |\beta_2| \] である。
これを使って損失関数を \[ L(\boldsymbol{\beta}) = S(\boldsymbol{\beta}) + \lambda \|(\boldsymbol{\beta})\|_1 \] と、する。
\(\lambda\)がある程度大きければ、\(L\)を最小化するために\(\|\boldsymbol{\beta}\|_1\) を小さくすることが必要になる。
\(\|\boldsymbol{\beta}\|_1\) を小さくするという制約は、どのようなことを意味するのだろうか。たとえば、 \(\|\boldsymbol{\beta}\|_1 \leq 1\) というのはどのような制約だろうか。 \[ \|\boldsymbol{\beta}\|_1 = |\beta_1| + |\beta_2| \leq 1 \] なので、パラメタ空間が点\((1, 0), (0, 1), (-1, 0), (0, -1)\)を通る正方形の内側（境界線を含む）に制限されるということである。
この状況を図に表す。
L1 <- tibble(x = seq(-1, 1, length.out = 100),
y_p = 1 - abs(x),
y_n = -1 * y_p) %>%
pivot_longer(cols = y_p:y_n,
names_to = "sign",
values_to = "y")
L1_1 <- cont1 +
geom_path(data = L1, aes(x = x, y = y, group = sign)) +
geom_vline(xintercept = 0,
linetype = "dotted",
color = "black") +
geom_hline(yintercept = 0,
linetype = "dotted",
color = "black") +
coord_fixed()
plot(L1_1)
この例では、\(\|\boldsymbol{\beta}\|_1 \leq 1\) という制約により、（おそらく）\(\hat{\boldsymbol{\beta}} = (1, 0)\) になるだろうということがわかる。 このように、\(L_1\)制約を課すと、パラメタの一部がゼロになる。
実際の分析では、\(\lambda\)の値によって制約の強さと調整する。 また、パラメタの次元は2よりも大きいことがほとんである。
制約付き最小二乗学習の実行
ISLR パッケージに含まれるCredit データを使って実習を行う。
data(Credit)
summary(Credit)## ID Income Limit Rating
## Min. : 1.0 Min. : 10.35 Min. : 855 Min. : 93.0
## 1st Qu.:100.8 1st Qu.: 21.01 1st Qu.: 3088 1st Qu.:247.2
## Median :200.5 Median : 33.12 Median : 4622 Median :344.0
## Mean :200.5 Mean : 45.22 Mean : 4736 Mean :354.9
## 3rd Qu.:300.2 3rd Qu.: 57.47 3rd Qu.: 5873 3rd Qu.:437.2
## Max. :400.0 Max. :186.63 Max. :13913 Max. :982.0
## Cards Age Education Gender Student
## Min. :1.000 Min. :23.00 Min. : 5.00 Male :193 No :360
## 1st Qu.:2.000 1st Qu.:41.75 1st Qu.:11.00 Female:207 Yes: 40
## Median :3.000 Median :56.00 Median :14.00
## Mean :2.958 Mean :55.67 Mean :13.45
## 3rd Qu.:4.000 3rd Qu.:70.00 3rd Qu.:16.00
## Max. :9.000 Max. :98.00 Max. :20.00
## Married Ethnicity Balance
## No :155 African American: 99 Min. : 0.00
## Yes:245 Asian :102 1st Qu.: 68.75
## Caucasian :199 Median : 459.50
## Mean : 520.01
## 3rd Qu.: 863.00
## Max. :1999.00これらの変数のうち、応答変数はBalance（利用残高）であり、個体を識別する変数であるIDを除く変数が特徴量である。
このデータを訓練データ (training data) と検証データ (test data) に分ける。そのために、rsample::initial_split() を使う。 全体の 3/4 の観測点を訓練用に使うことにする。
set.seed(2021-10-14)
Credit_split <- initial_split(Credit,
prop = 3/4)
Credit_train <- training(Credit_split) # 訓練（学習）用
Credit_test <- testing(Credit_split) # 検証（テスト）用nrow(Credit)## [1] 400nrow(Credit_train)## [1] 300nrow(Credit_test)## [1] 100訓練データCredit_trainを使い、応答変数と説明変数をそれぞれ標準化し、別々のオブジェクトとして保存しよう。 変数の標準化には、scale() 関数が使える。
y <- Credit_train %>%
pull(Balance) %>%
scale()特徴量には、質的変数がいくつか含まれている。 Gender、Student、Married は二値のカテゴリ変数なので、ダミー変数に変換する。 Ethnicity は3つのカテゴリをもつカテゴリ変数なので、2つのダミー変数を作る。
X <- Credit_train %>%
select(Income:Ethnicity) %>%
mutate(across(where(is.numeric), scale, .names = NULL)) %>%
mutate(Female = ifelse(Gender == "Female", 1, 0),
Student = ifelse(Student == "Yes", 1, 0),
Married = ifelse(Married == "Yes", 1, 0),
African = ifelse(Ethnicity == "African American", 1, 0),
Asian = ifelse(Ethnicity == "Asian", 1, 0)) %>%
select(!c(Gender, Ethnicity)) %>%
as.matrix()検証用データも同じ形式に揃えておく。
y_test <- Credit_test %>%
pull(Balance) %>%
scale()
X_test <- Credit_test %>%
select(Income:Ethnicity) %>%
mutate(across(where(is.numeric), scale, .names = NULL)) %>%
mutate(Female = ifelse(Gender == "Female", 1, 0),
Student = ifelse(Student == "Yes", 1, 0),
Married = ifelse(Married == "Yes", 1, 0),
African = ifelse(Ethnicity == "African American", 1, 0),
Asian = ifelse(Ethnicity == "Asian", 1, 0)) %>%
select(!c(Gender, Ethnicity)) %>%
as.matrix()制約なしの最小二乗学習
ols <- lm(y ~ -1 + X)
tidy(ols) %>%
select(term, estimate)## # A tibble: 11 × 2
## term estimate
## <chr> <dbl>
## 1 XIncome -0.593
## 2 XLimit 0.950
## 3 XRating 0.388
## 4 XCards 0.0572
## 5 XAge -0.0402
## 6 XEducation -0.0178
## 7 XStudent 0.929
## 8 XMarried -0.0642
## 9 XFemale -0.0460
## 10 XAfrican -0.0496
## 11 XAsian -0.00583RSSは、
sum(ols$residuals^2)## [1] 13.53664である。
決定係数は、
summary(ols)$r.squared## [1] 0.9547269ここで学習した\(\hat{f}_{\mathrm{o}}\) を利用して、検証データの\(y\)を予測してみる。まず、予測値を求める。
pred_ols <- X_test %*% as.matrix(coef(ols), ncol = 1)予測誤差 \(y - \hat{y}\) の二乗の平均値 (MSE) を計算する。
(mse_ols <- mean((y_test - pred_ols)^2))## [1] 0.0490556制約なしの最小二乗学習の平均二乗誤差は、0.0490556 である。
\(L_2\)制約付き最小二乗学習（リッジ回帰）
リッジ回帰を実行するためにglmnet::glmnet() を使う。glmnet() は正則化回帰を行うための関数で、リッジ回帰（\(L_2\)制約付き回帰）を実行するには引数 alpha = 0 を指定する。特徴量は既に標準化済みなので standerdize = FALSEに、切片はないので（応答変数も標準化済みなので） intercetp = FALSE にする。
ridge <- glmnet(x = X,
y = y,
alpha = 0,
intercept = FALSE,
standardize = FALSE)glmnet()は、正則化パラメタ\(\lambda\) (lambda) がさまざまな値をとる場合について、推定対象であるパラメタ\(\boldsymbol{\beta}\)の推定を実行してくれる。推定結果を図示してみよう。
plot(ridge, xvar = "lambda", label = TRUE)
この図の各線が、lambda の値に応じてそれぞれのパラメタ\(\beta_d\) (\(d = 1, 2, \dots, 11\)) の推定値がどのように変化するかを示している。たとえば、「1」というラベルがついた黒い線は、上で作った特徴量行列Xの第1列であるIncome（所得）の係数の推定値である。 それぞれの線について、lambdaの値が0のときには、制約なしの最小二乗学習の結果に一致する。lambdaの値が大きくなるにつれて、\(\beta_d\)は0に近づく。
glmnet() が利用したlambdaの値を確認してみよう。
ridge$lambda## [1] 856.34662643 780.27115239 710.95401379 647.79481873 590.24651249
## [6] 537.81063915 490.03302428 446.49984103 406.83402579 370.69201225
## [11] 337.76075559 307.75502099 280.41491316 255.50362515 232.80538731
## [16] 212.12359836 193.27912254 176.10873801 160.46372311 146.20856822
## [21] 133.21980200 121.38492197 110.60142006 100.77589474 91.82324201
## [26] 83.66591827 76.23326869 69.46091522 63.29019896 57.66767212
## [31] 52.54463507 47.87671452 43.62347916 39.74808950 36.21697878
## [36] 32.99956221 30.06797206 27.39681629 24.96295864 22.74531820
## [41] 20.72468682 18.88356277 17.20599910 15.67746556 14.28472273
## [46] 13.01570733 11.85942777 10.80586889 9.84590528 8.97122220
## [51] 8.17424354 7.44806628 6.78640060 6.18351547 5.63418900
## [56] 5.13366318 4.67760269 4.26205736 3.88342793 3.53843490
## [61] 3.22409009 2.93767082 2.67669624 2.43890592 2.22224024
## [66] 2.02482254 1.84494289 1.68104325 1.53170400 1.39563164
## [71] 1.27164757 1.15867790 1.05574415 0.96195475 0.87649734
## [76] 0.79863174 0.72768349 0.66303810 0.60413562 0.55046588
## [81] 0.50156401 0.45700644 0.41640725 0.37941478 0.34570861
## [86] 0.31499681 0.28701336 0.26151589 0.23828354 0.21711508
## [91] 0.19782718 0.18025276 0.16423960 0.14964901 0.13635460
## [96] 0.12424123 0.11320399 0.10314726 0.09398394 0.08563466このように、推定結果として得られたridgeのなかに、全部で100個のlambdaの値が大きい順に保存されていることがわかる。
lambdaの値が最も小さい（\(\lambda\) = 0.0856347）ときのパラメタの推定値を確認しよう。
coef(ridge)[, 100]## (Intercept) Income Limit Rating Cards Age
## 0.000000000 -0.388885573 0.565717389 0.561347602 0.050374228 -0.053147370
## Education Student Married Female African Asian
## -0.005389409 0.480420539 -0.045802632 -0.012373962 -0.022872021 0.006398403制約なしの最小二乗学習の結果よりも、それぞれの係数の推定値が0に近づいている（縮小されている）ことがわかる。 このときの、\(L_2\)ノルムを確認する。
coef(ridge)[, 100]^2 %>%
sum() %>%
sqrt()## [1] 1.012615lambdaの値が最も大きい（\(\lambda\) = 856.3466264）ときのパラメタの推定値を確認しよう。
coef(ridge)[, 1]## (Intercept) Income Limit Rating Cards
## 0.000000e+00 4.776402e-37 8.617265e-37 8.649966e-37 1.107144e-37
## Age Education Student Married Female
## -8.467948e-39 -2.817205e-38 9.051341e-38 3.415049e-39 8.674827e-39
## African Asian
## 1.270435e-38 3.062470e-40すべての係数の推定値がほぼ0になっていることがわかる。このときの\(L_2\)ノルムを確認する。
coef(ridge)[, 1]^2 %>%
sum() %>%
sqrt()## [1] 1.319279e-36lambda の値が大きくなると\(L_2\)ノルムが小さくなることがわかる。
predict() を使えば、lambda の値を指定して推定値を得ることができる。 lambda = 50 の場合の推定値を求めてみよう。
predict(ridge, s = 50, type = "coefficients")## 12 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) .
## Income 8.761427e-03
## Limit 1.627683e-02
## Rating 1.633978e-02
## Cards 2.109893e-03
## Age -2.646197e-04
## Education -5.130363e-04
## Student 1.789117e-03
## Married 5.385195e-05
## Female 1.675670e-04
## African 2.304232e-04
## Asian 1.042415e-05lambdaの値が大きいときにはパラメタの推定値が0に近づき、学習の分散が小さくなることが期待できる。代わりに、バイアスは大きくなる。バイアスと分散のバランスをとり、MSE（あるいは他の指標）が小さくなるようにlambda の値を選んで推定を行う必要がある。
交差検証 (cross validation) と呼ばれる手法（詳細は今後の授業で説明する）を利用して、MSE（平均二乗誤差）が最小になる lambda を見つける。
ridge_cv <- cv.glmnet(x = X,
y = y,
alpha = 0,
intercept = FALSE,
standardize = FALSE)
plot(ridge_cv)
MSE の最小値。
min(ridge_cv$cvm)## [1] 0.08606744MSEを最小にするlambdaの値。
ridge_cv$lambda.min## [1] 0.08563466このlambdaの値で\(L_2\)制約を課した最小二乗推定値は、
predict(ridge, s = ridge_cv$lambda.min, type = "coefficients")## 12 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) .
## Income -0.388885573
## Limit 0.565717389
## Rating 0.561347602
## Cards 0.050374228
## Age -0.053147370
## Education -0.005389409
## Student 0.480420539
## Married -0.045802632
## Female -0.012373962
## African -0.022872021
## Asian 0.006398403である。
ここで学習した\(\hat{f}_{\mathrm{r}}\) を利用して、検証データの\(y\)を予測してみる。まず、予測値を求める。
pred_ridge <- predict(ridge, s = ridge_cv$lambda.min, newx = X_test)予測誤差 \(y - \hat{y}\) の二乗の平均値 (MSE) を計算する。
(mse_ridge <- mean((y_test - pred_ridge)^2))## [1] 0.07501783\(L_2\)制約付き最小二乗学習の平均二乗誤差は、0.0750178 である。
\(L_1\)制約付き最小二乗学習（ラッソ）
ラッソの実行にもglmnet::glmnet() を使う。ラッソでは alpha = 1 を指定する。
lasso <- glmnet(x = X,
y = y,
alpha = 1,
intercept = FALSE,
standardize = FALSE)推定結果を図示してみよう。
plot(lasso, xvar = "lambda", label = TRUE)
この図の各線が、lambda の値に応じてそれぞれのパラメタ\(\beta_d\) (\(d = 1, 2, \dots, 11\)) の推定値がどのように変化するかを示している。たとえば、「1」というラベルがついた黒い線は、上で作った特徴量行列Xの第1列であるIncome（所得）の係数の推定値である。 それぞれの線について、lambdaの値が0のときには、制約なしの最小二乗学習の結果に一致する。lambdaの値が大きくなると、いくつかの\(\beta_d\)がぴったり0になることがわかる。
glmnet() が利用したlambdaの値を確認してみよう。
lasso$lambda## [1] 0.856346626 0.780271152 0.710954014 0.647794819 0.590246512 0.537810639
## [7] 0.490033024 0.446499841 0.406834026 0.370692012 0.337760756 0.307755021
## [13] 0.280414913 0.255503625 0.232805387 0.212123598 0.193279123 0.176108738
## [19] 0.160463723 0.146208568 0.133219802 0.121384922 0.110601420 0.100775895
## [25] 0.091823242 0.083665918 0.076233269 0.069460915 0.063290199 0.057667672
## [31] 0.052544635 0.047876715 0.043623479 0.039748090 0.036216979 0.032999562
## [37] 0.030067972 0.027396816 0.024962959 0.022745318 0.020724687 0.018883563
## [43] 0.017205999 0.015677466 0.014284723 0.013015707 0.011859428 0.010805869
## [49] 0.009845905 0.008971222 0.008174244 0.007448066 0.006786401 0.006183515
## [55] 0.005634189 0.005133663 0.004677603 0.004262057 0.003883428 0.003538435
## [61] 0.003224090 0.002937671 0.002676696 0.002438906 0.002222240 0.002024823
## [67] 0.001844943 0.001681043 0.001531704 0.001395632 0.001271648 0.001158678
## [73] 0.001055744このように、推定結果として得られたlassoのなかに、全部で73個のlambdaの値が大きい順に保存されていることがわかる。
lambdaの値が最も小さい（\(\lambda\) = 1.0557442）ときのパラメタの推定値を確認しよう。
coef(lasso)[, 73]## (Intercept) Income Limit Rating Cards Age
## 0.00000000 -0.58805365 0.86146327 0.47165992 0.05243384 -0.03944465
## Education Student Married Female African Asian
## -0.01586052 0.91327095 -0.06494059 -0.04476163 -0.04427810 0.00000000制約なしの最小二乗学習の結果よりも、推定値が0になったパラメタがあることがわかる。
lambdaの値がそこそこ大きい（\(\lambda\) = 232.8053873）ときのパラメタの推定値を確認しよう。
coef(lasso)[, 15]## (Intercept) Income Limit Rating Cards Age
## 0.0000000 0.0000000 0.0000000 0.6256267 0.0000000 0.0000000
## Education Student Married Female African Asian
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000Rating 以外の推定値がぴったり0になっていることがわかる。また、Ratingの推定値の絶対値は、制約なしの場合の推定値である0.3879896 よりも大きいことがわかる。この結果から、ラッソは説明変数の選択を行っていると考えることができる。
lambdaの値が最も大きい（\(\lambda\) = 856.3466264）ときのパラメタの推定値を確認しよう。
coef(lasso)[, 1]## (Intercept) Income Limit Rating Cards Age
## 0 0 0 0 0 0
## Education Student Married Female African Asian
## 0 0 0 0 0 0すべての係数の推定値がぴったり0になっていることがわかる。
lambdaの値が大きいときにはパラメタの推定値が0に近づき、学習の分散が小さくなることが期待できる。代わりに、バイアスは大きくなる。バイアスと分散のバランスをとり、MSE（あるいは他の指標）が小さくなるようにlambda の値を選んで推定を行う必要がある。
lasso_cv <- cv.glmnet(x = X,
y = y,
alpha = 1,
intercept = FALSE,
standardize = FALSE)
plot(lasso_cv)
MSEの最小値。
min(lasso_cv$cvm)## [1] 0.04878917MSEを最小にするlambdaの値。
lasso_cv$lambda.min## [1] 0.001531704このlambdaの値を使って\(L_1\)制約を課した最小二乗学習による推定値は、
predict(lasso, s = lasso_cv$lambda.min, type = "coefficients")## 12 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) .
## Income -0.58581081
## Limit 0.85045185
## Rating 0.48038868
## Cards 0.05160024
## Age -0.03914746
## Education -0.01518219
## Student 0.90704914
## Married -0.06444820
## Female -0.04398252
## African -0.04249131
## Asian .である。
ここで学習した\(\hat{f}_{\mathrm{l}}\) を利用して、検証データの\(y\)を予測してみよう。まず、予測値を求める。
pred_lasso <- predict(lasso, s = lasso_cv$lambda.min, newx = X_test)予測誤差 \(y - \hat{y}\) の二乗の平均値 (MSE) を計算する。
(mse_lasso <- mean((y_test - pred_lasso)^2))## [1] 0.04848514\(L_1\)制約付き最小二乗学習の平均二乗誤差は、0.0484851 である。
最小二乗学習の比較
ここで実行した3つの最小二乗学習について、予測の平均二乗誤差をもう1度確認しよう。
cbind(mse_ols, mse_ridge, mse_lasso)## mse_ols mse_ridge mse_lasso
## [1,] 0.0490556 0.07501783 0.04848514この分析例では、ラッソの予測が最も正確であることがわかる。しかし、その結果は、制約なしの最小二乗学習とあまり変わらない。今回の例では、制約を課さなくても、大きな過剰適合は起きていないようである。リッジ回帰の予測力は他の学習法を使った場合に比べて低い。制約を課せば予測性能が必ず良くなるというわけではないことがわかる。