【ロバスト回帰】最小二乗法と最小絶対値法

ロバスト回帰
単回帰
確率モデル
- ラプラス分布とガウス分布
  - 最尤推定でやってみる

ロバスト回帰

ロバスト回帰とは外れ値に頑強な回帰です。 Mean Squared Errorを最小化する、いわゆる最小二乗法の場合は外れ値が生じた場合に文字通り2乗で誤差が聞いてくるため、フィッティングが大きく外れ値側に寄ってしまいます。

これに対して、外れ値をデータに含んでいたとしても、それに大きく引っ張られることのない回帰の仕方をロバスト回帰と呼びます。今回はMean Absolute Errorを用いてロバスト回帰（の一種である最小絶対値法）を見ましょう。

単回帰

最小二乗法

今回は単回帰問題を扱います。データが

$$ D = \{ (x _ 1, y _ 1), ..., (x _ N, y _ N)\} $$

であるとして、モデルを

$$ y = a x + b $$

とします。ここで、実際には上記のモデルに完全には合致しないため各々のデータ $(x _ i , y _ i)$ は

$$ e _ i =y _ i - (a x _ i + b) $$

というズレを持っていることでしょう（もしもデータがモデルで表される直線に乗っているのであれば、 $e _ i = 0$ である）。さて、このズレは正にも負にもなりうるので、とにかくズレの大きさを表すために $e _ i ^ 2$ を考えて、この全データにおける平均

$$ {\rm MSE } =\frac{1}{N} \sum _ i ^ N e _ i ^ 2 = \frac{1}{N}\sum _ i ^ N \{y _ i - (a x _ i + b)\} ^ 2 $$

が小さくなるように $a , b$ を調整します。

最小絶対値法

モデルと実データのズレである

$$ e _ i =y _ i - (a x _ i + b) $$

を定式化するところまでは同じです。ズレ具合の大きさを表現する時に、二乗を使うのではなく絶対値を使うようにします。

$$ {\rm MAE } =\frac{1}{N} \sum _ i ^ N |e _ i| = \frac{1}{N}\sum _ i ^ N |y _ i - (a x _ i + b)| $$

二乗を使う場合と絶対値を使う場合では、モデルと実データがずれていた時に、パラメータ推定に与える影響力が変わります。二乗の方が大きなズレに敏感であるため、そのズレを無くそうと必死にフィッティングするようになります。

（※聞けば最小絶対値法の方が良さそうな気もしますが、実際は場合によります。場合によるというのは、扱っているデータの性質による（評価してみないとわからない）という意味でもありますし、最小二乗法はきれいに数式で解けたりしますが、絶対値は場合分けが必要だったり、結構面倒です）

MSEとMAEの比較

さて、あとはそれぞれのMSEとMAEの最適化を解くためにPyTorchを使います。PyTorchを使うのは良い選択ではありません（scipyを使いましょう）が、両方共、勾配法で同じ条件で殴ってみるということで比較してみようくらいのモチベーションです。

ライブラリをインポート

import torch

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("seaborn")

扱うデータの可視化

扱うデータの真の直線と、その直線にノイズが乗ったデータを10点程観測してみます。ただし、観測データのうち1点だけ、異様に値が大きくなっているというノイズを混入させます。

x_indexpoints = np.linspace(0, 10, 1000)
y_true = 2*x_indexpoints + 1

x = np.random.uniform(0, 10, 10)
y = 2*x + 1 + np.random.randn(len(x))

y[0] += 30

plt.scatter(x, y)
plt.plot(x_indexpoints, y_true)

f:id:s0sem0y:20190713162824p:plain — 異常値を含むデータ

モデル構築

1入力1出力の線形層を利用すればお望みの単回帰モデルが作れます。損失関数をそれぞれ、モデルの予測値の実際の値の誤差を使って

$$ \begin{align} {\rm MSE }& =\frac{1}{N} \sum _ i ^ N e _ i =\frac{1}{N} \sum _ i ^ N \{y _ i - (a x _ i + b)\} ^ 2 \\ {\rm MAE }& =\frac{1}{N} \sum _ i ^ N e _ i =\frac{1}{N} \sum _ i ^ N |y _ i - (a x _ i + b)| \end{align} $$

と書けるので、それも実装しておきましょう。最適化には勾配法を使います。

model_mse = torch.nn.Sequential(
    torch.nn.Linear(1, 1)
)
model_mae = torch.nn.Sequential(
    torch.nn.Linear(1, 1)
)


loss_mse = lambda y_pre, y_true: ((y_pre - y_true)**2).mean()
loss_mae = lambda y_pre, y_true: torch.abs(y_pre - y_true).mean()

