線形回帰超入門

こんにちは，ロジカル・アーツの笹原です．

会社から業務時間の一部を頂いて機械学習の勉強をしています．

機械学習は今年度から勉強を始めたのですが，学ぶことが多くてなかなか大変です．学習アルゴリズム一つとっても膨大な種類がありますし，そのアルゴリズムをどういった方法で実装するかを勉強していくのにも慣れないうちは時間がかかります．

何分私も初学者なので機械学習について説明できることは多くないのですが，自身の整理も兼ねて，ひとつ記事を書いてみようと思います．今回は，機械学習について学んでいく上でおそらく一番取っ付き易いであろう，線形回帰モデルをテーマにしました．

学習には，scikit-learn というライブラリに用意されているボストン市の住宅価格のデータセットを用います．住宅価格を予測する線形回帰モデルを決定することが今回の目標です．

実行環境
線形回帰とは
学習データについて
線形回帰モデルの学習について
- 学習（単回帰）
  - 学習結果の確認
- 学習（重回帰）
  - 学習結果の確認
おわりに
参考サイト

実行環境

環境は Google Colaboratory を使用します．Google アカウントさえあれば，だれでも使うことができます．時間的な制約はありますが，GPUが無料で使用できる，機械学習に用いられる代表的なライブラリが予め用意されている等，非常に有用です．

詳しいことは以下の記事にまとまっています．

qiita.com

また，言語は Python3 系を使用しています．

線形回帰とは

機械学習で用いられる学習アルゴリズムのうちの一つです．いわゆる教師あり学習 *1で用いられます．線形回帰自体は，統計学の回帰分析と呼ばれる手法の一部です．

線形回帰は，予測したい値が連続値を取る場合（今回の住宅価格も連続値になっています）に使用されます．この予測したい値を目的変数（または従属変数）と言い，目的変数を予測するのに用いる数値データを説明変数（または独立変数，特徴量）と言います．

線形回帰モデルとは，目的変数が $y$ ，説明変数が $x_1, x_2, \ldots, x_d$ の $d$ 個のとき，次の形で表されるモデルのことを言います： \begin{align} y = w_0 + w_1x_1 + w_2x_2 + \cdots + w_dx_d. \end{align}

ここで， $w_0$ はバイアス（または切片）， $w_1, w_2, \ldots, w_d$ は各変数 $x_i$ の重み（または回帰係数）を表します．特に $d=1$ のとき単回帰と言い， $d\geq2$ のときは重回帰と言います．

学習するとは，上記のパラメータ $w_i \ (0\leq i \leq d)$ を所定のアルゴリズムに従って決定することを指します．詳しくは学習の項目で述べます．

学習データについて

機械学習はデータありきです．学習に用いるデータが "良いもの*2" でなければ，そのデータから学習して得られるモデルが "良いもの" であることは期待できないでしょう．そのため機械学習では一般的に，データの前処理という作業が必要になります．様々な前処理の手法があるのですが，本記事の範囲を大きく超えるためここでは扱いません．今回用いる学習データセットはある程度整っているものなので，前処理はせず，いくつかの統計量を概観してみるだけにします．

まずはデータセットがどのようなものかを見ていきましょう．

ボストン市の住宅価格データセット

506 件のデータからなるデータセットで，各データは下の表にある 14 の属性を持っています．表の一番下にある属性 MEDV が今回予測する住宅価格になります．したがって，目的変数は MEDV，それ以外の 13 の属性が説明変数になります．なお，このデータセットには欠損値はありません．

属性名	説明
CRIM	町ごとの人口一人当たりの犯罪発生率
ZN	25,000 平方フィート以上の住宅区画の割合
INDUS	町ごとの小売業以外の事業の占める土地の割合
CHAS	チャールズ川沿いかどうか（川沿いなら 1，そうでなければ 0）
NOX	窒素酸化物の濃度（一千万分率）
RM	住居の平均部屋数
AGE	1940 年より前に建てられた居住者のいる物件の割合
DIS	5 つのボストン雇用センターからの重み付き距離
RAD	放射状高速道路へのアクセスのし易さの指標
TAX	$10,000 あたりの総額固定資産税率
PTRATIO	町ごとの教師一人当たりの生徒数
B	黒人の比率を Bk としたときの次の値；1000(Bk - 0.63)²
LSTAT	社会的地位の低い者の割合（%）
MEDV	居住者のいる住宅価格の中央値（$1000 単位)

