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

NumPyで行列 回帰分析 その1|Pythonで数学を学ぼう! 第22回

線形代数はタイヘン

前回前々回、ベクトルを取り上げました。ベクトルの次ということで行列をはじめていきます。ベクトルと行列は線形代数学という大きな数学へと結実します。

1冊にまとめられた線形代数のテキストを自力で終わらせるには、相当な根性と実力が必要になります。

前回取り上げたコサイン類似度とは、データをベクトル化することでn次元空間における2つのn次元ベクトルのなす角(実際にはそのコサイン)を計算することで比較するというアイディアです。

ベクトルは3次元が4次元に次元が大きくなってもさほど計算が面倒なことはありません。行列の場合、次元が1つ大きくなるだけで行列式・逆行列の計算はタイヘンになります。スラスラと手計算できるのは2次元どまりです。3次元、4次元の行列の計算ですら容易に手計算できません。

手計算ができないまま行列式・逆行列・連立分程式の仕組みをものにすることはできません。計算量の膨大さこそが行列にとって最大の壁です。

そこでPythonの登場です。

ベクトルと同様にNumPyライブラリを用いて、行列は容易に定義でき、NumPyに用意された関数を用いることで行列式・逆行列・連立分程式といった行列の様々な計算を1行のコードで済ませることができます。

これまでの連載同様にPythonで線形代数を学ぶことはもちろん可能です。しかし、それでも線形代数のテキストの章立ての長さは長いまま変わることはありません。線形代数とは、連立方程式の解法をいかにシステマティックに行うかという問題が原点にあります。

はたして、行列という新しい概念が誕生し、そこに見事な美しい法則が発見されていくことになったのです。行列の美しい法則の結晶が線形代数・線形代数学・線型代数・線型代数学(用語の組み合わせは4通り)です。

「NumPyで行列」として行列の計算を紹介していきます。行列は面白そうと思ってもらえるプログラムを用意しました。それが回帰分析です。多数のデータをもっともうまく近似する直線(曲線)を求めることで、予測シミュレーションを可能にする回帰分析はデータサイエンスの基本です。

Pythonによる回帰分析

さっそく、Pythonによる回帰分析のコードをみていただきましょう。

用意したデータは10人分の身長(cm)と体重(kg)です。まず、これらをnp.array()関数を用いて説明変数ベクトル(配列)xと目的変数(体重)ベクトルyを定義します。

次に、説明変数(身長)から説明変数行列Xを定義します。はたして、行列Xに対して転置行列、逆行列の計算をして最後に目的変数(体重)ベクトルyをかけることで回帰式の係数(回帰係数)が得られます。この回帰係数の計算コードはたった1行です。

最後に計算結果を出力します。

(身長と体重)回帰分析プログラム 「regression1.py」

>>> import numpy as np
>>>
>>>  # データ定義
>>>  x = np.array([170, 166, 157, 164, 157, 160, 163, 160, 153, 165]) # 説明変数 身長
>>>  y = np.array([57, 53, 47, 52, 47, 53, 48, 48, 43, 57]) # 目的変数 体重
>>>
>>>  # 説明変数行列Xの定義
>>>  ones = np.ones(len(x))
>>>  X = np.array([ones, x]).T
>>>
>>>  # 回帰係数の計算
>>>  theta = np.linalg.inv(X.T @ X) @ X.T @ y
>>>
>>>  print(f'説明変数行列 X =\n {X}')
>>>  print(f'回帰係数 θ_0 = {theta[0]}')
>>>  print(f'回帰係数 θ_1 = {theta[1]}')
>>>  print(f'回帰式 y = {theta[0]} + {theta[1]} * x')

このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1U1m17D--NUpoZ7wTijzbXlHCPfLb6Ibv/view?usp=sharing

実行結果

次がregression1.pyの実行結果です。最後に回帰式がxの1次式として出力されています。

これより身長が175cmであれば体重はいくつになるか計算できます。回帰式のxに175を代入して計算してみます。結果は約61.3kgです。

>>>  -78.77006507593245 + 0.8004338394793968 * 175
>>>  61.30585683296201

この例のように、説明変数が1つの場合には回帰式はxの1次式となるのでグラフは直線です。これを回帰直線といいます。

10人分のデータの散布図と回帰直線を描くプログラムは回帰分析プログラム「regression1.py」にmatplotlibライブラリを用いたグラフ描画部分を付け加えたのが回帰分析プログラム「regression2.py」です。

(身長と体重)回帰分析プログラム「regression2.py」

