ロバスト統計:外れ値を含むデータの扱い方

ロバスト統計とは

以下のような10個の観測値を得たとする。

2.773, 3.183, 2.969, 2.883, 3.229, 3.080, 3.204, 3.171, 2.798, 2.900

これらは標本平均 = 3.019、標本標準偏差 = 0.175であり、正規分布を仮定すれば95%信頼区間は

[3.019 - 1.96 × 0.175, 3.019 + 1.96 × 0.175] = [2.676, 3.362]

となる。実際このデータは平均 = 3、標準偏差 = 0.2 の正規乱数であるが、各統計量は良好な値を示している。しかしデータ収集時の記載ミスなどにより最後の値だけ小数点の位置がずれてしまったとすると、

2.773, 3.183, 2.969, 2.883, 3.229, 3.080, 3.204, 3.171, 2.798, 29.00

標本平均 = 5.629、標本標準偏差 = 8.213、95%信頼区間は [-10.468, 21.726]と、 一つの 外れ値 が含まれただけで途端に粗悪な推定結果となってしまった。 これは標本平均の性質として

ということがあり、一つの値に引っ張られやすいことが原因である。

ロバスト推定・統計 (robust estimation or statistics) とは、 外れ値に影響されにくい推定・統計的手法のことである。 最も単純なアイディアは、「直接眺めて外れ値を取り除く」ことである。 これは盲点となりやすいが可能であるならば最良と思われる。 本稿ではそれが困難である場合に用いるための計算手段を紹介する。

簡単なロバスト推定

外れ値を一つ含むデータ:

X* = { 2.773, 3.183, 2.969, 2.883, 3.229, 3.080, 3.204, 3.171, 2.798, 29.00 }

外れ値を取り除いたデータ:

X0 = { 2.773, 3.183, 2.969, 2.883, 3.229, 3.080, 3.204, 2.798, 3.171 }

外れ値を二つ含むデータ:

X** = { 2.773, 3.183, 2.969, 2.883, 3.229, 3.080, 3.204, 3.171, 27.98, 29.00 }

とする。以下ではデータの背後にある平均μや標準偏差σを 外れ値を含んだ X* あるいは X** からうまく推定する方法を考える。母集団分布は平均μを中心として左右対称であることは既知とする。

平均μの簡単なロバスト推定

 中央値

外れ値に影響を受けにくい推定量として広く知られているのは 中央値(median) である。 データX、標本数nに対して中央値は以下のように定義される。

Med(X*) = 3.126 と標本平均の値よりは大幅に真の値に近い。文字通り中央の値のみによるため外れ値の影響を直接には受けない。 ただし 外れ値が一方に偏る ような状況では注意が必要である。すなわち外れ値が大きい側(小さい側)に偏ることで、中央の位置が期待するよりも大きい側(小さい側)にずれてしまう。 中央の値しか見ていないためこのずれの影響が端的に表れてしまうのである。 実際、外れ値を二つ含むX** の中央値は:Med(X**) = 3.177 である。

 刈り込み平均

データに含まれる外れ値の割合はそこまで多くないものと思われるため、中央の値だけでなく中央の近傍を見て、もう少し多くのデータを参照しようというのが刈り込み平均(trimmed mean)である。 刈り込み平均とは、ある割合をαとすると、データの上側100α %と下側100α %を用いない標本平均のことである。式にすれば

となる。例えばαとして0.1とした場合、X* の刈り込み平均は3.065となり、外れ値を含まないデータの標本平均に近い。X**でもαを0.2とすれば刈り込み平均は3.139となる。この推定量の問題点はαを事前に設定しなければならない点にある。しかしながらこのαの値はある程度適当に決めても、 推定量の分散は小さくなる ことが知られており、 使いやすい推定量 である。 推定量の分散が小さくなることについてもう少し詳しく述べておく。表1は最尤推定量に対する刈り込み平均の漸近相対効率1である。母集団分布はt分布であり、νはt分布の自由度を表す。

表1:漸近相対効率
α
ν00.050.10.250.5
100.230.420.790.81
30.50.850.930.980.81
100.950.990.990.920.72
10.970.940.840.64

