一日坊主

雰囲気でやっている

PRML上巻 P10-11

昨日の続き.

多項式曲線フィッティングの線形方程式verを作成した.

動作確認のため,図1.6を再描画してみる.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.polynomial.polynomial as poly
np.random.seed(42)
# test data set (N=15)
# t = sin(2 \pi x) + gaussian noise
x4 = np.linspace(0, 1, 15)
t4 = np.sin(2 * np.pi * x4) + np.random.normal(0, 0.2, 15)
# test data set (N=100)
# t = sin(2 \pi x) + gaussian noise
x1 = np.linspace(0, 1, 100)
t3 = np.sin(2 * np.pi * x1) + np.random.normal(0, 0.2, 100)
# figure 1.6

def polyfit_linear_equation(x: np.ndarray, t: np.ndarray, M: int) -> np.ndarray:
    """solve linear equation
    
    Args:
        x (np.ndarray): explanatory variable (N dims)
        t (np.ndarray): response variable (N dims)
        M (int): order
    
    Returns:
        np.ndarray: coefficient (M dims)
    """
    # add intercept term
    M += 1
    # A: M x M matrix
    A = np.zeros(shape=(M, M))
    for i in range(M):
        for j in range(M):
            A[i][j] = np.sum(np.power(x, i + j))
    # T: 1 x M vector
    T = np.zeros(M)
    for i in range(M):
        T[i] = np.dot(np.power(x, i), t)
    # fitting with linear equations
    # FIXME: the result is flipped. To make the result in proper order, we must flip it.
    return np.linalg.solve(A, T)[::-1]

# fitting with polynomial functions
w9_15 = polyfit_linear_equation(x4, t4, 9)
w9_100 = polyfit_linear_equation(x1, t3, 9)

# visualize
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
axs[0].set_title('N = 15')
axs[1].set_title('N = 100')

for ax in axs.flat:
    ax.plot(x1, t1, color='green')
    ax.set(xlabel='x', ylabel='t')

axs[0].scatter(x4, t4, facecolor='None', edgecolor='blue')
axs[0].plot(x1, np.poly1d(w9_15)(x1), color='red')
axs[1].scatter(x1, t3, facecolor='None', edgecolor='blue')
axs[1].plot(x1, np.poly1d(w9_100)(x1), color='red')

plt.tight_layout()

f:id:twakamori:20210223225533p:plain
figure 1.6 (linear equation)

正しく動作しているようだ.

あらためて,線形方程式verを用いて図1.7を描画してみる.

なお,正則化した線形方程式には,演習問題1.2の解答である次式を用いる.

\displaystyle{
\sum_{j=0}^{M}(A_{ij}-\lambda I_{ij})w_j=T_i
}
# figure 1.7

def polyfit_linear_equation_regularize(x: np.ndarray, t: np.ndarray, M: int, alpha: float) -> np.ndarray:
    """solve linear equation
    
    Args:
        x (np.ndarray): explanatory variable (N dims)
        t (np.ndarray): response variable (N dims)
        M (int): order
        alpha (float): regularization parameter
    
    Returns:
        np.ndarray: coefficient (M dims)
    """
    # add intercept term
    M += 1
    # A: M x M matrix
    A = np.zeros(shape=(M, M))
    for i in range(M):
        for j in range(M):
            A[i][j] = np.sum(np.power(x, i + j))
    # I: M x M identity matrix
    I = np.identity(M)
    # T: 1 x M vector
    T = np.zeros(M)
    for i in range(M):
        T[i] = np.dot(np.power(x, i), t)
    A_reg = A - alpha * I
    # fitting with linear equations
    # FIXME: the result is flipped. To make the result in proper order, we must flip it.
    return np.linalg.solve(A_reg, T)[::-1]
w9_l18 = polyfit_linear_equation_regularize(x4, t4, 9, np.exp(-18))
w9_l0 = polyfit_linear_equation_regularize(x4, t4, 9, np.exp(0))
# visualize
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
axs[0].set_title(r'$\ln{\lambda} = -18$')
axs[1].set_title(r'$\ln{\lambda} = 0$')

for ax in axs.flat:
    ax.plot(x1, t1, color='green')
    ax.set(xlabel='x', ylabel='t')

