« 戻る次へ »

線形回帰モデル

モデル2：2次式を加える

モデル

前のページの単回帰分析では

\begin{equation} y = a + bx \end{equation}

という単純な一次式を用いて，回帰直線を求めた．ここでは2次の項を加えた

\begin{equation} y = w_0 + w_1x + w_2 x^2 \end{equation}

を用いて回帰式を求めてみよう．この場合は得られる回帰式が直線にはならない．しかしながら，\(f_1(x) = x\)，\(f_2(x) = x^2\) とおくと

\begin{equation} y = w_0 + w_1f_1(x) + w_2f_2(x) \end{equation}

と書くことができる．つまり，これは重み \(w_0, w_1, w_2\) の線形モデルになっているので，（得られる回帰式は曲線であっても）線形回帰モデルと呼ばれる．

目次に戻る

関数を定義する

まず，必要なライブラリを読み込んだ後に回帰式と残差2乗和を取得する関数を定義する．回帰式（15行目）に2次の項が増えただけで，残差2乗和のコードは1次式の場合と全く同じである．

import numpy as np
import pandas as pd
import scipy.optimize as optimize # 最適化
import matplotlib.pyplot as plt

# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')

"""
関数の定義
"""
def my_func(w, x):
    y = w[0] + w[1] * x + w[2] * x ** 2
    return y

"""
残差2乗和を求める関数
Residual Sum-of-Squares
"""
def get_rss(w, x, y):
    y_pred = my_func(w, x)
    error = (y - y_pred)**2
    return np.sum(error)

目次に戻る

CSVファイルを開いて最適化を行う

1次式の場合とほぼ同じ方法で最適化を行うことができる．僅かな違いは推定するパラメータ数が2から3に変化したことである．

# CSV ファイルを読み込む
url = "https://github.com/rinsaka/sample-data-sets/blob/master/lr.csv?raw=true"
df = pd.read_csv(url)

# NumPy 配列に変換する
x_data = df.loc[:, 'x'].values
y_data = df.loc[:, 'y'].values

# ネルダーミード法による最適化を行う
w = np.array([0.1, 0.1, 0.1]) # 初期値を設定
results_nm = optimize.minimize(get_rss, w, args=(x_data, y_data), method='Nelder-Mead')
print(results_nm)

 final_simplex: (array([[ 2.23792585,  1.99417077, -0.25169293],
       [ 2.23800909,  1.99414947, -0.2516921 ],
       [ 2.23786423,  1.99423825, -0.25169995],
       [ 2.23800163,  1.99413505, -0.25168841]]), array([148.915214  , 148.91521402, 148.91521403, 148.91521405]))
           fun: 148.91521400211883
       message: 'Optimization terminated successfully.'
          nfev: 263
           nit: 143
        status: 0
       success: True
             x: array([ 2.23792585,  1.99417077, -0.25169293])

上の結果を確認すると，目的関数すなわち残差2乗和が 148.92 となり，1次式の 244.84 から大きく改善されたことがわかる．また推定パラメータ数が2から3に増えたことで，問題が僅かに難しくなったことから，最適化での繰り返し回数 nit (Number of iterations performed by the optimizer) が 62 から 143 に増加していることもわかる．

目次に戻る

グラフを描く

最後に回帰式を散布図に重ねて描いてみる．下のコードは1次式の場合から変更の必要はありません．

# 最適解を使って回帰直線のデータを作成する
x_plot = np.linspace(0, 10, 100)
y_pred = my_func(results_nm["x"], x_plot)

# グラフを描く
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data, label="data")
ax.plot(x_plot, y_pred, label='model 2')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.legend()
# plt.savefig('lr-model2.png', dpi=300, facecolor='white')
plt.show()

目次に戻る

« 戻る次へ »