以下で実際にデータセットを読み込んで統計量を見ていきますが，その前に必要なライブラリや関数等をインポートします．

必要なライブラリのインポート

データセットに関するものだけでなく，本記事で必要になるライブラリ等を全てインポートしておきます．ライブラリの詳細については省略しますが，ここでは

scikit-learn (sklearn) : 機械学習
pandas : データ解析
matplotlib : グラフ描画

のために用いています．

import pandas as pd
from matplotlib import pyplot as plt
from sklearn import linear_model
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split

データセットの読み込み

まず，load_boston 関数を用いてボストン市の住宅価格のデータセットを読み込みます．引数に何も指定しなければ，Bunch オブジェクトとしてデータが得られます．以下のコードでは boston という名前の変数に代入しています：

boston = load_boston()

Bunchオブジェクトは辞書のような（dictionary-like）オブジェクトです．通常の辞書型とは異なり，属性参照もできます．

キー	値
data	データセット（目的変数以外）
target	目的変数
feature_names	説明変数名
DESCR	データセットの詳細
filename	データセットの場所

詳細は下記を参照してください：
scikit-learn.org

データセットの確認

データセットの概観を知りたいときは pandas ライブラリが便利です．describe メソッドでデータの統計量を概観することができるのですが，対象のデータが Series オブジェクト，または DataFrame オブジェクトである必要があります．今回は 2 次元の表形式データ構造をもつ DataFrame オブジェクトに変換します．

次のコードでは説明変数のデータ boston.data を DataFrame オブジェクトとして取得しています．変数名は pd_boston とし，引数 columns では列名を指定しています．

pd_boston = pd.DataFrame(boston.data, columns=boston.feature_names)

describe メソッドでは，デフォルトで以下の表にある統計量が各列に対して得られます．

項目名	内容
count	データ数
mean	平均値
std	標準偏差
min	最小値
25%	25 パーセンタイル
50%	50 パーセンタイル（中央値）
75%	75 パーセンタイル
max	最大値

詳しくは以下を参照してください：
pandas.pydata.org

説明変数のデータについて統計量を見てみます．

pd_boston.describe()

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT
count	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000	506.000000
mean	3.613524	11.363636	11.136779	0.069170	0.554695	6.284634	68.574901	3.795043	9.549407	408.237154	18.455534	356.674032	12.653063
std	8.601545	23.322453	6.860353	0.253994	0.115878	0.702617	28.148861	2.105710	8.707259	168.537116	2.164946	91.294864	7.141062
min	0.006320	0.000000	0.460000	0.000000	0.385000	3.561000	2.900000	1.129600	1.000000	187.000000	12.600000	0.320000	1.730000
25%	0.082045	0.000000	5.190000	0.000000	0.449000	5.885500	45.025000	2.100175	4.000000	279.000000	17.400000	375.377500	6.950000
50%	0.256510	0.000000	9.690000	0.000000	0.538000	6.208500	77.500000	3.207450	5.000000	330.000000	19.050000	391.440000	11.360000
75%	3.677083	12.500000	18.100000	0.000000	0.624000	6.623500	94.075000	5.188425	24.000000	666.000000	20.200000	396.225000	16.955000
max	88.976200	100.000000	27.740000	1.000000	0.871000	8.780000	100.000000	12.126500	24.000000	711.000000	22.000000	396.900000	37.970000

同様に，目的変数についても統計量を見てみましょう．columns には単一の値であってもタプルやリスト等の形式で代入しなければならないことに注意してください．

pd.DataFrame(boston.target, columns=['MEDV']).describe()

	MEDV
count	506.000000
mean	22.532806
std	9.197104
min	5.000000
25%	17.025000
50%	21.200000
75%	25.000000
max	50.000000

線形回帰モデルの学習について

線形回帰モデルの定義を思い出しましょう．次の数式で定義されるものでした:

\begin{align} y = w_0 + w_1x_1 + w_2x_2 + \cdots + w_dx_d. \tag{lr} \end{align}

ここで， $y$ は目的変数， $x_i \ (0\leq i \leq d)$ は説明変数， $w_0$ はバイアス， $w_i \ (1\leq i \leq d)$ は各変数 $x_i$ の重みを表します．