axs[0].scatter(x2, t2, facecolor='None', edgecolor='blue')
axs[0].plot(x1, np.poly1d(w9_l18.flatten())(x1), color='red')
axs[1].scatter(x2, t2, facecolor='None', edgecolor='blue')
axs[1].plot(x1, np.poly1d(w9_l0.flatten())(x1), color='red')

plt.tight_layout()

f:id:twakamori:20210223225737p:plain
figure 1.7

(元々のデータが書籍記載のものと異なるためか,書籍記載の図を再現できていない.\lambdaの値を変化させて確認したところ,正則化自体はうまく動作しているようだった.よってこの図はこのままとする.)

図から,\lambdaによって過学習が抑制され,背後の関数\sin(2\pi x)に近い表現が得られることがわかる.ただし,\lambdaの値を大きくしすぎると図1.7の右側に示すように当てはまりは再び悪くなる.

正則化が期待通り係数が大きい値を取らないようにしているかを見るために,当てはめた多項式の係数を次表に示す.(コードは省略)

\ln{\lambda}=-\infty \ln{\lambda}=-18 \ln{\lambda}=0
w_ 0^ * 20812.84 -2005.15 -0.95
w_ 1^ * -92408.71 3634.86 -0.85
w_ 2^ * 170879.63 614.50 -0.71
w_ 3^ * -170219.15 -3401.47 -0.54
w_ 4^ * 98565.04 -1228.53 -0.31
w_ 5^ * -33325.18 4548.16 -0.02
w_ 6^ * 6273.03 -2740.68 0.37
w_ 7^ * -609.66 618.99 0.84
w_ 8^ * 32.65 -40.79 1.19
w_ 9^ * -0.38 0.20 -0.38

表より,\lambdaの値を増やすと共にほとんどの係数は小さな値を取るようになることがわかる.

正則化項が汎化誤差に与える影響は,訓練集合とテスト集合の両方に対するRMS誤差(1.3)の値を\ln{\lambda}に対してプロットしてみれば良い.これを図1.8に示す.

# figure 1.8

RMSE_reg_train = []
RMSE_reg_test = []
w_values_reg = []
for alpha in range(-380, -199):
    w = polyfit_linear_equation_regularize(x2, t2, 9, np.exp(alpha / 10))
    w_values_reg.append(w.flatten())
    y_train = np.poly1d(w.flatten())(x2)
    y_test = np.poly1d(w.flatten())(x1)
    E_train = 0.5 * np.sum(np.square(y_train - t2))
    E_test = 0.5 * np.sum(np.square(y_test - t3))
    RMSE_reg_train.append(np.sqrt(2 * E_train / t2.size))
    RMSE_reg_test.append(np.sqrt(2 * E_test / t3.size))
ln_alpha = np.linspace(-38, -19, 181)
p1 = plt.plot(ln_alpha, RMSE_reg_train, color='blue')
p2 = plt.plot(ln_alpha, RMSE_reg_test, color='red')
plt.xlabel(r'$\ln{\lambda}$')
plt.ylabel(r'$E_{RMS}$')
plt.ylim(0, 1)
plt.legend([p1[0], p2[0]], ['訓練', 'テスト'])

f:id:twakamori:20210223230131p:plain
figure 1.8

むむ.書籍と随分違うグラフになってしまった. このグラフから何か読み取るのは困難だが,本来は,テスト誤差が最小となるような正則化パラメータ\lambdaを探索すればよい,ということになるはずである.

モデルの複雑さの問題は重要で,誤差関数を最小にするようなアプローチで実際の応用問題を解こうとする際には,モデルの複雑さを適切に決める方法を見つけなければならない. そのためには,得られたデータを,係数\mathbf{w}を決めるために使われる訓練集合と,それとは別の確認用集合(検証用集合; validation set)に分ける方法がある.確認用集合はホールドアウト集合(hold-out set)とも呼ばれ,モデルの複雑さ(Mまたは\lambda)を最適化するのに使われる.

これまでの多項式曲線フィッティングの議論は,主として直感に訴えるものであった. ここからは,確率論的な議論の枠組みを使い,パターン認識における問題を解くためのより原理的なアプローチを探っていく.

以上.