はじめに
最小二乗法による回帰分析では、観測値と推定値の差の2乗が最小になるようにフィッティングを行うため、外れ値の影響が大きい。
このため外れ値を含んだデータに対しては、適切な前処理を行わないと誤った、また直感とは反した結果を得ることとなる。
この外れ値の影響力を小さくすることにより、正しい結果に近づけようという回帰分析をロバスト回帰という。
ロバスト回帰には、例えば最小絶対値法やHuber回帰、Tukey回帰などがあるらしい。
しかし、それらの方法では微分が不連続で計算に不便な絶対値や、場合分けが登場し、何だか煩雑な感じである。
そこで今回、残差を代入したガウス型関数を最大化することにより、通常の最小二乗法と近い計算でありながら、外れ値の影響力を抑える方法を提案する。
とは言え、シンプルな発想ゆえ恐らく既出の手法である（名前をご存じの方がいれば教えていただけると嬉しい）。
方法
以下の式を最小化する。
S = -\sum_i\exp{\left[-\left(\frac{y_i-f(x_i)}{\sigma}\right)^2\right]} \equiv \sum_i s_i
例えば、
f(x) = ax+b
の場合、
\begin{align}
\frac{\partial S}{\partial a} &= \sigma^{-2}\sum_i (2ax_i^2-2(y_i-b)x_i)\exp{\left[-\left(\frac{y_i-f(x_i)}{\sigma}\right)^2\right]}\\
&= \sigma^{-2}\sum_i (2ax_i^2-2(y_i-b)x_i)s_i\\
&=0\\
\frac{\partial S}{\partial b} &= \sigma^{-2}\sum_i (2ax_i-2(y_i-b))\exp{\left[-\left(\frac{y_i-f(x_i)}{\sigma}\right)^2\right]}\\
&= \sigma^{-2}\sum_i (2ax_i-2(y_i-b))s_i\\
&=0
\end{align}
を満たす $a,b$ を求めれば良い。
しかし、この式は $\exp$ の項に $a,b$ が含まれているため非線形であり、逆行列を用いた係数の計算ができない。
そこで、係数を1度の計算で求めるのではなく逐次更新することとし、$\exp$ の項における $a,b$ に更新前の値を使用することで定数係数として扱う。
これを明記すると、
\begin{align}
\frac{\partial S}{\partial a_{n+1}} &= \sigma_n^{-2}\sum_i (2a_{n+1}x_i^2-2(y_i-b_{n+1})x_i)\exp{\left[-\left(\frac{y_i-(a_nx_i+b_n)}{\sigma_n}\right)^2\right]}\\
&= \sigma_n^{-2}\sum_i (2a_{n+1}x_i^2-2(y_i-b_{n+1})x_i)s_{i, n}\\
&=0\\
\frac{\partial S}{\partial b_{n+1}} &= \sigma_n^{-2}\sum_i (2a_{n+1}x_i-2(y_i-b_{n+1}))\exp{\left[-\left(\frac{y_i-(a_nx_i+b_n)}{\sigma_n}\right)^2\right]}\\
&= \sigma_n^{-2}\sum_i (2a_{n+1}x_i-2(y_i-b_{n+1}))s_{i, n}\\
&=0
\end{align}
であり、初期値 $a_0,b_0$ には今回通常の最小二乗法による係数を採用する。
また、$\sigma_n^2$ はハイパーパラメータとして扱っても良いが、今回は残差の分散を使用する。
$\sigma_n^2 \rightarrow \infty$ でフィッティング結果は最小二乗法と一致する。
非線形項が定数に置き換わってしまえば、上式は単なる重み付き最小二乗法であるので、以下の計算で機械的に求めることができる。
{\left(
\begin{array}\\
\sum_i x_i^2s_{i,n} & \sum_i x_is_{i,n}\\
\sum_i x_is_{i,n} & \sum_i s_{i,n}
\end{array}
\right)
\left(
\begin{array}\\
a_{n+1}\\
b_{n+1}
\end{array}
\right)
=
\left(
\begin{array}\\
\sum_i x_i y_is_{i,n} \\
\sum_i y_is_{i,n}
\end{array}
\right)
}\\
\rightarrow\boldsymbol{X_n}\boldsymbol{A_{n+1}}=\boldsymbol{Y_n}\\
\boldsymbol{X_n^T}\boldsymbol{X_n}\boldsymbol{A_{n+1}}=\boldsymbol{X_n^T}\boldsymbol{Y_n}\\
\boldsymbol{A_{n+1}}=(\boldsymbol{X_n^T}\boldsymbol{X_n})^{-1}\boldsymbol{X_n^T}\boldsymbol{Y_n}.
この重みは、誤差が真の回帰曲線を中心とした正規分布に従うと仮定したときに、観測点がその時点で想定している回帰曲線を中心とした正規分布に属している確率と解釈できる。
つまりその時点での観測点の信頼度により重み付けを行い、再度フィッティングする操作を繰り返す形だ。
Pythonによる実装
テスト用に色々準備
import numpy as np
import matplotlib.pyplot as plt
n = 1000
d = 1 # d次関数でテスト
x = np.linspace(0, 5, n)
c = np.random.randn(d+1) * 3 # ランダムな係数
y = np.poly1d(c)(x) + np.random.randn(n)/2 # ノイズを付加
y[np.random.randint(n/1.2, n, n//20)] *= -1.3 # 外れ値を付加
plt.scatter(x, y)
plt.show()
フィッティング
x_ = (x ** np.c_[:d+1])[::-1]
for i in range(20): # 何回か更新
if i==0: # 1周目は通常の最小二乗法
X = x_ @ x_.T
coef = np.linalg.pinv(X) @ (y * x_).sum(1)
else: # 2周目以降は今回の手法
e = y - np.poly1d(coef)(x) # 残差
S = np.exp(-e**2 / e.var()) # ガウシアン重み
X = (x_ * S) @ x_.T # 上式のXnに対応
coef = np.linalg.pinv(X) @ (y * x_ * S).sum(1)
プロット
plt.plot(x, np.poly1d(np.polyfit(x, y, d))(x), c='r', label='LSM') # 最小二乗法
plt.plot(x, np.poly1d(c)(x), c='orange', label='ground truth') # 真値
plt.plot(x, np.poly1d(coef)(x), c='k', ls='--', label='Gaussian-Weighted') # 本手法
plt.scatter(x, y, s=1)
plt.legend(frameon=0)
plt.show()
外れ値の影響が大きい最小二乗法の結果（赤）と比較して、かなり真値に近い回帰結果が得られている。
より高次の関数も試してみよう。
より高次の関数でのテスト
5次関数
外れ値の影響はほぼ完全に排除できているように見える。
さいごに
理論・実装ともに簡単ながら、かなり強力に外れ値の影響を抑えられるのでぜひ一度お試しください。