sonickun.log

備忘録

最小二乗法のロバスト推定についてまとめた

統計学

　最近，研究とはほとんど関係ないところで統計学にハマっているので自己満でロバスト推定をやってみた．それからRのggplotモジュールで描けるグラフの美しさをみなさんに知ってほしい．

最小二乗法

　最小二乗法とは，測定で得られた数値に組を適当なモデルから想定される関数(ここでは1次関数)を用いて近似するときに，残差の二乗和を最小とするような係数を決定する方法である．残差とは，二変量の観測値（Xi, Yi）に対してその回帰式が y ＝ax + b で与えられるときの Yi - (aXi + b) の値をいう．
　近似直線の式を"y = ax + b"とした時のa, bは以下の行列計算によって求めることができる．

\begin{pmatrix} \Sigma x_i^2 & \Sigma x_i \\ \Sigma x_i & \Sigma 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \Sigma x_i y_i \\ \Sigma y_i \end{pmatrix}

\begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \Sigma x_i^2 & \Sigma x_i \\ \Sigma x_i & \Sigma 1 \end{pmatrix}^{-1} \begin{pmatrix} \Sigma x_i y_i \\ \Sigma y_i \end{pmatrix}

　この辺は前提知識として扱うので，詳しくは「最小二乗法」でググって一番上に出てきたサイトを参照されたい．
最小二乗法について

　例えば，以下の様な分散図に対して最小二乗法で近似した直線が描ける．
f:id:sonickun:20150119215129p:plain

　ところが，計測の際に何らかの拍子でプロットが1つだけ大きく外れた値(外れ値)を取ると，以下の青色のグラフように近似が大きくずれる．
f:id:sonickun:20150119215615p:plain

　グレーの網掛けになっている部分は回帰モデルに対する95%信頼区間を表している(グレーの領域に95%の確率でプロットが存在するという意味)．外れ値がある方は分散が大きく，95%信頼区間も大きく広がってしまう．

　このような外れ値が最小二乗法に大きく影響を与えるという問題点を解決する方法として，ロバスト推定法がある．

ロバスト推定

　ロバスト推定とは，誤差があるデータに対してその誤差の影響を最小にすることを目的とした理論である．ここではロバスト推定法の中の「M推定法」を扱う(他にはL推定やR推定がある)．
　M推定法とは，二乗誤差基準では誤差が大きいほど値が大きくなるが，ある一定以上の誤差で値の上昇を抑えることで，外れ値による影響を小さくする方法である．最小二乗法のロバスト推定における誤差とは，近似した直線と各プロットとの距離のことをいう．また重み付けにはTukeyのBiweight推定法という手法を採用する．

　点(Xi,Yi)と近似直線との誤差をd = Yi - (aXi + b)，Wを誤差の許容範囲としたとき，誤差が大きいほど最小二乗に与える影響力を小さくするように，以下のような式を用いて重みw(d)を計算する．

\begin{equation} w(d) = \left \{ \begin{array}{l} \Bigl\{1- (\cfrac{d}{W})^2 \Bigr\}^2　　　(　|d| \leq W　)\\ 0　　　　　　　　　　(　|d| > W　) \end{array} \right. \end{equation}

　この重みの関数w(d)をグラフで表すと，以下のようになる．
f:id:sonickun:20150119225459p:plain
(引用:ロバスト推定法（Ｍ推定法） 画像処理ソリューション)

　この重みWiを各近似データ(Xi、Yi)に関して計算し、Wiを付加した最小二乗法を再度行い，ロバスト推定による近似式y = a'x + b'のa', b'を求める．

\begin{pmatrix} \Sigma w_i x_i^2 & \Sigma w_i x_i \\ \Sigma w_i x_i & \Sigma w_i \end{pmatrix} \begin{pmatrix} a' \\ b' \end{pmatrix} = \begin{pmatrix} \Sigma w_i x_i y_i \\ \Sigma w_i y_i \end{pmatrix}

\begin{pmatrix} a' \\ b' \end{pmatrix} = \begin{pmatrix} \Sigma w_i x_i^2 & \Sigma w_i x_i \\ \Sigma w_i x_i & \Sigma w_i \end{pmatrix}^{-1} \begin{pmatrix} \Sigma w_i x_i y_i \\ \Sigma w_i y_i \end{pmatrix}


　ロバスト推定後のグラフは以下のようになる．
f:id:sonickun:20150119230933p:plain

　青いグラフのように外れ値による影響が抑えられ，測定データに近い直線が描けた．

　最終的に作成したRスクリプトは以下のとおり．

library(ggplot2)
library(MASS)

pdf("robust.pdf", width=8, height=4)

set.seed(123)          # allow reproducible random numbers
orig <- within(data.frame(x=1:10),
               {
                 type <- "orig"
                 y <- rnorm(x, mean=x)
               }
               )
outlier <- orig
outlier$y[2] <- 20
outlier$type <- "outlier"


theme_set(theme_grey(base_size = 16))  # increase default font etc. size

#p <- ggplot(orig,aes(x,y)) + geom_point(colour="#F8766D")
#p <- p + geom_smooth(method="lm", se=FALSE, colour="#F8766D")

p <- ggplot(data = rbind(outlier, orig), aes(x, y, colour=type, linetype=type)) + geom_point() # "base" plot, with points only
p <- p + geom_smooth(method = "rlm")         #  fit & plot *robust* model + envelope

print(p)
dev.off()

参考