演習:多項式回帰

この演習では,パターン認識と機械学習 上 (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の最後のところ実は答えが書いてある. 変数名はなんでも良いがxyは後で使いたいので,xaytを使うこと.xaは要素数100程度あれば良い.

xa =  ??? # ここにコードを書く    # 要素数100の1次元配列
xa    # jupyterで表示するため
yt = ??? # ここにコードを書く    # 要素数100の1次元配列
yt
plt.plot(xa, yt)    # matplotlibでプロット

演習2:ガウスノイズ

下記の手順に従って,sin関数にガウスノイズを加えた10個のデータ点を作成してプロットせよ.

  1. 演習1と同じxの範囲で,ランダムに10点とった1次元のndarrayを作る.
  2. 10個のデータ$x$に対して,$y = \sin (2\pi x)$を計算して同サイズのndarrayを作る.
  3. 平均0,標準偏差0.1のガウス分布に従う乱数を10点生成して,同サイズのndarrayを作る.
  4. 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点のデータ$x_i$,$t_i$をこのモデルに当てはめて,具体的に書くと

$$ \begin{pmatrix} t_0 \\ t_1 \\ \vdots \\ t_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{Aw} = \mathbf{b}$の形の連立方程式を計算できるnp.linalg.solve()を用いて $\mathbf{w}$を求める.そのために

$$\mathbf{t} = \mathbf{\Phi} \mathbf{w}$$

$$\mathbf{Aw} = \mathbf{b}$$

の形に変形する.ここで$\mathbf{A}$は正方行列,$\mathbf{b}$はベクトルであり,

$$ \mathbf{A} = \Phi^\mathrm{T} \Phi $$ $$ \mathbf{b} = \Phi^\mathrm{T} \mathbf{t} $$ である.$\mathbf{A}$と$\mathbf{b}$を計算するプログラムを作成せよ.

# Aを計算
A = ????    # ここにコードを書く
A
# bを計算
b = ????    # ここにコードを書く
b

np.linalg.solve()を使って$\mathbf{w}$を求めよ.

# Aw = bを解く
w = np.linalg.solve(A, b)
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次で同様に回帰を行ない, 全て同時にプロットすることで比較する.

上記と同じxtxaを利用するのでまず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 = ???
    A = ???
    b = ???
    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)
前へ
次へ