PyLearnST 07: 線形回帰モデル¶
線形回帰モデルとは、$n$個の独立変数$x_i(i=1, \cdots, n)$から従属変数$\hat{y}$を導くモデルで、
$$ \hat{y}=\sum_{i=1}^{n} a_i x_i + b, $$と定式化されます。ここで、$a_i(i=1,\cdots,n)$と$b$は最適化すべきパラメータとなります。 $\hat{y}$は観測された実際のデータ$y$の推定値です。 そのため、分析者が対象$y$を予測・推定したいときに使用することが普通です。 良い予測とは、予測対象である$y$とその推定値$\hat{y}$のずれが小さいような予測モデルですので、パラメータ$a_1, \cdots, a_n, b$の最適な値$a_1^{\rm opt}, \cdots, a_n^{\rm opt}, b^{\rm opt}$は、
$$ (a_1^{\rm opt}, \cdots, a_n^{\rm opt}, b^{\rm opt}) = {\rm Argmin}_{a_1, \cdots, a_n, b} |\hat{y}-y|, $$によって定義されることが一般的です。と言っても少しややこしいと思うので、簡単化のため、独立変数の数を2個に設定して見ます。すなわち、$n=2$であり、
$$ \hat{y}=a_1 x_1 + a_2 x_2 +b, $$が線形回帰モデルとなります。ここで、
- $\hat{y}$: 試験得点の推定値
- $x_1$: 試験前の勉強時間
- $x_2$: 理解傾向
と定義すると、この式は「勉強時間と理解傾向から、試験得点を予測しようとしたモデルである」ということはわかるでしょうか。では、5時間勉強して、理解傾向が3のとき、試験得点は何点予測されるでしょうか、これは、
$$ \hat{y}=a_1 \times 5 + a_2 \times 3 +b, $$となります。$a_1, a_2, b$が未定なので、具体的な数値が求まりません。ここでわかるように、線形回帰モデルによりなんらかの値を予測するには、$a_1, a_2, b$の具体的な値を決めてやらなければいけないのです。これを決めるためには、$y, x_1, x_2$のデータをとにかくたくさん集める必要があります。そして、推定値$\hat{y}$と実測値$y$の差が最小化するような$a_1, a_2, b$を決めてあげるのです。ここでは、これを求める方法について説明します。
データセットの用意¶
本来は実際に、「試験得点」、「勉強時間」、「理解傾向」をアンケートなどでデータを集めなくてはいけませんが、練習ですので、いつも通りダミーデータを作ってみます。ここは今回の本質とは関係ないので、理解せず、コピペで実行して構いません。
import numpy as np
np.random.seed(4) # 擬似乱数シード： 毎回同じ乱数を出す
# 平均
mu = [50.0, 5.0, 5.0] # 試験得点、勉強時間、理解傾向の平均
# 分散共分散行列
# memo: 分散共分散行列は、対称行列、半正定値行列以外のとき、警告が出ます。
s_11 = 70 # 試験得点のばらつき
s_12 = 8 # 試験得点と勉強時間の関連
s_13 = 10 # 勉強時間と理解傾向の関連
s_21 = s_12 # 勉強時間と試験得点の関連
s_22 = 5 # 勉強時間のばらつき
s_23 = 0 # 勉強時間と理解傾向の関連
s_31 = s_13 # 理解傾向と試験得点の関連
s_32 = s_23 # 理解傾向と勉強時間の関連
s_33 = 5 # 理解傾向のばらつき
sigma = [[s_11, s_12, s_13], [s_21, s_22, s_23], [s_31, s_32, s_33]]
# 多変量正規分布による乱数生成
Dat = np.random.multivariate_normal(mu, sigma, 60)
# 小数点以下除去
Dat = np.round(Dat, 0)
# ダミーデータ
Y = Dat[:,0]
X1 = Dat[:,1]
X2 = Dat[:,2]
データを作成しました。せっかくなので、Y:試験得点、X1:勉強時間、X2:理解傾向のデータを可視化してみます。青の濃度が試験得点の推定値となります。勉強して、理解傾向も高いと、試験得点が高そうなのが、なんとなくみて取れます。
import matplotlib.pyplot as plt
plt.figure(figsize=(5, 4))
plt.scatter(X1, X2, c=Y, cmap="Blues", marker="o")
plt.grid(True)
plt.xlim(0, 11)
plt.ylim(0, 11)
plt.xlabel("X1: Studing time")
plt.ylabel("X2: Understanding")
plt.colorbar()
線形回帰モデルのパラメータ獲得¶
Y:試験得点、X1:勉強時間、X2:理解傾向のデータを作成してみました。これを実際に取得できたデータと仮定して、精度の高い試験得点の推定値$\hat{y}$を得るためのパラメータ$a_1, a_2, b$を求めてみます。
from sklearn.linear_model import LinearRegression
# X1とX2の結合
X1 = np.reshape(X1, (60, 1))
X2 = np.reshape(X2, (60, 1))
X=np.concatenate([X1, X2], 1)
ModelLR = LinearRegression() # 線形回帰モデルのセット
ModelLR.fit(X, Y) # パラメータ獲得
これでパラメータ獲得が終わりました。下記の変数に獲得された値が格納されています。
print('a1 = ', ModelLR.coef_[0]) # a1
print('a2 = ', ModelLR.coef_[1]) # a2
print('b = ', ModelLR.intercept_) # b
print('推定モデル: y^ = ', round(ModelLR.coef_[0], 2), ' x1 + ', round(ModelLR.coef_[1], 3) , ' x2 + ', round(ModelLR.intercept_, 3))
係数を見ると、X1:勉強時間よりも、X2:理解傾向の前についている重みの方が、あたいが大きいことが確認できます。すなわち、勉強時間よりも理解傾向の強さの方が、試験得点に大きな影響を与えることがわかります。このように、線形回帰モデルで得たパラメータになんらかの意味付けを行い、仮説検証的に分析する場合もあります。
予測値の取得方法¶
それでは、構築したモデルで予測値を得てみます。上に書いた式を実際に書いてもいいですが、predictを用いることで推定結果を得ることができます。試しに、5時間勉強して、理解傾向が3のときの試験得点の推定値を得てみます。
# 直書きで算出
Res01 = 5 * ModelLR.coef_[0] + 3 * ModelLR.coef_[1] + ModelLR.intercept_
print('試験得点の推定値: ', Res01)
# predictを使用して算出
Res02 = ModelLR.predict([[5, 3]])
print('試験得点の推定値: ', Res02)
小数点以下の解釈が違いますが、同じ値が得られています。どちらでも構いませんが、特に理由があるときを除いて、predictを使うようにすればokです。
切片無しの線形回帰モデルの構築方法¶
妥当かはともかくとして、「勉強時間が0で理解傾向が0のとき、0点を取るに決まっている！」といった仮説を立てる分析者もいるかもしれません。先ほど作成したモデルは、切片が30くらいあるので、最低30点は取れるという推定になってしまいます。入力が0のとき、推定が0になるという仮定に妥当性がある場合には、切片を強制的に0にしたモデルを構築することも、ありといえばありでしょう（例えば、釣りに行った時間を入力、つった魚の数を出力としたモデルは、切片が0にならないとおかしいですよね）。このようなときは、LinearRegressionのfit_intercept引数に、Falseを与えることで、切片0のモデルを実現できます。
ModelLR_noB = LinearRegression(fit_intercept=False) # 線形回帰モデルのセット
ModelLR_noB.fit(X, Y) # パラメータ獲得
print('a1 = ', ModelLR_noB.coef_[0]) # a1
print('a2 = ', ModelLR_noB.coef_[1]) # a2
print('b = ', ModelLR_noB.intercept_) # b
print(' *** 切片が0になりました *** ')
平均二乗誤差（Mean squared error; MSE）によるモデル評価¶
構築した回帰モデルの誤差評価を、平均二乗誤差により行います。実測と推定の差を二乗し、全データの平均をとったものです。0に近いほど良い性能ということになります。簡単な式ですので自分で書いてもいいですが、関数も用意されています。
from sklearn.metrics import mean_squared_error
# 切片ありモデルの誤差評価
Y_pred_01 = ModelLR.predict(X)
Mse_01 = mean_squared_error(Y, Y_pred_01)
print('切片ありモデル MSE = ', Mse_01)
# 切片なしモデルの誤差評価
Y_pred_02 = ModelLR_noB.predict(X)
Mse_02 = mean_squared_error(Y, Y_pred_02)
print('切片なしモデル MSE = ', Mse_02)
切片ありモデルの方が優れています。「勉強せず、理解していない人は、試験で0点だ！」という仮定にはちょっと無理があったみたいです。
平均二乗誤差平方根（Root mean squared error; RMSE）によるモデル評価¶
誤差の二乗の平均値では、二乗されているので若干わかりにくいです。そのため、MSEの平方根を取るという誤差評価方法もあります。これも、0に近いほど性能が良いことを示しています。
Rmse_01 = np.sqrt(Mse_01) # 平方根を算出
print('切片ありモデル RMSE = ', Rmse_01)
Rmse_02 = np.sqrt(Mse_02) # 平方根を算出
print('切片なしモデル RMSE = ', Rmse_02)
当たり前ですが、平方根をとっても大小関係は変わりません。誤差の二乗に対し、その二乗をキャンセルしたわけですから、「1データあたりの平均的な予測のずれ」を示しています。切片ありモデルは平均5点くらい予測を外す、切片なしモデルは平均12点も予測を外してしまうことが確認できます。
決定係数によるモデル評価¶
予測と実測の相関係数を加工したモデル評価方法です。予測と実測が完全に一致している場合、予測と実測の相関係数が1になるという性質を利用しています。1に近いほど、性能が良いことを意味します。この算出には、r2_scoreを使用します。引数に、実測と推定を入力すればokです。純粋に相関係数の2乗をとっている訳ではないので、マイナスになることもあります（マイナスは相当ダメなモデルであることを意味します）。
from sklearn.metrics import r2_score
print('切片ありモデルのR2: ', r2_score(Y, Y_pred_01))
print('切片なしモデルのR2: ', r2_score(Y, Y_pred_02))
散布図による誤差評価¶
色々とモデル評価の指標を説明しましたが、もうちょっと素直に行う方法もあります。横軸を実測、縦軸を予測として、散布図を見てみようというものです。実際にこれはかなり参考になります。
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 4)) # figureの縦横の大きさ
# 切片ありモデル
plt.subplot(1,2,1)
plt.scatter(Y, Y_pred_01, # 横軸、縦軸の値
s=30, # マーカサイズ
c='orange', # マーカの色（red, blue, pink, yellow, orange, blackなど）
marker='s', #マーカの形（.osなど：下記参照）
alpha=0.8, # 塗りつぶしの透明度（0：完全に透明〜１：透明にしない）
linewidths=0.5, # マーカの線の太さ
edgecolors='black') # マーカの線の太さ（マーカの色と同じ。不要の場合は''なしでNone）
plt.title("Intercept Model") # タイトル
plt.xlabel("Real value: y") # 軸名
plt.ylabel("Estimate value: yhat") #軸名
plt.grid(True) #グリッド線（True:引く、False:引かない）
plt.xlim(10, 80) # 横軸最小最大
plt.ylim(10, 80) # 縦軸最小最大
# 切片なしモデル
plt.subplot(1,2,2)
plt.scatter(Y, Y_pred_02, # 横軸、縦軸の値
s=30, # マーカサイズ
c='orange', # マーカの色（red, blue, pink, yellow, orange, blackなど）
marker='s', #マーカの形（.osなど：下記参照）
alpha=0.8, # 塗りつぶしの透明度（0：完全に透明〜１：透明にしない）
linewidths=0.5, # マーカの線の太さ
edgecolors='black') # マーカの線の太さ（マーカの色と同じ。不要の場合は''なしでNone）
plt.title("No Intercept Model") # タイトル
plt.xlabel("Real value: y") # 軸名
plt.ylabel("Estimate value: yhat") #軸名
plt.grid(True) #グリッド線（True:引く、False:引かない）
plt.xlim(10, 80) # 横軸最小最大
plt.ylim(10, 80) # 縦軸最小最大
plt.show() # グラフの表示
左が切片あり、右が切片なしです。横軸が実測、縦軸が推定です。実際の30点を30点と予測できるとベストな訳ですから、いいモデルは横軸と縦軸の値が一致するはずです。左は一致している傾向が見て取れますが、右は30点を20点、70点を60点と予測していたり、若干はずれが大きい傾向が見て取れます。このように、MSEやRMSEに加工するとわからなかった具体的な傾向を、散布図を利用すると見ることができるようになります。
MSEなど、1次元のスカラー値に加工された指標というものは、解釈を容易にするので大変便利です。しかし、大量の情報を有するデータを1次元に圧縮する行為は、確実になんらかの情報を損失されます（例えば、$2 \times 3 = 6$という2次元を1次元に圧縮する変換は一意に定まりますが、その逆変換は一意に定まりません。すなわち、元に戻せないのです。このように単純な変換ですら致命的な情報損失を引き起こします）。これにより、気づかないうちに大事な事実を見落としてしまう場合があります。高次元のデータ解析の場合には指標のみに頼らざるを得ない面も数多くありますが、もし可能なのであれば、指標に頼りつつも全体の傾向を見る癖をつけておいてください。