ほぼPython

Not技術ブログBut勉強ブログ 内容には誤りがあることが多いです

あらゆる次元に対応かつ比較的短いコードで線形回帰する

はじめに

前々回のブログでは最急降下法(勾配法)で直線モデルを求める - ほぼPythonという感じで、勾配法を使って2変数の場合の回帰をしました。

そして、前回のブログでは目で見て楽しむ重回帰分析 - ほぼPythonという感じで偏微分方程式を解くことで3変数の場合の回帰をしました。

今回はあらゆる次元に対応する線形回帰のコードを書いていきたいと思います。

数式

y軸の測定値(実際のデータ)をtとしておきます。また、データ数はNとします。このとき、予測されるモデルを
{}

\[ y = w_0x_0 + w_1x_1 + \cdots + w_{D-1}x_{D-1}+C\]


として、測定値との差の二乗の和の平均を各wで偏微分し、それらがすべて0と等しいという方程式を解くと結果的に、以下のようになります。

\[
\left(
\begin{array}{c}
w_0 \\
w_1 \\
\vdots \\
w_{D-1} \\
C\\
\end{array}
\right)
=
\left(
A^{\mathrm{T}}A \right)^{-1}A
\left(
\begin{array}{c}
t[0] \\
t[1] \\
\vdots \\
t[N-1] \\
\end{array}
\right)
\]

ただし、Aは以下のようにあらわされます。
行はデータ数、列は次元です。

\[
{ \scriptsize
A =
\left(
\begin{array}{lllll}
x_0[0] &x_1[0]&\cdots &x_{ \scriptsize D-1}[0] & 1\\
x_0[1] & x_1[1]&\cdots &x_ { \scriptsize D-1}[1] & 1\\
\vdots & \vdots &\vdots & \vdots&\vdots\\
x_0{ \scriptsize[N-1] }& x_1 { \scriptsize[N-1] }&\cdots &x_ { \scriptsize D-1}{ \scriptsize[N-1] }& 1\\
\end{array}
\right)}
\]

コード

上記の式を使えば、これまですごく長くてわかりにくかったコードをかなりすっきりと書くことが出来ます。
しかも、説明変数が2つでも100個でも、そしてデータ数が10個でも1000個でも、使うことが出来ます。(統計学的な問題は置いといて)

今回は実際に、前々回使ったデータを使って計算します。

import numpy as np

x0=np.array([i0 for i0 in range(0,10)]) #X0軸
x1=np.array([1 for i1 in range(0,10)]) #定数項のためにひとつ追加
X=np.matrix([x0,x1]) #説明変数を行列にする
t=np.matrix([5,1,6,8,10,12,17,16,20,25]) #y軸

Xt = X.T #転置行列を作る(縦ならび)
tt = t.T #転置行列を作る(縦ならび)

X_Xt = np.dot(X,Xt).I
X_tt = np.dot(X,tt)
w0,w1 = np.dot(X_Xt,X_tt)
print(w0,w1)

結果は[2.37575758, 1.30909091]となり、確かに正しいだろうと思われます。
前回までのコードよりかなり短くわかりやすくなりました。数学って素晴らしいですね・・・。

(ただ、数式上のAはプログラム上ではXtとなっており、縦横の関係が逆になってしまっています。行列を作る際に1つの命令のみで縦に並べる方法がわからなかったため、妥協しました・・・。わかる方いたら教えてください)