さて，学習とは "パラメータ $w_i \ (0\leq i \leq d)$ を決定すること" と述べましたが，当然何でもいいわけではありません．学習したと言うからには，少なくとも学習済みのデータには精度の良い結果を返すことが期待されます．精度はモデルが予測した値と実際の値との "ずれ" を評価することで判断することができます．したがって，学習したデータセットに対する "ずれ" が最小になるようにパラメータ $w_i \ (0\leq i \leq d)$ を決定すればよいことが分かります．線形回帰モデルでは多くの場合，"ずれ" は平均二乗誤差で表されます．

説明のために，いくつか記号を準備します．

まず，データセットは $n$ 個のデータからなるとし， $i$ 番目のデータ*3を $(x_1^{(i)}, \ldots, x_d^{(i)}, t_i)$ で表すことにします． $x_1^{(i)}, \ldots, x_d^{(i)}$ は説明変数の値， $t_i$ は目的変数の値に該当します．このとき，上記の線形回帰モデル $(\mathrm{lr})$ による $i$ 番目のデータの予測値を $\hat{y}_i$ とすると，平均二乗誤差 $E$ は

\begin{align} E &= \frac{1}{n}\{(t_1 -\hat{y}_1)^{2} + (t_2 - \hat{y}_2)^{2} + \cdots + (t_n - \hat{y}_n)^{2}\} \\ &= \frac{1}{n}\sum_{i=1}^{n}(t_i - \hat{y}_i)^{2} \end{align}

で表されます．ここで $\hat{y}_i = w_0 + w_1x_1^{(i)} + w_2x_2^{(i)} + \cdots + w_dx_d^{(i)}$ ．

平均二乗誤差は，各データの実際の値と予測値の差を二乗したものを足し合わせ，データセット全体で平均を取ったものになっています．直感的には，この誤差が小さくなればなるほどデータセットに対する線形回帰モデルの精度が向上すると考えられます．実際その通りで，線形回帰では学習によって平均二乗誤差を最小にするようなパラメータ $w_i \ (0\leq i \leq d)$ を求めます．

ところで，上記の方法では学習に用いたデータセットに対するモデルの精度は測れますが，未知のデータに対するモデルの精度がどの程度のものであるかはどのようにして測ればよいでしょうか？結論から言えば，データセットを訓練データとテストデータに分割します．訓練データを学習用のデータセットとして用い，テストデータを未知のデータとして扱うわけです．このとき，テストデータに対する線形回帰モデルの精度を未知のデータに対するモデルの精度とみなします．ここでは決定係数と呼ばれるものを扱います．決定係数 $R^{2}$ は次の式で定義されます：

\begin{align} R^{2} = 1 - \frac{\sum_{i=1}^{n}(t_i - \hat{y}_i)^{2}}{\sum_{i=1}^n(t_i - \bar{t})^{2}}. \end{align}

ここで， $\bar{t} = \frac{1}{n}\sum_{i=1}^nt_i$ は目的値のデータセットにおける平均を表します．詳しい説明は省略しますが，決定係数は説明変数が目的変数をどの程度説明できているかを表す指標で，値が 1 に近いほどモデルの精度が高いと見なされます．

説明が長くなりましたが，背後の処理はライブラリに任せてしまって，いよいよ学習をさせてみましょう．まずは単純な単回帰から試していきます．説明変数としては，住居の平均部屋数（RM）を選択します．

学習（単回帰）

単回帰の場合，線形回帰モデルは次のように表されました：

\begin{align} y = w_0 + w_1x_1. \tag{lr_simple} \end{align}

ここで，目的変数 $y$ は MEDV，説明変数 $x_1$ は RM です．

これから，訓練データに対して線形回帰モデル $(\mathrm{lr\_simple})$ から計算される平均二乗誤差を最小化するパラメータ $w_0, w_1$ を決定します．

まずは boston.data から学習に用いる説明変数 RM のデータを抜き出します．

X_rm = boston.data[:, 5:6]

次に，train_test_split 関数で説明変数のデータ X_rm と目的変数のデータ boston.target をそれぞれ訓練データとテストデータに分割します．変数 X_* が説明変数のデータ，y_* が目的変数のデータになります．テストデータの割合は引数 test_size で指定します．今回は 25% をテストデータとしました．

