【機械学習】線形基底関数モデル(線形回帰モデル)を理解して実装する

Posted: , Category: 回帰 , 機械学習 , 統計学

線形基底関数モデルについて理解しているでしょうか?

このモデルは線形回帰モデルとも呼ばれており、かの有名なPRMLでも序盤に登場するほど良く知られており、その考え方はガウス過程回帰など、多くの回帰のアルゴリズムの基礎となっているため、機械学習の学習を進める上で、理解をしておきたい回帰モデルです。また実用上もよく用いられています。

線形基底関数モデルは線形回帰モデルと呼ばれていますが、一般的に線形回帰とは違うことに注意してください。線形回帰は、基底関数が全て直線であったのに対し、線形回帰モデルは、基底関数が線形でない場合も扱うことが多くあります。と言っても、この違いについてもこの記事で解説していきます。

今回は、線形基底関数モデル(線形回帰モデル)について解説し、その後実際にPythonを用いて実装することを通して、線形基底関数モデルについて理解します。

本記事の内容
  • 線形基底関数モデルの概要
  • 線形基底関数モデルのパラメータを数式で導出
  • 線形基底関数モデルを実装

線形基底関数モデルとは

線形基底関数とは、説明変数の線形和で目的関数を予測する単回帰や重回帰などと違い、説明変数の多項式で目的変数を表現する手法です。線形基底関数モデルを理解するためには、まず単回帰や重回帰について理解する必要があるため、これらに不安がある人は、まず下記の記事をご覧ください。

【Python】単回帰分析の理論と実装をわかりやすく解説
単回帰分析は、機械学習の1つの手法である、直線近似の手法です。統計学や実務のデータサイエンス分野の教科書でも、一番最初に登場するくらい基本的な分析手法です。 単回帰分析を実際に行う場合は、Pythonやエクセル、Rなどの […]
【Python実装】重回帰分析の理論と実装を理解する
多変量解析の中で、重回帰分析は非常によく用いられている手法です。この記事を読んでいる人の中でも、その名前くらいは耳にしたことも多いのではないでしょうか。 重回帰分析は、説明変数が一般に多次元あるような場合の回帰問題を解く […]

単回帰分析では、説明変数を$x$としたとき、パラメータ$w_0, w_1$を用いて、

\begin{equation}
\hat{y} = w_0 + w_1 x
\end{equation}

の回帰式で、目的変数を予測するモデルでした。

例えば、体重(weight)と身長(height)のデータであれば、短回帰分析によりデータからこれらの回帰直線を引くことで、体重から身長を予測できるます。

これらのデータの場合、うまくデータを表現できている気がしますね。しかし、単回帰分析が常にうまく行くとは限りません。例えば次のようなデータではどうでしょうか。

このようなデータ点があると、直線では表現できないことがわかると思います。おそらく、このデータてんを表現するためには、こんな曲線を引いて回帰する必要があると思います。

このような関係性を表現するには、(1)式のような表現では表現力が不足しており、例えば

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

のような関数を仮定する必要があります。一般的に、線形規定関数モデルでは、説明変数の関数$\phi(x)$を任意の数$H$個用いて、次のような回帰式を用います。

線形基底関数モデルの回帰式
\begin{equation}
\begin{split}
\hat{y} &= w_0 +  w_1 \phi_1(x) + w_2\phi_2(x) + w_3\phi(x) + \cdots +  w_H \phi(x) \\
&= \sum_{j=0}^{H} w_j \phi_j(\bm{x}) \\
& = \bm{w}^T \bm{\phi({\bm{x}})}
\end{split}
\end{equation}

ここで、$\bm{\phi({\bm{x}})}$は$x$の関数$\phi(x)$のベクトルです。

線形基底関数モデル(線形回帰モデル)の回帰式は(3)のようになります。

ここで注意して欲しいのが、単回帰や重回帰の場合は、説明変数の次元の数 + 1だけ重み$w$が存在しましたが、線形基底関数モデルの$\phi(x)$は説明変数の数にはよらないことに注意してください。一般的に$H$は、線形基底関数モデルを設計する人が自由に選べるパラメータとなっています。

ここで、単回帰や重回帰で考えたように、全てのデータに対して、行列形式で(3)の回帰式を表現すると、次のようになります。

