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

Rとブートストラップ

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

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


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

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

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

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

パラメトリック・ブートストラップ法
(1)標本サイズがnである標本データ{x1,x2,…,xi,…,xn}の平均平均x、標準偏差sを計算する。
(2)n個の正規乱数z1,z2,…,zi,…,znを生成し、xi*=平均x+zisで新しい標本{x1*,x2*,…,xi*,…,xn*}を生成する。この標本による推定値を推定値 (例えば、平均平均x1)とする。
ノンパラメトリック・ブートストラップ法
(1)区間(0,1)をn等分した各区間の値を標本データ{x1,x2,…,xi,…,xn}に1対1で対応させる。
(2)n個の一様乱数u1,u2,…,ui,…,unを生成し、uiの値が含まれる区間に対応するxkxi*とし、新しい標本データ{x1*,x2*,…,xi*,…,xn*}を生成する。この標本から得られた推定値を推定値とする。

 両方法ともステップ(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={y1,y2,…,yn}、説明変数をX={x1,x2,…,xm}としたブートストラップ法による回帰分析の例を紹介する。そのアルゴリズムを次に示す。

(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)
summary(car.boot)の結果

 パッケージ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

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