ν = 1 のときコーシー分布であり、νが大きくなるにつれて裾の軽い分布(大きな値が出にくい分布)となり、ν = ∞ のとき標準正規分布である。α = 0 は標本平均、α = 0.5 は中央値に相当する。 表1を見ると、αとして 0.1 から 0.25 の間に設定しておけば、外れ値が多いデータから少ないデータどれにおいても、良好な値を示すことができる。一方で標本平均はやはり外れ値に弱く、中央値は分散が大きい推定量であることがわかる。

標準偏差σの簡単なロバスト推定

平均μをロバストに推定する方法は上に示した方法の他にも重み付き平均などがあるが、 標準偏差σについては妥当に推定することが困難であり、次に示すMADN2という方法が最も頻繁に使われる方法である。

 中央絶対偏差(MAD)

標本分散の考え方の順を追っていくと、

  1. 平均μを標本平均によって推定
  2. 各標本のばらつき度合いを標本平均との二乗誤差によって評価
  3. ばらつき度合いの標本平均をとる

となる。式は

標本標準偏差はこの値の平方根である。これが外れ値に脆弱であることは、標本平均がそうであることからも明らかである。 実際にX0の標本標準偏差は 0.181 であるのに対して、X* の標本標準偏差は 8.213 とかけ離れている。

そこで上記手順における標本平均の代わりに中央値を用いた統計量として

がある。これは中央絶対偏差(Median Absolute Deviation)と呼ばれる。 MAD でも外れ値に対してロバストなバラツキの値となるのだが、 これは標準偏差の不偏推定量ではない。これを以下のように正規化した MADN が実際には用いられる。

データが正規分布に従うとき に、これは標準偏差σの不偏推定量になる。 MADN(X*) = 0.1927と標本標準偏差よりも明らかに良い値が得られる。 ただ中央値を用いているので、やはり 外れ値が一方に偏る 状況では注意が必要である。MADN(X**) = 0.2261と少し大き目の値になっている。

外れ値のロバストな同定

データX* から外れ値を同定することを考える。よく行われるのは、データが正規分布に従っていると考え、その信頼区間を推定する方法である。冒頭でも示したがX* の95%信頼区間は [-10.468, 21.726] であるので、X* の外れ値29.00 はこれによって同定することができる。一見うまくいっているように思われるが、これは結果論であり信頼区間の推定としては非常に粗いものである。 同じ方法でX** の95%信頼区間を算出すると [-12.874, 29.169] となり、二つの外れ値は同定できなくなってしまう(マスク効果という)。 そこでロバストな推定量を用いて、平均μとして中央値Med、標準偏差σとしてMADNとすると

  • X* の95%信頼区間 = [ 2.748, 3.503 ]
  • X** の95%信頼区間 = [ 2.734, 3.620 ]

となる。外れ値を同定することができ、推定量としても尤もらしい値となっている。

線形回帰モデルのロバスト推定法

次に線形回帰モデルに対するロバストな推定法を紹介する。 サンプルデータとしては図1のような二次元のデータを用いる。

  • 図1の左図:Phonesデータセット
    1950 - 73 年までのベルギーにおける電話回数。 1964 - 69年に外れ値が見られるが、これらは電話時間が記録されてしまっている。
  • 図1の右図:正規乱数

    として発生させた乱数。ただし赤点はデータを上書きして混入させた外れ値である。

    ./robust-sample-1.png

    図1:サンプルデータ

線形回帰モデル

のパラメータβは最小二乗法(LSE)によって推定されるのが一般的である。 最小二乗法とは、誤差を

とすると、以下のように二乗和誤差関数を最小にするようなβを、その勾配を0と置くことによって求めるものである。

式を見ると、標本平均と同様に個々のデータの誤差を等価に扱っているので外れ値に脆弱であることがわかる。

LMS・LTS

誤差列に対して二乗和をとるのではなく、 中央値や刈り込み平均 を用いる方法として以下がある3

  • 二乗誤差中央値最小化(Least Median of Squares; LMS)
  • 二乗誤差刈り込み平均最小化(Least Trimmed Squares; LTS)

    誤差は昇順にソートされている。qの値としては、pを説明変数の次元とすると ( n + p + 1 ) / 2 がよく使われる。

図2はLSEとLTSをサンプルデータで試した結果である。 LSE-outlier は外れ値なしで最小二乗法を行った結果であり、括弧内の値は各推定結果における傾きの値 (β_1) である。

./robust-sample-2.png