optim_mse = torch.optim.SGD(model_mse.parameters() ,1e-3)
optim_mae = torch.optim.SGD(model_mae.parameters() ,1e-3)


def train(model, optimizer, loss, x, y_true):
    y_pre = model(x)
    loss_value = loss(y_pre, y_true)
    optimizer.zero_grad()
    loss_value.backward()
    optimizer.step()
    
    return loss_value.item()

学習

いざ学習します。 numpy.arrayをtorch.tensorにするのをお忘れなく。(下記コード程epochを回す必要はありませんが、損失関数が単峰なのでゴリ押しで行きます。単峰というのは語弊がありMAEは中心で微分不可能である。峰というよりはトンガリである。)

x_th = torch.tensor(x).float().reshape(-1, 1)
y_th = torch.tensor(y).float().reshape(-1, 1)

for epoch in range(10000):
    loss_value_mse = train(model_mse, optim_mse, loss_mse, x_th, y_th)
    loss_value_mae = train(model_mae, optim_mae, loss_mae, x_th, y_th)
    
    if (epoch + 1) % 1000 == 0:
        print("epoch {}  MSE LOSS {:.3f}, MAE LOSS {:.3f}".format(
            epoch + 1, loss_value_mse, loss_value_mae
        ))

結果の可視化

さて、出来上がったモデルによって引かれる直線を見てみましょう。

x_indexpoints_th = torch.tensor(x_indexpoints).float().reshape(-1, 1)
with torch.no_grad():
    y_pre_mse = model_mse(x_indexpoints_th).numpy()
    y_pre_mae = model_mae(x_indexpoints_th).numpy()

plt.figure(figsize=(12, 8), dpi=100)
plt.plot(x, y, "o")
plt.plot(x_indexpoints, y_true)
plt.plot(x_indexpoints, y_pre_mse)
plt.plot(x_indexpoints, y_pre_mae)
plt.legend(["obs data points", "true line", "regression on MSE", "regression on MAE"])

f:id:s0sem0y:20190713163602p:plain — 最小二乗法と最小絶対値法の比較

確かに、最小絶対値法は異常値に大きく引きずられること無くいい感じに直線を引いてくれていますね。

今回「異常値」だとか「外れ値」だとかいう言葉を使っていますが、そういうデータが実際に取れたということは、今後もそうなる可能性があるわけであり、本当に明確に混入した原因がわかっており、それが本来採取されるべきでなかったデータであると言えない限り、データから異常値と思われるものを削除してしまうことは、私は真面目なデータとの向き合い方だとは思いません。

なので、データを削除してきれいに整えてから、最小二乗法してしまえ、なんてことはすべきではないでしょう（削除するための基準として異常値であることを判定するモデリングができている→本来採取すべきでなかったデータと言える状況なら話は別）。

最小絶対値法なんてしていいの？

そもそも最小絶対値法を選択する時点で、それを頭のどこかで外れ値だと思っているのではないか？思っているのであれば削除して最小二乗法で良いのではないか？と思えてくることが時々あります。

しかし、問題は、そのような外れ値・異常値っぽいものが今後のデータ採取の中で「稀ではあるものの再現される」のかどうかであると思います。もしも滅多に起こらないけども、たまに生ずるということが繰り返されるのであれば、それは「異常値・外れ値」というの人間の感覚の問題であって、データを司る森羅万象からしたら「単に確率の低い、実在する現象」であるはずです。

外れ値を削除するというのは、「確率の低い、実在する現象」ということ自体を否定して、「そもそも存在しなかったはずのものである」と扱うことです。もしも、削除すべき外れ値が生ずるようなデータの採取・計測の仕方をしたならば、そのデータ収集自体をやり直すべきでしょう。ところが、そのデータ採取の1つの失敗が、今回集めた他のデータに影響していないとわかった場合どうでしょう。「データすべてを捨てるのではなく、失敗の影響を受けていなかったデータを残す」という姿勢になるはずです。すなわち外れ値を削除するというよりは、失敗の影響を受けていないデータを「残す」という姿勢です。

安易に外れ値だと言って削除してはなりません。それは「データ収集自体が何らかの要因で失敗した」というところから出発するはずです。その結果、失敗の影響を受けているデータを消して、影響を受けなかったデータだけを（もったいないから）残しているのです。外れ値削除というより生き残り確保です（持論です）。もはやそこそこ信用ならんデータなのです。