>>>  import numpy as np
>>>  import matplotlib.pyplot as plt
>>>
>>>
>>>  # データ定義
>>>  x = np.array([170, 166, 157, 164, 157, 160, 163, 160, 153, 165]) # 説明変数 身長
>>>  y = np.array([57, 53, 47, 52, 47, 53, 48, 48, 43, 57]) # 目的変数 体重
>>>
>>>
>>>  # 説明変数行列Xの定義
>>>  ones = np.ones(len(x))
>>>  X = np.array([ones, x]).T
>>>
>>>  # 回帰係数の計算
>>>  theta = np.linalg.inv(X.T @ X) @ X.T @ y
>>>
>>>  print(f'説明変数行列 X =\n {X}')
>>>  print(f'回帰係数 θ_0 = {theta[0]}')
>>>  print(f'回帰係数 θ_1 = {theta[1]}')
>>>  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 = θ_0 + θ_1 * xx # 回帰直線 グラフ用変数:xx,yy
>>>  yy = + theta[0] + theta[1] * xx # 回帰直線 グラフ用変数:xx,yy
>>>
>>>  plt.scatter(x, y, c='Blue') # データを青丸でプロット
>>>  plt.plot(xx, yy, c='Red') # 回帰直線を赤線でプロット
>>>  plt.xlabel('x:height(cm)')
>>>  plt.ylabel('y:weight(kg)')
>>>  plt.show()

このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1ZE5B6X7Sp5ugeBgvb7DZgyEsjsZIE_mg/view?usp=sharing

実行結果

numpy.polyfit関数を使った回帰分析

NumPyライブラリには回帰分析のためのpolyfit関数が用意されています。これを用いればコードはもっと簡単になります。
>>>  theta_1 , theta_0 = np.polyfit(x, y, 1)

第3引数を1とすることで回帰式を1次式と指示します。
回帰係数はxの係数(theta_1)、切片(theta_0)の順序で出力されます。

(身長と体重)回帰分析プログラム「regression3.py」

>>>  import matplotlib.pyplot as plt
>>>
>>>  # データ定義
>>>  x = np.array([170, 166, 157, 164, 157, 160, 163, 160, 153, 165]) # 説明変数 身長
>>>  y = np.array([57, 53, 47, 52, 47, 53, 48, 48, 43, 57]) # 目的変数 体重
>>>
>>>  # 回帰係数の計算
>>>  theta_1 , theta_0 = np.polyfit(x, y, 1)
>>>
>>>  print(f'回帰係数 θ_0 = {theta_0}')
>>>  print(f'回帰係数 θ_1 = {theta_1}')
>>>  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 = θ_0 + θ_1 * xx # 回帰直線 グラフ用変数:xx,yy
>>>  yy = + theta_0 + theta_1 * xx # 回帰直線 グラフ用変数:xx,yy
>>>
>>>  plt.scatter(x, y, c='Blue') # データを青丸でプロット
>>>  plt.plot(xx, yy, c='Red') # 回帰直線を赤線でプロット
>>>  plt.xlabel('x:height(cm)')
>>>  plt.ylabel('y:weight(kg)')
>>>  plt.show()

このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1X5KCOubqdvUyl5NPIBMthXbqdfTFfjwj/view?usp=sharing

実行結果

回帰係数は先の結果とほぼ同じ結果になることが確認できます。

回帰分析の鍵「最小2乗法」

いかがでしたでしょうか。身長、体重のデータを近似する回帰直線を求めるコードはこのように数行でできてしまいます。それもこれも次の1行のおかげです。

>>>  # 回帰係数の計算
>>>  theta = np.linalg.inv(X.T @ X) @ X.T @ y

なぜこのような行列の計算(行列の積、転置行列、逆行列)によって回帰係数が得られるのか。今後の連載でじっくりみていきたいと思います。そこで展開されるのは驚きの連続です。

時代は今から二百年前のフランスに遡ります。1789年7月14日、バスティーユ襲撃が起きます。所謂、フランス革命です。フランス革命の中で生まれたのがメートルです。

登場人物は全員数学者、ラグランジュ、ラプラス、ドランブルそしてガウス。メートルという統一単位をつくるために数学は総動員されました。はたして、彼らの闘いの中から新しい数学──最小2乗法が誕生しました。最小2乗法のおかげで回帰係数の計算が可能になります。

上記の1行の数式は、行列とメートルとフランス革命という壮大な物語から生みだされました。次回以降をご期待ください。

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

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

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

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

© 2022 株式会社インフォマティクス