機械学習ノート: 線形回帰分析part1(参考:The element of statistical learning等)

#1. 要約
　本記事では線形回帰分析を取り上げ，係数の推定量の導出，プログラムの実装例を提示しました．回帰分析は機械学習か?という疑問を持ったあなた!気持ちはわかります．
　本記事は私の勉強まとめです．おそらく探せばもっとわかりやすい記事があると思うので，色々サーチしてみてください．
#2. 線形回帰分析とは
#2-1. 何に使われる？
　回帰分析は実用面でも広く利用されている解析手法で，スイーツなど(商品ならだいたい何でもいいです)の売り上げの予測に使われたり，現在持っているデータの説明のために用いられたりします．
　例えば以下のような散布図が得られたとします．
このようなデータを説明するときに，一つ一つのデータ点を列挙するよりもデータを要約するようななんらかの指標(e.g.平均，分散等)を用いた方が効果的かつ楽です．回帰問題では独立変数が変化したとき従属変数はどのように変化するかを関数(今回は直線の関数．曲線を用いてフィッティングする方法もあります)を用いて表すことを目的としています．

しかしながら，直線は以下のようにいくらでも好きに引けます(緑が最小二乗法による回帰直線)．
　そこでデータと直線の当てはまりを表す指標を定義して，当てはまりが最もよくなる直線を最適解として選択します．

#3. 導出
　データと直線の当てはまりを表す指標として，最小二乗基準や最尤法などがよく利用されています．以下では最小二乗基準を最適化する直線としての回帰直線の導出を行います．
###3-1. データが1次元の場合
　独立変数を$x_i$，従属変数を$y_i$とすると，従属変数と独立変数の関係(モデル式)は以下のように表されます．$$y_i = a+bx_i+e_i$$ここで$a$は直線の切片を，$b$は回帰係数を，$e_i$はi番目の個体(先ほどの散布図のデータ点だと思ってください)の回帰係数と独立変数の積と従属変数によって説明できなかった部分(=誤差)をそれぞれ表しています．最小二乗法はこの誤差をの二乗和を__モデルとデータの当てはまりの悪さ__の基準として利用し，これを最小化することで最適なモデルを求めます．今回の例でいうと，回帰係数$b$と切片$a$を求めることがこれに対応します．誤差$e_i$を赤線として先ほどの散布図で視覚化したのが以下の図です．
このように，データ点からy軸に平行に回帰直線に下ろした線が誤差として定義されます．最小二乗基準では，この誤差の二乗和を最小化するので，目的関数は$$LS(a,b)=\sum_{i=1}^ne_i^2=\sum_{i=1}^n(y_i-a-bx_i)^2$$で定義されます．これをa,bについて微分して0と置くことでa,bの最適値を計算することができます．それぞれを実際に微分してみると，

\begin{eqnarray}
\frac{\partial}{\partial a}LS(a,b) &=& -2\sum_{i=1}^n(y_i-a-bx_i)=0\tag{1}\\
\frac{\partial}{\partial b}LS(a,b) &=& -2\sum_{i=1}^nx_i(y_i-a-bx_i)=0\tag{2}
\end{eqnarray}

となり，(1)から，$$\widehat{a} = \bar{y}-\widehat{b}\bar{x}$$となり，(2)から

\begin{eqnarray}
\widehat{b}\sum_{i=1}^nx_i^2&=&\sum_{i=1}^nx_iy_i - \widehat{a}\sum_{i=1}^n x_i\\
\Leftrightarrow  \widehat{b}\sum_{i=1}^nx_i^2&=&\sum_{i=1}^nx_iy_i - \widehat{a}n\bar{x}\\
\Leftrightarrow  \widehat{b}\sum_{i=1}^nx_i^2&=&\sum_{i=1}^nx_iy_i - n\bar{x}(\bar{y}-\widehat{b}\bar{x})\\
\Leftrightarrow \widehat{b}(\sum_{i=1}^nx_i^2-n\bar{x}^2)&=&\sum_{i=1}^nx_iy_i-n\bar{x}\bar{y}\\
\Leftrightarrow \widehat{b}s_{x}^2&=&s_{xy}^2\\
\Leftrightarrow \widehat{b} &=& \frac{s_{xy}^2}{s_x^2}
\end{eqnarray}

が得られる．以上を整理すると，