これに対して外れ値（と言われるもの）を残しながらフィッティングするのは、外れ値と呼びつつも、それは「極端に確率が低い事象のデータである」と認めて、モデルの中に組み込んで推定を行っているということです。最小絶対値法は、もっとも期待値の高いデータの中心から異様に離れた場所にデータが観測される確率を、最小二乗法より大きめに見積もったモデルになります。なので、フィッティングした直線自体が外れ値に引きずられることはなく、「ああ、そういうこともたまに起こるんだったな」という姿勢でフィッティングを行います。

確率モデル

ラプラス分布とガウス分布

最小二乗法がガウス分布に従うノイズが乗っているとみなした場合の最尤推定であるというのは有名な話です。一方、最小絶対値法はラプラス分布に従うノイズが乗っているとみなした場合の最尤推定と一致します。

ガウス分布に比べ、ラプラス分布は裾野が広く、遠くに飛んでしまういわゆる異常値・外れ値みたいなものがデータの計測上生じうることを受け入れた分布となっています。ちなみにそれぞれ

$$ \begin{align} {\rm Normal}(\mu, \sigma) &= \frac{1}{\sqrt {2\pi \sigma ^ 2} } \exp \left[ - \frac{(x-\mu) ^ 2}{2\sigma ^2}\right] \\ {\rm Laplace}(\mu, \phi) &= \frac{1}{ 2 \phi } \exp \left[ - \frac{|x-\mu|}{2\phi}\right] \end{align} $$

となっています。ここでラプラス分布の対数を取ってみましょう。

$$ \begin{align} \log [{\rm Laplace}(\mu, \phi)] &=\log \left[ \frac{1}{ 2 \phi } \exp \left[ - \frac{|x-\mu|}{2\phi}\right] \right]\\ & = \log \frac{1}{ 2 \phi } - \frac{|x-\mu|}{2\phi} & = - \log 2\phi - \frac{|x-\mu|}{2\phi} \end{align} $$

つまり、この符号を反転させたもの（すなわち負の対数確率）を見ると

$$ \log 2\phi + \frac{|x-\mu|}{2\phi} $$

となっています。このとき $\phi$ を適当な定数でおいてしまえば、 $\mu$ だけがパラメータとなります。最尤推定で最小化の目的関数となる負の対数尤度は、個々のデータの負の対数確率の和ですから、これは結局最小絶対値法と $\mu$ の推定に限ってみれば一致してしまうという理屈になります。

最尤推定でやってみる

最小絶対値法と最小二乗法を見たときと同じようにPyTorchで最尤推定を実施してみます。

model_gauss = torch.nn.Sequential(
    torch.nn.Linear(1, 1)
)
model_laplace = torch.nn.Sequential(
    torch.nn.Linear(1, 1)
)


def negative_gaussian_likelihood(y_pre, y_true):
    dist = torch.distributions.Normal(loc=y_pre, scale=torch.tensor(1.))
    return - dist.log_prob(y_true).mean()

def negative_laplace_likelihood(y_pre, y_true):
    dist = torch.distributions.Laplace(loc=y_pre, scale=torch.tensor(1.))
    return - dist.log_prob(y_true).mean()
    

optim_gauss = torch.optim.SGD(model_gauss.parameters() ,1e-3)
optim_laplace = torch.optim.SGD(model_laplace.parameters() ,1e-3)

for epoch in range(10000):
    loss_value_gauss = train(model_gauss, optim_gauss, 
                             negative_gaussian_likelihood, x_th, y_th)
    loss_value_laplace = train(model_laplace, optim_laplace, 
                             negative_laplace_likelihood, x_th, y_th)
    
    if (epoch + 1) % 1000 == 0:
        print("epoch {}  Gauss LOSS {:.3f}, Laplace LOSS {:.3f}".format(
            epoch + 1, loss_value_gauss, loss_value_laplace
        ))

with torch.no_grad():
    y_pre_gauss = model_gauss(x_indexpoints_th).numpy()
    y_pre_laplace = model_laplace(x_indexpoints_th).numpy()

plt.figure(figsize=(12, 8), dpi=100)
plt.plot(x, y, "o")
plt.plot(x_indexpoints, y_true)
plt.plot(x_indexpoints, y_pre_gauss)
plt.plot(x_indexpoints, y_pre_laplace)
plt.legend(["obs data points", "true line", "regression on Gauss", "regression on Laplace"])

f:id:s0sem0y:20190713172833p:plain — ガウス分布とラプラス分布での推定

最小二乗法と最小絶対値法を見たときと同じような結果が得られました。確率分布を仮定すれば、自分がばらつきをどのように考えているかを明示的に指定できるため、非常に明快です。最終的に最尤推定は損失関数の最小化問題に帰着されますが、損失関数を意味もわからずごちゃごちゃ選んだりいじったいるするのに疲れた方は、是非確率モデルを勉強してみてください。