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

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

NumPyのpolyfit関数で回帰分析

数値計算ライブラリ「NumPy」は回帰分析のためのprofit関数を用意しています。

さっそく使ってみましょう。次の表1の10組のデータについて回帰分析を行います。

表1:身長と体重のデータ

  • STEP1. データを1次元配列として格納する。
  • STEP2. STEP1の1次元配列をprofit関数に渡すことで回帰係数が得られる。
  • STEP3. 回帰直線のグラフを描画する。

プログラム「kaiki1.py」回帰分析(1次式)プログラム

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>>
>>> # STEP1. データを1次元配列として格納する
>>> x1_data = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185])
>>> y_data = np.array([ 53, 59, 64, 75, 72, 74, 83, 87, 84, 90])
>>>
>>> # STEP2. STEP1の1次元配列をprofit関数に渡すことで回帰係数が得られる
>>> a, b = np.polyfit(x1_data, y_data, 1)
>>> print("回帰式:y =", a, "x1 +", b)
>>>
>>> # STEP3. 回帰直線のグラフを描画
>>> margin = 10
>>> x_min, x_max = x1_data.min() - margin, x1_data.max() + margin
>>>
>>> # 回帰式の直線を描くためのデータを作る
>>> X1 = np.arange(x_min, x_max, 1)
>>> Y = a * X1 + b
>>>
>>> plt.figure(dpi=200)
>>> plt.scatter(x1_data, y_data, c='Blue') # 表1のデータをプロット
>>> plt.plot(X1, Y, c='Red') # 回帰直線を描画
>>> plt.xlabel('x : Height (cm)')
>>> plt.ylabel('y : Weight (kg)')
>>> plt.show()

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

回帰直線y=ax+bの回帰係数aとbを求めるコードはSTEP2の1行だけです。

>>> a, b = np.polyfit(x1_data, y_data, 1)

このあとデータ・プロットと回帰直線を描くプログラムが続きます。体裁良くグラフを描くために少し長くなっていますが、回帰直線だけを求めるのであればSTEP2までの数行で済みます。

STEP1.  データを1次元配列として格納する

NumPyのarray関数を用いて表1のデータを1次元配列として定義します。

>>> import numpy as np
>>>
>>> x1_data = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185])
>>> y_data = np.array([ 53, 59, 64, 75, 72, 74, 83, 87, 84, 90])

x1_dataが身長、y_dataが体重のデータです。

STEP2. STEP1の1次元配列をprofit関数に渡すことで回帰係数が得られる

次のようにpolyfit関数を実行することで回帰係数のaとbが得られます。

>>> a, b = np.polyfit(x1_data, y_data, 1)

第3引数の1は回帰式を1次式にする指定です。

STEP3. 回帰直線のグラフを描画する

グラフ描画のためのライブラリ「Matplotlib」を用いて、得られた回帰式を可視化します。ここでのポイントは、回帰直線のxの範囲をうまく決めてあげることです。

描画する直線を10組のデータの両端のデータを結ぶ線分とするよりも、両端にマージンをつけることで回帰直線らしくなります。

x1_dataの最小値と最大値を取得しそれらにマージンをつけることで、10個の点の領域よりも大きな範囲の回帰直線を描きます。

>>> margin = 10
>>> x_min, x_max = x1_data.min() - margin, x1_data.max() + margin

変数marginの値を変えることで伸び幅の調整ができます。

プログラム実行結果

kaiki1.pyを実行してみましょう。次のように回帰式がコンソール出力されます。

回帰式:y = 1.546155406776225 x1 + -197.868736051938

この回帰式を使って、身長が160cmの体重を推測できます。x1に160を代入して計算してみます。

>>> 1.546155406776225 * 160 -197.868736051938

結果は「49.51612903225802」となり、約50kgであると推測できます。回帰分析を使うことで推測・予測ができることがデータサイエンスの現場で使われている理由です。

そして図1のグラフも合わせて出力されます。これを見ると、回帰直線が10個の点の中央あたりを通ることがわかります。実にうまい具合に回帰式が得られていることが直感的にわかります。

この後、最小2乗法によりこのような直線が求められることをじっくり見ていきます。

図1:「kaiki1.py」の実行結果

polyfit関数で2次式の回帰式を求める

polyfit関数の第3引数を2として回帰式を2次式 y = a*x1² + b*x1 + c
としてみましょう。

x1_dataの2乗に誤差を加えた値をy_dataとする実験データを20組準備します。NumPyライブラリに用意されたnp.random.randn関数が便利です。平均0、分散1(標準偏差1)の正規分布(標準正規分布)に従う乱数を返す関数です。

>>> x1_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19])
>>> y_data = np.round(x1_data ** 2 + np.random.randn(20) * 10)

polyfit関数の第3引数を2とすると、回帰係数は3つ返されるので次のようにします。

>>> a, b, c = np.polyfit(x1_data, y_data, 2)

それに合わせて回帰式の出力も次のように変更します。

>>> print("回帰式:y =", a, "x1² +", b ,"x1 +", c)

プログラム「kaiki2.py」回帰分析(2次関数)プログラム

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>>
>>> x1_data = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19])
>>> y_data = np.round(x1_data ** 2 + np.random.randn(20) * 10)
>>> print(y_data)
>>>
>>> # STEP2. STEP1の1次元配列をprofit関数に渡すことで回帰係数が得られる
>>> a, b, c = np.polyfit(x1_data, y_data, 2)
>>> print("回帰式:y =", a, "x1² +", b,"x1 +", c)
>>>
>>> # STEP3. 回帰式のグラフを描画
>>> margin = 2
>>> x_min, x_max = x1_data.min() - margin, x1_data.max() + margin
>>>
>>> # 回帰式を描くためのデータを作る
>>> X1 = np.arange(x_min, x_max, 1)
>>> Y = a * X1**2 + b * X1 + c
>>>
>>> plt.figure(dpi=200)
>>> plt.scatter(x1_data, y_data, c='Blue') # データをプロット
>>> plt.plot(X1, Y, c='Red') # 回帰式を描画
>>> plt.xlabel('x1')
>>> plt.ylabel('y')
>>> plt.show()

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

プログラム実行結果

kaiki2.pyを実行してみましょう。次のようにy_dataと回帰式がコンソール出力されます。y_dataはx1_dataの値を2乗した値をもとにしたものなのでx1²の係数aは1に近い0.91です。

[ 3. -9. 1. 3. 27. 27. 44. 48. 54. 89. 108. 147. 146. 168. 198. 219. 256. 297. 319. 364.]
回帰式:y = 0.9175780359990887 x1² + 1.7171451355661895 x1 + -4.183766233766231

そして図2のグラフも合わせて出力されます。これを見ると、回帰式の2次式(放物線)は20個のデータの点にうまくフィットしていることが見て取れます。

図2:2次式の回帰式

以上が説明変数1個の場合の回帰分析で単回帰分析とも呼ばれます。次回は説明変数2個の場合の重回帰分析へと進みます。

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

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

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

あわせて読みたい

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