X_train, X_test, y_train, y_test = train_test_split(X_rm, boston.target, test_size=0.25)

学習は， LinearRegression インスタンスに対して fit メソッドを呼び出すことでできます．ここではインスタンスを変数 clf_simple に代入しています．fit メソッドは引数に訓練データを与えることで，そのデータについて学習します．

clf_simple = linear_model.LinearRegression()
clf_simple.fit(X_train, y_train)

出力結果：

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

学習結果の確認

インスタンスに対して coef_ 属性を参照することで，重み $w_1$ が得られます．

clf_simple.coef_

出力結果：

array([8.58873838])

同様に intercept_ 属性を参照して，バイアス $w_0$ が得られます．

clf_simple.intercept_

出力結果：

-31.586881645149255

したがって，学習して得られた線形モデルは次のようになります：

\begin{align} y = -31.586881645149255 + 8.58873838\times x_1. \end{align}

今度は決定係数を見てみましょう．score メソッドにテストデータを渡すことで得られます．

clf_simple.score(X_test, y_test)

出力結果：

0.5243770944009284

あまり高くありませんね．．．

最後に結果を図示してみましょう．まずは predict メソッドで線形回帰モデルの予測結果を取得します．ここではデータ全体にわたって取得します．

predictions = clf_simple.predict(X_rm)

下の図は横軸を住居の平均部屋数（RM），縦軸を住宅価格の中央値（MEDV）とし，データセットの値を青丸で，得られた線形回帰モデルを赤い直線で表しています．

一応，密集しているところを通ってはいそうですね．

plt.scatter(boston.data[:,5], boston.target, color='blue')
plt.plot(boston.data[:,5], predictions, color='red')
plt.xlabel('RM')
plt.ylabel('MEDV')
plt.show()

f:id:logicalarts:20200115153756p:plain

学習（重回帰）

単回帰の場合と同様のことをしていきます．唯一異なるのは，学習に用いる説明変数の数が 1 から 13 になったことです．

なお，説明変数の順序は上記の表の上からの順序通りとなっています．

データセットの分割

X_train, X_test, y_train, y_test = train_test_split(boston.data, boston.target, test_size=0.25)

fit メソッドによる学習

clf = linear_model.LinearRegression()
clf.fit(X_train, y_train)

出力結果：

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

学習結果の確認

重み $w_i \ (1\leq i \leq 13)$ の表示

clf.coef_

出力結果：

array([-8.82485842e-02,  5.50257058e-02,  2.48789873e-02,  2.86344620e+00,
           -2.00802142e+01,  3.42428640e+00,  1.26609585e-02, -1.51068084e+00,
            3.27813800e-01, -1.14496429e-02, -9.71262354e-01,  8.77940578e-03,
           -5.87441678e-01])

バイアス $w_0$ の表示

clf.intercept_

出力結果：

40.0202774509255

決定係数の表示

clf.score(X_test, y_test)

出力結果：

0.7497022307896131

ぐっと上がりましたね．

残念ながら，単回帰のようにモデルを図示するには視認できる空間の次元が不足しています（14 次元必要です）．重回帰の場合は，各重み $w_i \ (1\leq i \leq 13)$ の大きさを図示してみましょう．重みの絶対値が大きいほど，その説明変数が目的変数に与える影響が大きいということになります．

plt.scatter(boston.feature_names, clf.coef_, color='blue')
plt.xticks(rotation=45)
plt.show()

f:id:logicalarts:20200115153803p:plain

図から CHAS，NOX，RM が特に影響が強いことが分かります．CHAS，RM が大きいほど MEDV も大きくなり，NOX が大きいほど MEDV は小さくなります．

おわりに

お疲れ様でした．本記事では線形回帰モデルをテーマに扱いましたが，データセットを訓練データとテストデータに分割することや "ずれ" を最小化するようにパラメータを決定するといった手法・考え方は機械学習全般に共通するところでもあります．線形回帰モデルについてはごく基本的なところは押さえたつもりではありますが，ここでは扱わなかったことも多くあります．いずれはそこも踏まえた内容の記事を書こうと目論んではいるのですが，いつになることやら．．．

少しでもお役に立てば幸いです．