#freeze
[[FrontPage]]
#contents
2011/06/05からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://sage.math.canterbury.ac.nz/home/pub/96/
また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* 第3章-予測分布をSageで試す [#k3f9ca83]
[[sage/PRML-逐次ベイズ学習]]では、
重みwの分布について解いていましたが、
実際には新しいｘの値たいするtを予測することが目的となります。
ここではPRML3.3.2の予測分布、図3.8
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.8a.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.8b.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.8c.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.8d.jpg,center,30%);
をSageを使って試してみます。
*** 準備 [#u42a23ce]
[[sage/PRML-線形回帰]]
で使ったデータと同じものを座標Xと目的値tにセットし、関数Φを定義します。
sageへの入力:
#pre{{
# 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への入力:
#pre{{
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
data_plt = list_plot(zip(X, t))
}}
*** ガウス基底関数の定義 [#y1690db6]
[[sage/PRML- エビデンス近似]]
のようにガウス基底関数を定義しますが、今回は0から1の範囲で定義します。
近似に使う場合には、_phiのようにj=0が1となるような項を追加します。
sageへの入力:
#pre{{
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への入力:
#pre{{
# ガウス基底関数で三角関数の例題を近似
# j=0に対応するため_phiを以下のように定義
def _phi(x, j):
if j == 0:
return 1
else:
return _phi_gauss(x, j-1)
# 初期化
alpha = 2
beta = 25
}}
** 予測分布 [#dee55f44]
予測分布は、式(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への入力:
#pre{{
# 予測分布をプロットする
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)
}}
*** 予測分布のプロット [#h0cffc72]
以下にM=1, 2, 4, 10の予測分布を示します。パラメータは、
- \(\alpha = 2\)
- \(\beta = 25\)
としました。
sageへの入力:
#pre{{
# 1点の場合
X_ = vector([X[4]])
t_ = vector([t[4]])
_plot_predict(X_, t_)
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# 2点の場合
X_ = vector([X[4], X[6]])
t_ = vector([t[4], t[6]])
_plot_predict(X_, t_)
}}
&ref(sage0-1.png);
sageへの入力:
#pre{{
# 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_)
}}
&ref(sage0-2.png);
sageへの入力:
#pre{{
# 10点の場合
_plot_predict(X, t)
}}
&ref(sage0-3.png);
** コメント [#ndf1cb8f]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[1])
#vote(おもしろかった[3],そうでもない[0],わかりずらい[1])
皆様のご意見、ご希望をお待ちしております。
- 自分としては、PRMLの図3.8をある程度再現できた思っていたのですが、分かりづらいとのコメントはちょっときつかったです。 -- [[竹本 浩]] &new{2013-06-25 (火) 20:49:17};
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
