コラム 人と星とともにある数学 数学

NumPyで行列 回帰分析 その8 最小2乗法|Pythonで数学を学ぼう! 第29回

重回帰分析 回帰係数の算出

回帰分析シリーズは最終回です。これで重回帰分析を計算する準備ができました。さっそく、身長と足の長さと体重についての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

  • この記事を書いた人
  • 最新記事

桜井進(さくらいすすむ)様

1968年山形県生まれ。 サイエンスナビゲーター®。株式会社sakurAi Science Factory 代表取締役CEO。 (略歴) 東京工業大学理学部数学科卒、同大学大学院院社会理工学研究科博士課程中退。 東京理科大学大学院非常勤講師。 理数教育研究所Rimse「算数・数学の自由研究」中央審査委員。 高校数学教科書「数学活用」(啓林館)著者。 公益財団法人 中央教育研究所 理事。 国土地理院研究評価委員会委員。 2000年にサイエンスナビゲーターを名乗り、数学の驚きと感動を伝える講演活動をスタート。東京工業大学世界文明センターフェローを経て現在に至る。 子どもから大人までを対象とした講演会は年間70回以上。 全国で反響を呼び、テレビ・新聞・雑誌など様々なメディアに出演。 著書に『感動する!数学』『わくわく数の世界の大冒険』『面白くて眠れなくなる数学』など50冊以上。 サイエンスナビゲーターは株式会社sakurAi Science Factoryの登録商標です。

あわせて読みたい

-コラム, 人と星とともにある数学, 数学
-,