[連載] フリーソフトによるデータ解析・マイニング 第44回

Rとブートストラップ

1.ブートストラップとは

  統計学の主な目的の1つは、標本データを用いて母集団の性質を推測することである。同じ母集団から抽出した標本であっても、無作為であるため標本を構成する要素、標本のサイズが異なると、それらの統計量(比率、平均、分散など)は異なる。従って、標本データを用いて母集団の性質を推測する際には常に誤差が伴う。
  正規分布N (μ , σ2) の母集団から抽出した大きさn の無作為標本の平均はN ( μ , σ2 / n ) に従うことが知られている。σ は一定の条件のもとでは標本の不偏標準偏差を用いることも可能である。このように正規分布、t 分布、x 2 分布などの確率分布を用いて母数やモデルの推定およびその推定の誤差を計算することができる。しかし、問題によっては確率分布を仮定できないケースも少なくない。そこで、1970年代にエフロン (Efron) は確率分布の性質に頼らないブートストラップ (bootstrap) という方法を提唱した。ブートストラップの語源に関しては、インターネットでも検索できる。
  データサイエンスの分野では、1つの標本から復元抽出を繰り返して大量の標本を生成し、それらの標本から推定値 を計算し、母集団の性質やモデルの推測の誤差などを分析する方法をブートストラップ法と呼ぶ。
  ブートストラップ法では母数 の推定量は標本から生成したブートストラップ標本の推定量 を用いて推定する。1つの標本からリサンプリングを繰り返して生成する標本をブートストラップ標本と呼ぶ。図1にブートストラップ法のイメージを示す。

図1 ブートストラップ法のイメージ

図1 ブートストラップ法のイメージ

2.ブートストラップ標本の生成

  ブートストラップ標本の生成には幾つかの方法が提案されているが、確率分布型を仮定するパラメトリック・ブートストラップ法と確率分布型を仮定しないノンパラメトリック・ブートストラップ法に大別される。そのアルゴリズムの例を次に示す。

 (1) パラメトリック・ブートストラップ法

 @ 標本サイズがnである標本データ { x 1 , x 2 , … , x i , … , x n } の平均 平均x、標準偏差s を計算する。
 A n 個の正規乱数 z 1 , z 2 , … , z i , … , z n を生成し、xi*=平均x+zis で新しい標本 { x 1* , x 2* , … , x i* , … , x n* } を生成する。この標本による推定値を推定値 (例えば、平均 平均x1 ) とする。

 (2) ノンパラメトリック・ブートストラップ法

 @ 区間 (0,1) をn 等分した各区間の値を標本データ { x 1 , x 2 , … , x i , … , x n } に1対1で対応させる。
 A n 個の一様乱数 u 1 , u 2 , … , u i , … , u n を生成し、u i の値が含まれる区間に対応する x kx i* とし、新しい標本データ { x 1* , x 2* , … , x i* , … , x n* } を生成する。この標本から得られた推定値を 推定値 とする。

  両方法ともステップ (2) を B 回繰り返し、B 個の標本の推定値 推定値 を求める。その推定値、標準偏差、バイアスはそれぞれ次の式で求める。

推定値、標準偏差、バイアスそれぞれ次の式

  また、確率分布関数は

確率分布関数
により推定できる。B 個の推測値を大小順に並べたB ×α 番目の値を100α%点とする。
繰り返しの回数B については、推定値の標準誤差を求める場合は約100〜200回、確率分布関数の値や100α%点を求める場合は1000〜2000回が必要であるとされている (小西 (1988, 2004))

3.区間推定

  ブートストラップ標本を生成し、平均の信頼区間を推定する例を次に示す。用いるデータは平均170、標準偏差6である正規分布の母集団から抽出したサイズが20の標本であるとする。標本データ生成のコマンドを示す。

> set.seed(20)
> sam<-rnorm(20,170,6)
> round(sam[1:5],2)
[1] 165.75 166.71 170.16 169.54 165.03

  パラメトリック・ブートストラップ法による2000個のブートストラップ標本を生成し、その平均を求めるコマンドを次に示す。ただし、標準偏差は標本の標準偏差を用いる。

> tt<-numeric(0)
> ME<-mean(sam); SD<-sd(sam)
> for(i in 1:2000){z<-rnorm(20,0,1);bx<-ME+z*SD; tt<-cbind(tt,mean(bx))}

  求めた2000個の標本平均の平均と95%の信頼区間を100 %点で求める例を示す。100 %点は、関数 quantile を用いて求めることができる。

> mean(tt)
[1] 170.2390
> quantile(tt,p=c(0.025,0.975))
  2.5%  97.5%
168.1689 172.3980

  ノンパラメトリックのブートストラップ標本は、関数 sample を用いて生成することができる。関数 sample を用いて生成したブートストラップ標本を用いて区間推定を行うコマンドを次に示す。

> tt<-numeric(0)
> for(i in 1:2000){bs<- sample(sam,20,replace = TRUE);me<-mean(bs);tt<-cbind(tt,me) }
> mean(tt)
[1] 170.1977
> quantile(tt,p=c(0.025,0.975))
  2.5%  97.5%
