【機械学習】単回帰分析をpythonで実装してみましょう
今回は、教師あり学習の基礎中の基礎である「単回帰分析」を実装します。
教師あり学習とは、説明変数（インプット）から目的変数（アウトプット）を予測するモデルを求める手法です。
訓練データには目的変数や説明変数があり、あらかじめ作ったモデルに訓練データの説明変数を入力し、そのモデルからの出力が訓練データの目的変数に近づくようにモデルのパラメータを調整することで学習していきます。
単回帰分析
教師あり学習は目的変数のデータ形式によって、いくつかの種類に分類できます。
目的変数が株価など数値を取る場合を回帰(regression)、「男性・女性」「幼児・小学生・学生・大人」などのカテゴリになる場合を分類（classification）といいます。
教師あり学習のアルゴリズム（手法）には、重回帰（multiple linear regression）、ロジスティック回帰（logistic regression）、k近傍法（k-Nearest Neighbors）、決定木（Decision Tree）、サポートベクターマシン（Support Vector Machine）、ランダムフォレスト（Random Forest）、勾配ブースティング（Gradient Boosting）等があります。
これらの手法は、回帰で使われるときもあれば、分類で使われるときもあるので、注意しましょう。
単回帰分析は、以下のように目的変数1つが説明変数1つ(+定数項と誤差項)からなる線形の式で表し、予測を行う手法です。
推定するxの変数β1や、定数項β0のことを「パラメーター」と呼びます。
重みをもとに目的変数に対する説明変数の寄与を考えられる+線形モデルのため、根拠が明確で扱いやすいなどの利点があります。
別のコンテンツで解説しますが、説明変数を2個以上にすると重回帰分析になり、目的変数に対して非線形活性化関数(例えばシグモイド関数)を加えると単純パーセプトロン、何層も非活性化関数を重ねたものが多層パーセプトロン(深層学習)になります。
補足：数学的背景
では、線形代数を使って単回帰分析の数学的背景をご説明します。
$$f(x)=w_{1}x+w_{0}$$
単回帰では上のような直線を想定します。
xについては特徴ベクトル、wについては重みベクトルなので、以下のように整理できます。
$$\boldsymbol{x}=[x,1]^T$$
$$\boldsymbol{w}=[w_{1},w_{0}]^T$$
よって下のようなベクトルの内積の形で直線を表すことができます。
重回帰分析なら高次元での平面または超平面ですね。
ちなみに\(R^n\) において次元がちょうど 1 だけ小さい \(n − 1\) 次元 の空間のことを超平面 (hyperplane) といいます。
$$\begin{equation} \mathbf{y} = \mathbf{X}\mathbf{w} + \boldsymbol{\epsilon} \end{equation}$$
さて、線型回帰分析では予測値\(\hat{y}\)と教師データ\(y\)の二乗誤差を最小化するような重みベクトル\(\boldsymbol{w}\)を探したいです。
この損失関数を残差平方和と呼びます。
$$\sum_{i}^N||y_{i}-\hat{y}_{i}||^2$$
手法によって損失関数の種類は違ってくるし、目的によっても変わってきます。
RSSは全ての予測に対する誤差の合計を示すため、データの数が多いほど大きな値になる傾向があります。
平均二乗誤差と残差平方和は違うのか
平均二乗誤差（MSE）: MSEはRSSをデータの数で割ったもので、予測誤差の平均を示します。
MSEはデータの数によらず誤差の度合いを示すため、データセットの大きさが異なる場合でもモデルの性能を比較しやすいです。
定義式は以下の通りです。
$$\begin{equation} MSE = \frac{1}{n}\sum_{i=1}^{n} (y_i – \hat{y}_i)^2 \end{equation}$$
主な違いはRSSが誤差の総量を、MSEが誤差の平均を表している点です。
詳しくは以下のコンテンツをご覧ください
【MSEを最小化】ガウス・マルコフの定理と最良線形不偏推定量について
この重みベクトルは解析的に求めることができます。
まず損失関数であるRSS(残差平方和)を表します。
$$\begin{equation} \text{RSS}(\mathbf{w}) = (\mathbf{y} – \mathbf{X}\mathbf{w})^T(\mathbf{y} – \mathbf{X}\mathbf{w}) \end{equation}$$
次に重みベクトルwについて微分してイコール０にします。
$$\begin{equation} \frac{\partial \text{RSS}(\mathbf{w})}{\partial \mathbf{w}} = -2\mathbf{X}^T(\mathbf{y} – \mathbf{X}\mathbf{w}) = 0 \end{equation}$$
$$\begin{equation} \mathbf{X}^T\mathbf{y} – \mathbf{X}^T\mathbf{X}\mathbf{w} = 0 \end{equation}$$
$$\begin{equation} \mathbf{w} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \end{equation}$$
このように解析的に解を求めることができました。
ただし、\(\mathbf{X}^T\mathbf{X}\)が正則な場合に限ります。
では、pythonで実装していきます。
CODE
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# irisデータセットのロード
from sklearn.datasets import load_iris
iris_dataset = load_iris()
df = pd.DataFrame(iris_dataset.data, columns=iris_dataset.feature_names)
df.head()scikitlearnのirisデータセットを使います。
このままでは扱いづらいので上のようにデータフレーム型に変換しています。
df.head()で数行だけデータフレームを見てみましょう。
| sepal length (cm) | sepal width (cm) | petal length (cm) | petal width (cm) | |
| 0 | 5.1 | 3.5 | 1.4 | 0.2 |
| 1 | 4.9 | 3.0 | 1.4 | 0.2 |
| 2 | 4.7 | 3.2 | 1.3 | 0.2 |
| 3 | 4.6 | 3.1 | 1.5 | 0.2 |
| 4 | 5.0 | 3.6 | 1.4 | 0.2 |
一旦欠損値があるかどうかを確認します。
df.isnull().values.sum()こちら結果は0でした、欠損値はないようです。
次は相関行列を見て、各カラム同士の関係を見てみましょう。
corr_mat = df.corr(method='pearson')
corr_mat| sepal length (cm) | sepal width (cm) | petal length (cm) | petal width (cm) | |
| sepal length (cm) | 1.000000 | -0.117570 | 0.871754 | 0.817941 |
| sepal width (cm) | -0.117570 | 1.000000 | -0.428440 | -0.366126 |
| petal length (cm) | 0.871754 | -0.428440 | 1.000000 | 0.962865 |
| petal width (cm) | 0.817941 | -0.366126 | 0.962865 | 1.000000 |
df.corrの引数にはどの方法で相関係数を求めるかを入力することができます。
今回は一般的なピアソンの相関係数にしましたが、スピアマン(df.corr(“spearman”))の順位相関係数やケンドールの順位相関係数(df.corr(“kendall”))にすることもできます。
順位相関係数について詳しく知りたい方は、【外れ値に対処】順位相関係数と相関係数の違いについて | pythonをご覧ください。
さて、相関行列を見ると、sepal lengthとpetal lengthには強い相関があるようにみられます。
今回は、sepal lengthを説明変数にし、petal lengthを目的変数とした単回帰分析を行ってみましょう。
from sklearn import linear_model
X = df[["sepal length (cm)"]]
Y = df[["petal length (cm)"]]
model_lr = linear_model.LinearRegression()
model_lr.fit(X, Y)
#決定係数を算出する
my_result = model_lr.score(X, Y)
my_result決定係数は、0.759945‥となりました。まずまずのようです。
決定係数については、決定係数とは？説明変数の確らしさを図る指標の一つ。をご覧ください。
最後に、係数と定数項を確認しましょう。
w = model.coef_
b = model.intercept_
print(w)
print(b)\(w = 1.85843298, b = -7.10144337\)となりました。
これは、\(petal length = -7.10144337 + 1.85843298*sepal length +u_{i}\)という線形の式ができるということを意味します。
他の教師あり学習モデルについては、以下をご覧ください。
【XGB】交差検証法を使った勾配ブースティング決定木の実装|python
【python】Ridge(リッジ)回帰で多重共線性を解決する話
【ランダムフォレスト】ブートストラップ法を決定木に応用|python
【判別問題】サポートベクトルマシン(SVM)の仕組み|python
【kaggle】ベイズ最適化とXGBでtitanicの予測問題を解く|python