\begin{equation}
\begin{split}
\begin{pmatrix}
\hat{y_1} \\
\hat{y_2} \\
\vdots  \\
\hat{y_D} \\
\end{pmatrix} 
& = 
\begin{pmatrix}
1 & \phi_1(x_{1}) &  \phi_2(x_{1}) & \cdots & \phi_H(x_1) \\
1 &  \phi_1(x_{2}) & \phi_2(x_{2}) & \cdots & \phi_H(x_2) \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
1 &  \phi_1(x_{N}) & \phi_2(x_{N}) & \cdots &  \phi_H(x_N) \\
\end{pmatrix}
\begin{pmatrix}
w_0 \\
w_1 \\
\vdots  \\
w_H \\
\end{pmatrix} \\
&= \bm{\Phi} \bm{w}
\end{split}
\end{equation}

と表現することができます。

では、この$\bm{w}$をどのように決めるのでしょうか、となりますが、(4)式は重回帰分析と同じ形式をしているため、重回帰分析で行ったのと同様に、パラメータベクトル$\bm{w}$で偏微分をすることで、最適なパラメータを決定することができます。

ここでは結論だけ記載します。

線形基底関数モデルのパラメータの解
\begin{equation}
\begin{split}
\bm{w} = (\bm{\Phi^T\Phi})^{-1}\bm{\Phi}^T\bm{y}
\end{split}
\end{equation}

線形基底関数モデルの実装

理論だけ勉強しても、技術の実態を理解することは難しいので、実際に実装してみましょう。

線形基底関数モデルは、説明変数を与えた際に、目的変数の値を予測する回帰モデルになります。

今回、線形基底関数モデルで予測したいデータセット

今回、線形基底関数モデルで予測したいデータセットは次のように人工データを用意します。

本当は、実データを用いて、線形基底関数モデルを試しかったのですが、良いデータセットが見つからなかったので、今回は人工データセットで試します。

人工データセットの真のモデルは、次のような関数系になるとします。

今回の回帰で使うデータの真のモデル
\begin{equation}
y = x^4 + 4cos (10x) + 6sin(3x)
\end{equation}

この関数を図示すると、このようになります。今回は$-3~3$の定義域で考えることにします。

2時間数でもプロットできるような感じですが、$-2 ~ 2$あたりが三角関数の項が効いてきて、グネグネ動いていますね。この関数は中々回帰するのが難しそうです。

import numpy as np
import matplotlib.pyplot as plt

def true_model(x):
    return x ** 4 + 4 * np.cos(10 * x) + 6 * np.sin(3*x)

xs = np.linspace(-3, 3, 100)
y = true_model(xs)

fig, ax = plt.subplots(figsize=(10, 2), dpi=100)

ax.plot(xs, y, c="mediumaquamarine")
ax.set_title(r"$True Model: y=x^4+4cos(10x)+6sin(3x)$")
ax.legend()

続いて、この関数からデータをサンプリングして、データセットを作ってみます。

(2)の関数からそのままデータを取得してもいいのですが、実際のデータは何かしらのノイズが必ず入るので、今回もノイズを加えてデータを取得してみます。ノイズは、平均が0、分散が$\sigma^2$の正規分布とします。

$\sigma$の値もなんでも良いのですが、取り合えす今回は、$\sigma = 5$としてデータを生成してみます。データ数は$N=500$で生成してみます。

N = 500
noice_sigma = 5

rand_x = np.random.uniform(low=np.min(xs), high=np.max(xs), size=N)
rand_y = true_model(rand_x) + np.random.normal(loc=0, scale=noice_sigma, size=N)

fig, ax = plt.subplots(figsize=(10, 2), dpi=100)

ax.plot(xs, y, c="mediumaquamarine")

ax.set_title(r"$True Model: y=x^4+4cos(10x)+6sin(3x)$")
ax.fill_between(x=xs, y1=y-2.0*noice_sigma, y2=y+2.0*noice_sigma,  color='mediumspringgreen', alpha=0.2, linestyle='--')

ax.scatter(rand_x, rand_y, color="coral", s=2)

ax.legend()

薄いライムグリーンの箇所は、2$\sigma$区間を示しています。2$\sigma$区間なので、およそ95%のデータはこの区間に入っていますが、稀にこの区間から外れるようなデータも出てきています。

データだけをプロットするとこのようになります。

これで正解データは作れました。このデータをもとに、線形基底関数モデルを当て込んでいきます。

デザイン行列(計画行列)を実装する

【広告】
統計学的にあなたの悩みを解決します。
仕事やプライベートでお悩みの方は、ベテラン占い師 蓮若菜にご相談ください。

機械学習と情報技術