ブートストラップ法とそのロバスト推定への応用
M2011MM051 室梅秀平
指導教員：木村美善

1

はじめに2.2

ロバスト推定量とは, 観測したデータに外れ値が混入していることが予測される場合に有効な推定量である. 通常, 最小二乗法などの一般の推定量においては, 誤差は互
いに独立で同一分布に従う対称な確率変数であると仮定して推定を行なう. しかし, その仮定に従わない観測値がデータに含まれていた場合, これらの推定量は良い推定量である根拠を失ってしまう. このように, 仮定したモデルから逸脱する観測値がデータに含まれていた場合でも,その影響を抑え, 良い推定が得られることを目的とした推定量がロバスト推定量である.
ブートストラップ法とは, 観測された限られたデータを元に新たな標本をリサンプリングする手法である. 新しい標本をいくらでもリサンプリングできるので, それらを用いて推定量の偏りや分散, 分布などの統計的性質を大量のデータを用いた反復計算に置き換えて推定することが可能である. 計算機が必要となるが, 複雑な数式を知らなくてもモンテカルロ式に推定を行なえることが大きな魅力である.
ロバスト推定とブートストラップ法はどちらも優れた手法であるが, この二つの手法は相性が良くなく, そのまま組み合わせて使用するといくつかの問題が発生する. これらの問題について考察を行い, その改善案を模索, 検証することが本研究の主な目的である.
また, この問題に対する既知の改善案である SalibianBarrera and Zamar [2] の Fast Robust Bootstrap(FRB)法についても考察を行なう.

2

定義

2.1

