はじめに #
この演習では，パターン認識と機械学習 上 (PRML 上)の第1章に出てくる多項式回帰を再現することで，線形回帰の基礎を学習する．
データはsin関数にガウスノイズを乗せたものを10点用意し，0次，1次，3次および9次の多項式で回帰する．
正規方程式を NumPy のライブラリを使って解く．
ある程度コードの骨格は載せているが，??????の部分などを自分で考えて書くこと．
NumPy とグラフの描画に Matplotlib を用いるので先に下記のページを見てそれぞれの基本を学んでおくこと．
Jupyter Notebookでの演習を想定．下記のようにライブラリのインポートも必要． ここに出てくるグラフ出力例は Matplotlib > よく使う設定と保存 の設定も用いている．
import numpy as np
import matplotlib.pyplot as plt演習1：sin関数のプロット #
この演習での目的関数は
$$ y = \sin(2\pi x) $$ここで， \(x\) の範囲は \(0 \le x < 1 \) を考える．
sin関数のプロットは NumPy > 数学関数 のページに書いてあるので，そのまま用いる．ただし， 変数x，yは後で別のもので使いたいので，それぞれ x_denseおよびy_sinに変更してある．
# ---------- x, y
x_dense = np.linspace(0, 1, 100)
y_sin = np.sin(2 * np.pi * x_dense)
# ---------- plot
fig, ax = plt.subplots(1, 1)
ax.plot(x_dense, y_sin, label=r'$\sin(2\pi x)$')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
自分でプロットして確認すること．
演習2：訓練データ #
訓練データとして，sin関数にガウスノイズを加えた10個のデータを作成する． 下記の手順にしたがって訓練データを生成し，sin関数とともにプロットすること． 一様乱数やガウス分布に従う乱数生成については NumPy > 乱数 を参考．
- \(0 \le x < 1 \) の範囲で，ランダムに10点とった1次元配列（ndarray）
x_trainを作る． x_trainに対して，\(y_\mathrm{train} = \sin(2\pi x_\mathrm{train})\) を計算し，同サイズの1次元配列y_trainを作る．- 平均0，標準偏差0.1のガウス分布に従う乱数を10点生成して，同サイズの1次元配列
noiseを作る． y_trainにnoiseを加えて更新して，sin関数とともにプロットする
rng = np.random.default_rng(123)のように何か適当な数値を seed として入れておくとよい
rng = np.random.default_rng()
x_train = ??????
x_train # jupyterで出力用y_train = ?????? # sin
y_train # jupyterで出力用noise = ?????? # ガウスノイズ
noise # jupyterで出力用y_train = y_train + noise # ガウスノイズを加えて更新
# ---------- plot
fig, ax = plt.subplots(1, 1)
ax.plot(x_dense, y_sin, label=r'$\sin(2\pi x)$')
ax.plot(x_train, y_train, 'o', markeredgecolor='black', label='Training data')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
演習3：計画行列 #
最終的には以下の4つのモデルで，演習2で生成したデータ10点の線形回帰を行いたい．
- 0次：\( y = w_0 \)
- 1次：\( y = w_0 + w_1 x \)
- 3次：\( y = w_0 + w_1 x + w_2 x^2 + w_3 x^3 \)
- 9次：\( y = \sum_{j=0}^9 w_j x^j \)
そのためにここではまず計画行列 \(\Phi\) を計算する関数design_matrix(x, M)を作る．
ここで，\(M\) はモデルの次数． ちなみに定数項を入れると項の数は \( (M+1) \) になる．
線形回帰モデルは
$$ \mathbf{y} = \Phi \mathbf{w} $$と表され，作成した10点の訓練データをこのモデルに当てはめると
$$ \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_9 \end{pmatrix} = \begin{pmatrix} x_0^0 & x_0^1 & \cdots & x_0^{M} \\ x_1^0 & x_1^1 & \cdots & x_1^{M} \\ \vdots & \vdots & \ddots & \vdots \\ x_9^0 & x_9^1 & \cdots & x_9^{M} \\ \end{pmatrix} \begin{pmatrix} w_0 \\ \vdots \\ w_{M} \end{pmatrix} $$となる．ここで
$$ \Phi = \begin{pmatrix} x_0^0 & x_0^1 & \cdots & x_0^{M} \\ x_1^0 & x_1^1 & \cdots & x_1^{M} \\ \vdots & \vdots & \ddots & \vdots \\ x_9^0 & x_9^1 & \cdots & x_9^{M} \end{pmatrix} $$が計画行列である．この計画行列を返す関数は色々な書き方があるが，例えば以下のようになる．
def design_matrix(x, M):
x_mat = x.reshape(-1, 1) # 10行1列の2次元配列に変換
p = np.arange(M + 1) # 0, 1, 2, ..., M
dmat = np.power(x_mat, p)
return dmatnp.arange()やnp.power()については下記ページを参考にすること．
チェックのために，\(M=3\) の場合の計画行列を画面出力してプログラムが合っているかどうか確認すること．
M = 3
phi = design_matrix(x_train, M)
phi # jupyterで出力用例えば下記は一例．
array([[1.00000000e+00, 5.35370744e-01, 2.86621833e-01, 1.53448944e-01],
[1.00000000e+00, 2.51244518e-01, 6.31238078e-02, 1.58595107e-02],
[1.00000000e+00, 3.59168825e-01, 1.29002245e-01, 4.63335848e-02],
[1.00000000e+00, 5.44489549e-01, 2.96468869e-01, 1.61424200e-01],
[1.00000000e+00, 1.68456471e-01, 2.83775827e-02, 4.78038745e-03],
[1.00000000e+00, 4.44173207e-02, 1.97289838e-03, 8.76308599e-05],
[1.00000000e+00, 7.67506529e-01, 5.89066272e-01, 4.52112210e-01],
[1.00000000e+00, 3.48437878e-01, 1.21408955e-01, 4.23034785e-02],
[1.00000000e+00, 8.79593377e-01, 7.73684508e-01, 6.80527769e-01],
[1.00000000e+00, 6.41433853e-01, 4.11437388e-01, 2.63909869e-01]])演習4：フィッティング #
ここでは，\(M=3\) で試しにフィッティングしてみる．
計画行列が計算できれば，正規方程式よりが計算できる． 正規方程式の導出などは 最小二乗法 のPDFを見ること．
PDFにも書いてあるが，ランク落ちや数値計算上不安定になる場合に備えて，
この演習では，\( \Phi \mathbf{w} = \mathbf{y} \) の形の最小二乗解を計算できるnp.linalg.lstsq(phi, y)を用いて \( \mathbf{w} \) を求める．
np.linalg.lstsq(phi, y)の返り値はタプルで
(w, Sums of squared residuals, rank of phi, singular values of phi)のようになっているので，はじめのインデックスを取り出せばそれが \( \mathbf{w} \) になる．
w, residuals, rank, s = np.linalg.lstsq(phi, y_train)
w # jupyterで出力用例えば \( \mathbf{w} \) は以下のように出力される．
array([ -0.09749882, 10.87685823, -31.30137186, 20.41339222])\( \mathbf{w} \) が求まったので，\(x\) を指定すればその多項式のモデル
$$ y = w_0 + w_1 x + w_2 x^2 + w_3 x^3 $$
が計算できる．sin関数を描くときに使ったx_denseを用いて，それに対応するy_dense_predictを計算すれば，連続的なモデルを描画できる．
\( \mathbf{y} = \Phi \mathbf{w} \) であったことを思い出すと，10点の訓練データx_trainの代わりに，x_denseの計画行列phi_denseを計算して，\( \mathbf{w} \) との積を取ればy_dense_predictは簡単に計算できる．
phi_dense = design_matrix(x_dense, M)
y_dense_predict = ??????y_dense_predictが計算できたら，sin関数と訓練データとともにプロットすること．
例えば下図のようなものが得られる．
fig, ax = plt.subplots(1, 1)
ax.plot(x_dense, y_sin, label=r'$\sin(2\pi x)$')
ax.plot(x_train, y_train, 'o', markeredgecolor='black', label='Training data')
ax.plot(x_dense, y_dense_predict, label='Prediction')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
演習5：0次，1次，3次，9次の多項式回帰 #
ここまでで，3次多項式による回帰ができたはずなので，0次，1次，3次，9次で同様に回帰を行ない， 全て同時にプロットすることで比較する．
M の値を変化させるためのリスト用意する．
mlist = [0, 1, 3, 9]あとはMをfor文で回しながら，3次の回帰の時と同様に計算すれば良い．
fig, ax = plt.subplots(4, 1, figsize=(8, 12))
for i, M in enumerate(mlist):
phi = ??????
w, residuals, rank, s = ??????
phi_dense = ??????
y_dense_predict = ??????
ax[i].plot(x_dense, y_sin, label=r'$\sin(2\pi x)$')
ax[i].plot(x_train, y_train, 'o', markeredgecolor='black', label='Training data')
ax[i].plot(x_dense, y_dense_predict, label='Prediction')
ax[i].set_ylim(-2, 2) # yの描画範囲を設定
ax[i].set_ylabel(r'$y$')
ax[i].text( # M = 0などのテキストの表示
0.05, 0.90,
rf'$M = {M}$',
transform=ax[i].transAxes,
va='top',
bbox=dict(facecolor='white'),
fontsize=14
)
ax[-1].set_xlabel(r'$x$') # x軸ラベルは一番下のグラフだけに表示
ax[0].legend(fontsize=12) # 一番上のグラフだけに凡例を表示．少しフォントを小さくしているfig, ax = plt.subplots(4, 1, figsize=(8, 12))のようにして４つのグラフを縦に並べる．figsizeは自由に調節して良い．ax.plot()ではなく，ax[0].plot(),ax[1].plot()… のように使う
enumerate()の使い方は各自で調べること．
下記のような結果が得られればうまく回帰できているので確認すること．
