目次
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個の場合の重回帰分析へと進みます。