線形回帰モデルロバスト・ブートストラップ法を適用するロバスト回帰推定量として, MM 回帰推定量 (Yohai [7]) を用いる.
MM 回帰推定量は, まず尺度推定値 σ̂n の推定を行い, そMMβ n を推定する. それぞの推定値を用いて回帰推定値 β̂れの推定の式には損失関数 ρ0 , ρ1 が組み込まれ, この関数が外れ値の重みを低く抑えることでロバストな推定量になっている.
β ) の内, 最尺度の推定値 σ̂n は次の式 (3) を満たす σ̂n (β小のものである.
)(0n1∑yi − xiβρ0= b.
(3)β)
nσ̂n (βi=1β ).
σ̂n = min σ̂n (β(yi , x1i , . . . , xpi ) , i = 1, . . . , n は一般的な分布 H に従う
0独立な確率ベクトルとする. また, xi = (1, x1i , . . . , xpi ) ∈Rp+1 とする. xi によって yi を説明する次のような線形回帰モデルを考える.
0yi = xiβ + σεi ,i = 1, ..., n.

(1)

β = (β0 , β1 , . . . , βp ) は回帰パラメータ, εi ∼ N (0, 1) は誤差, σ はその尺度である.
外れ値の発生やモデル分布からのずれを考慮に入れるため, 分布 H は興味がある分布 H0 の -汚染近傍 H に属すると仮定する.
H = {H = (1 − )H0 + H ∗ |H ∗ ∈ M}.

(2)

 は 0 ≤  ≤ 1/2 であり, M は確率分布の全体からなる集合である.

(4)βb は EΦ [ρ0 ] = b(Φ は標準正規分布) となるような定数β ) となるような β を Sである. なお, このとき σ̂n = σ̂n (βSβ n とする.
回帰推定値と呼び, β̂MMβn
次に, σ̂n を用いて回帰推定値 β̂を求める. MM 回MMβ n は次の式 (5) を満たすものである.
帰推定値 β̂1∑ 0ρn i=1 1n(0MM )βnyi − xiβ̂σ̂nxi = 0.

(5)損失関数 ρ0 , ρ1 は連続微分可能な偶関数であり, ρ0 の取り方は漸近効率を, ρ1 の取り方は破綻点をそれぞれ左右する (Huber [1] を参照).

2.30MM 回帰推定量

ブートストラップ法

ここでは, 経験分布を用いたノンパラメトリックなブートストラップ法について定義する. 観測したデータに基づく経験分布を H1 とする. この経験分布 H1 とは, n 個の観測値の各点で確率 1/n を持つ分布である. H1 から無作為復元抽出された大きさ n の標本をブートストラップ0標本と呼び, (yi∗ , x∗1i , . . . , x∗pi ) , i = 1, ..., n と表す. また,x∗i = (1, x∗1i , . . . , x∗pi )0 , i = 1, ..., n とする.
ブートストラップ標本を元に, 線形モデルの係数 β を∗βn と推定した推定量をブートストラップ推定量と呼び β̂∗β n は次の式を満た表す. 前述の MM 回帰推定量の場合, β̂す値となる.
)(0 ∗nβn1 ∑ 0 yi∗ − x∗i β̂x∗i = 0.
(6)ρ1nσ̂n∗i=1同じく, 尺度のブートストラップ推定量 σ̂n∗ は, 次の式をβ ) の中で最小のものである.
満たす σ̂n∗ (β)(0n1∑yi∗ − x∗i β= b.
(7)ρ0β)
nσ̂n∗ (βi=1β ).
σ̂n∗ = min σ̂n∗ (β(8)βブートストラップ標本は既知の分布 H1 から無作為に生成されるため, 何度でも取り直すことが出来る. そのため,ブートストラップ標本を大量に生成してブートストラップ推定量を大量に用意すると, モンテカルロ法の要領で推定量の分布を推定することも出来る (汪他 [6] を参照).

3

ロバスト推定量への適用の際の問題点

このようなブートストラップ法をロバスト推定量に適用する場合, Salibian-Barrera and Zamar [2] によると次のような問題が発生する.

3.1

計算コスト

ブートストラップ法では, 何千何万ものブートストラップ標本を生成し, そのそれぞれに対して推定量の計算を行なう. しかし, 表 1 にあるようにロバスト推定量の計算コストは最小二乗法などと比べて非常に大きく, 大量に反復計算を行なうと計算コストが膨大なものとなってしまう.

MMβn
図 1 外れ値が混入したデータへの MM 推定量 β̂4
表 1 サンプル数 n, 次元数 p の場合の主なロバスト推定量の計算時間 (秒)(n, p)(50, 1)(100, 1)(50, 3)(100, 3)LSMSMMτ

0.00008
0.032
0.42
0.45
0.45

0.00009
0.055
0.92
0.97
0.94

0.00008
0.11
1.04
1.15
1.31

0.00009
0.21
2.21
2.43
1.88

問題点への改善案

これらの問題点に対する改善案を模索したところ, ロバ
スト推定量のブートストラップ推定量を算出する際の計算アルゴリズムを一部簡略化することで問題点を改善できることが分かった.
ロバスト推定量を計算する際の一般的な計算アルゴリ
ズムは, Salibian-Barrera, Willems and Zamar [4] によると, 大まかには次のようになる.

β n の候補となる β を N 個用意する.
1. 解となる β̂
2. 用意した候補のそれぞれに対し, その候補の近くにある局所最適解を探し出す.

3.2

ブートストラップ標本の外れ値の割合

ブートストラップ標本は, 元の標本の経験分布から無作為復元抽出されたものである. ロバスト推定を適用する
からには元の標本には外れ値が含まれていることが予想されるが, 元の標本に含まれるどの観測値も平等に抽出され得るため, 外れ値も正常な値もどちらも平等に選ばれ得る. また, 復元抽出であるため 1 つのブートストラップ標
本に同じ外れ値が複数抽出されることも考えられる. その結果, ブートストラップ標本に含まれる外れ値の割合は一定にならず, 外れ値が多く含まれ過ぎているブートストラップ標本も生成されてしまう.
ロバスト推定は外れ値を含むデータに対しても良い推定を行なう方法であるが, あまりに外れ値が多く含まれているとロバスト推定であっても良い推定にはならない. 外
れ値が多く含まれ過ぎてしまったブートストラップ標本から算出されたブートストラップ推定値は信頼できない値となってしまう. 図 1 は外れ値を含んだサンプルデー
タのロバスト推定量に対してブートストラップ法を適用した結果であるが, 良いと思われる推定値の付近 (A) から離れたところ (B) にも多くの推定値が流れてしまっている. これらは, 外れ値が含まれ過ぎたブートストラップ標本によるものである.

3. 見つけた局所最適解をそれぞれ比較し, 最も良いものを解となる推定量とする.
手順 1 で複数の候補を用意する理由は, 非凸の最適化問題だからである. 非凸の最適化問題の場合, 最適解に辿り着くためには最適解にある程度近い位置の初期解が必要となるが, それが何処にあるかが分からないため複数用意して虱潰しにする必要がある. 用意すべき候補の数 N はモデルの次数と, 「最適解に辿り着く初期解が N に含まれる確率」をどれだけ必要とするかで次のように変わってくる.
表 2 モデルの次数 p に対して用意すべき係数の候補の数

p+12351095%99%

101622349414530664713

手順 2, 3 は解の候補の数 N の分だけ計算が必要となるため, N の増加と比例して計算コストは増大する.
そこで, 考案したのが, N の数を減らす方法である. Nが複数必要なのは最適解が何処にあるかが分からないた∗
β n は元のめであるが, ブートストラップ推定量の最適解 β̂β n の近くにあると考えデータに対する推定量の最適解 β̂られる. それならば, ブートストラップ推定量の計算の際β n 一つだけで良いのではないかと考えた.には N は β̂また, このように初期解を一つに絞ることはもう 1 つの問題点に対する対策にもなる. 図 2 は, 外れ値を含むデータに対してロバスト推定を行なう際に, β の変動に対して目的関数がどう変動するかを調べたものである.

トストラップ法の場合 517 秒かかっていたのに対し, こちらの手法では 55 秒と 1/10 程度で済んでいる.

図 4 改良した手法によるブートストラップ推定量5
図 2 β に対する目的関数の変動
最適解は A であるが, B にも局所最適解が出来ている.
次に, このデータに対するブートストラップ標本で, 外れ値が含まれ過ぎてしまったデータを用意して同じように変動を調べたのが図 3 である.

Fast Robust Bootstrap

MM 推定量に対してブートストラップ法を適用する際に,前述の問題点を解決する手法として Salibian-Barrera andZamar [2] で提案された Fast Robust Bootstrap(FRB) 法がある
この手法では, 先程提案した改善案よりもさらに少ない計算コストでブートストラップを行なうことが可能である. (棚瀬 [5] 参照)
R∗

β n では, 先に挙げた問題点を解決するた
FRB 推定量 β̂
めに, ブートストラップ推定量の算出過程を簡略化している. まず, 次のような重み関数を定義する.
0

ωi =

ρ1 (ri /σ̂n )ri(9)0

なお, ri = yi − xiβ である. これを用いて MM 推定量の式 (5) を整理すると次のようになる.
MMβnβ̂
=[ n∑

ωi xi xii=1

図 3 β に対する目的関数の変動最適解であったはずの A が B より大きくなってしまっている. この時, 通常ならば B が最適解と判断されるが,提案した方法の場合 A の付近の局所最適解を最適解と判断するため, この場合でも A が推定量として採用される.
図 1 で使ったサンプルデータに対してこの方法でブートストラップを行なうと図 4 のようになった. 図 1 では右下にあった, 外れ値の影響を強く受けたブートストラップ推定量の群がなくなっている. 計算時間も通常のブー0]−1 n∑ωi xi yi(10)i=1

これは ωi を重み関数とした重み付き最小二乗法の形であり, 非凸の線形最適化問題であった MM 推定量の算出はただの線形問題になる. その結果, 計算コストは大幅に軽減されることになる.
しかし, 実際には ωi の中に変数 β が含まれているため,この問題は線形問題になってはいない. そこで, FRB 法
では ωi の中の β を変数ではなく定数に置き換えることでこの問題を線形問題とする.
∗β を繰り返し計
具体的には, ブートストラップ推定量 β̂
算する過程において, ωi の中の β を元の観測値から得ら
MMβnれた推定量 β̂∗βn =β̂で代用する.

[ n∑i=10
ωi∗ x∗i x∗i]−1 n∑i=1ωi∗ x∗i yi∗

(11)

0MM

β n である.
ここで, ri∗ = yi∗ − x∗i β̂
これにより, 計算コストを大幅に削減することが出来るが, 一部の変数を置き換えたことによる誤差が生じる. それを是正するため, Kn を定義する.
[ n]−1 n∑ 00∑00Kn = σ̂nρ1 (ri /σ̂n , xi )xi xi
ωi xi xi
i=1

それより後に発表された論文である Salibian-Barrera andAelst [3] における FRB 推定値の導出式には違いがあり,改良がなされていた. ここで紹介した導出式は新しいほうのものである.

6(12)

i=1

この Kn は元のデータに対してのみ計算し, ブートストラップ推定量を算出する過程で再計算を行う必要はない.
R∗β n は次のようにこれを用いて, 最終的な FRB 推定量 β̂なる.
R∗MM∗MM
β n = β̂
β n + Kn (β̂β n − β̂βn )β̂(13)
このように計算式を簡略化することで非凸の最適化問
題をただの線形問題へと変換し, 計算コストを大幅に削減することが出来る. また, これによりもう 1 つの問題点も同時に改善される. 重み関数 ωi∗ の中で, β は元のデータに対して算出されたものをそのまま使っている. そのため, ブートストラップ標本においてリサンプリングされた一つ一つの観測値に対し, 重み関数 ωi∗ が与えるウエイトは, 元のデータにおけるそれと同じになる. すなわち, 元のデータに対する分析において外れ値と判断された観測値はブートストラップ標本においても常に外れ値として扱われ, 推定量の算出に対してほとんど影響を及ぼすことが出来ない. 従って, ブートストラップ標本に外れ値が多R∗

β n は良く含まれすぎてしまった場合でも, FRB 推定量 β̂
い推定であり続ける.
図 1, 図 4 と同じサンプルデータへの MM 推定量に対して, FRB 法によるブートストラップを行なった結果が
図 5 である. 外れ値の影響を受けすぎた推定量は見られない. 計算時間については 2 秒程度しか要しておらず, 提案した改善案と比較しても段違いに速い.

おわりに

本研究では, ロバスト推定量に対してブートストラップ法を適用する際の問題点について考察し, その改善案についても模索した. 拙いながらも自身による改善案を提案したが, 既に提案されている FRB 法を超えるものとは言い難いものとなった. 提案したブートストラップは FRB法とは異なるアプローチに基づいたものであったが, この方向よりも FRB 法の方向へと研究を深めていくべきと感じた.
また, 他に改善案を考えるとするならば, ブートストラップ法自体を工夫する手段もあると考えている. 本研究では, ブートストラップ標本をリサンプリングする際には最も一般的な手法を土台にして研究を行なっていたが, ブートストラップ法には様々なバリエーションが存在する. 例
えば, 元のデータに対して最初にロバスト推定を行なった段階で, 外れ値と見做される観測値はブートストラップ標
本のリサンプリング時には除かれるようにすれば, ブートストラップ標本に含まれる外れ値の割合など気にしないで済むだろう.
研究結果としては少々中途半端な結果となったことが心残りではあるが, また機会があればもう少し研究を進めてみたい.

参考文献[1] Huber, P.J.(1981). Robust Statistics, Wiley, NewYork.
[2] Salibian-Barrera, M. and Zamar, R.H.(2002). Bootstrapping robust estimates of regression, The Annals of Statistics, 30, 556-582.
[3] Salibian-Barrera, M. and Aelst, S.V.(2008). Robustmodel selection using fast and robust bootstrap,
Computational Statistics and Data Analysis, 52,5121-5135.
[4] Salibian-Barrera, M., Willems, G. and Zamar,
R.(2008). The Fast-τ Estimator for Regression,
Journal of Computational and Graphical Statistics,
17, 659-682.
[5] 棚瀬暁俊 (2008). 線形回帰における MM-回帰推定量とロバスト・ブートストラップ法, 南山大学大学院
数理情報研究科 修士論文.
[6] 汪金芳・大内俊二・景平・田栗正章 (1992). ブート
ストラップ法 -最近までの発展と今後の展望-, 行動
計量学, 19-2, 50-81.

図 5 FRB 法によるブートストラップ推定量
なお, FRB 法の初出は Salibian-Barrera and Zamar [2]であるが, この時に提案された FRB 推定値の導出式と,[7] Yohai, V.J.(1987). High breakdown-point and higheﬃciency robust estimates for regression, The Annals of Statistics, 15, 642-656.