予測性の良い

モデルを構築する1：

罰則項付き最小二乗法

汎化性能が高いモデルが欲しい

　我々は説明変数をモデルに組み込むことで、その説明変数を使って興味ある現象を予測しようとする。予測の言うからには、今手元にあるデータだけでなく、将来取れるデータでも同じような結果になっていないと意味がないだろう。予測精度の良いモデルとは、どのように構築すればよいのだろう。素朴な感覚としては、今手元のデータで構築したモデルの誤差が小さいものが予測精度がよさそうに思える。ところが、この誤差は説明変数を増やせば増やすほど、言い換えればモデルを複雑にするほど、小さくできる。例え、その説明変数が予測にあまり寄与しなくても、小さくできるのである。もっと厄介なことに、モデルが過剰に複雑だと、新たに得られたデータに関しての予測には全く対応できない状態になる。このような状態は、元のデータだけを過剰に学習した状態、過学習overfittingと呼ばれる。また、傾向として過学習状態では、説明変数の絶対値が大きくなる傾向がある。一方で、説明変数が少なすぎてデータの予測ができない状態を未学習underfittingと呼ぶ。過学習も未学習も予測精度を下げるため、統計モデリングを行う上で避けなければならない。手元のデータだけでなく、将来取れるデータなどでも予測性を担保することを汎化generalizationという。つまり、私たちが最も望んでいるのは、汎化性能の高いモデルなのだ。では、どうやって説明変数は選んでいけばよいのだろう。人間が選べば、そこに恣意的な力が働く可能性も捨てきれない。本稿では、このような変数選択の術として、罰則項付き最小二乗法（正則化最小二乗法）regularized least squares methodについて紹介する。

説明変数が多すぎるとモデルの予測性が下がる

　とはいえ、上記の現象は本当なのか、気になるところだろう。そこでまずは上記の現象の確認を行う。以下のようなデータを解析しよう。

------------------------------------------------------

library(plotn)

x1 <- c(22.3, 18.9, 0.4, 5.5, -1.3, 3.8, 22.1, 12.1, 6.2, 13.9, 6.1,

24, 13.2, 18.3, 14.2, -5.4, 9.3, -3.6, -7.5, -8, 7.8, 4.5, 9,

11.1, -0.5, 1.9, 19.8, 17.9, 14.1, -4.1, -4.7, -4.4, 5.5, 26.8,

15.6, 9.5, 29.2, 13.7, -5.7, -8.9, -9, 9.3, 12.2, 16.3, 21.4,

-4, 0.3, 26.6, 14.7, 2, 13.9, -4.7, 27.8, -3.3, -8.5, 21.5, 20.4)

y <- c(-51.3, -39.6, -7.2, -28.9, -53.4, 11, -74.5, -44.1, 1.6, -2.3,

-8.4, -58.8, -41.5, -25.8, -31.9, -0.9, 11.4, -2.3, -1.8, 29.7,

-15.7, -17.4, -21.8, -26, 16.7, 14.7, -33.2, -47.4, 6.8, 47.6,

39.3, 18.5, -26.1, -27, -28.9, -11.9, -38.1, 3.1, -5.3, 14.3,

-9.9, -7.9, -30.9, 10.8, -46.7, 17.6, 45.1, -62.8, -24.1, -5.7,

-10, 8.1, -53.6, -7.5, 21.6, -32.2, -32.2)

d <- data.frame(x1 = x1, y = y, group = 1)

head(d)

## x1 y group

## 1 22.3 -51.3 1

## 2 18.9 -39.6 1

## 3 0.4 -7.2 1

## 4 5.5 -28.9 1

## 5 -1.3 -53.4 1

## 6 3.8 11.0 1

plotn(y ~ x1, data = d) #図1の描画

------------------------------------------------------

図1 データの様子

　さあ、このデータ、皆さんはどんな回帰を考えるだろう。1次回帰？、ちょっと曲がってる気もするし2次回帰？　正直、データ数も少ないし、決めるのは難しい盤面だ。けれども、生態学ではこういうデータを解析する必要が出てくるものだ。そこで、とりあえずx1の多項式で表される信じて、4次式までの多項式回帰を試してみよう。

