目次
重回帰分析 回帰係数の算出
回帰分析シリーズは今回が最終回です。これで重回帰分析を計算する準備ができました。さっそく身長と足の長さと体重についての10人分のデータに対して回帰係数を計算してみましょう。
説明変数は2つ、身長x1と足の長さx2、目的変数は体重yです。説明変数行列Xをつくっていきます。
Xの1列目は、データ数のn個だけ1となるので、ones = np.ones(len(x1))として成分が1である1次元配列をつくっておきます。1次元配列ones、x1、x2を並べて行列をつくりますが、この行列は型が図11にあるような説明変数行列Xではありません。これを転置することで説明変数行列Xが完成します。
念のためprint(X)として出力して確かめてみます。あとは回帰係数ベクトルθを求める前回28回の公式(4)をpythonコードに置き換えていきます。逆行列はnp.array.inv()関数を用います。
プログラム「multiple_regression_coefficient.py」｜重回帰分析回帰係数算出プログラム
>>> # multiple_regression _coefficient.py
>>> import numpy as np
>>>
>>> x1 = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185]) # 身長
>>> x2 = np.array([23, 24, 25, 27, 25, 24, 28, 29, 28, 29]) # 足の長さ
>>> y = np.array([ 53, 59, 64, 75, 72, 74, 83, 87, 84, 90]) # 体重
>>>
>>> ones = np.ones(len(x1))
>>> X = np.array([ones, x1, x2]).T
>>> print(X)
>>>
>>> theta = np.linalg.inv(X.T @ X) @ X.T @ y
>>>
>>> print(f'θ0 = {theta[0]}')
>>> print(f'θ1 = {theta[1]}')
>>> print(f'θ2 = {theta[2]}')
プログラム実行結果
「python multiple_regression _coefficient.py」と入力することで実行します。 次のように説明変数行列Xと3つの回帰係数が出力されます。
このプログラムは単回帰分析にも使うことができます。説明変数x2と回帰係数θ2部分を削除する（3箇所変更）
実行結果は第25回の回帰分析の結果
回帰式：y = 1.546155406776225 x1 + -197.868736051938
と同じであることが確認できます。
matplotlibによるデータと結果の可視化
最後に回帰分析と重回帰分析についてデータと結果をグラフに表してみましょう。
グラフによりデータの傾向を概観することができます。これまでのプログラムにグラフ描画のコードを付け足したのがプログラム「regression_analysis.py」です。
matplotlibライブラリとmatplotlib日本語表示ライブラリをインポートします。データのプロットはplt.scatter(x, y, c='Blue')だけで済みます。回帰直線を描くためにいくつか準備します。
まず説明変数xのデータからxの範囲を決めます。回帰直線を描くための変数が必要です。xとyはデータ定義に使われているのでxxとyyとします。
プログラム「regression_analysis.py」回帰分析｜回帰直線表示プログラム
>>> # regression_analysis.py
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import japanize_matplotlib
>>>
>>>
>>> x = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185]) # 身長
>>> y = np.array([52, 60, 63, 65, 78, 70, 72, 85, 90, 94]) # 体重
>>>
>>> theta_1 = np.cov(x, y, bias=True)[0][1]/np.var(x)
>>> theta_0 = np.mean(y) - theta_1 * np.mean(x)
>>>
>>> print(f'回帰直線の傾きθ_1 = {theta_1}')
>>> print(f'回帰直線の切片θ_0 = {theta_0}')
>>> print(f'回帰直線 y = {theta_0} + {theta_1}x')
>>>
>>>
>>> # データと回帰直線 グラフ描画
>>> margin = 1
>>> x_min, x_max = x.min()-margin, x.max()+margin # xの範囲
>>> xx = np.arange(x_min, x_max, 1) # 回帰直線グラフ用変数:xx
>>>
>>> yy = + theta_0 + theta_1 * xx # 回帰直線 グラフ用変数:xx,yy
>>>
>>> plt.scatter(x, y, c='Blue') # データを青丸でプロット
>>> plt.plot(xx, yy, c='Red') # 回帰直線を赤線でプロット
>>> plt.xlabel('x:身長（cm）')
>>> plt.ylabel('y:体重（kg）')
>>> plt.show()
プログラム実行結果
「python regression_analysis.py」と入力することで実行します。次のように回帰係数と回帰直線の式、グラフにデータ（青丸）と回帰直線（赤線）が表示されます。
3Dを表示するためにはmatplotlibライブラリの他にmpl_toolkits.mplot3dライブラリをインポートします。matplotlib日本語表示ライブラリもインポートします。データのプロットはax.scatter3D(x1, x2, y, color='Blue')とします。
回帰平面を描く場合も回帰直線の場合と同じようにいくつか準備が必要です。まず説明変数x1とx2のデータからx1とx2の範囲を決めます。回帰平面を描くための変数が必要です。x1とx2とyはデータ定義に使われているのでxx1とxx2とyyとします。
プログラム「multiple_regression_analysis.py」重回帰分析｜データ・回帰直線表示プログラム
>>> # multiple_regression_analysis.py
>>> import numpy as np
>>> import japanize_matplotlib
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D #3Dplot
>>>
>>> # データ定義
>>> x1 = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185]) # 説明変数 身長
>>> x2 = np.array([23, 24, 25, 27, 25, 24, 28, 29, 28, 29]) # 説明変数 足の長さ
>>> y = np.array([52, 60, 63, 65, 78, 70, 72, 85, 90, 94]) # 目的変数 体重
>>>
>>> ones = np.ones(len(x1))
>>> X = np.array([ones, x1, x2]).T # 説明変数行列
>>> print(f'説明変数行列X=\n{X}')
>>>
>>> # 回帰係数ベクトルの計算
>>> theta = np.linalg.inv(X.T @ X) @ X.T @ y
>>>
>>> print(f'回帰係数 θ_0 = {theta[0]}')
>>> print(f'回帰係数 θ_1 = {theta[1]}')
>>> print(f'回帰係数 θ_2 = {theta[2]}')
>>> print(f'回帰平面 y = {theta[0]} + {theta[1]}*x1 + {theta[2]}*x2')
>>>
>>>
>>> # データと回帰平面 グラフ描画
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111, projection='3d')
>>>
>>> ax.scatter3D(x1, x2, y, color='Blue') # データを青丸でプロット
>>> ax.set_xlabel("x1:身長（cm）")
>>> ax.set_ylabel("x2:足の長さ（cm）")
>>> ax.set_zlabel("y:体重（kg）")
>>>
>>> # x1とx2の範囲
>>> mesh_size = 1
>>> margin = 0.01
>>> x1_min, x1_max = x1.min()-margin, x1.max()+margin
>>> x2_min, x2_max = x2.min()-margin, x2.max()+margin
>>> x1_range = np.arange(x1_min, x1_max, mesh_size)
>>> x2_range = np.arange(x2_min, x2_max, mesh_size)
>>> xx1, xx2 = np.meshgrid(x1_range, x2_range)
>>>
>>> yy = (theta[0] + theta[1] * xx1 + theta[2] * xx2) # 回帰平面
>>>
>>> ax.plot_surface(xx1, xx2, yy, color='red', alpha=0.5) # alphaは透過度
プログラム実行結果
「python multiple_regression_analysis.py」と入力することで実行します。 次のように説明変数行列と回帰係数と回帰平面の式、グラフにデータ（青丸）と回帰平面（赤線）が出力されます。
説明変数が3つ以上ある重回帰分析
以上の重回帰分析は説明変数2つの場合です。結果が回帰平面として3D可視化できるため重回帰分析の理解にはうってつけです。
説明変数が3つ以上ある重回帰分析では回帰式を可視化することはできませんが、リスト4の説明変数をx3、x4、…と増やすことで回帰式を得ることができます。身長、足の長さに加えて年齢を説明変数x3としてみます。
変更する箇所は次の4行です。
>>> x3 = np.array([43, 20, 54, 27, 33, 19, 28, 40, 29, 50]) # 説明変数 年齢
>>> X = np.array([ones, x1, x2, x3]).T # 説明変数行列
>>> print(f'回帰係数 θ_3 = {theta[3]}')
>>> print(f'回帰式 y = {theta[0]} + {theta[1]}*x1 + {theta[2]}*x2 + {theta[3]}*x3')
プログラム「multiple_regression_analysis_3.py」｜重回帰分析（3変数）プログラム
>>> # multiple_regression_analysis_3.py
>>> import numpy as np
>>> import japanize_matplotlib
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D #3Dplot
>>>
>>> # データ定義
>>> x1 = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185]) # 説明変数 身長
>>> x2 = np.array([23, 24, 25, 27, 25, 24, 28, 29, 28, 29]) # 説明変数 足の長さ
>>> x3 = np.array([43, 20, 54, 27, 33, 19, 28, 40, 29, 50]) # 説明変数 年齢
>>> y = np.array([52, 60, 63, 65, 78, 70, 72, 85, 90, 94]) # 目的変数 体重
>>>
>>> ones = np.ones(len(x1))
>>> X = np.array([ones, x1, x2, x3]).T # 説明変数行列
>>> print(f'説明変数行列X=\n{X}')
>>>
>>> # 回帰係数ベクトルの計算
>>> theta = np.linalg.inv(X.T @ X) @ X.T @ y
>>>
>>> print(f'回帰係数 θ_0 = {theta[0]}')
>>> print(f'回帰係数 θ_1 = {theta[1]}')
>>> print(f'回帰係数 θ_2 = {theta[2]}')
>>> print(f'回帰係数 θ_3 = {theta[3]}')
>>> print(f'回帰式 y = {theta[0]} + {theta[1]}*x1 + {theta[2]}*x2 + {theta[3]}*x3')
プログラム実行結果
「python multiple_regression_analysis_3.py」と入力することで実行します。 次のように説明変数行列、回帰係数、回帰式が出力されます。
以上のPythonコードは「Pythonで数学を学ぼう！第29回.ipynb」としてまとめてあります。
このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1VUVzP1OpWdniGqeW-Dkiif3NLn4soktp/view?usp=share_link
