標本データにおける外れ値の考え方②

こんにちは、データ分析の際に厄介者として扱われがちな外れ値について考えてきたいと思います。

外れ値とは、「データのメインボディから外れている値」のことをいいます。例えば、以下のようなデータがあるとき、25.5は怪しいと考えるのは自然だと思います。

$$
\mathcal{X} = \{5.6, 5.7, 5.4, 5.5, 5.2, 5.3, 5.8, 5.4, 5.6, 25.5 \}
$$

このような、値が現れたときにどのように考えていくかを書いていきます。

※こちら外れ値については2部作で書いていきます。この記事は2ndで、外れ値の対処法やロバスト推定の初歩の部分を書いていきます。

前回の記事はこちらです。

外れ値による影響

前回は外れ値を簡単に外れ値と決めつけるな！と主張してきましたが、今回の記事では外れ値っぽいものは全て外れ値として話を進めます。

標本データに外れ値が混入している場合、様々な推定量が正しく計算できないことがあります。

例えば、$${\mathcal{X}}$$の平均値は以下のようになります。

$$
E(\mathcal{X}) = \frac{5.6 + 5.7 + 5.4 + 5.5 + 5.2 + 5.3 + 5.8 + 5.4 + 5.6 + 25.5}{10} = 7.5
$$

一方で、$${\mathcal{X}}$$の外れ値である25.5を除外した場合の平均値は以下のようになります。

$$
E(\mathcal{X}) = \frac{5.6 + 5.7 + 5.4 + 5.5 + 5.2 + 5.3 + 5.8 + 5.4 + 5.6}{9} = 5.5
$$

標本データを見る限りだと、後者の方が正しい平均値を推定できています。外れ値が混入することにより、母集団分布に対する推定量がずれてしまうということが分かります。

平均値に限らず、以下の推定量なども外れ値の影響を受けてしまいます。

• 分散：大きくなる

• 相関係数：0に近づく

• 回帰係数：ブレが大きくなる

また、信頼区間も広がってしまい、それが外れ値同定を難しくしている場合もあります。$${\mathcal{X}}$$の平均値の95%信頼区間はおおよそ次のようになります。外れ値がない標本データを考えると、かなり広い範囲を取っているようにみえます。

$$
(-4.901292, 19.901292)
$$

それでは、外れ値がある上で正しい推定量を導くためにはどうしたらよいでしょうか。そのために、外れ値の混入を考慮した解析方法としてロバスト推定があります。

ロバスト推定の基礎

外れ値が混入している場合の推定として、分布の平均値を頑強に推定できる推定量中心に何種類か紹介します。

中央値

外れ値に頑強な推定量として真っ先に思い浮かべるのは中央値だと思います。中央値は順序統計値の内の一つで、値の大小の順番の情報のみを使って推定量を算出します。なので、値の絶対量の大きさはあまり影響しません。

実際に$${\mathcal{X}}$$の中央値は以下のようになり、外れ値の影響を排除できています。

> median(x)
[1] 5.55

刈り込み平均

以下の推定量、刈り込み平均を定義します。

$$
\hat{\mu}_{\alpha} = \frac{1}{n - 2m} \sum_{i = m + 1}^{n - m}{x_{[i]}}, m = \lfloor x \alpha \rfloor, 0 \leqq \alpha \leqq 1
$$

※見慣れない記号が沢山ありそうなので、補足します。

• $${x_{[i]}}$$は標本データを小さい順から並べたデータのi番目の要素

• $${\lfloor x \alpha \rfloor}$$は$${x \alpha}$$を超えない最大の整数

即ち、刈り込み平均は上下$${100 * \alpha}$$%のデータを除外したデータの平均値となります。

$${\alpha = 0.1}$$の場合の$${\mathcal{X}}$$の刈り込み平均は以下のようになります。確かに外れ値の影響を無視できています。

> mean(x, trim=0.1)
[1] 5.5375

しかし、前もって外れ値の割合をあらかじめ見積もる必要があります。（一応、大雑把に決めても妥当な数値になると言われています。）

重み付き平均

重み付き平均は、仮定した母集団分布から離れたデータの影響をほぼ無視するために定義された推定量です。母数$${\theta}$$を固定した重み関数を$${w(x_{i}: \theta), \sum{w(x_{i}: \theta)} = 1}$$とすると、重み付き平均は以下のようになります。