------------------------------------------------------

fit1 <- lm(y ~ 1, data = d) #0次回帰

fit2 <- summary(fit1)

fit2

## Call:

## lm(formula = y ~ 1, data = d)

## Residuals:

## Min 1Q Median 3Q Max

## -60.867 -18.567 3.633 20.433 61.233

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) -13.633 3.657 -3.728 0.000452 ***

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Residual standard error: 27.61 on 56 degrees of freedom

fit3 <- lm(y ~ x1, data = d) #1次回帰

fit4 <- summary(fit3)

fit4

## Call:

## lm(formula = y ~ x1, data = d)

## Residuals:

## Min 1Q Median 3Q Max

## -57.418 -12.369 -2.204 13.313 43.990

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 1.655 3.230 0.512 0.61

## x1 -1.817 0.235 -7.733 2.38e-10 ***

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Residual standard error: 19.29 on 55 degrees of freedom

## Multiple R-squared: 0.5209, Adjusted R-squared: 0.5122

## F-statistic: 59.79 on 1 and 55 DF, p-value: 2.383e-10

fit5 <- lm(y ~ x1 + I(x1^2), data = d) #2次回帰

fit6 <- summary(fit5)

fit6

## Call:

## lm(formula = y ~ x1 + I(x1^2), data = d)

## Residuals:

## Min 1Q Median 3Q Max

## -57.728 -12.189 -0.869 12.052 43.178

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 2.38015 3.38912 0.702 0.48551

## x1 -1.52074 0.46615 -3.262 0.00192 **

## I(x1^2) -0.01705 0.02310 -0.738 0.46369

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Residual standard error: 19.37 on 54 degrees of freedom

## Multiple R-squared: 0.5257, Adjusted R-squared: 0.5081

## F-statistic: 29.92 on 2 and 54 DF, p-value: 1.797e-09

fit7 <- lm(y ~ x1 + I(x1^2) + I(x1^3), data = d) #3次回帰

fit8 <- summary(fit7)

fit8

## Call:

## lm(formula = y ~ x1 + I(x1^2) + I(x1^3), data = d)

## Residuals:

## Min 1Q Median 3Q Max

## -59.594 -11.504 -1.864 11.566 41.107

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 4.429609 4.679656 0.947 0.3482

## x1 -1.438476 0.486083 -2.959 0.0046 **

## I(x1^2) -0.060399 0.071703 -0.842 0.4034

## I(x1^3) 0.001495 0.002339 0.639 0.5256

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Residual standard error: 19.47 on 53 degrees of freedom

## Multiple R-squared: 0.5293, Adjusted R-squared: 0.5026

## F-statistic: 19.86 on 3 and 53 DF, p-value: 9.268e-09

fit9 <- lm(y ~ x1 + I(x1^2) + I(x1^3) + I(x1^4), data = d) #4次回帰

fit10 <- summary(fit9)

fit10

## Call:

## lm(formula = y ~ x1 + I(x1^2) + I(x1^3) + I(x1^4), data = d)

## Residuals:

## Min 1Q Median 3Q Max

## -58.265 -12.548 -1.948 10.393 42.026

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 3.4219934 5.1516507 0.664 0.509

## x1 -1.1492949 0.7717053 -1.489 0.142

## I(x1^2) -0.0339957 0.0904585 -0.376 0.709

## I(x1^3) -0.0028603 0.0092863 -0.308 0.759

## I(x1^4) 0.0001129 0.0002329 0.485 0.630

## Residual standard error: 19.62 on 52 degrees of freedom

## Multiple R-squared: 0.5314, Adjusted R-squared: 0.4954

## F-statistic: 14.74 on 4 and 52 DF, p-value: 4.089e-08

------------------------------------------------------

さて、負の相関があるのはほぼ確実なので、1次回帰の結果は納得だろう。問題は2次以上の項を考慮するかである。さて、当てはまりの良さの指標であるR二乗値を確認すると、モデルが複雑にするほど大きな値を示していることがわかるだろう。これらの関係を図示してみよう。

------------------------------------------------------

xx <- seq(min(d$x1), max(d$x1), length = 200)

coefs1 <- coef(fit2)[,1]