図2:線形回帰モデル推定結果(LSE・LTS)

  • 図2:左図(Phones)について
    63年および70年のデータについて記録ミスという情報はないとは言え、前後と比較すると少し不自然であることからこれらも記録ミスに類するものと言えるのではないだろうか。その観点からすれば、LSE-outlierよりもLSTの方がより好ましく、刈り込みが功を奏したと言えよう。
  • 図2:右図(正規乱数)について
    LTSの結果はLSEよりもデータのメインボディに沿っており良好である。だが外れ値ではない左下の数点も刈り込んでしまっているため、LSE-outlierよりはフィットできていない。

M推定

ここまでに示した推定量は、外れ値が混入していた場合はそれを無視して算出しようというものである。以下では無視するのではなく、データの外れ度合いに応じて重みを設定する推定方法について述べる。そのような方法を M推定 4という。

最小二乗法における推定方程式の誤差を以下のように関数ψを介して表現する。

ψを恒等関数とすれば最小二乗法の推定方程式となるが、先述のように目的変数yに外れ値が含まれているとそれに対応する誤差rの絶対値は大きくなり、推定結果に多大な影響を及ぼしてしまうことになる。そこで誤差が大きい場合は小さい重みになるような関数ψをうまく選んでやることで、外れ値の影響を減じようというのがこの方法である。 ψの候補としては複数あるが、典型的な二つを以下に挙げる。

  • Huber型
    ./robust-uber.png
  • Tuckey's biweight型
    ./robust-tuckey.png

どちらの型でも c を与える必要があり、その役割は外れ値の目安といったものである。 cを決める簡単な方法としては、データが正規分布に従っていると仮定し、その標準偏差σを MADN を用いて推定するという方法がある。95%信頼区間を用いるならば c = 1.96 σ とすればよい。その他の方法やψの種類などについては別途資料 3 4 を参照されたい。 図3はサンプルデータ正規乱数に対して、Huber型およびTuckey's biweight型(以下 Tuckey型)を用いたM推定による推定結果である。 Tuckey型の推定結果はLSE-outlierとほぼ同じ結果となっており、この分布とは相性が良いようである。 何回かシミュレーションをやり直して同様に推定を行ったが、Tuckey型は常にこのような良い推定値を出力し続けた。

./robust-sample-3.png

図3:線形回帰モデル推定結果(M推定)

まとめ

ロバスト統計におけるいくつかの計算手段を紹介した。この他にも順位相関係数のような順位統計量を用いるR推定や、ランダムにサンプルを抽出しLMSを繰り返すRANSAC 5 といった方法も提案されている。最適な手法やパラメータ、外れ値の定義などは当然データによって異なると思われ、十分にデータを吟味した上で用いたい。線形回帰モデルで見たようにモデルにうまく適用できれば、現実のデータとのギャップを解消することができると思われる。 本稿における結果はRおよびそのパッケージMASS(robustbase)を使用して算出した。

参考資料

  1. 藤澤洋徳 「統計数理研究所公開講座 ロバスト統計 -外れ値への対処の仕方-」(2013)
  2. Peter J. Rousseeuw and Annick M. Leroy “Robust regression and outlier detection” (1987)
  3. Maronna, Martin and Yohai “Robust Statistics: Theory and Methods” (2006)
  4. Huber and Ronchetti “Robust Statistics 2nd Edition” (2008)
  5. Martin A. Fischler and Robert C. Bolles “Random Sample Consensus” (1980)

Footnotes:

1 推定量の漸近分散の比率。漸近分散とは標本数が十分たくさんとれるときの推定量の分散。表中の値は最尤推定量の漸近分散を刈り込み平均による推定量の漸近分散で割ったもの。最尤推定量は、データが十分たくさんとれるとき、真の値に近づき分散はあらゆる推定量の中で最も小さくなる(漸近正規性およびCramer-Raoの不等式より)。よって最尤推定量の漸近分散との比を見ることで、刈り込み平均の分散を評価する。

2 四分位範囲も一般的であるが、ロバスト推定においてはMADNほぼ一択のようである。

3 解法については別途資料 2 を参照されたい。

4 なぜ「M」なのかは明確ではないようである。

Author: 数理システム知識工学部

Date: 2014-04-09

Generated by Org version 7.8.11 with Emacs version 24

© Mathematical Systems Inc. 2014