最小2乗法の復習

昨日のTokyo.Rで、判別分析とかロバスト推定とかの話を聞いて、久しぶりにその辺りの復習をしようと思い立ちました。GW中に、線形回帰からロジスティック回帰まで復習できればいいなぁと思っています。

さて、線形回帰を行うにあたって、重要な要素である、最小2乗法について復習します。
最小2乗法の理解を深めるために実際に手を動かしてみます。

最小2乗法については、最小二乗法についてなどを見て頂くとして、今回は

これなら分かる最適化数学―基礎原理から計算手法まで
これなら分かる最適化数学―基礎原理から計算手法まで金谷 健一

共立出版 2005-09
売り上げランキング : 40962


Amazonで詳しく見る
by G-Tools
のP106にあるで連立1次方程式による表現を利用します。(こっちの方が拡張性がありますので)。

\begin{pmatrix}            \sum_{\alpha=1}^{N} x_{\alpha}^2 & \sum_{\alpha=1}^{N} x_{\alpha} \\ \sum_{\alpha=1}^{N} x_{\alpha} & \sum_{\alpha=1}^{N} 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{\alpha=1}^{N} x_{\alpha} y_{\alpha} \\ \sum_{\alpha=1}^{N} y_{\alpha} \end{pmatrix}
まで落とし込んだ後、データが与えられて実際に直線を引くまでを行います。
ここまでの式変形が一つの肝だとは思いますが、今回はそれが分かってるというところからのスターですね。
線形の単回帰分析なので、そこまで難しくはありません。

y=ax+bの直線について最小2乗法適用

参考文献の例題を解いてみます。

鉱石の密度x(g/cm^3)と鉄含有量y(%)の間の関係

x 2.8 2.9 3.0 3.1 3.2 3,2 3.2 3.3 3.4
y 30 26 33 31 33 35 37 36 33

このデータについて、y=ax+bの直線を求めてみます。

# XとYの値
Xs <- c(2.8, 2.9, 3.0, 3.1, 3.2, 3.2, 3.2, 3.3, 3.4)
Ys <- c(30, 26, 33, 31, 33, 35, 37, 36, 33)

# とりあえず点のplot
plot(Xs, Ys)

#X^2とxyを求める
Xs2 <- sum(Xs^2)
XY <- sum(Xs * Ys)

# aとbについての連立一次方程式を解く
est <- solve(matrix(c(Xs2, sum(Xs), sum(Xs), length(Xs)), 2, 2),
             matrix(c(XY, sum(Ys)), 2, 1))

# 与えられた点と一緒にプロット
plot(Xs, Ys)
abline(b=est[1,1], a=est[2,1], col="red")

# R組み込みのlm関数の結果と比較
test.frame <- data.frame(cbind(Xs, Ys))
test.lm <- lm(Ys~Xs, data=test.frame)
abline(test.lm, col="blue")

結果

(当たり前ですが、双方ほぼ同じ結果が得られました)

> est
          [,1]
[1,] 12.067669
[2,] -5.011278
>
>
> test.lm

Call:
lm(formula = Ys ~ Xs, data = test.frame)

Coefficients:
(Intercept)           Xs 
     -5.011       12.068

結果の図

こんな感じで直感通りの結果になっています。

lm関数の結果とほぼ一致していて自前実装の推定結果は重なって見えません。
なかなかいい感じです。

ちなみに、外れ値が入ると悲しい結果になります。
極端に変な値には対応していないからですね。
ここら辺がロバスト性と関係してくるところです。

多項式への当てはめ、へと続けていきたいと思います。