Leave one out法とBootstrap法でモデル選択してみる



正しいモデル選択をするためにモデル評価が大事になってきます。いくつか方法がありますが、代表的なものを列挙します。

* Hold Out法
* Cross Validation法
* Leave One Out法 
* Bootstrap法
* WAIC
* WBIC

leave one outはcross validationにおいて、データ数とグループ数が等しい場合を指します。まずはこの方法を使って、性能評価を行い、正しくモデル選択ができるか検証してみます。

## テストデータの作成
下の因果グラフに基づき、テストデータを生成します。誤差項はすべて、正規分布に従うとします。1回のサンプリング数は10とします。
<img width="300" alt="スクリーンショット 2018-08-18 17.17.21.png" src="https://qiita-image-store.s3.amazonaws.com/0/55421/6fcde0ef-5154-380a-6a08-daa9384f3df1.png">

```R
n= 10
mu = c(0,0)
sd = c(1,1)
cor =.2
cov = sd %*% t(sd) * cor
diag(cov) = sd^2
generate_data = function(){
    mvrnorm(n,mu,cov) %>% as.tibble() %>% mutate(V3 = rnorm(n,0,1),y=3 * V1 + rnorm(n,0,.5))
}
```
作成されたデータの概観をみてみましょう。
![download.png](https://qiita-image-store.s3.amazonaws.com/0/55421/b51f11fc-20e1-6845-c7e5-0aee4aeee5bc.png)

モデルの候補
1. Y ~ V1 (正しいモデル)
1. Y ~ V1+V2
1. Y ~ V1+V2+V3
1. Y ~ (V1+V2)^2 (相互作用も考慮したモデル)
1. Y ~ (V1+V2+V3)^2

## Leave one out法
では、Leave one outに基づいてモデルを選択してみます。

```R
models = list(y~V1, y~V1+V2, y~V1+V2+V3)
result = models %>% map(~{
    model = .x
    1:nrow(df) %>% map(
        ~(lm(model, data=df[-.x,]) %>% predict(df[.x,]) - df$y[.x])^2
    ) 
})
result %>% simplify_all() %>% set_names(str_c("model_",1:3)) %>% as.tibble %>% summarise_all(mean) %>% which.min 
# >>model_1
```
上記のデータセットでは正しくモデル1が選ばれたようです。
次に上記の流れ(データ生成+モデル選択)を1000回繰り返します。果たして正しいモデルを選択できたのでしょうか?

```R
models = list(y~V1, y~V1+V2, y~V1+V2+V3, y~(V1+V2)^2, y~(V1+V2+V3)^2)
1:100 %>% map(~{
    df = generate_data()
    result = models %>% map(~{
        model = .x
        1:nrow(df) %>% map(
            ~(lm(model, data=df[-.x,]) %>% predict(df[.x,]) - df$y[.x])^2
        ) 
    })
    result %>% simplify_all() %>% set_names(1:length(models)) %>% as.tibble %>% summarise_all(mean) %>% which.min 
})

result %>% unlist %>% as_tibble %>% 
ggplot(aes(value)) + geom_bar() + labs(x="モデル",y="カウント数")
```
さて、気になる結果です。一番選択されやすいモデルは正しいモデル①であったものの、それ以外のモデルも少なからず選択されています。

| Model| Count|
|-----:|-----:|
|     1|   639|
|     2|   130|
|     3|   101|
|     4|    96|
|     5|    34|

![download-1.png](https://qiita-image-store.s3.amazonaws.com/0/55421/ea06f0a7-1787-577e-ea0b-e0e4fb261ba5.png)

次にサンプルサイズを100に増やして、同じように実験してみます。
結果はこちら。

| Model| Count|
|-----:|-----:|
|     1|   701|
|     2|    95|
|     3|    89|
|     4|    76|
|     5|    39|

![download-1のコピー.png](https://qiita-image-store.s3.amazonaws.com/0/55421/0c7447fe-b54f-a0c8-6506-3b0701fc15c4.png)

想像どおり、サンプル数を増やすことによって、より正しくモデル選択ができるようになりました。
しかし、それでも正答率は70%ほど。データのみから正しいモデルを選択することはかなり困難な印象です。Leave one out以外にWAICも汎化誤差の近似として使いますが、性能は似たようなものらしいです(うろ覚え)。
今回は正規分布かつ線形モデルでテストデータを作成しました。
しかし、現実は背景因子の確率分布も複雑であり、因果関係も線形構造ではないことが多いでしょう。
やはり、その道の専門家と手を取り合い、少しでも多くのドメイン知識に基づいて、モデルを作成していく&候補を絞っていく作業が重要になってくるのではないでしょうか。

## Bootstrap法(追記)

1つの標本からリサンプリングを繰り返して生成する標本をブートストラップ標本と呼びます。
手元のデータから最尤推定で求めたパラメタでMSEを計算しても、真のデータ分布を使ったMSEはそれより大きくなる。それをブートストラップ標本を使って修正する方法です。詳しくは「はじパタ第2章」参照してください。

```R
result = 1:100 %>% map(~{
    df = generate_data()
    models %>% map(~{
        model=.x
        expected_bias = 1:100 %>% map_dbl(~{
            sample = df %>% sample_frac(replace=T) 
            fit = lm(model,data=sample)
            fitted_MSE = fit$residuals^2 %>% mean 
            predicted_MSE = (predict(fit, df) - df$y)^2 %>% mean
            bias = predicted_MSE-fitted_MSE
            bias
        }) %>% mean
        lm(model,df)$residuals %>% .^2 %>% mean + expected_bias
    }) %>% which.min   
})
```
気になる結果は。。。

| Model| Count|
|-----:|-----:|
|     1|   622|
|     2|   164|
|     3|   134|
|     4|    72|
|     5|     8|

![download-1のコピー3.png](https://qiita-image-store.s3.amazonaws.com/0/55421/70951914-26af-57ff-d7aa-1ac033eb7b3d.png)

今回のテストデータだけからは結論付けれませんが、LOOVと比較すると性能が悪いような気がします。