はじめに
『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでのスクラッチ実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。
この記事は、3.1.1項「最尤推定と最小二乗法」の内容です。線形基底関数モデルの最尤推定をPythonで実装します。
【数式読解編】
【前節の内容】
【他の節一覧】
【この節の内容】
3.1.1 線形基底関数モデル
線形基底関数モデルの最尤推定を実装します。
利用するライブラリを読み込みます。
# 3.1.1項で利用するライブラリ import numpy as np import matplotlib.pyplot as plt
・基底関数の作成
まずは、3種類の基底関数(多項式基底関数・ガウス基底関数・シグモイド基底関数)を実装します。基底関数については「【Python】3.1.0：基底関数の作図【PRMLのノート】 - からっぽのしょこ」を参照してください。
多項式基底関数を作成します。
# 多項式基底関数を作成 def phi_polynomial(x, j): return x**j
多項式$\sum_{m=0}^{M-1} w_m \phi_m(x)$について、$m$番目の基底関数を次の式とします。
phi_polynomial()を使って、多項式基底関数の計算を行う計画行列を作成します。
# 多項式基底関数の計画行列を作成:(M > 2) def Phi_polynomial(x_n, M, _x_n=None): # 変数を初期化 phi_x_nm = np.ones((len(x_n), M)) # 列ごとに多項式基底関数による変換 for m in range(1, M): phi_x_nm[:, m] = phi_polynomial(x_n, m) return phi_x_nm
計画行列$\boldsymbol{\Phi}$は、$\boldsymbol{\Phi}_{m,n} = \phi_m(x_n)$である行列です。
$N$次元ベクトルの入力変数$\mathbf{x}$を、非線形関数(基底関数)により値を変換し、$N \times M$の行列に変換します。
ただし、0番目の基底関数(行列の0列目)については、$\phi_0(x) = 1$とします。そのため、phi_x_nmの全ての要素を1として作成しておき、(0から数えて)1列目以降を順番に計算して代入します。多項式基底関数においては$x^0 = 1$なので、0列目から計算しても問題ありません。
第3引数の_x_nは、他の関数と対応させるためのダミーです。
続いて、ガウス基底関数を作成します。
# ガウス基底関数を作成 def phi_gauss(x, mu, s): a_n = -0.5 * (x - mu)**2 / s**2 return np.exp(a_n)
ガウス基底関数は、次の式で定義されます。
基底関数ごとに、パラメータ$\mu_m, s_m$を持ちます。
phi_gauss()を使って、計画行列を作成します。
# ガウス基底関数の計画行列を作成:(M > 2) def Phi_gauss(x_n, M, _x_n=None): # パラメータ設定用の入力値を設定 if _x_n is None: _x_n = x_n.copy() # M-1個のパラメータを作成:(M > 1) s = np.std(_x_n) # 標準偏差で固定 if M == 2: # 入力値の平均 mu_m = np.array([np.mean(_x_n)]) elif M == 3: # 調整幅を指定 sgm = s * 1.0 # 平均を中心に標準偏差のn倍の範囲 mu_m = np.array([np.mean(_x_n) - sgm, np.mean(_x_n) + sgm]) elif M > 3: # 調整幅を指定 sgm = s * 0.25 # 最小値から最大値を等分 mu_m = np.linspace(np.min(_x_n) + sgm, np.max(_x_n) - sgm, num=M-1) # 変数を初期化 phi_x_nm = np.ones((len(x_n), M)) # 列ごとにガウス基底関数による変換 for m in range(1, M): phi_x_nm[:, m] = phi_gauss(x_n, mu_m[m-1], s) return phi_x_nm
0番目の$\phi_0(x)$を除くと、基底関数は$\phi_1(x), \cdots, \phi_{M-1}(x)$の$M - 1$個です。対応する$M - 1$個のパラメータを$\boldsymbol{\mu} = (\mu_1, \cdots, \mu_{M-1})$で表します。$s_m$については、この例では全ての基底関数で1つの値$s$とします。
関数の実行時に、観測データ$\mathbf{x}$から$\boldsymbol{\mu}, s$の値を設定します。
$s$は、$\mathbf{x}$の標準偏差とします。標準偏差はnp.std()で計算できます。
$\boldsymbol{\mu}$は、$\mathbf{x}$の最小値から最大値までを$M - 1$個に等分した値とします。ただし、$M = 2, 3$のときは$\boldsymbol{\mu}$が$\mathbf{x}$の最小値と最大値になってしまうので、条件分岐で設定方法を変えることにします。
$M = 2$のときは、基底関数のパラメータは1つなので、$\mathbf{x}$の平均値とします。$M = 3$のときは、$\mathbf{x}$の平均値から標準偏差$s$の定数倍を引いた値と足した値とします。
また、作図時には観測データ$\mathbf{x}$とは別の値$x_{*}$に対して非線形変換を行います。その場合、基底関数のパラメータ$\boldsymbol{\mu}, s$が変わってしまいます。そこで、第3引数に観測データ$\mathbf{x}$を渡すことで、学習時と同じ値を設定できるように実装しています。
(もっとスマートで適切な設定方法があれば教えて下さい。)
# ロジスティックシグモイド関数を作成 def sigmoid(x): return 1.0 / (1.0 + np.exp(-x)) # シグモイド基底関数を作成 def phi_sigmoid(x, mu, s): a_n = (x - mu) / s return sigmoid(a_n)
シグモイド基底関数は、ロジスティックシグモイド関数
を用いて、次の式で定義されます。
Phi_sigmoid()を使って、計画行列を作成します。
# シグモイド基底関数の計画行列を作成:(M > 2) def Phi_sigmoid(x_n, M, _x_n=None): # パラメータ設定用の入力値を設定 if _x_n is None: _x_n = x_n.copy() # M-1個のパラメータを作成:(M > 1) s = np.std(_x_n) * 0.5 # 標準偏差で固定 if M == 2: # 入力値の平均 mu_m = np.array([np.mean(_x_n)]) elif M == 3: # 調整幅を指定 sgm = s * 1.0 # 平均を中心に標準偏差のn倍の範囲 mu_m = np.array([np.mean(_x_n) - sgm, np.mean(_x_n) + sgm]) elif M > 3: # 調整幅を指定 sgm = s * 0.25 # 最小値から最大値を等分 mu_m = np.linspace(np.min(_x_n) + sgm, np.max(_x_n) - sgm, num=M-1) # 変数を初期化 phi_x_nm = np.ones((len(x_n), M)) # 列ごとにシグモイド基底関数による変換 for m in range(1, M): phi_x_nm[:, m] = phi_sigmoid(x_n, mu_m[m-1], s) return phi_x_nm
ガウス基底関数と同様に$\boldsymbol{\mu}, s$を設定します。
・モデルの設定
次に、データを生成する関数を設定します。
データ生成関数(推定対象の関数)を指定します。
# 真の関数を作成 def y_true(x): # 計算式を指定 return np.sin(2.0 * np.pi * x)
この例では、$y = \sin(2 \pi x)$とします。
作図用にx軸の点を作成します。
# 作図用のxの値を指定 x_vals = np.linspace(0.0, 1.0, num=101)
グラフとして表示する範囲を指定して、x軸の点を作成します。
真のパラメータを指定します。
# 真の精度パラメータを指定 beta_true = 3.0 # 真の標準偏差を計算 sigma_true = np.sqrt(1.0 / beta_true) sigma_true print(sigma_true)
0.5773502691896257
ガウスノイズ$\epsilon \sim \mathcal{N}(\epsilon | 0, \beta^{-1})$の精度(分散の逆数)$\beta$をbeta_trueとして指定します。
また、標準偏差$\sigma = \sqrt{\frac{1}{\beta}}$を計算して、sigma_trueとします。(標準偏差を直接指定しても問題ありません。)
真の関数の出力を計算します。
# 真のモデルを計算 y_true_vals = y_true(x_vals) print(y_true_vals[:5])
[0. 0.06279052 0.12533323 0.18738131 0.24868989]
始めに作成した関数y_ture()で計算します。
真の関数のグラフを作成します。
# 真のモデルを作図 plt.figure(figsize=(12, 8)) plt.plot(x_vals, y_true_vals, color='darkturquoise', label='true model') # 真のモデル plt.fill_between(x=x_vals, y1=y_true_vals-2.0*sigma_true, y2=y_true_vals+2.0*sigma_true, color='darkturquoise', alpha=0.2, linestyle='--') # 真のノイズ範囲 plt.xlabel('x') plt.ylabel('t') plt.suptitle('$t = \\sin(2 \pi x)$', fontsize=20) plt.title('$\\beta=' + str(beta_true) + '$', loc='left') plt.grid() plt.legend() plt.show()
ノイズの範囲の目安として、期待値$y$から標準偏差の2倍の範囲を計算して、fill_between()で描画します。
この曲線を推定(近似)するのが目標です。
・データの生成
学習に利用するデータを人工的に作成します。
観測データを作成します。
# データ数を指定 N = 100 # データを生成 x_n = np.random.uniform(low=np.min(x_vals), high=np.max(x_vals), size=N) t_n = y_true(x_n) + np.random.normal(loc=0.0, scale=sigma_true, size=N)
x_valsの最小値から最大値までの一様乱数に従って入力データ$\mathbf{x} = \{x_1, \cdots, x_N\}$を生成します。
真の関数y_true()の出力に平均0・精度$\beta$のガウス乱数$\epsilon_n$を加えて、目標値$\mathbf{t} = \{t_1, \cdots, t_N\}$を作成します。
観測データの散布図を真のモデルに重ねて描画します。
# 観測データの散布図を作成 plt.figure(figsize=(12, 8)) plt.plot(x_vals, y_true_vals, color='darkturquoise', label='true model') # 真のモデル plt.fill_between(x=x_vals, y1=y_true_vals-2.0*sigma_true, y2=y_true_vals+2.0*sigma_true, color='darkturquoise', alpha=0.2, linestyle='--') # 真のノイズ範囲 plt.scatter(x_n, t_n, color='orange') # 観測データ plt.xlabel('x') plt.ylabel('t') plt.suptitle('$t_n = \\sin(2 \pi x_n) + \epsilon_n$', fontsize=20) plt.title('$N=' + str(N) + ', \\beta=' + str(beta_true) + '$', loc='left') plt.grid() plt.legend() plt.show()
x軸の各点においてy軸方向に、平均$\sin(2 \pi x)$・精度$\beta$のガウス分布$\mathcal{N}(t | \sin(2 \pi x), \beta^{-1})$になっています。
各データが概ね$2 \sigma$の範囲に収まっているのを確認できます。
この観測データ$\mathbf{x}, \mathbf{t}$からモデルを推定します。
・最尤推定
3つの基底関数を用いて、線形基底関数モデルの最尤推定を行います。
・多項式基底関数の場合
多項式基底関数を用いて推定します。
利用する基底関数を設定します。
# 基底関数を指定
phi = Phi_polynomial
関数phi()として、多項式基底関数Phi_polynomial()を複製します。
パラメータの最尤推定量を計算します。
# 基底関数の数を指定 M = 6 # 基底関数により入力を変換 phi_x_nm = phi(x_n, M) # 重みパラメータの最尤解を計算 w_ml_m = np.linalg.inv(np.dot(phi_x_nm.T, phi_x_nm)).dot(phi_x_nm.T).dot(t_n.reshape(-1, 1)).flatten() print(w_ml_m) # 分散パラメータの最尤解を計算 sigma2_ml = np.sum((t_n - np.dot(phi_x_nm, w_ml_m.reshape(-1, 1)).flatten())**2) / N print(sigma2_ml) # 精度パラメータの最尤解を計算 beta_ml = 1.0 / sigma2_ml print(beta_ml)
[ 0.2264529 5.66131787 -7.90246035 -27.17944274 44.18510743
-15.02524282]
0.28994217378387593
3.448963587978767
尤度関数を最大化する重みパラメータ$\mathbf{w}_{\mathrm{ML}}$と分散パラメータ$\frac{1}{\beta_{\mathrm{ML}}}$を次の式で計算します。
ただし、$t_n - \mathbf{w}^{\top} \phi(\mathbf{x}_n)$の計算を$N$個同時に行うため、$\mathbf{t} - \boldsymbol{\Phi}^{\top} \mathbf{w}$で計算しています。
sigma2_hatの逆数が精度パラメータの最尤推定量$\beta_{\mathrm{ML}}$になります。これは計算する必要はありません。
推定したパラメータを用いて、モデルの出力を計算します。
# 作図用のxを変換 phi_x_valsm = phi(x_vals, M) # 回帰曲線を計算 y_ml_vals = np.dot(phi_x_valsm, w_ml_m.reshape(-1, 1)).flatten()
推定した重みパラメータ$\mathbf{w}_{\mathrm{ML}}$を使って、未知の入力値(x軸の値)$x_{*}$に対する出力値(y軸の値)$t_{*}$を次の式で計算します。
下の式でx_valsの全ての値を同時に計算します。
回帰曲線(近似曲線)を真のモデルと重ねて描画します。
# 回帰曲線を作図 plt.figure(figsize=(12, 8)) plt.plot(x_vals, y_true_vals, color='darkturquoise', label='true model') # 真のモデル plt.fill_between(x=x_vals, y1=y_true_vals-2.0*sigma_true, y2=y_true_vals+2.0*sigma_true, color='darkturquoise', alpha=0.2, linestyle=':') # 真のノイズ範囲 plt.plot(x_vals, y_ml_vals, color='blue', label='$w^T \phi(x)$') # 推定したモデル plt.fill_between(x=x_vals, y1=y_ml_vals-2.0*np.sqrt(sigma2_ml), y2=y_ml_vals+2.0*np.sqrt(sigma2_ml), color='blue', alpha=0.2, linestyle=':') # 推定したノイズ範囲 plt.scatter(x_n, t_n, color='orange') # 観測データ plt.xlabel('x') plt.ylabel('t') plt.suptitle('Linear Basis Function Model : Polynomial', fontsize=20) plt.title('$N=' + str(N) + ', M=' + str(M) + ', w=(' + ', '.join([str(w) for w in np.round(w_ml_m, 2)]) + ')' + ', \\beta=' + str(np.round(beta_ml, 2)) + '$', loc='left') plt.legend() plt.grid() plt.show()
推定したモデル(青色の曲線)がデータに当て嵌まり、また真のモデル(青緑色の曲線)を近似できているのが分かります。真のモデルと同様に、推定した分散パラメータsigma2_mlを用いて精度を可視化しています。
以上で、多項式基底関数を用いた線形基底関数モデルの最尤推定を実装できました。
最後に、基底関数と重みパラメータと近似曲線の関係を確認します。
$M$個の基底関数$\phi_m(x)$をそれぞれ作図します。
# 基底関数を作図 plt.figure(figsize=(12, 8)) for m in range(M): plt.plot(x_vals, phi_x_valsm[:, m], label='$\phi_' + str(m) + '(x)$') # 基底関数 plt.xlabel('$x$') plt.ylabel('$\phi_m(x)$') plt.suptitle('Polynomial Basis Function', fontsize=20) plt.legend() plt.grid() plt.show()
これは、3.1.0項で確認した多項式基底関数のグラフです。
$M$個の基底関数$\phi_m(x)$をそれぞれ対応する重み$w_m$と掛け合わせて、近似曲線$\mathbf{w}_{\mathrm{ML}}^{\top} \boldsymbol{\phi}(x)$と重ねて描画します。
# 重み付き基底関数を作図 plt.figure(figsize=(12, 8)) for m in range(M): plt.plot(x_vals, w_ml_m[m] * phi_x_valsm[:, m], linestyle='--', label='$w_' + str(m) + '\phi_' + str(m) + '(x)$') # 重み付き基底関数 plt.plot(x_vals, y_ml_vals, color='blue', label='$w^T \phi(x)$') plt.xlabel('x') plt.ylabel('t') plt.suptitle('Polynomial Basis Function', fontsize=20) plt.title('w=(' + ', '.join([str(w) for w in np.round(w_ml_m, 2)]) + ')', loc='left') plt.legend() plt.grid() plt.ylim(-5.1, 5.1) plt.show()
$M$個の基底関数(破線)をy軸方向に和をとった値が近似曲線(青色の実線)になります。重みパラメータの大小によって各基底関数の影響度合い(グラフの形状)が変わっているのを確認できます。また、負の値であればy軸方向に反転します。
１つの直線(バイアス)$w_0 = w_0 \phi_0(x)$と$M - 1$個の曲線$w_1 \phi_1(x), \cdots, w_{M-1} \phi_{M-1}(x)$によって、近似曲線(青色の実線)が構成されているのが分かります。
・ガウス基底関数の場合
続いて、ガウス基底関数を用いて行います。
ガウス基底関数の計画行列Phi_gauss()をphi()として設定します。
# 基底関数を指定 phi = Phi_gauss # 基底関数の数を指定 M = 6 # 基底関数により入力を変換 phi_x_nm = phi(x_n, M) # 重みパラメータの最尤解を計算 w_ml_m = np.linalg.inv(np.dot(phi_x_nm.T, phi_x_nm)).dot(phi_x_nm.T).dot(t_n.reshape(-1, 1)).flatten() print(w_ml_m) # 分散パラメータの最尤解を計算 sigma2_ml = np.sum((t_n - np.dot(phi_x_nm, w_ml_m.reshape(-1, 1)).flatten())**2) / N print(sigma2_ml) # 精度パラメータの最尤解を計算 beta_ml = 1.0 / sigma2_ml print(beta_ml) # 作図用のxを変換 phi_x_valsm = phi(x_vals, M, x_n) # 回帰曲線を計算 y_ml_vals = np.dot(phi_x_valsm, w_ml_m.reshape(-1, 1)).flatten()
[ -1.48467994 3.1923181 -4.69238328 9.88594362 -10.44907788 5.84362878] 0.293458414521164 3.407637847535228
phi()でx軸の値x_valsを変換する際に、第3引数にx_nを指定する必要があります。
他の処理は先ほどと同じです。
・シグモイド基底関数の場合
最後に、シグモイド基底関数を用いて行います。
シグモイド基底関数の計画行列Phi_sigmoid()をphi()として設定します。
# 基底関数を指定 phi = Phi_sigmoid # 基底関数の数を指定 M = 6 # 基底関数により入力を変換 phi_x_nm = phi(x_n, M) # 重みパラメータの最尤解を計算 w_ml_m = np.linalg.inv(np.dot(phi_x_nm.T, phi_x_nm)).dot(phi_x_nm.T).dot(t_n.reshape(-1, 1)).flatten() print(w_ml_m) # 分散パラメータの最尤解を計算 sigma2_ml = np.sum((t_n - np.dot(phi_x_nm, w_ml_m.reshape(-1, 1)).flatten())**2) / N print(sigma2_ml) # 精度パラメータの最尤解を計算 beta_ml = 1.0 / sigma2_ml print(beta_ml) # 作図用のxを変換 phi_x_valsm = phi(x_vals, M, x_n) # ガウス基底関数・シグモイド基底関数 # 回帰曲線を計算 y_ml_vals = np.dot(phi_x_valsm, w_ml_m.reshape(-1, 1)).flatten()
[-1.16587199 4.96087243 -4.10636291 0.66709391 -5.4687201 7.5356283 ]
0.2892052141418903
3.457752319463294
こちらもphi()の第3引数にx_nを設定します。
以上で、3種類の基底関数による線形基底関数モデルの最尤推定を実装できました。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.
おわりに
Rで実装したときから色々アップデートしたので、Rの方を書き直したくなってきた。
【次節の内容】