yy1 <- rep(coefs1[1], length(xx))

coefs2 <- coef(fit4)[,1]

yy2 <- coefs2[1] + coefs2[2]*xx

coefs3 <- coef(fit6)[,1]

yy3 <- coefs3[1] + coefs3[2]*xx + coefs3[3]*(xx)^2

coefs4 <- coef(fit8)[,1]

yy4 <- coefs4[1] + coefs4[2]*xx + coefs4[3]*(xx)^2 + coefs4[4]*(xx)^3

coefs5 <- coef(fit10)[,1]

yy5 <- coefs5[1] + coefs5[2]*xx + coefs5[3]*(xx)^2 + coefs5[4]*(xx)^3 + coefs5[5]*(xx)^4

plotn(y ~ x1, data = d, legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1,

pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図2の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

points(xx, yy2, col = "red", type = "l"),

points(xx, yy3, col = "blue", type = "l"),

points(xx, yy4, col = "darkgreen", type = "l"),

points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

図2 データの様子とそれに当てはめた多項式回帰、0～4は最大の次数を示す

上の図はどう映るだろうか。確かに3次や4次はデータに合わせて、曲線をフィットしているので、データの平均値を通っている気もする。けれども一方で、本当に真の構造はこんなにうねうねしてるのかという気持ちにもなる。次に、考えられるシンプルな変数選択は、有意な効果があればそれを選択する方法である。ところが、1次の項ですら、モデルによって有意になったりならなかったりだ。果たして、この中から予測性のいいモデルを選べるだろうか。また、確認しておくと、係数の絶対値の和はモデルが複雑になるほど、基本的に大きな値を示す。

------------------------------------------------------

#係数の絶対値の和

sum(abs(coefs1))

## [1] 13.63333

sum(abs(coefs2))

## [1] 3.472766

sum(abs(coefs3))

## [1] 3.917942

sum(abs(coefs4))

## [1] 5.92998

sum(abs(coefs5))

## [1] 4.608257

------------------------------------------------------

　そこで、始めに得たデータと同じ条件で、新たに次のようなデータが取れたとする。このデータも図示してみよう。

------------------------------------------------------

x1_2 <- c(-4.7, 16.8, 28.1, 9.7, 25.5, -4.6, 1.3, 13.7, 5.4, 12.7,

12.3, 25.6, 17.8, 24.7, -8.3, 22.4, -5.5, 18.9, 3.3, 26.5,

27.4, 8.1, 4.5, 5, 19.3, 20.5, 14.7, -8.9, 29.7, 27, -3.8,

11.7, 0.8, 3.2, 26.1, 8.9, 25.3, -7.7, 0.7, -8.4, 28.9, 24.1,

-6.4, 15.8, 13.1, 16.8, 8.2, 28.3, 22.9, 1.1, 22.3, 23.6,

-7.9, 4.7, -6, 3.8, 10.1)

y_2 <- c(0.1, -11.6, -76, -24.6, -29.4, 21.6, 22.2, -78.4, -11.8, -17.5,

-21.9, -47.2, -44.3, -26.9, 15.6, -36.1, 12, -25.1, -19.3, -2.2,

-46.7, -11.7, 3.5, 0.6, -35.1, -24.4, -32.2, 19, -80.3, -24, 9.5,

-22.6, -20.9, -15.8, -36.3, 1.5, 0.4, -9.2, 17.5, -3.1, -76.9,

-52.4, 23.4, -64.2, 1.9, -30.7, -20.6, -46.3, -5.4, 0, -63.6,

-30, 3.6, -27.8, 50.6, -20.7, 14.7)

d2 <- data.frame(x1 = x1_2, y = y_2, group = 2)

d <- rbind(d,d2)

plotn(y ~ x1, data = d, mode = "m", group = "group", pch = c(16, 2),

col.dot = c("black","red"),

legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1,

pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図3の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

points(xx, yy2, col = "red", type = "l"),

points(xx, yy3, col = "blue", type = "l"),

points(xx, yy4, col = "darkgreen", type = "l"),

points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

図3 データの様子とそれに当てはめた多項式回帰、0～4は最大の次数を示す。赤い三角は新たに得られたデータ。

どうだろう。なんだか、両端部分は結構怪しいんじゃないか、そんな気がしてくるだろう。では、予測精度の指標として、モデルから予測される値と実際のデータの偏差の二乗和を確認してみよう。group 1のデータで行った多項式回帰の係数を使って以下のように計算する。

------------------------------------------------------

#データ1のyの推定値

xe <- d$x1[d$group == 1]

ye1_1 <- coefs1[1]

ye1_2 <- coefs2[1] + coefs2[2]*xe

ye1_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye1_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye1_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4

error1 <- c(sum((d$y[d$group == 1] - ye1_1)^2),

sum((d$y[d$group == 1] - ye1_2)^2),

sum((d$y[d$group == 1] - ye1_3)^2),

sum((d$y[d$group == 1] - ye1_4)^2),

sum((d$y[d$group == 1] - ye1_5)^2))

error1

## [1] 42695.57 20456.67 20252.38 20097.53 20007.10

#データ2のyの推定値

xe <- d$x1[d$group == 2]

ye2_1 <- coefs1[1]

ye2_2 <- coefs2[1] + coefs2[2]*xe

ye2_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye2_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye2_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4

error2 <- c(sum((d$y[d$group == 2] - ye2_1)^2),

sum((d$y[d$group == 2] - ye2_2)^2),

sum((d$y[d$group == 2] - ye2_3)^2),

sum((d$y[d$group == 2] - ye2_4)^2),

sum((d$y[d$group == 2] - ye2_5)^2))

error2

## [1] 46232.83 21896.20 22444.60 22598.95 23676.39

barplotn(rbind(error1, error2), beside = T, ylab = "error", names = 0:4,

legend = T, leg.lab = c("group1", "group2"), leg.sp = 4) #図4の描画

------------------------------------------------------

図4 group1データで行った各多項式回帰の係数を使った誤差の大きさ。

group 1のデータの誤差は次数が大きくなるほど小さくなる。これはR二乗値からも確認できている。ところが、group 2のデータに対する誤差は次数が上がるほど大きくなることがわかる。つまり、高次項の多項式モデルは、group 1を過剰に学習しており、group 2のデータの予測はうまく行ってないのである。この状態がまさに過学習である。

　実は、これらのデータは以下のような設定で、1次式を基盤として作成している。新たにデータを生成して、その時の様子を確認しよう。

------------------------------------------------------

b1 <- 5

b2 <- -2

n <- 57

x1_3 <- runif(n, -10, 30)

y_3 <- rnorm(n = n, mean = b1 + b2*x1_3, sd = 20)

d3 <- data.frame(x1 = x1_3, y = y_3, group = 3)

d <- rbind(d,d3)

plotn(y ~ x1, data = d, mode = "m", group = "group", pch = c(16, 2, 4),

col.dot = c("black","red", "red"),

legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1,

pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図5の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

points(xx, yy2, col = "red", type = "l"),

points(xx, yy3, col = "blue", type = "l"),

points(xx, yy4, col = "darkgreen", type = "l"),

points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

図5 データの様子とそれに当てはめた多項式回帰、0～4は最大の次数を示す。赤い三角と〆は新たに得られたデータ。

いよいよ怪しい気分だろう。誤差を確認しよう。

------------------------------------------------------

#データ3のyの推定値

xe <- d$x1[d$group == 3]

ye3_1 <- coefs1[1]

ye3_2 <- coefs2[1] + coefs2[2]*xe

ye3_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye3_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye3_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4

error3 <- c(sum((d$y[d$group == 3] - ye3_1)^2),

sum((d$y[d$group == 3] - ye3_2)^2),

sum((d$y[d$group == 3] - ye3_3)^2),

sum((d$y[d$group == 3] - ye3_4)^2),

sum((d$y[d$group == 3] - ye3_5)^2))

error3

## [1] 48109.64 18615.40 19748.39 20001.61 20587.06

barplotn(rbind(error1, error2, error3), beside = T, ylab = "error", names = 0:4,

legend = T, leg.lab = c("group1", "group2", "group3"), leg.sp = 4) #図6の描画

------------------------------------------------------

図6 group1データで行った各多項式回帰の係数を使った誤差の大きさ。

すると、1次回帰が最もマシな成績をおさめていることがわかる。こんな感じで、過剰に複雑なモデルは必ずしもデータ一般に対して予測性は高いわけではないのである。

罰則項付き最小二乗法

　そこで汎化性能を高めるために考えられた手法が、罰則項付き最小二乗法regularized least squares methodである。正則化最小二乗法とも呼ばれ、英語の直訳はこちらである。ただ、罰則項付き最小二乗法のほうが、解析の様子が非常にイメージしやすい。罰則項を課すことを正則化するregularization、と表現することもある。最小二乗法は、モデルとデータの偏差二乗和を最小にするパラメータを選択するのであった。ところがこれだけだと、モデルをどんどん複雑にすればより小さな偏差平方和を得ることができるため、あたかも過剰に複雑なモデルが予測性が良いかのようになってしまう。そこで、以下の式を最小にするようなパラメータを探索する。

式の前半は今までの最小二乗法でも登場した偏差平方和そのものである。f(xi)と一般の形で表示しているが、わかりにくければb1 + b2 x1,iという形を代入して想像してもらってもよい。後半は罰則項（ペナルティー項、正則化項）と呼ばれる部分で、まさに式の形は罰則項付き最小二乗法である。罰則項は係数の和をλだけ補正して足している。つまり、もし偏差平方和部分で絶対値の大きな係数が多く推定されてしまうと、その分だけRは大きな値を示す、ことになる。ということは、式全体として最小化されるということは、極端に大きな係数が出ないように、場合によっては係数の値を0に制御しつつも、誤差を最小にするようなパラメータを推定していることになる。λは罰則項による制御の強さを表しており、大きくなればその分だけ罰則が強くなるので、より強く大きな係数が出ないような制御が働くことになる。逆にλが小さい、もしくは0であれば、それは通常の最小二乗法そのものである。λがこの解析の肝であるが、この値は人間が決めるハイパーパラメータである。どうやって決めるかは、あとで紹介することとして、上記の式を最小にする、とはどういうことかをもう少し掘り下げよう。

ラグランジュの未定乗数法に帰着して考える

　上記の罰則項について、以下の制約を満たすようにパラメータを求めることを考える。

この後に登場するのは上記の式についてp = 1 or 2の時なので、これらの時にg(b1……bm)が示す領域を考える。簡単のため2パラメータで考えよう。以下の図でグレーで示される領域がg(b1, b2) ≦ 0である。

このようなパラメータに制約があるときの最小値を求めるためには、ラグランジュの未定乗数法を活用するとよいのであった。そこでラグランジュ乗数λLを導入して、ラグランジュ関数Lを以下のように考える。

では、ラグランジュの未定乗数法の手続きにのっとると、

ここでベクトル∇h（ナブラh）と∇gは関数hや関数gの勾配ベクトルであり、勾配ベクトルはその関数の値が上昇する方向に向いているのであった。これを図示すると以下のようになる。

制約条件gはその内側はg < 0であり、外側はg > 0であるから、∇gは領域の外を向くようなベクトルである。一方、hはもともと偏差平方和の式だったので、だいたいは下に凸のおわん型の関数である。上図のようにhの最小値がg < 0の領域の外に位置する時、∇h = ∇gを満たす点における関数hの勾配ベクトルは∇gとは逆方向を向いていることになる。つまり、式に立ち戻ればλL < 0であることがわかる。そこで、λ'L = - λLと置きなおせば、ラグランジュ関数は罰則項付き二乗和の式を使って、以下のように書きなおせる。

Rを最小にするパラメータの条件はRをパラメータで偏微分し、偏微分係数が0になる点だから、上記の式で考えると、

以上から、Rの最小値を求める問題は、制約条件（罰則項）g ≦ 0の下での偏差平方和の最小値を求める問題に帰着できるということである。λ = λ'Lμであり、λ'Lはラグランジュの未定乗数法により求まる値だから、λの設定は、μの設定に依存していることになる。μと制約条件gの関係を図で確認すると、以下のようになる。

制約条件gの「円」の半径はμの逆数に依存しているから、μが小さいほど大きな円、μが大きいほど小さな円の制約条件となる。μとλは比例しているから、λが大きいときは円が小さくなる。逆に、λが小さいほど円が大きくなり、λ = 0とは無限大の半径を持つ円、すなわち、制約がない最小二乗法となる。当然、円の半径が大きいほど、偏差平方和の関数の真の最小値が制約条件g < 0に含まれる可能性が高くなるので、制約なしの最小二乗法と一致する。これらの点から、λ（μ）がまさに制約の強さである。

LASSO回帰/Ridge回帰

　では、罰則項の次数によって、パラメータ推定はどのように変わるのか、概要を見てみよう。特にp = 1の時をLASSO（ラッソ、Least Absolute Shrinkage and Selection Operator）回帰、p = 2の時をRidge（リッジ）回帰と呼ぶ。

上図の左がLASSO回帰、右がRidge回帰である。LASSO回帰では制約条件が四角形となっているため、鋭い頂点が存在する。そのため、通常の最小二乗法の最適パラメータと比較して、パラメータ全体が小さな値をとる、だけでなく効果の弱いパラメータが0になりやすくなっている。一方、Ridge回帰は効果の弱いパラメータを0にするだけの力はなく、全体的にパラメータを小さくするだけである。全体的にパラメータを小さくすることで、モデルは極端なブレを生み出すことが少なくなり、汎化性能が高まるのである。また、詳細には立ち入らないが、以下のようなLASSO回帰とRidge回帰をα:(1-α)の割合で混ぜた罰則を科すことで、これら2つの回帰のいいとこどりを目指すのが、エラスティック・ネットElastic netと呼ばれる方法である。

制約の強さλをどう決める？ーk-fold交差検証法

　制約の強さλは、k-fold交差検証法 k-fold cross-validationと呼ばれる方法で決定する。交差検証法をクロス・ヴァリデーションよ呼ぶこともある。上記の例では、データの生成過程がわかっていたので、新たにデータを生成してゆき、未知のデータでの当てはまりが悪くなることを示した。しかし、実際は1つのデータから予測性の高いモデルを構築せねばならない。そこで、手元のデータをk分割することを考える。分割したデータに個別にD1、D2……Dkと名前を付けることとして、1回目はD1を除いたデータを学習用として回帰を行い、パラメータを推定する。その後、推定したパラメータを使ったモデルとD1検証用データの偏差平方和（平均二乗誤差）を求める。次に、D2を除いたデータでパラメータを推定し、……とこれをDkまで繰り返す。例えば、k = 5であれば以下のようになる。最終検証用データに関しては後述する。

するとk個分の平均二乗誤差が得られる。平均的に小さな平均二乗誤差を得るには、学習用データに合わせるのではなくパラメータ推定に使わなかった検証用データとの平均二乗誤差が小さくならなければならない。このように、手元のデータを疑似的に既知データ（学習用データ）と未知データ（検証用データ）に分けることで、汎化を達成するのである。この手続きを様々なλで実行し、平均的に小さな平均二乗誤差が得られたλをハイパーパラメータとして用いる。

　なお、今回はそこまでの検証をしないが、実はこれだけでは汎化は不十分である。k個分割したとは言えども、すべてのデータを最終的には学習に使ってしまっているので、このk-foldの中で過学習しているかもしれないからである。したがって、厳密な評価には、学習に一度も使わないデータを最終検証用としてより分けておいて、出来上がったモデルと最終検証用データとの誤差をモデルの性能として評価するべきである。

Rで行う罰則項付き最小二乗法（LASSO回帰）

　では、Rで罰則項付き最小二乗法を実行してみよう。特に、効果の弱い項を削除するモデル選択に興味があるので、LASSO回帰を行う。glmnetパッケージが必要なので、あらかじめインストールしておこう。解析のデータは上記の例で出したものとする。glmnet関数の引数、alphaによって実行する回帰が異なり、alpha = 1ならLASSO回帰、0 < alpha < 1ならエラスティック・ネット、alpha = 0ならRidge回帰である（上記のエラスティック・ネットの式参照）。glmnet関数をそのまま実行すれば、様々なλでの解析結果を返してくれる。

------------------------------------------------------

library(glmnet)