R言語第八講評価モデルの交差検証法分析事例
タイトル
複数の線形モデルに対してAutoデータセットでフィットしたテストエラー率を評価します.Autoデータセットは、ISLRパッケージの1つであるオートバイ関連のデータが存在するデータセットであり、読者は自分でISLRパッケージをダウンロードし、Autoデータセットをロードすることができる.
関連資料
クロス検証は,機械学習においてモデルの構築とモデルパラメータの検証によく用いられる方法である.クロス検証は、その名の通り、データを繰り返し使用し、得られたサンプルデータを分割し、異なる訓練セットとテストセットに組み合わせ、訓練セットでモデルを訓練し、テストセットでモデル予測の良し悪しを評価することである.これに基づいて、複数の異なるトレーニングセットとテストセットが得られ、あるトレーニングセットのあるサンプルが次回、テストセットのサンプル、すなわち「交差」になる可能性がある.
では、クロス検証はいつ必要ですか?クロス検証はデータが十分でない場合に使用されます.例えば、私の日常プロジェクトでは、一般的な適当な問題に対して、データサンプル量が1万本未満であれば、クロス検証を採用して最適化選択モデルを訓練します.サンプルが1万個より大きい場合、私たちは一般的にランダムにデータを3つに分け、1つはトレーニングセット(Training Set)、1つは検証セット(Validation Set)、最後の1つはテストセット(Test Set)に分けます.モデルを訓練セットで訓練し,モデル予測の良し悪しと選択モデルとその対応するパラメータを検証セットで評価した.最終的に得られたモデルをテストセットに再利用し,最終的にどのモデルと対応パラメータを使用するかを決定する.
クロス検証に戻ります.分割方法によって、クロス検証は次の3つに分けられます.
1つ目は単純交差検証であり,いわゆる単純であり,他の交差検証方法とは対照的である.まず、サンプルデータをランダムに2つの部分(例えば、70%のトレーニングセット、30%のテストセット)に分け、トレーニングセットでモデルをトレーニングし、テストセットでモデルとパラメータを検証します.次に,サンプルを乱し,トレーニングセットとテストセットを再選択し,トレーニングデータと検証モデルを継続した.最後に,損失関数を選択して最適モデルとパラメータを評価した.
2つ目はSターンクロス検証(S-Folder Cross Validation)です.1つ目の方法とは異なり,Sターンクロス検証ではサンプルデータをランダムにS部に分け,毎回ランダムにS−1部をトレーニングセットとして選択し,残りの1部をテストセットとする.このラウンドが完了すると、S-1部をランダムに選択してデータを訓練します.いくつかのホイール(Sより小さい)の後,最適モデルとパラメータを評価するために損失関数を選択した.
3つ目は、第2のケースの特例であるLeave-one-out Cross Validation(Leave-one-out Cross Validation)を残すことであり、このときSはサンプル数Nに等しく、N個のサンプルに対して、毎回N-1個のサンプルを選択してデータを訓練し、モデル予測の良し悪しを検証するためにサンプルを残す.この方法は主にサンプル量が非常に少ない場合に用いられる.例えば、通常の適当な問題に対して、Nが50未満の場合、私は一般的にクロス検証を採用する.
得られたモデルの良し悪しを損失関数で測定する繰り返し交差検証により,最終的にはより良いモデルを得ることができる.では、この3つの状況では、いったい私たちはどの方法を選ぶべきなのでしょうか.一言でまとめると、データに対して初歩的なモデルを構築するだけで、深く分析するのではなく、簡単に交差して検証すればいいのです.そうでない場合はS折れ交差で検証します.サンプル量が少ない場合は、Sターンクロス検証の特例を用いてクロス検証を残します.
また,比較的特殊なクロス検証方式もあり,サンプル量が少ない場合に用いられる.セルフサービス法(bootstrapping)と呼ばれています.例えば、m個のサンプル(mが小さい)があり、このm個のサンプルの中でランダムに1個のサンプルを採取し、訓練セットに入れ、サンプリングが終わったらサンプルを戻します.このようにm回繰り返し採取し,m個の試料からなる訓練セットを得た.もちろん,このmサンプルには重複するサンプルデータがある可能性が高い.また、サンプリングされていないサンプルをテストセットとして使用します.これにより、クロス検証が行われる.我々のトレーニングセットには重複データがあるため、データの分布が変化するため、トレーニング結果には推定偏差があるため、データ量が本当に少なく、例えば20個未満でない限り、この方法はあまり一般的ではありません.
ランダムサンプリング関数sampleの説明:
> x=1:10
> sample(x=x)
[1] 3 5 9 6 10 7 2 1 8 4第1行のコードはxベクトルxに値を付与し、第2行のコードはxベクトルをランダムにサンプリングすることを示す.結果出力はサンプリング毎に抽出された結果であり,このサンプリングは無戻しサンプリング−最大n回,nはベクトル中の要素の個数であることが分かった.
このベクトルで抽出する要素の数を指定するには、パラメータsizeを追加する必要があります.
> x=1:1000
> sample(x=x,size=20)
[1] 66 891 606 924 871 374 879 573 284 305 914 792 398 497 721 897 324 437
[19] 901 33これは1~1000の正の整数でサンプリングされ、sizeがサンプリングの回数を指定して20回抽出したもので、結果は上記の通りである.
これらはすべて無戻しサンプリングです.無戻しサンプリングとは、ある要素が選択されると、その全体にその要素は存在しないということである.戻りサンプルがある場合は、パラメータrepalce=Tを追加する必要があります.
> x=1:10
> sample(x=x,size=5,replace=T)
[1] 4 7 2 4 8「replace」は繰り返しの意味です.すなわち、元素のサンプリングを繰り返すことができ、いわゆるリターンサンプリングがある.上記の結果を見て,元素4は5回のランダムサンプリングの過程で2回抽出された.
コードの位置が関数内のパラメータの位置に対応している場合は、次のような関数のパラメータを書かなくてもいいです.
> x=1:10
> sample(x,20,T)
[1] 1 2 2 1 5 5 5 9 9 5 2 9 8 3 4 8 8 8 1 1簡単な運用:サイコロを投げたり、コインを投げたり(これはサンプリングを紹介する上で必ず紹介する内容かもしれません)、すべてサンプリングに戻したりします.
ここで、sample関数の場合、パラメータxは数値であっても文字であってもよく、実際にはパラメータxは任意のベクトルを表す.
> a=c("A","B")
> sample(x=a,size=10,replace=T)
[1] "B" "A" "A" "A" "B" "A" "A" "B" "A" "A"上記のコードはコインを投げて10回投げたと理解でき、そのうち正面(A)と裏面(B)が現れる回数は繰り返すことができる.
上記のサンプリングプロセスは、各要素が抽出される確率が等しく、ランダムサンプリングと呼ばれる.
抽出要素の確率が必ずしも等しいとは限らない場合があります(一般的な2つの分布確率問題など)、パラメータprob、すなわち「probability」(確率)の略語を追加する必要があります.ある医者が患者にある手術に成功する確率が80%だと仮定すると、今彼が20例の患者に手術をしたのは、何回成功したのだろうか.コードは次のとおりです.
> x=c("S","F")
> sample(x,size=20,replace=T,prob=c(0.8,0.2))
[1] "F" "S" "S" "S" "S" "S" "S" "S" "S" "S" "S" "S" "F" "S" "S" "F" "S" "S"
[19] "F" "S"ここで、「S」は成功を表し、「F」は失敗を表す.
> x=c(1,3,5,7)
> sample(x,size=20,replace=T,prob=c(0.1,0.2,0.3,0.9))
[1] 3 5 7 3 7 3 7 5 3 7 7 7 1 5 7 5 7 7 3 7これらのコードは、各要素に対して1つの確率を与えることができ、各確率は独立していることを示しています.すなわち、パラメータprobでは、必ずしもすべての要素の確率を加算して1に等しくなく、ある要素が抽出される確率だけを表しています.
sample関数の場合、そのパラメータxは、上述した文字のサンプリングのようなR内の任意のオブジェクトであってもよい.同じ機能の関数sampleもあります.int,「int」は「intger」の略,すなわち「整数」である.パラメータnは正の整数でなければなりません.
> x=-10.5:7.5
> sample(x=x,size=3)
> sample.int(n=x,size=3)
[1] -5.5 -7.5 0.5
Error in sample.int(x, size = 3) : invalid first argument第1行のコードは-10.5から7.5の等差数列を生成し、結果出力の第1行はsampleの結果である.2行目はsampleです.intの結果、正の整数ではないため、「最初の引数は無効です」というエラーが表示されます.残りの使い方はampleと同じです.
じっけん
始める前にsetを使います.seed()関数は、Rの乱数生成器にシードを設定.これにより,以下に示すものと全く同じ結果が得られる.一般に,クロス検証法のようにランダム性を含む解析法を用いる場合,ランダムシードを設定することで,次回と全く同じ結果を得ることができる.
まずsample()関数で観測セットを半分に分け,元の392個の観測から196個の観測のあるサブセットをランダムに選択し,訓練セットとした.
> library(ISLR)
> set.seed(1)
# 1-392 196 。
> train=sample(392,196)
[1] 324 167 129 299 270 187 307 85 277 362 330 263 329 79 213 37 105
[18] 217 366 165 290 383 89 289 340 326 382 42 111 20 44 343 70 121
[35] 40 172 25 248 198 39 298 280 160 14 130 45 22 206 230 193 104
[52] 367 255 341 342 103 331 13 296 375 176 279 110 84 29 141 252 221
[69] 108 304 33 347 149 287 102 145 118 323 107 64 224 337 51 325 372
[86] 138 390 389 282 143 285 170 48 204 295 24 181 214 225 163 43 1
[103] 328 78 284 116 233 61 86 374 49 242 246 247 239 219 135 364 363
[120] 310 53 348 65 376 124 77 218 98 194 19 31 174 237 75 16 358
[137] 9 50 92 122 152 386 207 244 229 350 355 391 223 373 309 140 126
[154] 349 344 319 258 15 271 388 195 201 318 17 212 127 133 41 384 392
[171] 159 117 72 36 315 294 157 378 313 306 272 106 185 88 281 228 238
[188] 368 80 30 93 234 220 240 369 164lm()関数のsubsetオプションを用いて,訓練セットの観測だけで線形回帰モデルをフィッティングした.
> lm.fit=lm(mpg~horsepower,data=Auto,subset=train)ここでpredict()関数を用いて全392観測の応答変数を推定し,検証セット196観測の平均二乗誤差をmean()関数を用いて計算した.なお、以下の-train指標は、訓練集中していない観測のみを選択することを意味する.(説明:'-'はsubset)
> attach(Auto)
> mean((mpg-predict(lm.fit,Auto))[-train]^2)
[1] 23.26601以上のデータから,線形回帰フィッティングモデルを用いて生成した試験の平均二乗誤差推定は23.26601であることが分かった.次にpoly()関数を用いて,二次および三次多項式回帰による試験誤差を推定する.
> lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)#poly
> mean((mpg-predict(lm.fit2,Auto))[-train]^2)
[1] 18.71646
> lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)#
> mean((mpg-predict(lm.fit3,Auto))[-train]^2)
[1] 18.79401この2つの誤り率はそれぞれ18.71646と18.79401であった.異なるトレーニングセットを選択すると、検証セットに異なる誤差が得られます.
> set.seed(2)
> train=sample(392,196)
> lm.fit=lm(mpg~horsepower,subset=train)
> mean((mpg-predict(lm.fit,Auto))[-train]^2)
[1] 25.72651
> lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
> mean((mpg-predict(lm.fit2,Auto))[-train]^2)
[1] 20.43036
> lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
> mean((mpg-predict(lm.fit3,Auto))[-train]^2)
[1] 20.38533ケーススタディが終わり、仲間たちは学んだかな(*^^*)