Rのoptim関数でロバスト回帰（LMSとLAV）

Rのoptim関数でロバスト回帰（LMSとLAV）

• 1.
  2012.06 作成 optim 関数による OLS, LAV 及び LMS 回帰 OLS による通常の回帰分析 以下のような回帰モデルを考える。 y βX ε (1) ここで、各変数は以下のように定義する。 y1 1 x11  x p1 0 1 y  , X     , β 1 , ε   yn 1 x1 n  x pn n p y は n 個の説明変数ベクトル、X は p 次元 n 個の独立変数のデータ行列の頭に要素が全て 1 の列ベ クトルを付与した、(p+1)×n の行列。 は回帰係数ベクトルで、 β その要素のうち 0 は切片にあたる。 ε は誤差項。 β ˆ β このとき、OLS による回帰係数の不偏推定量 は、以下のような行列計算により得られる。 ˆ β t [ XX ] 1 t Xy (2)
• 2.
  ◆ R による乱数データを用いたデモ(1) set.seed(204) theta1 <- c(5, 2, 3) # 回帰パラメータ x1 <- matrix(runif(200, min=1, max=10), nc=2, nr=100) # x の値 x2 <- cbind(rep(1, dim((x1)[1]), x1) y1 <- x2 %*% theta1 + rnorm(100) # 推計値に正規残差付与 # プロット library(scatterplot3d) scatterplot3d(cbind(x1, y1), highlight.3d=TRUE, type="h", pch=16) scatterplot3d(cbind(x1, y1), highlight.3d=TRUE, type="h", pch=16, angle=145) # lm 関数による回帰係数の OLS 推計 lm1 <- lm(y1~x1) lm1$coeff # lm 関数による回帰係数の推計結果 # 行列計算による回帰係数の OLS 推計 solve(t(x2) %*% x2) %*% t(x2) %*% y1
• 3.
  [angle=145] [angle=90 (default)] > lm1$coeff# lm 関数による回帰係数の推計結果 (Intercept) x11 x12 5.092885 1.9718333.025425 > # 行列計算による回帰係数の OLS 推計 > solve(t(x2) %*% x2) %*% t(x2) %*% y1 [,1] [1,] 5.092885 [2,] 1.971833 [3,] 3.025425
• 4.
  OLS と残差の関係 t ˆ ˆ 得られた回帰係数 β を用いた y の推計値を y 、残差を r r1 ,  , rn とすると、OLS は残差 r の二乗和を最小化しているので、以下のように表現することができる。 n 2 minimize ˆ ri (3) β i 1 OLS が回帰係数の最適化の指標として、残差の二乗の和を使用しているために、残差の大きいデー タが二乗のオーダーで影響を持つことになり、そのために OLS は大きな残差を持つ外れ値の影響を 非常に大きく受けることになる。 LAV(Least Absolute Value)回帰 上述のような大きな残差の影響を弱めるために、Edgeworth(1887) が考案したのが、残差の絶対 値の和を最小化する LAV 回帰。これは L1 回帰と呼ばれる(L2 は最小二乗法)。 n minimize ri (4) ˆ β i 1 この改良により、回帰推計への外れ値の影響を弱めることはできるが、ただしもはやこの回帰係数 は、(2)のような行列計算により求めることはできず、線形計画法を必要とする。梃子比の大きな外 れ値に弱い。
• 5.
  LMS(Least Median ofSquares)回帰 OLS の場合、最小化するのは残差の二乗和だが、これを残差の二乗の中央値に変えたもの。 Hampel(1975) のアイデアに基づいて Rousseeuw(1984) が提案した。 2 minimize med ri (5) ˆ β i LMS は break down point(破綻点)が 50%で説明変数の外れ値にも強いが、計算効率が悪いのが 難点とされる。 ◆ R による乱数データを用いたデモ(2) optim 関数を用いて、上述の OLS, LAV、及び LMS 回帰を適用する。 # 乱数データ作成 [モデル: y1 = 5 + 2×x1 + 3×x2] set.seed(204) theta1 <- c(5, 2, 3) # 回帰パラメータ x1 <- matrix(runif(200, min=1, max=10), nc=2, nr=100) # x の値 x2 <- cbind(rep(1, length(y1)), x1) y1 <- x2 %*% theta1 + rnorm(100) # 推計値に正規残差付与 # OLS による重回帰 lm1 <- lm(y1~x1) ##### ここまでデモ(1)と同じ
• 6.
  # optim の最適化計算に用いる初期値の設定。ここでは適当な値をセット。 int1 <- c(4, 1, 2) # 初期値 1 モデルのパラメータにやや近い int2<- c(0, 0, 0) # 初期値 2 いいかげんな値 # optim 関数に入れるための、残差の二乗和を最小化する OLS 関数 myLM <- function(int1){ t(y1 - x2 %*% int1) %*% (y1 - x2 %*% int1) } # optim 関数に入れるための、残差の絶対値を最小化する LAD 関数 myLAV<- function(int1){ sum(abs(y1 - x2 %*% int1)) } # optim 関数に入れるための、残差の中央値を最小化する LMS 関数 myLMS <- function(int1){ median(t(y1 - x2 %*% int1)^2) } # 係数比較 lm1$coeff optim(int1, myLM) $par # 初期値 1 optim(int2, myLM) $par # 初期値 2 optim(int1, myLAV) $par # 初期値 1 optim(int2, myLAV) $par # 初期値 2 optim(int1, myLMS) $par # 初期値 1 optim(int2, myLMS) $par # 初期値 2
• 7.
  結果の概要 [1-1] lm関数の OLS(全データ) [1-2] optim 関数を用いた OLS(全データ) 、初期値 1 [1-3] optim 関数を用いた OLS(全データ) 、初期値 2 [2-1] optim 関数を用いた LAV、初期値 1 [2-2] optim 関数を用いた LAV、初期値 2 [3-1] optim 関数を用いた LMS、初期値 1 [3-2] optim 関数を用いた LMS、初期値 2 (Intercept) x11 x12 本来の値 5 2 3 [1-1] 5.092885 1.971833 3.025425 [1-2] 5.094411 1.971991 3.025025 [1-3] 5.090921 1.972390 3.025526 [2-1] 4.946654 2.013080 3.026913 [2-2] 4.939468 2.014040 3.026640 [3-1] 4.103185 2.097506 3.086221 [3-2] 3.713079 1.806043 3.350095
• 8.
  結果からわかること =>optimを使った LMS の場合、初期値は結果への影響が大きい。 Reference Hampel, F. R. (1975) Beyond Location Parameters: Robust Concepts and Methods (with Discussion), Bulletin of the ISI, Vol.46, pp.375-391 Rousseeuw, P. J. (1984) Least median of Squares Regression, Journal of the American Statistical Association, Vol.79, pp.871-880 E 入手済［RW-04 の Ref.］ Rousseeuw& Leroy(1987), Robust Regression and Outlier Detection, John Wiley & Sons, Inc.