| 表記 | 内容 | |
|---|---|---|
| $y_i$ | 目標変数の観測値 | |
| $x_{i,j}$ | 入力変数の観測値 | |
| ${\bf x}_i$ | $=(x_{i,1},\cdots,x_{i,p})^T$:入力変数のベクトルの観測値 | |
| $N$ | 観測データ数($i=1,\cdots,N$) | |
| $p$ | 入力変数の次元数($j=1,\cdots,p$) | |
| $D$ | $={({\bf x}_i,y_i) | i=1,\cdots,N}$:有限データ集合 |
| $y$ | 目標変数(target)、目的変数、非説明変数、従属変数、応答変数 | |
| $x_j$ | 入力変数(input)、説明変数、独立変数、特徴量(feature)、predictor | |
| ${\bf x}$ | 入力変数のベクトル |
＊目標変数や入力変数は様々な名称で呼ばれています。
既知のパラメータ$\beta_j$とし、パラメータのベクトルを${\bf \beta}=(\beta_1,\cdots,\beta_p)^T$とします。線形回帰モデルでは、入力${\bf x}$、既知のパラメータ${\bf \beta}$を用いて、目標変数$y$を回帰超平面（回帰直線）で予想します。 \begin{eqnarray} y({\bf x})&=&\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \\ &=&\sum_{j=0}^{p} \beta_j x_j \\ &=&{\bf x}^T {\bf \beta} \end{eqnarray} ＊ダミー変数$x_0=1$としています。
例として、Boston house-prices (ボストン市の住宅価格)のデータを用います。$i$を1から5まで出力します（ただし、Pythonのインデックスは0からはじまります）。
In [6]:
%pwd
Out[6]:
In [12]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
%matplotlib inline
In [2]:
boston = datasets.load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df.head(5)
Out[2]:
観測データ数$N$と入力変数と目標変数を合わせた次元数$p+1$をshapeメソッドで出力します。
In [13]:
df.shape
Out[13]:
データの各属性は以下の値が格納されています。
In [1]:
#from IPython.display import HTML
| 属性名 | 内容 |
|---|---|
| CRIM | 人口1人当たりの犯罪発生数 |
| ZN | 25,000平方フィート以上の住居区画の占める割合 |
| INDUS | 小売業以外の商業が占める面積の割合 |
| CHAS | チャールズ川によるダミー変数(1:川の周辺,0:それ以外) |
| NOX | NOxの濃度 |
| RM | 住居の平均部屋数 |
| AGE | 1940年より前に建てられた物件の割合 |
| DIS | 5つのボストン市の雇用施設からの距離(重み付け済) |
| RAD | 環状高速道路へのアクセスしやすさ |
| TAX | $10,000ドルあたりの不動産税率の総計 |
| PTRATIO | 町毎の児童と教師の比率 |
| B | 町毎の黒人(Bk)の比率を次の式で表したもの。1000(Bk–0.63)^2 |
| LSTAT | 給与の低い職業に従事する人口の割合(％) |
In [15]:
df.corr()
Out[15]:
In [13]:
sns.heatmap(df.corr())
Out[13]:
各データをユークリッド空間$\mathbb{R}^{p+1}$上の点とし、２次元のみ可視化します。
In [3]:
from highcharts import Highchart
In [22]:
x_label = 'AGE'
y_label = 'NOX'
H = Highchart(width=500, height=500)
data = list(zip(df[x_label], df[y_label]))
H.add_data_set(data, 'scatter', '観測値')
H.set_options('xAxis', {'title': {'text': x_label}})
H.set_options('yAxis', {'title': {'text': y_label}, 'lineWidth': 2})
H.set_options('title', {'text': '散布図'})
H
Out[22]:
In [23]:
import webbrowser
open('regression_scatter.html', 'w').write(H.buildhtml())
Out[23]:
In [24]:
webbrowser.open('regression_scatter.html')
Out[24]:
線形回帰モデルでは、 \begin{eqnarray} y({\bf x})&=&\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \\ &=&\sum_{j=0}^{p} \beta_j x_j \\ &=&{\bf x}^T {\bf \beta} \end{eqnarray} とし、パラメータ${\bf \beta}$を既知としました。実際には、${\bf \beta}$はわからないため、推定する必要があります。
どんな直線がいいかいろいろ試してみましょう。
In [25]:
from ipywidgets import interact
In [45]:
b0_min = int(df[y_label].min()*1000)/1000
b0_max = int(df[y_label].max()*1000)/1000
b0_step = int((df[y_label].max()-df[y_label].min())/100*1000)/1000
b_min = int(df[x_label].min()*1000)/1000
b_max = int(df[x_label].max()*1000)/1000
b1_step = int((b0_max-b0_min)/(b_max-b_min)*10000)/100000
b1_min = b1_step-b1_step*10
b1_max = b1_step+b1_step*10
In [48]:
@interact
def plot(b0=(b0_min,b0_max,b0_step),b1=(b1_min,b1_max,b1_step)):
x_tmp = np.linspace(int(df[x_label].min()),int(df[x_label].max())+1,50)
y_tmp = b0 + b1 * x_tmp
fig = plt.figure(figsize=(7,7))
plt.xlim([int(df[x_label].min()+0.5),int(df[x_label].max())+1])
plt.ylim([int(df[y_label].min()+0.5),int(df[y_label].max())+1])
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.scatter(df[x_label],df[y_label])
plt.plot(x_tmp,y_tmp,color='r')
plt.grid()
plt.show()
スライダーを動かすことで、切片と傾きを調節できます。
最小２乗法は誤差の２乗和を最小にするようにパラメータ${\bf \beta}$を決定します。
観測誤差$\epsilon_i$を \begin{eqnarray} \epsilon_i &=& y_i - y({\bf x}_i) \\ &=& y_i - \sum_{j=0}^{p} \beta_j x_{i,j} \\ \end{eqnarray} とします。観測誤差のベクトルは \begin{equation} {\bf \epsilon}= \begin{bmatrix} y_1 - {\bf x}_1^T {\bf \beta} \\ \vdots \\ y_N - {\bf x}_N^T {\bf \beta} \end{bmatrix}= \begin{bmatrix} y_1 \\ \vdots \\ y_N \end{bmatrix}- \begin{bmatrix} {\bf x}_1^T {\bf \beta} \\ \vdots \\ {\bf x}_N^T {\bf \beta} \end{bmatrix}={\bf y}- \begin{bmatrix} {\bf x}_1^T \\ \vdots \\ {\bf x}_N^T \end{bmatrix} {\bf \beta}={\bf y}-{\bf X}{\bf \beta} \end{equation} となります。ただし、${\bf y}=(y_1,\cdots,y_N)^T$、${\bf X}=(x_{ij})$とします。
誤差の２乗和を最小化する以下の問題を解くことにより、${\bf \beta}$を求めます。 $$\min_{\bf \beta} \sum_{i=1}^N (y_i - {\bf x}_i^T {\bf \beta})^2$$
損失関数を \begin{eqnarray} L({\bf \beta})&=&\sum_{i=1}^N (y_i - {\bf x}_i^T {\bf \beta})^2 \\ &=&(y_1 - {\bf x}_1^T {\bf \beta})(y_1 - {\bf x}_1^T {\bf \beta})+\cdots +(y_N - {\bf x}_N^T {\bf \beta})(y_N - {\bf x}_N^T {\bf \beta})\\ &=&({\bf y}-{\bf X}{\bf \beta})^T ({\bf y}-{\bf X}{\bf \beta}) \\ &=&({\bf y}^T-{\bf \beta}^T {\bf X}^T)({\bf y}-{\bf X}{\bf \beta}) \\ &=&{\bf y}^T {\bf y}-{\bf y}^T {\bf X}{\bf \beta}-{\bf \beta}^T {\bf X}^T {\bf y}+{\bf \beta}^T {\bf X}^T {\bf X}{\bf \beta}\\ &=&{\bf y}^T {\bf y}-2{\bf \beta}^T {\bf X}^T {\bf y}+{\bf \beta}^T {\bf X}^T {\bf X}{\bf \beta} \end{eqnarray} とします。${\bf \beta}$に関する多変数の２次関数なので、勾配が０となる点${\bf \beta}$を求めます。 \begin{eqnarray} \nabla L({\bf \beta})&=&-2X^T Y+2X^T X{\bf \beta}\\ &=&-2(X^T Y-X^T X{\bf \beta})\\ &=&0 \end{eqnarray} これを変形すると、正規方程式 $$X^T X{\bf \beta}=X^T Y$$ が導けます。$X^T X$が正則行列（逆行列を持つ）ならば $$\hat{{\bf \beta}}=(X^T X)^{-1} X^T Y$$ と解が求まります。また、ヘッセ行列$\nabla^2 L({\bf \beta})=2X^T X$が半正定値行列（${\bf \beta}^T {\bf X}^T {\bf X}{\bf \beta}=||{\bf X}{\bf \beta}||^2 \geq 0$を満たす）なので、$L({\bf \beta})$は凸関数（下に凸な関数）です。よって、$\hat{{\bf \beta}}=(X^T X)^{-1} X^T Y$が誤差の２乗和が最小値をとる最適解となり、これを推定値とします。
実際に最小２乗法を計算し可視化してみます。
In [49]:
cols = list(df.columns)
cols.remove('RM')
In [50]:
X = df[cols]#入力変数
X['dammy'] = 1
y = df['RM']#目標変数
print(X.head(5))
print(y.head(5))
In [51]:
#Xのムーア・ペンローズの擬似逆行列
X_mp = np.linalg.inv(X.T.dot(X)).dot(X.T)
#正規方程式の解
beta_hat = X_mp.dot(y)
beta_hat_s = pd.Series(beta_hat,index=X.columns)
print('正規方程式の解')
print(beta_hat_s)
In [56]:
x_label = 'DIS'
y_label = 'RM'
x_tmp = np.linspace(int(df[x_label].min()),int(df[x_label].max())+1,50)
y_tmp = beta_hat_s['dammy'] + beta_hat_s[x_label] * x_tmp
fig = plt.figure(figsize=(7,7))
plt.xlim([int(df[x_label].min())-1,int(df[x_label].max())+1])
plt.ylim([int(df[y_label].min())-1,int(df[y_label].max())+1])
plt.xlabel(x_label)
plt.ylabel(y_label)
plt.scatter(df[x_label],df[y_label])
plt.plot(x_tmp,y_tmp,color='r')
plt.grid()
plt.show()
＊単回帰ではないので、２次元プロットにおいて、もっともらしい直線にならないことに留意してください。
データ集合$D$は母集団から抽出された標本であると仮定します。すなわち、観測値$y_i$は確率変数$$Y=y(X)+\epsilon$$の確率分布に従って、生成された値と考えます。ここで、$X$は入力変数の確率変数のベクトルです。また、$\epsilon$~ $N(0,\sigma^2)$（$\epsilon$は期待値0、分散$\sigma^2$の正規分布に従う）と仮定します。観測値$y_i$は有限集合$\{y_i|i=1,\cdots,N\}$の要素でしたが、$Y$は無限集合$\mathbb{R}$を考えています。 $Y$の期待値は \begin{eqnarray} E[Y]&=&E[y(X)+\epsilon]\\ &=&E[y(X)]+E[\epsilon]\\ &=&E[y(X)]\\ &=&E[X^T {\bf \beta}]\\ &=&E[X^T]{\bf \beta} \end{eqnarray} ですが、$X=x$のときの条件付き期待値は \begin{eqnarray} E[Y|X=x]&=&E[y(x)+\epsilon]\\ &=&y(x)+E[\epsilon]\\ &=&y(x)\\ &=&x^T {\bf \beta}\\ \end{eqnarray} となります。また、$X=x$のときの$Y$の条件付き分散は \begin{eqnarray} Var[Y|X=x]&=&E[(Y-E[Y])^2|X=x]\\ &=&E[(y(X)+\epsilon-E[y(X)])^2|X=x]\\ &=&E[(y(x)+\epsilon-E[y(x)])^2]\\ &=&E[(y(x)+\epsilon-y(x))^2]\\ &=&E[\epsilon^2]\\ &=&\sigma^2 \end{eqnarray} となります。よって、$X=x$のとき$Y$~$N(y(x),\sigma^2)$（$Y$は期待値$y(x)$、分散$\sigma^2$の正規分布に従う）となります。
同時分布を$X$で条件付けしている。
EPEを最小化するためには，それぞれの$x$について最小化を行う．$c=f(x)$とすると \begin{equation} f(x)=argmin_{c} E_{Y|X} [(Y-c)^2|X=x] \end{equation} と求められる．この解は条件付き期待値 \begin{equation} f(x)=E[Y|X=x] \end{equation}
パラメータ${\bf \beta}$のとき、データ${\bf x}_i,y_i$を観測する確率 $$P[{\bf x}_i,y_i|{\bf \beta}]$$ とする。各サンプルは独立に観測されると仮定し、尤度関数 $$\Pi_i P[{\bf x}_i,y_i|{\bf \beta}]$$ と定義する。 対数尤度関数は $$L({\bf \beta}) = \sum_{i=1}^{N} \log P[{\bf x}_i,y_i|{\bf \beta}]$$ となる。正規分布を仮定すると \begin{eqnarray} L({\bf \beta}) &=& \sum_{i=1}^{N} \log \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{(-\frac{(y_i - {\bf x}_i^T {\bf \beta})^2}{2\sigma^2})} \\ &=& \sum_{i=1}^{N}( \log(1) - \log (\sqrt{2 \pi \sigma^2}) -\frac{(y_i - {\bf x}_i^T {\bf \beta})^2}{2\sigma^2})) \end{eqnarray} パラメータは${\bf \beta}$のみなので対数尤度の最小化は $$\min_{\bf \beta} \sum_{i=1}^N (y_i - {\bf x}_i^T {\bf \beta})^2$$ に等しい。
In [ ]: