![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2011/06/05からのアクセス回数 5913
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://sage.math.canterbury.ac.nz/home/pub/96/
また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
sage/PRML-逐次ベイズ学習では、 重みwの分布について解いていましたが、 実際には新しいｘの値たいするtを予測することが目的となります。 ここではPRML3.3.2の予測分布、図3.8
をSageを使って試してみます。
sage/PRML-線形回帰 で使ったデータと同じものを座標Xと目的値tにセットし、関数Φを定義します。
sageへの入力:
# PRML fig.3.8の再現
# PRMLのsin曲線のデータ
data = matrix([
[0.000000, 0.349486],
[0.111111, 0.830839],
[0.222222, 1.007332],
[0.333333, 0.971507],
[0.444444, 0.133066],
[0.555556, 0.166823],
[0.666667, -0.848307],
[0.777778, -0.445686],
[0.888889, -0.563567],
[1.000000, 0.261502],
]);
X = data.column(0)
t = data.column(1)
# データを増やす場合
# N = 25
# X = vector([random() for i in range(25)])
# t = vector([(sin(2*pi*x) + +gauss(0, 0.2)).n() for x in X.list()]);
sageへの入力:
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
data_plt = list_plot(zip(X, t))
sage/PRML- エビデンス近似 のようにガウス基底関数を定義しますが、今回は0から1の範囲で定義します。
近似に使う場合には、_phiのようにj=0が1となるような項を追加します。
sageへの入力:
from pylab import linspace
# 定数をセット
M=9
mu = linspace(0, 1, M)
s = 0.1
s_sq = (s)^2
# ガウス基底関数
def _phi_gauss(x, j):
return e^(-1*(x - mu[j])^2/(2* s_sq))
sageへの入力:
# ガウス基底関数で三角関数の例題を近似
# j=0に対応するため_phiを以下のように定義
def _phi(x, j):
if j == 0:
return 1
else:
return _phi_gauss(x, j-1)
# 初期化
alpha = 2
beta = 25
予測分布は、式(3.58) $$ p(t|x, t, \alpha, \beta) = \mathcal{N}(t|m_N^T\phi(x), \sigma_N^2(x)) $$ から計算し、その分散$\sigma_N^2(x)$は、 $$ \sigma_N^2(x) = \frac{1}{\beta} + \phi(x)^T S_N \phi(x) $$ で与えられます。
sageへの入力:
# 予測分布をプロットする
def _plot_predict(X, t):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()])
Phi_t = Phi.transpose()
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi).inverse()*Phi_t;
# 平均の重み
Wml = beta*Phi_dag * t
f = lambda x : (sum((Wml[i]*_phi(x, i)) for i in range(0, (M+1))))
# 分散(標準偏差)
def s(x):
phi_x = vector([_phi(x, i).n() for i in range(M+1)])
S = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi).inverse()
s_sqr = 1/beta + phi_x * S * phi_x
return sqrt(s_sqr)
data_plt = list_plot(zip(X, t))
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcolor='grey')
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcolor='grey')
y_plt = plot(f, [x, 0, 1], rgbcolor='red')
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymin=-1.5, ymax=1.5)
以下にM=1, 2, 4, 10の予測分布を示します。パラメータは、
としました。
sageへの入力:
# 1点の場合 X_ = vector([X[4]]) t_ = vector([t[4]]) _plot_predict(X_, t_)
sageへの入力:
# 2点の場合 X_ = vector([X[4], X[6]]) t_ = vector([t[4], t[6]]) _plot_predict(X_, t_)
sageへの入力:
# 4点の場合 X_ = vector([X[0], X[4], X[6], X[9]]) t_ = vector([t[0], t[4], t[6], t[9]]) _plot_predict(X_, t_)
sageへの入力:
# 10点の場合 _plot_predict(X, t)
皆様のご意見、ご希望をお待ちしております。
