この演習では，パターン認識と機械学習 上 (PRML 上)の第1章に出てくる多項式回帰を再現することで，線形回帰の基礎を学習する． データはsin関数にガウスノイズを乗せたものを10点用意し，0次，1次，3次および9次の多項式で回帰する． 行列が正則かどうかなどの細かいところは今回は無視して，正規方程式をnumpyのライブラリを使って解く．
jupyter notebookやjupyter labでの演習を想定している．下記のようにライブラリのインポートも必要．
#%matplotlib inline # jupyter notebookでは必要．JupyterLabでは不要
import numpy as np
import matplotlib.pyplot as plt
演習1：sin関数のプロット
$$y = \sin (2\pi x)$$
をプロットせよ．ただし$x$の範囲は
$$0 \le x < 1$$
とする．
このサイトの研究ノートのNumPyの最後のところ実は答えが書いてある．
変数名はなんでも良いがxとyは後で使いたいので，xaとytを使うこと．xaは要素数100程度あれば良い．
xa = ??? # ここにコードを書く # 要素数100の1次元配列
xa # jupyterで表示するため
yt = ??? # ここにコードを書く # 要素数100の1次元配列
yt
plt.plot(xa, yt) # matplotlibでプロット
演習2：ガウスノイズ
下記の手順に従って，sin関数にガウスノイズを加えた10個のデータ点を作成してプロットせよ．
- 演習1と同じ
xの範囲で，ランダムに10点とった1次元のndarrayを作る． - 10個のデータ$x$に対して，$y = \sin (2\pi x)$を計算して同サイズのndarrayを作る．
- 平均0，標準偏差0.1のガウス分布に従う乱数を10点生成して，同サイズのndarrayを作る．
- 2と3で作ったものを足して，元のsin関数とともにプロットする．
x = ??? # ここにコードを書く # ランダム10点
x
tt = ??? # ここにコードを書く # 上記乱数xに対するsin(2 pi x)
tt
err = ??? # ここにコードを書く # 平均0，標準偏差0.1のガウス分布に従う乱数を10点
err
t = tt + err # sin関数にノイズを加える
plt.plot(xa, yt) # 元のsin関数をプロット
plt.plot(x, t, 'o') # ガウスノイズを加えた10点をプロット
演習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$
そのためにここではまず計画行列$\mathbf{\Phi}$を計算する関数design_matrix(x, M)
を作る．ここで，$M$はモデルの次数．
ちなみに定数項を入れると項の数は$(M+1)$になる．
線形回帰モデルは
$$\mathbf{y} = \mathbf{\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} $$
ここで
$$\mathbf{\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} $$
が計画行列である．この計画行列を返す関数を作成せよ．
色々な書き方があるが，
np.arrange(), np.power(), reshape()あたりを使うと
スッキリ書ける．（研究ノートのNumPyのページを参照）
# ---------- 計画行列dmatを返す関数
def design_matrix(x, M):
???? # ここにコードを書く
????
return dmat
チェックのために，$M = 3$を使って実際に計画行列を作成，出力し プログラムが合っているかどうか確認すること．
# ---------- 多項式の次数を指定
M = 3
# ---------- 計画行列を上で作った関数を使って計算
phi = design_matrix(x, M)
phi # jupyterで表示用
演習4：フィッティング
ここでは，演習3の$M=3$の結果を使用する．
計画行列が計算できれば，正規方程式より$w$が計算できる．
正規方程式の導出などは最小二乗法のノートを見ること．
ノートにも書いてあるように，ランク落ちや数値計算上不安定になる場合に備えて，
この演習では，$\mathbf{\Phi} \mathbf{w} = \mathbf{t} $の形の最小二乗解を計算できるnp.linalg.lstsq(phi, t)を用いて
$\mathbf{w}$を求める．
np.linalg.lstsq(phi, t)の返り値はタプルで(w, Sums of squared residuals, rank of phi, singular values of phi)のようになっているので，はじめのインデックスを取り出せばそれが$\mathbf{w}$になる
# wを計算する
w = np.linalg.lstsq(phi, t)[0]
w # w表示用
演習5：モデルの描画
$\mathbf{w}$が求まったので，$\mathbf{x}$を指定すればその多項式のモデルが計算できる
$y = w_0 + w_1 x^1 + w_2 x^2 + \cdots$
演習1でsin関数を描くときに使ったxaを使って，xaに対応するyを計算する．
演習4で得られた$\mathbf{w}$を用いたモデルをプロットせよ．
xa # jupyter上でxaを表示して確認
$\mathbf{y} = \mathbf{\Phi} \mathbf{w}$であったことを思い出すと，
10点のデータ点xの代わりにxaの計画行列を計算して，
$\mathbf{w}$との積をとれば$\mathbf{y}$は簡単に計算できる．
phi2 = design_matrix(xa, m) # xaの計画行列を計算
y = ???? # ここにコードを書く
# 元のsin関数，10個のデータ点と重ねて計算したモデルyを描画
plt.plot(xa, yt) # 元のsin関数
plt.plot(x, t, 'o') # 訓練データ
plt.plot(xa, y) # 予測モデル
演習6：0次，1次，3次，9次の多項式回帰
演習5までで，3次多項式による回帰ができたはずなので，0次，1次，3次，9次で同様に回帰を行ない， 全て同時にプロットすることで比較する．
上記と同じx，t，xaを利用するのでまずjupyter上で変数の確認をする．
x # 表示するだけ
t # 表示するだけ
xa # 表示するだけ
Mの値を変化させるためのリストを用意する．
mlist = [0, 1, 3, 9]
これまでの演習と同様に回帰してグラフをプロットする．
plt.subplotsについては各自で調べる．figsizeなどは自由に調節して良い．
enumerate()の使い方も覚えておくこと．
fig, ax = plt.subplots(4, 1, figsize=(8, 12))
for i, M in enumerate(mlist):
phi = ???
w = ???
phi2 = ??? # yを計算するためのxaの計画行列
y = ???
ax[i].set_ylim(-1.5, 1.5) # yの描画範囲を設定
ax[i].plot(xa, yt)
ax[i].plot(x, t, 'o')
ax[i].plot(xa, y)