目次
回帰直線を1行のコードで描く
前回、回帰直線を得るのにnumpy.profit関数を用いる方法を紹介しました。
>>> a, b = np.polyfit(x1_data, y_data, 1)
として、回帰直線の係数（回帰係数）を取り出した後に
>>> X1 = np.arange(x_min, x_max, 1)
>>> Y = a * X1 + b
>>> plt.plot(X1, Y, c='Red')
として回帰直線を描くという方法です。今回はこれを1行で済ませる方法から紹介していきます。
>>> plt.plot(x, np.poly1d(np.polyfit(x, y, 1))(x), label='d=1',c='Red')
これで上記4行のコードと同じように回帰直線を描くことができます。したがって、すべてのコードはさらに簡潔にできます。
前回と同じデータを用いて回帰分析を行います。
表1：身長と体重のデータ
プログラム「kaiki1-1py」回帰分析（1次式）プログラム
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>>
>>> # STEP1. データを1次元配列として格納する
>>> x = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185])
>>> y = np.array([ 53, 59, 64, 75, 72, 74, 83, 87, 84, 90])
>>>
>>>
>>> # plt.figure(dpi=200)
>>> # plt.scatter(x, y, c='Blue') # 表1のデータをプロット
>>> plt.plot(x, np.poly1d(np.polyfit(x, y, 1))(x), label='d=1',c='Red') # 回帰直線を描画
>>> # plt.xlabel('x : Height (cm)')
>>> # plt.ylabel('y : Weight (kg)')
>>> plt.show()
このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1WOduM-_33JycPdbSQTz6HwiG0LQi_JX3/view?usp=sharing
プログラム実行結果
kaiki1-1.pyを実行してみましょう。次のグラフが出力されます。
図1：「kaiki1-1.py」の実行結果
これならラクです。matplotlibを用いたグラフ描画は高機能の分、手間がかかります。特に、
>>> X1 = np.arange(x_min, x_max, 1)
のようにxの値をつくってあげるところが億劫です。そうしなくてもいいということです。
実際これでは回帰分析らしくないので、表のデータプロットと軸にラベルをつけてみます。kaiki1-1.pyの後半にある4つのコメントアウトを除きます。次がその結果です。
図2：「kaiki1-1.py」（コメントアウト除去）の実行結果
重回帰分析
いよいよ本丸の重回帰分析に進んでいきます。説明変数が2個以上のデータに対して回帰式を求めます。これまで回帰分析に用いた身長、体重に足の長さを説明変数に加えます。
すると身長x1と足の長さx2の説明変数と目的変数（体重y）の間に1次式（線形）の関係があると仮定して、回帰式を求めることができます。
回帰式を図示すると平面となります。回帰平面と呼ばれます。
図3：重回帰分析
Pythonで重回帰分析を行う方法はいくつもありますが、回帰分析（単回帰分析）のように簡単なコードを書くことは難しいです。
しかし次に紹介するプログラムはグラフ以外のライブラリは、NumPyだけで3行のコードで回帰係数を求めることができるものです。
プログラム「jukaiki1.py」重回帰分析プログラム
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D #3Dplot
>>> import japanize_matplotlib
>>>
>>> # データ定義
>>> 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(dpi=200)
>>> 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は透過度
>>> plt.show()
このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1xclNcNzlbMrQ7ExzFSkxtSqGGr8ioYY9/view?usp=sharing
プログラム実行結果
jukaiki1.pyを実行してみましょう。次の3Dグラフが出力されます。
回帰係数を求めるためのコードは次の3行です。
>>> ones = np.ones(len(x1))
>>> X = np.array([ones, x1, x2]).T # 説明変数行列
>>> theta = np.linalg.inv(X.T @ X) @ X.T @ y
回帰係数はベクトルデータthetaとして得られます。この3行のコードの中に含まれるのが最小2乗法と線形代数学です。
次回はいよいよ回帰分析の核心へと迫っていきます。