$$
\hat{\mu}_{w}(\mathcal{X}) = \sum^{n}_{i = 1}{w(x_{i}: \theta) x_{i}}
$$

実際に、N(5, 1)を仮定して$${\mathcal{X}}$$の重み付き平均を求めると以下のようになります。こちらも外れ値の影響をうまく無視できています。

> mean(dnorm(x, mean=5, sd=1) * x / sum(dnorm(x, mean=5, sd=1)))
[1] 0.5483626

この重み付き平均では、事前に重み関数およびその母数を仮定する必要があります。この問題についてはM推定の解を母数として利用することができます。

中央絶対偏差

平均値については、前述した推定量でロバストな推定が行えるようになりました。次に分散をロバストに推定する方法を考えます。

考え方としては、平均値の代わりに中央値を活用して、代表値からのズレの大きさを用います。そのようにして、定義された以下の数式を中央絶対偏差と呼びます。分母の0.675は正規化項です。
（0.675という数値は実際の標準偏差と比較したときにこのくらいの値になるからだそうです。）

$$
MADN(\mathcal{X}) = Med(\{ | x_{i} - Med(\mathcal{X}) | \}^{n}_{i = 1}) / 0.675
$$

実際に$${\mathcal{X}}$$の中央絶対偏差を求めると、以下のようになります。通常の標準偏差よりもかなり値を小さく抑えることができています。そして、25.5を除外した$${\mathcal{X}}$$を考えると、こちらの方が妥当な値とみることができます。

> mad(x) # 中央絶対偏差 
[1] 0.22239
> sd(x)  # 標本標準偏差
[1] 6.32719

しかし、この中央絶対偏差は正規分布に従わないデータに対しては不偏推定量ではありません。ここだけは注意が必要です。

実務では、標本標準偏差の代わりにこの中央絶対偏差で平均の信頼区間を計算することで外れ値を割り出すという手法が使われたりします。

より一般化されたロバスト推定

※ここから先はやや上級者向けになります。

推定量の推定を考える上で外せないのが、最尤推定量です。ようは、標本データに対して、確率関数が最大となる母数を探す方法で、以下の数式で定義されます。

$$
\hat{\theta} = \argmax_{\theta} \prod^{n}_{i = 1} f(x_{i}; \theta) = \argmin_{\theta} \sum^{n}_{i = 1} - \log f(x_{i}; \theta)
$$

この$${- \log f}$$をより抽象化した関数$${\rho}$$を用いた以下の方程式を解くことをM推定と呼びます。（抽象化しない$${- \log f}$$をそのまま使う場合は、最尤法になります。）

$$
\hat{\theta} = \argmin_{\theta} \sum^{n}_{i = 1} \rho(x_{i}; \theta)
$$

ここで、$${\psi(x) = \frac{\partial}{\partial x} \rho(x)}$$とすると、上の方程式は以下の形で単純化することができます。（一般的にはこちらの形を使うことが多いです。）

$$
\sum^{n}_{i = 1} \psi(x_{i}; \theta) = 0
$$

ここから、平均のM推定に絞ってどのようにこの方程式が利用されているかを見ていきます。

平均のM推定

母数に平均を用いるので、関数$${\psi(x; \mu) = \psi(x - \mu)}$$を考えます。ちなみに、標本平均もこの$${\psi}$$にあてはめることができます。

ここでのM推定は以下の方程式を解くことです。そして、平均のM推定を行うときによく利用される$${\psi(x - \mu)}$$は4つあります。

$$
\sum^{n}_{i = 1} \psi(x_{i} - \mu) = 0
$$

• 刈り込み型

$$
\psi(x - \mu) =
\begin{cases}
&x - \mu &(|x - \mu| \leqq c) \\
&0 &(|x - \mu| > c)
\end{cases},
c > 0
$$

この関数は、刈り込み平均と同じような意図で、事前に決めた定数$${c}$$よりも大きい値を外れ値とみなして、除外して平均を推定します。

• フーバー型

