Stanの開発者でもある統計学界の重鎮、Andrew Gelmanがこんなブログ記事をupしていました。
ちなみに@berobero11さんがこの件についてこんなコメントをされてました。
AkiらのPSIS-LOOがWAICより良いとする論文が出た。https://t.co/BWYNALp88K
— Kentaro Matsuura (@berobero11) October 21, 2016
渡辺先生の反論はこちら。https://t.co/MLQQvQuwM7
・WAICはCVではなく汎化誤差を小さくするもの
・MCMCを何度もやった場合の揺らぎを見よ
WAICとcross validationの関係については渡辺澄夫先生の本にも当然のように記述があって、そこで密林でポチって読んでみたのですが、いかんせん僕の数学力が弱すぎてまだ全然理解できてない有様ですorz
- 作者: 渡辺澄夫
- 出版社/メーカー: コロナ社
- 発売日: 2012/03
- メディア: 単行本
- 購入: 1人 クリック: 4回
- この商品を含むブログ (4件) を見る
というザマでこの件について自分では何も勉強出来てないんですが(汗)、そもそも論としてモデル選択においてAICとcross validationってどんな関係性にあるんだっけ?というのをふと思ったので、大雑把に実験してみようかと思います。
データセットの準備
既に挙動がはっきりしているものということで、テニス四大大会データセットを用います。
'men.txt'を学習データに、'women.txt'をholdoutデータに使います。
# Training > dm <- read.table('men.txt',sep='\t',header=T) > dm <- dm[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)] > dm$Result<-as.factor(dm$Result) # Holdout > dw <- read.table('women.txt',sep='\t',header=T) > dw <- dw[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]
何も考えずにロジスティック回帰モデルを推定して、AIC / LOO ACC / holdout ACCを出してみる
とりあえず素のロジスティック回帰をやった上で、AIC / LOO ACC / holdout ACCを出してみましょう。
> dm.glm <- glm(Result~., dm, family=binomial) > summary(dm.glm) Call: glm(formula = Result ~ ., family = binomial, data = dm) Deviance Residuals: Min 1Q Median 3Q Max -2.91243 -0.14289 -0.00641 0.08835 3.04139 Coefficients: (2 not defined because of singularities) Estimate Std. Error z value Pr(>|z|) (Intercept) -1.378486 5.064233 -0.272 0.785468 FSP.1 0.116125 0.065947 1.761 0.078258 . FSW.1 0.198896 0.043094 4.615 3.92e-06 *** SSP.1 NA NA NA NA SSW.1 0.290083 0.069708 4.161 3.16e-05 *** ACE.1 -0.040074 0.048374 -0.828 0.407430 DBF.1 -0.136684 0.082271 -1.661 0.096636 . WNR.1 0.074055 0.029644 2.498 0.012485 * UFE.1 -0.007523 0.030277 -0.248 0.803779 BPC.1 0.393523 0.102120 3.854 0.000116 *** BPW.1 0.313412 0.069018 4.541 5.60e-06 *** NPA.1 0.012960 0.042315 0.306 0.759398 NPW.1 -0.014493 0.036409 -0.398 0.690578 FSP.2 -0.065368 0.066044 -0.990 0.322288 FSW.2 -0.263955 0.047725 -5.531 3.19e-08 *** SSP.2 NA NA NA NA SSW.2 -0.247359 0.063257 -3.910 9.22e-05 *** ACE.2 0.026067 0.038823 0.671 0.501937 DBF.2 0.117291 0.084828 1.383 0.166759 WNR.2 -0.032345 0.028113 -1.151 0.249925 UFE.2 -0.024302 0.028906 -0.841 0.400490 BPC.2 -0.329055 0.111769 -2.944 0.003239 ** BPW.2 -0.351452 0.087206 -4.030 5.57e-05 *** NPA.2 -0.036414 0.059266 -0.614 0.538935 NPW.2 0.047514 0.042770 1.111 0.266603 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 680.42 on 490 degrees of freedom Residual deviance: 150.87 on 468 degrees of freedom AIC: 196.87 Number of Fisher Scoring iterations: 8 > res <- rep(0, nrow(dm)) > for (i in 1:nrow(dm)){ + train <- dm[-i,] + test <- dm[i,] + model <- glm(Result~., train, family=binomial) + pred <- round(predict(model, newdata=test, type='response'),0) + res[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0) + } > sum(res)/nrow(dm) [1] 0.9205703 > sum(diag(table(dw$Result, round(predict(dm.glm, newdata=dw[,-1], type='response'),0))))/nrow(dw) [1] 0.9269912
AIC = 196.87, LOO ACC = 0.921, holdout ACC = 0.927という結果になりました。
ステップワイズ法でAIC最小のモデルを選んで、AIC / LOO ACC / holdout ACCを出してみる
普通に{MASS}パッケージのstepAICを使うだけです。交互作用は想定していないのでご注意ください。
> library(MASS) > stepAIC(dm.glm) # 中略 # Df Deviance AIC <none> 153.86 185.86 - DBF.2 1 155.88 185.88 - FSP.2 1 156.11 186.11 - NPW.2 1 156.16 186.16 - DBF.1 1 156.25 186.25 - FSP.1 1 158.54 188.54 - UFE.2 1 161.87 191.87 - WNR.1 1 165.91 195.91 - SSW.2 1 176.82 206.82 - SSW.1 1 178.59 208.59 - BPW.2 1 180.71 210.71 - FSW.1 1 181.31 211.31 - BPC.2 1 182.61 212.61 - BPC.1 1 187.76 217.76 - BPW.1 1 189.80 219.80 - FSW.2 1 200.18 230.18 Call: glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2, family = binomial, data = dm) Coefficients: (Intercept) FSP.1 FSW.1 SSW.1 DBF.1 WNR.1 BPC.1 BPW.1 -0.67573 0.12962 0.18332 0.29546 -0.11832 0.05049 0.38807 0.32564 FSP.2 FSW.2 SSW.2 DBF.2 UFE.2 BPC.2 BPW.2 NPW.2 -0.09027 -0.25356 -0.26049 0.11584 -0.04035 -0.39068 -0.31849 0.02078 Degrees of Freedom: 490 Total (i.e. Null); 475 Residual Null Deviance: 680.4 Residual Deviance: 153.9 AIC: 185.9 > dm.glm.aic <- glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + + BPW.2 + NPW.2, family = binomial, data = dm) > summary(dm.glm.aic) Call: glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2, family = binomial, data = dm) Deviance Residuals: Min 1Q Median 3Q Max -2.8678 -0.1372 -0.0069 0.1008 3.0421 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.67573 5.02335 -0.135 0.89299 FSP.1 0.12962 0.06182 2.097 0.03602 * FSW.1 0.18332 0.03890 4.713 2.44e-06 *** SSW.1 0.29546 0.06685 4.420 9.87e-06 *** DBF.1 -0.11832 0.07914 -1.495 0.13490 WNR.1 0.05049 0.01536 3.287 0.00101 ** BPC.1 0.38807 0.08934 4.344 1.40e-05 *** BPW.1 0.32564 0.06178 5.271 1.36e-07 *** FSP.2 -0.09027 0.06117 -1.476 0.14003 FSW.2 -0.25356 0.04424 -5.732 9.94e-09 *** SSW.2 -0.26049 0.06132 -4.248 2.16e-05 *** DBF.2 0.11584 0.08149 1.422 0.15516 UFE.2 -0.04035 0.01475 -2.735 0.00624 ** BPC.2 -0.39068 0.09954 -3.925 8.67e-05 *** BPW.2 -0.31849 0.06928 -4.598 4.28e-06 *** NPW.2 0.02078 0.01396 1.488 0.13665 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 680.42 on 490 degrees of freedom Residual deviance: 153.86 on 475 degrees of freedom AIC: 185.86 Number of Fisher Scoring iterations: 8 > res.aic <- rep(0, nrow(dm)) > for (i in 1:nrow(dm)){ + train <- dm[-i,] + test <- dm[i,] + model <- glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + + BPW.2 + NPW.2, family = binomial, data = train) + pred <- round(predict(model, newdata=test, type='response'),0) + res.aic[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0) + } > sum(res.aic)/nrow(dm) [1] 0.9287169 > sum(diag(table(dw$Result, round(predict(dm.glm.aic, newdata=dw[,-1], type='response'),0))))/nrow(dw) [1] 0.9314159
AIC = 185.86, LOO ACC = 0.929, holdout ACC = 0.931という結果になりました。
L1正則化モデル
{glmnet}でalpha = 1としてL1正則化ロジスティック回帰でやってみます。ちなみにLOOは2つ考え方があるかと思うのですが、一応両論併記になるように「あらかじめ決め打ちでL1正則化モデルを固定して回す」のと「都度L1正則化モデルを10-folds CVでパラメータ決定して回す」のとを使い分けています。ただしholdout ACCは10-fold CVでパラメータ決定したものにしてあります。
> library(glmnet) > dm.cv.glmnet <- cv.glmnet(as.matrix(dm[,-1]), as.integer(dm[,1])-1, family="binomial", alpha=1) > plot(dm.cv.glmnet) > coef(dm.cv.glmnet, s=dm.cv.glmnet$lambda.min) 25 x 1 sparse Matrix of class "dgCMatrix" 1 (Intercept) 3.533402e-01 FSP.1 3.805604e-02 FSW.1 1.179697e-01 SSP.1 -3.275595e-05 SSW.1 1.475791e-01 ACE.1 . DBF.1 -8.934231e-02 WNR.1 3.628403e-02 UFE.1 -7.839983e-03 BPC.1 3.758665e-01 BPW.1 2.064167e-01 NPA.1 . NPW.1 . FSP.2 -2.924528e-02 FSW.2 -1.568441e-01 SSP.2 . SSW.2 -1.324209e-01 ACE.2 1.233763e-02 DBF.2 4.032510e-02 WNR.2 -2.071361e-02 UFE.2 -6.114823e-06 BPC.2 -3.648171e-01 BPW.2 -1.985184e-01 NPA.2 . NPW.2 1.340329e-02 # 一旦glmを使って直す > dm.glm.l1 <- glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + WNR.1 + + UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2, family = binomial, data = dm) > summary(dm.glm.l1) Call: glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + WNR.1 + UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2, family = binomial, data = dm) Deviance Residuals: Min 1Q Median 3Q Max -2.8984 -0.1487 -0.0076 0.1013 3.1776 Coefficients: (1 not defined because of singularities) Estimate Std. Error z value Pr(>|z|) (Intercept) -1.430594 5.007621 -0.286 0.775121 FSP.1 0.131931 0.063677 2.072 0.038276 * FSW.1 0.183218 0.039567 4.631 3.65e-06 *** SSP.1 NA NA NA NA SSW.1 0.290247 0.068486 4.238 2.25e-05 *** DBF.1 -0.129578 0.082095 -1.578 0.114479 WNR.1 0.067105 0.026901 2.495 0.012612 * UFE.1 -0.008339 0.028484 -0.293 0.769697 BPC.1 0.375908 0.088707 4.238 2.26e-05 *** BPW.1 0.313854 0.063923 4.910 9.11e-07 *** FSP.2 -0.080219 0.064170 -1.250 0.211257 FSW.2 -0.257784 0.046546 -5.538 3.05e-08 *** SSW.2 -0.250462 0.062454 -4.010 6.06e-05 *** ACE.2 0.031084 0.038078 0.816 0.414312 DBF.2 0.098090 0.082402 1.190 0.233895 WNR.2 -0.032804 0.027685 -1.185 0.236053 UFE.2 -0.017172 0.027139 -0.633 0.526898 BPC.2 -0.350231 0.100279 -3.493 0.000478 *** BPW.2 -0.318826 0.074215 -4.296 1.74e-05 *** NPW.2 0.020829 0.014036 1.484 0.137833 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 680.42 on 490 degrees of freedom Residual deviance: 152.04 on 472 degrees of freedom AIC: 190.04 Number of Fisher Scoring iterations: 8 > res.l1 <- rep(0, nrow(dm)) > for (i in 1:nrow(dm)){ + train <- dm[-i,] + test <- dm[i,] + model <- glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + WNR.1 + + UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 + + BPW.2 + NPW.2, family = binomial, data = train) + pred <- round(predict(model, newdata=test, type='response'),0) + res.l1[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0) + } > sum(res.l1)/nrow(dm) [1] 0.9205703 > sum(diag(table(dw$Result, round(predict(dm.glm.l1, newdata=dw[,-1], type='response'),0))))/nrow(dw) [1] 0.9269912 # {glmnet}縛り > res.l1 <- rep(0, nrow(dm)) > for (i in 1:nrow(dm)){ + train <- dm[-i,] + test <- dm[i,] + model <- cv.glmnet(as.matrix(train[,-1]),as.integer(train[,1])-1,family="binomial",alpha=1) + pred <- round(predict(model, newx=as.matrix(test[,-1]), type='response', s=model$lambda.min),0) + res.l1[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0) + } > sum(res.l1)/nrow(dm) [1] 0.9205703 > sum(diag(table(dw$Result, round(predict(dm.cv.glmnet, newx=as.matrix(dw[,-1]), s=dm.cv.glmnet$lambda.min, type='response'),0))))/nrow(dw) [1] 0.9336283
AIC = 190.04, LOO ACC = 0.921, holdout ACC = 0.934という結果になりました。
3つのモデルのパフォーマンスを比べてみる
単にプロットしてみると以下のような結果が得られます。
全体としてみればもちろんAICが小さければ小さいほどACCは2つとも改善しているのが分かるかと思いますが、LOOとholdoutとで微妙に結果が割れています。これについては@ibaibabaibai先生が本来のWAICとLOOとの比較についてコメントされているのと話が似ているのかなと。
MCMCで事後分布を計算するとき,1個のサンプルの部分の尤度(独立サンプルを過程)の逆数で重みをつけて計算すればleave-oneしたことになるから,いろいろ抜いたのが一回のRUNでできるわけ.しかしこれ,1個抜いてある程度大きく変わる場所を抜いたら不味いことになりそう. https://t.co/uLWMbeAAZy
— baibai (@ibaibabaibai) October 22, 2016
つまり、LOOというのは結局のところ学習データの内部に完結させているものなので、たまたま学習データ自体に変なばらつきがあって毎回抽出するごとに推定されたモデルと真のデータ生成過程とが大きく乖離しているようなケースではLOOは悪くなるけれども、holdoutのような大きめのデータセットに対しては「均される」のでパフォーマンスが良くなるみたいなことはあるのかもしれないな、とはちょっと思ってます。
とはいえ、本当は間を埋めるような微妙なAICの値のモデルを幾つか選んできてサンプルサイズを増やした方が良かったのでしょうが、とりあえずザッとやってみましたという程度なので今回はこの辺にしときます。。。
追記
何と@genkuroki先生から言及いただいてしまいました。以下そのまま貼ります。
Re:RT LOOCV(一個抜き交差検証)と他の規準を比較するときにはそれらが予測分布と真の分布のKullback-Leibler情報量(+定数)の推定値であることを忘れてはいけない。LOOCVがどんなにうまく近似計算できていてもそれだけだと肝腎なことがわからない。続く #数楽
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 続き。確率分布Pによる確率分布Qの予測誤差(シミュレーションの誤差)の大きさをKullback-Leibler情報量D(Q||P)=∫Q(x)log(Q(x)/P(x))dxで測るという考え方はSanovの定理(←決定的に重要)を知っていれば非常に納得できることです。続く
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 続き。だから、LOOCVの近似計算などを含む「予測誤差の大きさを比較するための情報量規準」どうしを正しく比較する方法を理解するためには、予測誤差がKL情報量で測られるということ(実質的にSanovの定理)を知っている必要があります。これが理解のための第一の関門。続く
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 続き。一個抜き交差検証LOOCVは直観的にわかり易い考え方なので、これこそが予測誤差の計測法の決定版だと誤解してしまう危険性があると思う。サンプルから1個抜いて(その1個はサンプル全体を動かす)残りを使って予測分布を作り、抜いた1個の予測精度を見るのがLOOCV。続く
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 続き。サンプルを生成した真の分布は未知のままなので、サンプルだけから予測分布による真の分布の予測精度を推定しなければいけない。サンプル内で予測精度を測ってみる一個抜き交差検証LOOCVは直観的に非常にわかり易いと思う。しかし、LOOCVも予測誤差の推定値でしかない。続く
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 続き。サンプルだけからは知ることができない本当の予測誤差は予測分布と真の分布のKL情報量で測られます。数値実験するときには真の分布をカンニングできるので真の予測誤差を数値計算できます。これとの比較が重要。その辺のことは常識になっているのかなあと思いました。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 予測誤差を測るための各種情報量規準の正しい比較の例が https://t.co/A9N31GtCH6 の注4にあります。ソースも公開で素晴らしい。しかしMatlabには一般人は高価なので手が出ない。RStanで検証するためのコードを誰か公開すると社会貢献になると思う。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 用語説明。未知である真の分布Qが生成したサンプルをもとに作った予測分布Pによる真の分布Qの予測誤差D(Q||P)とQのエントロピーの和(それは-∫Q(x)log P(x) dxに等しい)を汎化誤差(generalization error)と呼ぶようです。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 真の分布が未知であるからこそ、それを推定したいので、自前で数値実験するとき以外には、汎化誤差を推定するときに真の分布をカンニングすることはできません。しかし、サンプルだけから汎化誤差を推定する手段が複数開発されています。最初のAICのアイデアは画期的だったと思う。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 LOOCVも汎化誤差の推定量です。一個抜き交差検証そのものが特権的な重要性を持っているわけではない。交差検証のアイデアはわかりやすいので、この点は誤解し易い点なのかなと思いました。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 過剰学習による予測精度の悪化の問題を扱う場合には、まず、予測分布による真の分布のシミュレーションの誤差を表すKL情報量について学ぶ必要があると思いました。サンプルだけから予測誤差を見積もる方法はKL情報量+定数(=汎化誤差)の推定法になっています。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 私が書いたKullback-Leibler情報量とSanovの定理の易しい場合の解説はこちら→ https://t.co/40sUnHaUCD
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
KL情報量は相対エントロピーの-1倍と同じものです。統計力学を学んだ人は相対エントロピーについては知っていると思う。
#数楽 数値実験する場合には、AIC、DIC、LOOCV(の変種)、WAICなどだけではなく、真の分布をカンニングすることによって計算できる汎化誤差=KL情報量+エントロピーも計算して比較した方がよいと思う。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
#数楽 汎化誤差を推定する方法が複数あることは喜ばしいことだと思う。
— 黒木玄 Gen Kuroki (@genkuroki) 2016年10月27日
実用的に十分な推定が得られる方法の中で計算コストが最も低いものを使うこともできるし、比較によって信頼性を増すことにも利用できる。
こういう考え方なら喧嘩せずにみんなハッピーになれる。