168.1203 172.3522

  Rにはブートストラップに関連するパッケージ boot, simpleboot, bootstrap などがある。いずれも CRAN ミラーサイトからダウンロードできる。
  パッケージ boot の中にはブートストラップの推定量を求める関数 boot があるが引数が若干多いので、ここではパッケージ simpleboot の中の関数 one.boot を用いることにする。関数 one.boot の書き式を次に示す。ただし、引数に関しては最も基本的な項目のみを示す。

one.boot(data, FUN, R, ...)

  引数 data には用いる標本データ、引数FUNには推定する統計量の関数 (mean, median, あるいは自作関数) を指定する。引数Rには生成するブートストラップ標本の数を指定する。デフォルトの値のままブートストラップの推定値のベクトル 推定値のベクトル を生成するためには、これらの3つの引数のみを指定すればよい。
  標本データ sam を用いた10個のブートストラップ標本の平均値を求める例を次に示す。ブートストラップの平均ベクトルは$tに格納されている。生成されたブートストラップの統計量はノンパラメトリック法によるものである。

>library(simpleboot); set.seed(20)
> b.mean <- one.boot(sam, mean, 10)
> b.mean$t
    [,1]
[1,] 170.4848
<中略>
[10,] 169.1167

  パッケージ simpleboot の関数が返す結果のオブジェクト形式は、パッケージ boot を用いた結果と一致する。パッケージ boot には信頼区間を求める関数 boot.ci がある。関数 boot.ci の書き式を次に示す。ただし、引数に関しては最も基本的な項目のみを示す。

boot.ci(boot.out, conf = 0.95, type = "all", ...)

  引数 boot.out はパッケージ boot、あるいは simpleboot で生成したブートストラップの結果オブジェクトである。引数 conf は信頼区間であり、デフォルトには95%の信頼区間が指定されている。引数 type は信頼区間を求める方法であり、正規分布の近似 (normal) 法、basic 法、ブートストラップ t (studentized) 法、パーセンタイル (percentile) 法、BCa 法の5種類の方法を指定することができる。ただし、ブートストラップ t 法に関しては、ブートストラップ標本の推定値を生成する際にブートストラップ t 法を用いなければならない。関数 boot、one.boot には関連の引数がある。これらのアルゴリズムに関しては汪(2006)が詳しい。
  関数の boot.ci の引数 type をデフォルト値 "all" のままで実行するとブートストラップ t 法を除いた4種類の信頼区間を返す。
  関数 boot.ci に実装された5種類の区間推定法の中の正規分布の近似法を除いた4種類の境界値の計算法を表1に示す。

表1 4種類の区間推定の計算式
方法 有意水準 α の点
 Basic  確率分布関数
 Studentized  確率分布関数
 Percentile  確率分布関数
 BCa  確率分布関数

  表の中の θ は標本データの推定値、Fb はブートストラップの累積分布関数、z α はブートストラップ t の累積分布関数、Fs は有意水準 α の標準正規分布の値、a , c は定数である。定数 ac の決め方に関しては汪・田栗 (2003), 丹後 (2000) に紹介されている。ブートストラップの標本の数を2000とした標本平均の信頼区間を、関数 boot.ci を用いて求める例を次に示す。

> library(boot)
> b.mean <- one.boot(sam, mean, 2000)
> boot.ci(b.mean)

 boot.ci(b.mean)の結果

  返された結果 Intervals の下部の Normal は通常よく用いられている正規分布を近似した区間推定の結果である。Percentile の結果は推定値を小さいものから大きい順に並べた数列の100α%点である。
  異なる計算方法であるのにもかかわらず4種類の推定値が非常によく近似している。乱数を用いる計算結果は、同じ方法を繰り返してもその結果が同じになるとは限らない。勿論、同一の乱数シードを用いるなどの工夫を行うことで再現することは可能である。
  より安定した区間推定値を得るためにはブートストラップ標本の数を増やすことが1つの方法である。上記の例のような計算量であればブートストラップ標本を1万にしてもRでは直ちに計算結果が返される。

4.回帰分析とブートストラップ

  ブートストラップ法は多くの統計データ処理に用いられている。本項では目的変数を Y = { y 1 , y 2 , … , y n }、説明変数を X = { x 1 , x 2 , … , x m } としたブートストラップ法による回帰分析の例を紹介する。そのアルゴリズムを次に示す。

 (1) 観測の標本データを用いて回帰モデル Y=R(X) を作成する。
 (2) 残差 E=Y-Y を用いてノンパラメトリック・ブートストラップ標本 E * = { ε1*, ε2*, … , εn* } を生成する。
 (3) 回帰モデルの予測値に残差のブートストラップ標本を加え、新しい目的変数のブートストラップ標本 Y*=Y+E* を作成する。
 (4) X ,Y * を用いて回帰モデルを作成する。
 (5) ステップ(2)〜(4)をB 回繰り返す。

  回帰係数はB回のブートストラップ標本の回帰係数の平均、回帰係数の推定誤差はブートストラップ標本の回帰係数の標準偏差を用いて求める。
  Rの中のデータ cars を用いた例を次に示す。

