非線形関数の回帰 - 元データ分析の会社で働いていた人の四方山話

前回、 $y=ax + b$ の線形単回帰について、最小２乗法で確認しました。

今回は非線形の $y=ax^2 + bx + c$ について見ていきます。
a,b,cについては、前回の拡張を利用して
$\begin{pmatrix} \sum_{\alpha=1}^{N} x_{\alpha}^4 & \sum_{\alpha=1}^{N} x_{\alpha}^3 & \sum_{\alpha=1}^{N} x_{\alpha}^2 \\ \sum_{\alpha=1}^{N} x_{\alpha}^3 & \sum_{\alpha=1}^{N} x_{\alpha}^2 & \sum_{\alpha=1}^{N} x_{\alpha} \\ \sum_{\alpha=1}^{N} x_{\alpha}^2 & \sum_{\alpha=1}^{N} x_{\alpha} & \sum_{\alpha=1}^{N} \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} \sum_{\alpha=1}^{N} x_{\alpha}^2 y_{\alpha} \\ \sum_{\alpha=1}^{N} x_{\alpha} y_{\alpha} \\ \sum_{\alpha=1}^{N} y_{\alpha} \end{pmatrix}$

という形式で求められます。
これを解くと、係数のa,b,cが求まるので、当てはめる２次関数が求まるというわけですね。

今回も

これなら分かる最適化数学―基礎原理から計算手法まで
	金谷健一共立出版 2005-09 売り上げランキング : 40962 Amazonで詳しく見る by G-Tools

のP111、例題4.5の(-1,0), (0,-2), (0,-1), (1,0)に、 $y=ax^2 + bx + c$ を当てはめてみることを考えます。

Rのコードです。
比較の為に、Rの組み込みの非線形回帰の関数nlsを利用しています。

# XとYの値
Xs <- c(-1.0, 0, 0, 1.0)
Ys <- c(0, -2.0, -1.0, 0)

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

#行列の各要素を先に求めておく
Xs4 <- sum(Xs^4)
Xs3 <- sum(Xs^3)
Xs2 <- sum(Xs^2)
X2Y <- sum(Xs^2 * Ys)
XY <- sum(Xs * Ys)

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

# 与えられた点と一緒にプロット
#png()
plot(Xs, Ys)
curve(est[1,1]*(x^2) + est[2,1]*x + est[3,1], add=TRUE, col="red")

# R組み込みのnls関数の結果と比較
test.frame <- data.frame(cbind(Xs, Ys))
test.nls <- nls(Ys~a*Xs^2+b*Xs+c, start=c(a=1,b=1,c=-1), data=test.frame)
#dev.off()

結果です。

> est
[,1]
[1,] 1.5
[2,] 0.0
[3,] -1.5
> test.nls
Nonlinear regression model
model: Ys ~ a * Xs^2 + b * Xs + c
data: test.frame
a b c
1.5 0.0 -1.5
residual sum-of-squares: 0.5
Number of iterations to convergence: 1
Achieved convergence tolerance: 4.441e-16

図では

のようになりました。

はい、うまく動いているようです。
曲線も絶妙な場所を通っています。

で、ここまでで分かるように、この最小２乗法は一般化できます。
（数式を書くのはめんどくさいので割愛）
すなわち、任意のn次式の係数を推定できるということですね。
どんなに複雑なモデルの当てはめもできてしまいます。
それに意味があるかは別ですが・・・

しかし、この最小２乗法は、ノイズに異常に弱いので気をつけてる必要があります。
というのも、上記の例で言うと、(4,-10)とかの外れ値が加わると・・・

のように、２次関数の向きが逆さまになるくらい影響を受けてしまいます。
このような外れ値やノイズへの耐性を持った推定方法がロバスト推定になります。

次回は余裕があればロバスト推定を扱いたいところです。