外れ値を考慮したコーシー分布に基づく線形回帰~最尤法とベイズの両面から~

今まで何度も

普通の回帰分析は「外れ値に引っ張られる」から注意した方が良い

という記事を上げてきました。

ただそうは言っても、ちょっと他のデータから外れていそうなデータを無根拠に弾いていては、

それは単なる恣意的なデータ操作では？

と言われかねません。

単にデータ自体の分布の裾が重く、そもそもたまに他から外れた値が飛び出すデータ生成構造を持っている可能性だってあるわけです。

そんなデータに対して、それでも線形回帰を行いたいというときどうすればよいのでしょうというのが今回の話。

先に答えから。

コーシー分布を使いましょう。
コーシー分布ってなんだ？って人はこちらをご覧ください。
期待値の無いコーシー分布の平均を取ると何が起こるか - バナナでもわかる話

ここでコーシー分布の確率密度関数を考えておきます。
$f(y;\gamma,y_{mid})=\frac{1}{\pi}\frac{\gamma}{\gamma^2+(y-y_{mid})^2}$

$y_{mid}$ が中央値のパラメータです。

「コーシー分布は、期待値が無いから回帰分析は無理だろう」というような気がするのですが、中央値に線形の仮定を与えてやりましょう。

$y_{mid}=\alpha+\beta X$

この下で最尤推定量を求めてやれば線形仮定の下中央値に回帰する分析が出来ますね。
最尤推定量に関する説明はこちら
基礎からイメージで学ぶ統計学~最尤推定量編~ - バナナでもわかる話

今回は計算を簡略化するため次の仮定の下で考えることにします。
$y_{mid}=\beta X$

$f(y_i;\gamma,\beta,X_i)=\prod_{i=1}^n(\frac{1}{\pi}\frac{\gamma}{\gamma^2+(y-\beta X_i)^2})$

対数を取ります。
$l(\gamma,\beta,X_i;y_i)=log(f)=\sum_{i=1}^n(-log(\pi)+log(\gamma)-log(\gamma^2+(y_i-\beta X_i)^2))$

$\gamma$ で偏微分してやると
$\frac{n}{\gamma}-2\gamma \sum_{i=1}^n\frac{1}{(\gamma^2+(y_i-\beta X_i)^2)}$

$\beta$ で偏微分してやると
$2\sum_{i=1}^n X_i(y_i-\beta X_i)\frac{1}{\gamma^2+(y_i-\beta X_i)^2}$

それぞれイコール0と置いてやって方程式を解きます......!!!

色々いじってみてください。
解けません！

最尤推定量は往々にしてこういう問題が発生します。これを解析的に解が求まらないと言います。

そこでどうするかというと、逐次的に値を変えていって解を探索するという手法で求めることになります。

ベイズではどうするかというと
$\gamma>0$ に注意しつつ事前分布を設定し、先ほどと同じような要領で尤度を決めます。
そして事後分布を求めるという手順を踏みます。

Stanにはcauchyという関数があるので、モデルの記述はラクチンですね。

model{
 for(n in 1:N)
  Y[n]~cauchy(b*X[n],gamma);
}

裾が重い分布ならt分布だってアリでは？何故コーシー分布？って思う方もいるかもしれません。

実際その通りでt分布でも良いと思います。でもまあちょっとt分布の関数ってゴチャゴチャしてるし、texで打つのもしんどいので今回は考えないことにします。やったことないけど多分出来ると思います。ベイズだともっと簡単にモデリング出来そうですよね。