$$
\psi(x - \mu) =
\begin{cases}
&-c &(x - \mu < -c) \\
&x - \mu &(|x - \mu| \leqq c) \\
&c &(x - \mu > c)
\end{cases},
c > 0
$$

この関数は、刈り込み平均と似ています。事前に決めた定数$${c}$$よりも大きい値を外れ値とみなして、$${c}$$に置換した上で平均を推定します。刈り込み平均と異なるのは、外れ値の影響を全く無視していません。

※ある凸最適化問題では積極的に使われるらしいです。

• Bisquare型

$$
\psi(x - \mu) =
\begin{cases}
&(x - \mu) \{ 1 - (\frac{x - \mu}{c})^{2} \}^{2} &(|x - \mu| \leqq c) \\
&0 &(|x - \mu| > c)
\end{cases},
c > 0
$$

この関数は、やや複雑な形をしていますが、刈り込み平均を微分可能な形にしたものです。即ち、$${c}$$の前後が連続となるように値を調整したものです。

事前に決めた定数$${c}$$よりも大きい値の影響を無視するという点は共通しています。

• 重み付き型

重み付き平均と同じ発想で、平均から離れるほど重みが減衰するような関数を考えます。例えば、母集団分布の密度関数が$${\phi(x; \mu)}$$で与えられるとき、関数は以下のように与えられます。

$$
\psi(x_{i} - \mu) = \phi(x; \mu)^{\gamma} (x - \mu), \gamma > 0
$$

この関数をM推定の方程式に当てはめて変換すると、平均の推定値が求まります。

$$
\mu = \frac{\sum^{n}_{i = 1} \phi(x_{i}; \mu)^{\gamma} x_{i}}{\sum^{n}_{i = 1} \phi(x_{i}; \mu)^{\gamma}}
$$

再下降性

様々な$${\psi}$$を見てきましたが、関数は何でも良いわけではありません。特に外れ値の影響を除外する上では、$${|x| \rightarrow \infty}$$において発散しないような関数を使う必要があります。

特に、$${\psi}$$が以下の性質を満たす場合、再下降性を持つと呼び、外れ値を影響をほぼなくす効果があります。

$$
\lim_{|x| \rightarrow \infty} \psi(x) = 0
$$

また、裾が重い分布を使った最尤推定によるM推定も再下降性を持ち、外れ値に有効です。これは尤度型と呼びます。

M推定の解法

M推定の方程式は多くの場合解析的でないため、アルゴリズムを使って近似解を得ます。例えば、平均のM推定においては、以下のようなアルゴリズムによって母数$${\mu}$$を推定します。

$$
\begin{aligned}
0 &= \sum^{n}_{i = 1} \psi(x_{i} - \mu) = \sum^{n}_{i = 1} \frac{\psi(x_{i} - \mu)}{x_{i} - \mu} (x_{i} - \mu) \\
&= \sum^{n}_{i = 1} W(x_{i} - \mu) (x_{i} - \mu)
\end{aligned}
$$

ただし、

$$
W(x) =
\begin{cases}
\frac{\psi(x)}{x} &x \neq 0 \\
\psi'(x) &x = 0
\end{cases}
$$

この方程式を変形して、

$$
\mu = \frac{W(x_{i} - \mu) x_{i}}{W(x_{i} - \mu)}
$$

この、右辺の$${\mu}$$を使って左辺の$${\mu}$$を更新するような漸化式を作成することで最終的に収束した$${\mu}$$を推定値とします。$${\mu}$$の初期値は一般的に中央値を使います。

$$
\mu^{k + 1} = \frac{W(x_{i} - \mu^{k}) x_{i}}{W(x_{i} - \mu^{k})}
$$

分散などの他の推定量も、方程式を上手く変形することによって推定値を求めます。（このときにダイバージェンス理論が役に立ちます。）

また、信頼区間や検定などの他の推定理論は漸近理論から導くことができます。長くなるので今回は割愛しますが、ぜひ調べてみてください！！

ということで、外れ値の対処法やその理論について色々と書いてみました。実務ではM推定までは必要ないと思います。しかし、その理論立てもしっかり考えられていることが伝われば幸いです。

個人的には思いっきり数学が書けて楽しくなりましたw。趣味でやる分にはとても良いと思いました！

最後までお読みいただきありがとうございました。