> data(cars)
> car.lm<-lm(dist~speed,data=cars)
> car.boot<-lm.boot(car.lm, R = 2000)
> summary(car.boot)

BOOTSTRAP OF LINEAR MODEL (method = rows)
Original Model Fit
------------------

Call:
lm(formula = dist ~ speed, data = cars)

Coefficients:

(Intercept) speed
-17.579 3.932

Bootstrap SD's:

(Intercept) speed
5.8648169 0.4192157

  パッケージ simpleboot には関数 lm.boot のブートストラップ単回帰結果の作図関数 plot.lm.simpleboot があり、plot に略して用いる。この関数は散布図の回帰直線にブートストラップの2倍の標準誤差の区間を表示する。その例を次に示す。

> plot(car.boot,xlab="speed", ylab="dist")

図2 lm.bootの回帰図

図2 lm.boot の回帰図

  パッケージ simpleboot には、局所多項式回帰関数 loess による回帰結果についてブートストラップ回帰を行う関数 loess.boot がある。データ cars を用いた2000回のブートストラップの回帰結果のグラフを次に示す。

図3 関数loess.bootの回帰図

図3 関数 loess.boot の回帰図

  パッケージ boot, bootstrap にもブートストラップ回帰関数が用意されている。

5.クラスター分析とブートストラップ

  階層的クラスター分析で得られた樹形図がデータの構造を表す真の樹形図であるどうかは判断しがたい。従って、樹形図を用いたクラスターの結果を解釈する際には主観的になりがちである。
  そこで、得られた樹形図が真の樹形図であるかどうかをブートストラップ法によって評価する研究が行われている。
  クラスター分析におけるブートストラップ法の適応は、ブートストラップ標本が仮説を支持する相対頻度 (ブートストラップ確率) を仮説検定の確率値として用いる。下平 (2002) は、標本データの変数の次元と異なるブートストラップ標本を生成するマルチスケールブートストラップ法によるブートストラップ確率を用いて樹形図を評価する方法を提案した。そのアルゴリズムを含む解説論文の PDF ファイルがインターネット上で入手できるので、その説明は省略する。その理論に基づいたパッケージ pvclust は CRAN ミラーサイトからダウンロードできる。クラスターの生成とブートストラップ確率を計算する関数は pvclust である。その書き式を次に示す。

pvclust(data, method.hclust="average",method.dist="correlation",use.cor="pairwise.complete.obs",nboot=1000,…)

  引数 method.hclust は、関数 hclust の引数 method に対応する。引数 method.dist では距離関数 dist の引数および相関係数のような類似度関数を指定する。引数 use.cor は、相関係数行列を求める関数 cor の引数 use に相当する。引数 nboot はブートストラップ標本の数であり、1000以上が推奨されているが、データサイズが大きい場合は計算結果が返されるのを待つのに辛抱を要する。
  データ iris の中の異なる2種類のデータの一部分を用いた例を次に示す。ここではキャンベラ距離を用いる

> y<-t(iris[c(51:60,141:150),1:4])
> library(pvclust)
> irisy.pvwar <- pvclust(y, method.dist="can",nboot=5000, method.hclust="ward",r=seq(.5,1,by=.1))
> plot(irisy.pvwar,hang=-1,cex.pv=1, col.pv=c(4,2,8),float=0.02 )

図4 ブートストラップ確率と樹形図1

図4 ブートストラップ確率と樹形図1

  この結果の場合はブートストラップ確率 (ノードの右側の数値 bp) が5%以下のものはない。2つのクラスターに分ける際のそれぞれのクラスターの最上部のノードのブートストラップ確率は50%を超えている。都合よくこの2つのクラスターは真のクラスターと一致する。パッケージの中には確率値を用いてクラスターを見つける関数 pvpick がある。

> pvpick(irisy.pvwar, alpha=0.50, pv="bp")
<結果は省略>

  参考のため、同じキャンベラ距離に最遠隣法を用いたブートストラップ標本数を5000とした結果を図5に示す。このように、同じデータと距離に対して、異なるクラスター方法を用いると、そのブートストラップ確率も大きく異なる。

図5 ブートストラップ確率と樹形図2

図5 ブートストラップ確率と樹形図2

参考文献
[1] 小西 貞則(1988). ブートストラップ法による推定量の誤差評価、パソコンによるデータ解析(村上・田村編)、p.123-142, 朝倉書店
[2] 小西 貞則, 北川 源四郎(2004). 情報量基準、朝倉書店
[3] 丹後 俊郎(2000). 統計モデル入門, 朝倉書店
[4] 下平 英寿(2002). ブートストラップ法によるクラスター分析のバラツキ評価、統計数理、第50巻第1号、p.33-44. (http://www.ism.ac.jp/editsec/toukei/pdf/50-1-033.pdf)
[5] 汪 金芳・田栗 正章(2003). ブートストラップ法入門, 「計算統計T:統計科学のフロンティア11」, 岩波書