回帰③・基底関数

Image from Gyazo

基底関数

これまでは回帰式として線形関数(一次式)を採用してきました.

$$f(x) = w_0 + w_1 x $$

しかし,下記のようなデータを対象とする場合は,直線ではなく曲線を用いた方が,傾向を正しく表現できそうです. このデータは,愛知県衛生研究所が公開しているインフルエンザの発症数です. 変数$x$は 経過週 ,変数$y$は インフルエンザの報告数(定点当たり) を表しています.

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]) #週
y = np.array([0.26,0.25,0.19,0.34,0.89,1.78,2.51,4.81,10.97,18.02,19.28,32.42,62.12,62.18,62.60,49.18,32.57,25.93,21.17,14.68]) #報告数

Image from Gyazo

このような場合は,基底関数 を利用します. 線形関数の$x$を,任意の$\phi(x)$という基底関数に置き換えるという方法です. 基底関数には,どんな関数を用いても良いのですが,よく用いられるのは, 多項式基底ガウス基底 と呼ばれる関数です. 今回は,これら 基底関数 を利用した回帰について学びます.

ノートブックの作成

Anacondaを導入していない場合は,まず最初に数値解析のためのライブラリ SciPypipでインストールする必要があります. インストールしたフォルダで,PowerShellを起動し,下記のコマンドを実行してください.

> .\Scripts\pip install scipy

Jupyter Notebook を起動し,新規にノートブックを作成してください. ノートブックのタイトルは Notebook5 とします. ノートブックの作成方法は第1回の資料を参照してください. また,numpymatplotlib.pyplot を導入しておいてください.

import numpy as np
import matplotlib.pyplot as plt

多項式基底

多項式基底は下記の式で表されます. ここで,$i$は組み合わせる基底関数の番号を表します.

$$ \phi_i(x) = x^i $$

回帰式

4つの基底関数を組み合わせると,回帰式は下記のように表されます (2つの基底関数を組み合わせた場合は線形関数と一致する). この回帰式では,$w_0$,$w_1$,$w_2$,$w_3$の4つのパラメータの最適値を求めることが目的です.

$$ f(x) = w_0 \cdot \phi_0(x) + w_1 \cdot \phi_1(x) + w_2 \cdot \phi_2(x) + w_3 \cdot \phi_3(x) $$

$$ = w_0 \cdot x^0 + w_1 \cdot x^1 + w_2 \cdot x^2 + w_3 \cdot x^3 $$

$$ = w_0 + w_1 \cdot x + w_2 \cdot x^2 + w_3 \cdot x^3 $$

この回帰式を関数 fx として定義しておきます.

def fx(w, x): # 回帰式の関数定義
  y_ = w[0]  + w[1] * x + w[2] * x**2  + w[3] * x**3
  return y_

目的関数

目的関数$E$は,これまでと同じ2乗誤差で与えます.

$$E=\sum_{n=0}^{N-1} (y_n - f(x_n))^2$$

この目的関数を関数 gx として定義しておきます.

def gx(w, x, y): # 目的関数(2乗誤差)の関数定義
  y_ = fx(w, x)
  e = np.sum((y - y_) ** 2)
  return e

数値解の導出

数値解の導出には, SciPy という数値解析のためのライブラリを導入します. 同ライブラリの minimize 関数は,指定した関数を最小化するようなパラメータを導出してくれます. 引数には,目的関数 gx ,パラメータの初期値 w_init , 最適化の対象とならないパラメータ (x, y) , 最適化アルゴリズム powell を指定します.

[In:]

from scipy.optimize import minimize # ライブラリの導入

w_init = np.array([0, 0, 0, 0]) # パラメータの初期値
result = minimize(gx, w_init, args=(x, y), method="powell") # Powell法で最適化
print(result.x) # 導出されたパラメータ

[Out:]

[ 20.22701023 -14.62576848   2.32622593  -0.08134343]

この結果,$w_0 \simeq 20.23$,$w_1 \simeq -14.63$,$w_2 \simeq 2.32$,$w_3 \simeq 0.08$であることがわかりました. この値を基に回帰式をグラフ化すると,曲線を利用して対象のデータを近似していることがわかります.

plt.scatter(x, y)
plt.plot(x, fx(result.x, x))
plt.xlabel("week")
plt.ylabel("infuluenza")

Image from Gyazo

ガウス基底

ガウス基底は下記の式で表されます. ここで,$i$は組み合わせる基底関数の番号を表します. 関数の中心位置を表す$\mu_i$,関数の広がりを表す$s$は,事前に適切な値に設定します.

$$ \phi_i(x) = \exp \left( - \frac{ (x - \mu_i)^2}{2s^2} \right) $$

ガウス基底を関数 gauss として定義しておきます.

def gauss(x, mu, s):
  g = np.exp( -(x - mu) ** 2 / (2 * s**2))
  return g

ガウス基底の形状を確認しましょう. $\mu=0$,$s=3$のグラフを描いてみます. ガウス基底は上に凸の分布であり,正規分布(ガウス分布)と同じ形状であることが分かります. 正規分布の正規化(面積の総和を1にする)のための係数を取り除いた式がガウス基底です.

x_ = np.arange(-10, 10, 0.1)
y_ = gauss(x_, 0, 3)
plt.plot(x_, y_)

Image from Gyazo

回帰式

4つの基底関数を組み合わせることにします. このとき,$\phi_0(x)=1$とするため,実際は残りの3つのガウス基底を組み合わせることになります.

$$ f(x) = w_0 \cdot 1 + w_1 \cdot \phi_1(x) + w_2 \cdot \phi_2(x) + w_3 \cdot \phi_3(x) $$

この3つのガウス基底の平均$\mu$と間隔$s$は下記のように設定します.

$$ (\mu_1=5, s_1=5), (\mu_2=10, s_2=5), (\mu_3=15, s_3=5), $$

これは3つのガウス基底を均等な間隔で配置することを意味しています.

Image from Gyazo

この回帰式を関数 fx として定義しておきます.

def fx(w, x): # 回帰式の関数定義
  y_ = w[0]  + w[1] * gauss(x, 5, 5) + w[2] * gauss(x, 10, 5)  + w[3] * gauss(x, 15, 5)
  return y_

目的関数

目的関数$E$は,多項式基底と同じ2乗誤差で与えます.

数値解の導出

多項式基底と同様に SciPy で最適化します.

[In:]

w_init = np.array([0, 0, 0, 0]) # パラメータの初期値
result = minimize(gx, w_init, args=(x, y), method="powell") # Powell法で最適化
print(result.x) # 導出されたパラメータ

[Out:]

[-60.68720705  74.75886194 -44.73746539 130.07080173]

この結果,$w_0 \simeq -60.69$,$w_1 \simeq -74.76$,$w_2 \simeq -44.74$,$w_3 \simeq 130.07$であることがわかりました. この値を基に回帰式をグラフ化すると,多項式基底と同様に曲線で対象のデータを近似していることがわかります.

plt.scatter(x, y)
plt.plot(x, fx(result.x, x))
plt.xlabel("week")
plt.ylabel("infuluenza")

Image from Gyazo

参考書籍