Pythonで実装 PRML 第3章 ベイズ線形回帰
パターン認識と機械学習、第3章からベイズ線形回帰が予測分布を収束させていく様子を再現します。
やっていることは、1章でやったベイズ推定とほぼ一緒で、基底関数がガウス関数(3.4)になっているだけです。
ですので、コードもほぼほぼ流用してます。
1章と比較すると、3章では数式の表現がベクトルを扱うように書きなおされてますが、この辺何か意図あるのでしょうか。
(3.16)の特徴ベクトルの概念をここで紹介したからこの書き方なのかと思ってますが、どうなんでしょう。
実装の大まかな流れ
①求めたいのは予測分布(3.57)。
②基底関数は(3.4)の。
基底関数は9個と指定されているがで、配置の仕方に指示がないので均等に分布するように配置。
③予測分布の平均、分散を求める式は下記の3つなので、これをそれぞれ実装して予測分布を求める。
コード
図を描くにあたって、サンプルデータが増えるに従って収束する様子を@naoya_tさんの記事にあるもの以上に美しく描くアイデアを出す自信が無かったので、ほぼ写経な感じで参考にさせて頂きました。ありがとうございました。
import numpy as np from numpy.linalg import inv import pandas as pd from pylab import * import matplotlib.pyplot as plt %matplotlib inline def func(x_train, y_train): # (3.4) Gaussian basis function def phi(s, mu_i, M, x): return np.array([exp(-(x - mu)**2 / (2 * s**2)) for mu in mu_i]).reshape((M, 1)) #(3.53)' ((1.70)) Mean of predictive distribution def m(x, x_train, y_train, S): sum = np.array(zeros((M, 1))) for n in xrange(len(x_train)): sum += np.dot(phi(s, mu_i, M, x_train[n]), y_train[n]) return Beta * phi(s, mu_i, M, x).T.dot(S).dot(sum) #(3.59)'((1.71)) Variance of predictive distribution def s2(x, S): return 1.0/Beta + phi(s, mu_i, M, x).T.dot(S).dot(phi(s, mu_i, M, x)) #(3.53)' ((1.72)) def S(x_train, y_train): I = np.identity(M) Sigma = np.zeros((M, M)) for n in range(len(x_train)): Sigma += np.dot(phi(s, mu_i, M, x_train[n]), phi(s, mu_i, M, x_train[n]).T) S_inv = alpha*I + Beta*Sigma S = inv(S_inv) return S #params for prior probability alpha = 0.1 Beta = 9 s = 0.1 #use 9 gaussian basis functions M = 9 # Assign basis functions mu_i = np.linspace(0, 1, M) S = S(x_train, y_train) #Sine curve x_real = np.arange(0, 1, 0.01) y_real = np.sin(2*np.pi*x_real) #Seek predictive distribution corespponding to entire x mean = [m(x, x_train, y_train, S)[0,0] for x in x_real] variance = [s2(x, S)[0,0] for x in x_real] SD = np.sqrt(variance) upper = mean + SD lower = mean - SD plot(x_train, y_train, 'bo') plot(x_real, y_real, 'g-') plot(x_real, mean, 'r-') fill_between(x_real, upper, lower, color='pink') xlim(0.0, 1.0) ylim(-2, 2) title("Figure 3.8") show() if __name__ == "__main__": # Maximum observed data points (underlining function plus some noise) x_train = np.linspace(0, 1, 26) #Set "small level of random noise having a Gaussian distribution" loc = 0 scale = 0.3 y_train = np.sin(2*np.pi*x_train) + np.random.normal(loc,scale,26) #Sample data pick up def randidx(n, k): r = range(n) shuffle(r) return sort(r[0:k]) for k in (1, 2, 5, 26): indices = randidx(size(x_train), k) func(x_train[indices], y_train[indices])
結果