統計ソフトの多くに付属しているボストンデータを例に重回帰について試してみる。
ボストンデータとは、ボストン市内の住宅物件の価格と様々な社会的データをまとめたものであり、それぞれの社会的データから価格を予測(評価)する手法を考察しようとするものである。
以下のような項目(Attribute)と物件の価格のデータが500件あまり用意されている。
Attribute | description |
---|---|
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 | 給与の低い職業に従事する人口の割合 (%) |
これらの項目を説明変数(feature 特徴量)、住宅物件の価格を目的変数(target)として、両者が線形な関係にあると仮定して、その係数を求めることを線形重回帰と呼ぶ。
説明変数を$(x_1, x_2, \cdots, x_M)$、目的変数を$y$で表せば、
$$ y = w_1 x_1 + w_2 x_2 + \cdots + w_M x_M $$参考:"scikit-learn boston 重回帰"くらいで検索すると山のようにページがでてくるが、上位のものとしては
そのほか、ランダムフォレストを使った分析もある。(後の回にふれるかも)
重回帰は、複数の 以下では、最初にあげたページに沿って実行(訓練)し、係数を決めた後、テストデータを使って予測の良さを調べてみる。
from sklearn.datasets import load_boston
boston = load_boston() # データセットの読み込み
# Bostonデータの項目確認
print(boston.DESCR)
以下で、randome_state=0として分割の仕方を固定している。実際には、何種類もの分割を試して、モデルの良さを確かめる、交差検定を行う。
from sklearn.model_selection import train_test_split
# 入出力の切り分け
x = boston['data'] # 物件の情報
t = boston['target'] # 家賃
x_train, x_test, t_train, t_test = train_test_split(x, t, test_size=0.3, random_state=0) # データの70%を訓練用に、残りをテスト用に
print(x_train[:3,:]) # 試しに最初の3行を出力
順序は
その後、
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x_train, t_train)
model.score(x_train, t_train) # 決定係数R^2の出力
決定係数
テストデータの実測値を$(t_1, t_2, \cdots, t_N)$、予測(推定)値を$(y_1, y_2, \cdots, y_N)$としたとき $$ R^2 = 1 - \frac{\sum_1^N (y_i - t_i)^2}{ \sum_1^N (t_i - \bar{t})^2} $$
# テストデータの予測値と実測値比較
x_predict = model.predict(x_test) # テストデータ(物件)の予測値
import matplotlib.pyplot as plt
plt.plot(x_predict, t_test, "o")
plt.xlabel("予測値")
plt.ylabel("実測値")
授業では、もっぱら説明変数が1つで目的変数との関係が非線形な場合を扱う。
ここで線形重回帰に触れたのは、重回帰の場合
$$ y = w_1 x_1 + w_2 x_2 + \cdots + w_M x_M $$で使ったLinearRegressionは非線形な回帰で用いることができるという点である。
一つの説明変数$x$と目的変数$y$が非線形な関係にあると推量される場合、、例えばその関係が「べき関数」だとすれば、各項目(Attribute)を$x$のべき乗として $$ y = w_1 x + w_2 x^2 + \cdots w_M x^M $$
のようにかけば、係数$\{w_i\}$を線形重回帰で推定することと同じであうことがわかるだろう。
$N$個のデータがあるとき、説明変数としては、各行に$x^1, x^2,\cdots, x^{M-1}$の値をセットした$N$行の行列を与えることになる。
[補足] 一般に、1つの説明変数$x$、1つの目的変数$y$の間に非線形な関係がある場合、$x$に対するいくつかの関数を用意し、それの線形結合によって表す。
$$ y = w_1\phi_1 (x) + w_2 \phi_2 (x) + \cdots + w_M\phi_M (x) $$のように書いて、線形重回帰と同様のスキームで扱う。このような場合、 $\{\phi_i\}$のセットを「基底」と呼び、データの性質によって、べき以外に、ガウシアンやシグモイド関数($\tanh$関数)、フーリエ基底(三角関数)を用いることが多い。フーリエ基底の場合はwaveletと呼ばれる。
# 回帰係数の出力
print(model.coef_)
# 回帰係数とfeatureの名前を結合してみる
import numpy as np
coef_names = np.stack([boston.feature_names, model.coef_])
print(coef_names)
# 上記が見にくいので、pandasのDataFrameに直して出力してみる
import pandas as pd
result= pd.DataFrame(coef_names)
# NOXが負の大きな値に、RMが正の比較的大きな値になっていることがわかる。
result
# 行と列を入れ替え
result = pd.DataFrame(coef_names).T
# 絶対値で降順にソート
result['abs_coef'] = np.abs(result[1].astype(float))
result.sort_values('abs_coef', ascending=False)
係数の大きさは、単純には寄与度を表すが、注意が必要 ⇒ 一番下の記述を参照!
その他はWeb上のいろいろなページを参照願いたい。
# 元のデータをpandasのDataFrameに入れて処理
import pandas as pd
boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)# 説明変数(data)
boston_df['PRICE'] = boston.target # 目的変数(target)をDataFrameの列として追加
boston_df.head()
import matplotlib.pyplot as plt
# 全データ項目の散布図
pd.plotting.scatter_matrix(boston_df,figsize=(20,20))
plt.plot()
一番下の散布図に注目してみよう。
縦軸は、目的変数の「価格」である。
RMと相関、NOXとLSTATが逆相関になっていることが目立つだろう。
RMとNOXの重回帰係数の絶対値が大きくなっていることと符合する。が、LSTATの係数が小さいのはなぜだろうか? ⇒ データの「正規化」あるいは標準化が必要なことを示唆している。
(時間があれば実際に行ってみよう)
ボストンデータのような様々な属性のデータを使って、予測する場合には、木構造にデータを分割するランダムフォレストによる回帰分析が良く用いられるる。 たとえば、アルゴリズムの説明を抜きにシンプルなコードのみを説明したページとして
https://hinomaruc.hatenablog.com/entry/2019/11/14/200857
を上げて置く。