\begin{eqnarray}
\left\{\begin{array}
1\widehat{a} &=& \bar{y} - \frac{s_{xy}^2}{s_x^2}\bar{x}\\
\widehat{b} &=& \frac{s_{xy}^2}{s_x^2}
\end{array}\right.
\end{eqnarray}

となることがわかりました．データが１次元の場合は，x，yの平均と，xの標本分散，x，yの共分散を計算すれば，切片aと回帰係数bは計算することができます．

###3-2. データがp次元の場合
　次にデータがp次元の場合($p\ge 2$，n個の観測点に対して変数がp個ある場合)を考えます．この場合，線形回帰分析は重回帰分析と呼ばれます．独立変数のデータを$n\times p$行列X，従属変数を$n\times 1$ベクトルyで表すとすると，モデル式は$$y =\beta_01_p + X\beta+e$$のようにかけます．$\beta_0,\beta,e$はそれぞれ切片ベクトル，係数ベクトル，誤差ベクトルを表しています．一般的に切片項は$X=[1_p,X]$，$\beta = [\beta_0,\beta]^T$として$$y=X\beta+e$$のように表されます．最小二乗法による損失関数は，

\begin{eqnarray}
LS(\beta) &=& \sum_{i=1}^n(y_i-\sum_{j=0}^px_{ij}\beta_j)^2\\
&=&(y-X\beta)^T(y-X\beta)\\
&=&y^Ty-2\beta^TX^Ty + \beta^TX^TX\beta
\end{eqnarray}

これを$\beta$で微分すると，

\begin{eqnarray}
\frac{\partial LS(\beta)}{\partial \beta}=-2X^T(y-X\beta)\tag{3}\\
\frac{\partial ^2LS(\beta)}{\partial \beta \partial \beta^T}=2X^TX
\end{eqnarray}

となります．ここで，Xが列フルランクだと仮定すると，$X^TX$は正定値行列となるので，逆行列をもちます．よって(3)を0と置くことで，$\beta$の推定値は$$\widehat{\beta}=(X^TX)^{-1}X^Ty$$となることがわかります．先ほどのデータが１次元の場合のbの推定値と形が似ていますね．
　逆行列の計算を行うため，途中でXが列フルランクだと仮定しましたが，この仮定が成り立たない場合はリッジなどの正則化項を損失関数に加えることで逆行列の計算ができるようになります．正則化の話はまたの機会に解説します．
#4. 実装例
　最初の3枚の図と重回帰分析のプログラムです．重回帰分析はsklearn.datasetsのBoston住宅価格データセットに対して行っています．

{python3, mlinear_reg.py}

# -*- coding: utf-8 -*-

#ライブラリの読み込み
import numpy as np 
import matplotlib.pyplot as plt 
from sklearn import datasets as sd
import pandas as pd 


# 重回帰分析の関数

def mlinear_regression(data,target):
	X = data.astype("float64")
	n,p = X.shape
	ones = np.ones(n).reshape(n,1)
	X1 = np.hstack([ones,X])
	beta = np.dot(np.linalg.inv(np.dot(X1.T,X1)),np.dot(X1.T,target))
	yhat = np.dot(X1,beta)
	return beta, yhat


if __name__ == "__main__":

	np.random.seed(10)
	x = np.linspace(-5,5,50)
	y = 3 * x + np.random.normal(0,10,len(x))

	#散布図
	plt.scatter(x,y,color="black")
	#plt.savefig("適当なディレクトリ")
	plt.show()

	#散布図with回帰直線
	plt.scatter(x,y,color="black")
	plt.plot(x,3*x,color="green")
	#plt.savefig("適当なディレクトリ")
	plt.show()

	#散布図with直線いっぱい
	plt.scatter(x,y,color="black")
	plt.plot(x,3*x,color="green")
	plt.plot(x,x,color="red")
	plt.plot(x,5*x,color="blue")
	#plt.savefig("適当なディレクトリ")
	plt.show()


	#重回帰分析
	###データの読み込み
	boston = sd.load_boston()
	independent = boston.data
	dependent = boston.target
	names = np.insert(boston.feature_names,0,"beta0")

	beta, yhat = mlinear_regression(independent, dependent)

	beta_df = pd.DataFrame(beta)
	beta_df.index = list(names)

	print(beta_df)
	#print("\nyhat\n",yhat)

#参考文献
Hastie, Trevor, Tibshirani, Robert, Friedman, Jerome(2009).The elements of statistical learning
稲垣宣夫「数理統計学」