PRML上巻 P10-11
昨日の続き.
多項式曲線フィッティングの線形方程式verを作成した.
動作確認のため,図を再描画してみる.
%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()
正しく動作しているようだ.
あらためて,線形方程式verを用いて図を描画してみる.
なお,正則化した線形方程式には,演習問題1.2の解答である次式を用いる.
# 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()
(元々のデータが書籍記載のものと異なるためか,書籍記載の図を再現できていない.の値を変化させて確認したところ,正則化自体はうまく動作しているようだった.よってこの図はこのままとする.)
図から,によって過学習が抑制され,背後の関数に近い表現が得られることがわかる.ただし,の値を大きくしすぎると図の右側に示すように当てはまりは再び悪くなる.
正則化が期待通り係数が大きい値を取らないようにしているかを見るために,当てはめた多項式の係数を次表に示す.(コードは省略)
20812.84 | -2005.15 | -0.95 | |
-92408.71 | 3634.86 | -0.85 | |
170879.63 | 614.50 | -0.71 | |
-170219.15 | -3401.47 | -0.54 | |
98565.04 | -1228.53 | -0.31 | |
-33325.18 | 4548.16 | -0.02 | |
6273.03 | -2740.68 | 0.37 | |
-609.66 | 618.99 | 0.84 | |
32.65 | -40.79 | 1.19 | |
-0.38 | 0.20 | -0.38 |
表より,の値を増やすと共にほとんどの係数は小さな値を取るようになることがわかる.
正則化項が汎化誤差に与える影響は,訓練集合とテスト集合の両方に対するRMS誤差の値をに対してプロットしてみれば良い.これを図に示す.
# 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]], ['訓練', 'テスト'])
むむ.書籍と随分違うグラフになってしまった. このグラフから何か読み取るのは困難だが,本来は,テスト誤差が最小となるような正則化パラメータを探索すればよい,ということになるはずである.
モデルの複雑さの問題は重要で,誤差関数を最小にするようなアプローチで実際の応用問題を解こうとする際には,モデルの複雑さを適切に決める方法を見つけなければならない. そのためには,得られたデータを,係数を決めるために使われる訓練集合と,それとは別の確認用集合(検証用集合; validation set)に分ける方法がある.確認用集合はホールドアウト集合(hold-out set)とも呼ばれ,モデルの複雑さ(または)を最適化するのに使われる.
これまでの多項式曲線フィッティングの議論は,主として直感に訴えるものであった. ここからは,確率論的な議論の枠組みを使い,パターン認識における問題を解くためのより原理的なアプローチを探っていく.
以上.