sage/PRML- エビデンス近似
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/04からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/95/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第3章-エビデンス近似をSageで試す [#vecd04fc]
[[sage/PRML-線形回帰]]では、
与えられたモデルパラメータM, α、βに対して解いていましたが、
実際には自分でベストなモデルパラメータM, α、βを求める必要...
[[パターン認識と機械学習>http://www.amazon.co.jp/dp/44311...
の３章ではエビデンス関数を使って最適なパラメータを求める...
ここでは、Sageを使ってエビデンス関数の評価、最適なモデル...
** エビデンス関数の評価 [#xf7b5eab]
エビデンス関数の対数表現は、式(3.86)
$$
\ln p(t|\alpha, \beta) = \frac{M}{2} \ln \alpha + \frac{...
$$
で与えられます。
ここで、$E(m_N)$は、式(3.82)
$$
E(m_N) = \frac{\beta}{2}|| t = \Phi m_N ||^2 + \frac{\alp...
$$
$m_N$は、式(3.84)
$$
m_N = \beta A^{-1} \Phi^T t
$$
Aは、式(3.81)
$$
A = \alpha I + \beta \Phi^T \Phi
$$
目指すは、図3.14です。
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
*** 準備 [#ac2b67a1]
[[sage/PRML-線形回帰]]
で使ったデータと同じものを座標Xと目的値tにセットし、関数Φ...
sageへの入力:
#pre{{
# PRML fig.3.14の再現
# 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...
}}
sageへの入力:
#pre{{
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** エビデンスの計算 [#p34968e0]
まずは、α、βを固定値
- $\alpha = 5 \times 10^{-3}$
- $\beta = 11.1$
として、多項式の次数をMを0から9まで変えた値を計算し、プロ...
一番良いエビデンスは、M=4あたりになります。
理由は分かりませんが、図3.14のようなM=3以降の急激な下降は...
sageへの入力:
#pre{{
# α=5*10^-3, β=11.1を固定して計算すると
ln_p_List = []
N = len(t)
alpha = 5*10^-3
beta = 11.1
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta*A.inverse() * Phi_t * t
# エビデンスの対数
res = (t - Phi*m_N)
E_mN = beta/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/...
ln_p_List += [ln_p_t]
list_plot(ln_p_List, ymin = -26, ymax = -5).show()
# 結果は、M=3以降でfig.3.14のような急激な降下は見られない
}}
&ref(sage0.png);
*** β_MLを使ってエビデンスを計算 [#rf539535]
そこで、βを平均値の分散$\beta_{ML}$を使って計算してみまし...
$$
\frac{1}{\beta_{ML}} = \frac{1}{N} \sum (y(x_n, W_ML) - t...
$$
期待に反し、M=0でも高い値となり、あまり良くありません。
sageへの入力:
#pre{{
# 1/β_ML = 1/N Σ{y(x_n,W_ML) - t_n}^2から計算したβ_MLでエ...
beta_ML_List = []
ln_p_List = []
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# エビデンスの対数
res = (t - Phi*m_N)
beta_ML = N / (res*res)
E_mN = beta_ML/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta_ML) - E_mN -...
ln_p_List += [ln_p_t]
beta_ML_List += [beta_ML]
list_plot(ln_p_List, ymin = -26, ymax = 0).show()
#print beta_ML_List
}}
&ref(sage0-1.png);
*** どうして図3.14と合わないのか？ [#q2979465]
どうして図3.14のような結果にならないのかと調べたのですが...
方がいらして、私と同じ傾向になったとの記述がありました。
[[基底関数を色々変えて、線形回帰のエビデンスを計算してみ...
私の推測では、ln |A|の計算をA.norm()と間違えたのではない...
sageへの入力:
#pre{{
# 当初、式3.85, 3.86のミスプリと考えたのですが、
# →http://d.hatena.ne.jp/n_shuyo/20090715/evidence#c
# n_shuyoさんのコメントで 正則行列の行列式は、逆行列の...
# 式3.85, 3.86 は正しい！
# 私の推測：図3.14の計算で ln |A|とするところをln A.norm(...
ln_p_List = []
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# エビデンスの対数
res = (t - Phi*m_N)
E_mN = beta/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/...
ln_p_List += [ln_p_t]
list_plot(ln_p_List)
}}
&ref(sage0-2.png);
*** 最適なα、βを求める [#tf19c14f]
α、βを固定にした場合に、エビデンスが最大はM=4だったので、...
最適なα、βを求めてみます。
エビデンスを最大にするαは、式(3.92)
$$
\alpha = \frac{\gamma}{m_N^T m_N}
$$
ここで\(\gamma\)は、式(3.87)
$$
(\beta \Phi^T \Phi) u_i = \lambda_i u_i
$$
の固有値から、式(3.91)
$$
\gamma = \sum_{i} \frac{\lambda_i}{\alpha + \lambda_i}
$$
で求まります。
βについても、式(3.95)
$$
\frac{1}{\beta} = \frac{1}{N - \gamma} \sum_{n=1}^{N} ( t...
$$
\(\gamma\)は、\(\alpha\)に依存し、\(m_N\)も\(\alpha\)に依...
できません。
+ 適当な \(\alpha, \beta\)を初期値とする
+ \(m_N, \gamma\)を求め
+ 式(3.95)から\(\beta\)を求める
+ 式(3.92)から\(\alpha\)を求め、2番目からを繰り返す
以下の例では、
- $\alpha = 5 \times 10^{-3}$
- $\beta = 11.1$
から20回の計算を繰り返し、\(\gamma, \alpha, \beta\)の値を...
sageへの入力:
#pre{{
# M=4の場合のα、βを決める
M=4
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# Φ^T Φは固定なので、先に固有値を計算
B = Phi_t*Phi
lambs = B.eigenvalues()
# 初期化
alpha = 5*10^(-3)
beta = 11.1
gamma = M
for i in range(20):
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
gamma = sum((lamb*beta) /(alpha + (lamb*beta)) for la...
alpha = gamma/(m_N*m_N)
res = (t - Phi*m_N)
beta = (N - gamma) / (res*res)
print gamma, alpha, beta
}}
sageの出力：
#pre{{
4.02295319874 0.0160029572678 15.5703330431
3.86957616939 0.0192233025727 14.6079344569
3.81804019954 0.0208658867531 14.0162355337
3.79079625631 0.021876053777 13.6574090928
3.77410819657 0.0225480023858 13.4248103487
3.76308615077 0.0230152846069 13.2666278994
3.75546956914 0.0233495381967 13.1554382026
3.75004882783 0.0235932327887 13.0754487343
3.74611222743 0.0237732924978 13.0169432556
3.74321234677 0.0239076177514 12.973632867
3.7410540217 0.0240085312871 12.9412853436
3.73943541844 0.0240847392383 12.9169657574
3.7382147331 0.0241425141272 12.8985911035
3.73729026269 0.0241864427411 12.8846563016
3.73658790674 0.0242199173426 12.8740587014
3.73605302024 0.0242454685376 12.8659818046
3.73564493065 0.0242649966702 12.8598160009
3.73533314922 0.0242799360403 12.8551032449
3.73509469602 0.0242913734086 12.8514976848
3.73491217789 0.0243001346625 12.848737196
}}
*** γと2αEw(mN)のグラフ [#qa302662]
M=4の多項式回帰に対して、図3.16(a)
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
のように\(\gammaと2\alpha E_w(m_N)\)のグラフを\(\ln \alph...
形が図3.16とは異なり、曲線の交差する点もα=0.024だとln αの...
交差しています。（何かちがう！）
sageへの入力:
#pre{{
# γをln(a)の関数として定義
def g(ln_a):
return abs(sum((lamb*beta) /(e^ln_a + (lamb*beta)) fo...
a_plt = plot(g, [x,-6, 6], color='red')
}}
sageへの入力:
#pre{{
# 2αEw(m_N)をln(a)の関数として定義
def aEw(ln_a):
a = (e^ln_a).n()
A = a*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# m_Nいくつかの点が複素数になるためabs(a*m_N*m_N)とした
return abs(a*m_N*m_N)
Ew_plt = plot(aEw, [x,-6, 6], color='blue')
(a_plt + Ew_plt).show()
}}
&ref(sage0-3.png);
*** 残差と対数エビデンスのグラフ [#z3abc719]
同様に、M=4の多項式回帰に対して、図3.16(b)
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
のように残差と対数エビデンスを\(\ln \alpha\)の関係を図化...
やはり形が、図3.16とは異なります。
sageへの入力:
#pre{{
b_phi_2 = beta*Phi_t*Phi
Phi_t_t = Phi_t* t
def ln_pt(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
E_mX = beta/2 * res*res + a/2 * m_N*m_N
# E_mXいくつかの点が複素数になるためabs(E_mX)とした
return (M/2 * ln_a + N/2 * ln(beta) - abs(E_mX) -1/2*...
#
def res2(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
return abs(res*res)
}}
sageへの入力:
#pre{{
ln_pt_plt = plot(ln_pt, [x,-6, 6], color='red')
res_2_plt = plot(res2, [x,-6, 6])
(ln_pt_plt + res_2_plt).show()
}}
&ref(sage0-4.png);
*** M=4に対する最適な解 [#bff6d673]
M=4の多項式回帰に対する最適なα、βを使って\(sin(2 \pi x)\)...
sageへの入力:
#pre{{
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
data_plt = list_plot(zip(X, t))
}}
sageへの入力:
#pre{{
# 求まったα、βでベイズ的なフィッティングを再計算
alpha = 0.0243001346625
beta = 12.848737196
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)))
# 分散
def s(x):
phi_x = vector([x^i for i in range(M+1)])
S = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi)....
s_sqr = 1/beta + phi_x * S * phi_x
return sqrt(s_sqr)
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red')
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-5.png);
M=9の最適解（図1.17）
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
と比べると若干フィッティングが良くありませんが、与えられ...
** どうしてエビデンス関数が合わないのか [#ic963cc3]
つぎにどうして、図3.16と合わないのか、調べてみます。
図3.16の説明をよく読むと、基底関数が多項式ではなく、ガウ...
*** 基底関数のチェック [#z831c5cd]
ガウス基底関数がどのようなものなのか、図3.1を再現しながら...
以下に、
- 多項式の基底関数
- ガウス基底関数
$$
\phi_j(x) = exp \left \{ - \frac{(x - \mu_j)^2}{2 s^2} \r...
$$
- シグモイド基底関数
$$
\phi_j(x) = \sigma \left( \frac{x - \mu_j}{s} \right)
$$
$$
\sigma(a) = \frac{1}{1 + exp(-a)}
$$
を表示します。
ガウス基底関数、シグモイド基底関数では、μの値として-1から...
とり、それを基底関数の数で等分した値をとるみたいです。
[[線形回帰モデルとか>http://d.hatena.ne.jp/smly/20100630/...
sageへの入力:
#pre{{
# 基底関数のチェック
from pylab import linspace
M = 11
mu = linspace(-1, 1, M)
s = 0.1
s_sq = (s)^2
# 多項式基底関数定義
def _phi_poly(x, j):
return x^j
# ガウス基底関数
def _phi_gauss(x, j):
return e^(-1*(x - mu[j])^2/(2* s_sq))
# ロジスティック基底関数
def _logi_sig(x):
return 1/(1 + e^(-1*x))
def _phi_sigmoid(x, j):
return _logi_sig((x - mu[j])/s)
}}
sageへの入力:
#pre{{
x = var('x')
poly_plt =Graphics()
for j in range(1, M+1):
f = lambda x : _phi_poly(x, j)
poly_plt += plot(f, [x, -1, 1])
poly_plt.show()
}}
&ref(sage0-6.png);
sageへの入力:
#pre{{
gauss_plt =Graphics()
for j in range(M):
f = lambda x : _phi_gauss(x, j)
gauss_plt += plot(f, [x, -1, 1])
gauss_plt.show()
}}
&ref(sage0-7.png);
sageへの入力:
#pre{{
sigmoid_plt =Graphics()
for j in range(M):
f = lambda x : _phi_sigmoid(x, j)
sigmoid_plt += plot(f, [x, -1, 1])
sigmoid_plt.show()
}}
&ref(sage0-8.png);
*** ガウス基底関数を使ったγと2αEw(mN)のグラフ [#p159acdf]
M=9個のガウス基底関数を使ったγと2αEw(mN)のグラフを以下に...
関数の形状と交差点の位置が図3.16と合います。
sageへの入力:
#pre{{
# ガウス基底関数で三角関数の例題を近似
# j=0に対応するため_phiを以下のように定義
def _phi(x, j):
if j == 0:
return 1
else:
return _phi_gauss(x, j-1)
# 定数をセット
M=9
mu = linspace(0, 1, M)
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# Φ^T Φは固定なので、先に固有値を計算
B = Phi_t*Phi
lambs = B.eigenvalues()
# 初期化
alpha = 5*10^(-3)
beta = 11.1
gamma = M
for i in range(20):
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
gamma = sum((lamb*beta) /(alpha + (lamb*beta)) for la...
alpha = gamma/(m_N*m_N)
#res = (t - Phi*m_N)
#beta = (N - gamma) / (res*res)
print gamma, alpha, beta
}}
#pre{{
8.97526712356 2.74314585233 11.1000000000000
5.9805253866 5.87614916325 11.1000000000000
5.02659378174 6.10302152117 11.1000000000000
4.97666250731 6.12118986054 11.1000000000000
4.97273606757 6.12265708879 11.1000000000000
4.97241943407 6.1227756571 11.1000000000000
4.97239384954 6.12278523925 11.1000000000000
4.97239178194 6.12278601364 11.1000000000000
4.97239161484 6.12278607622 11.1000000000000
4.97239160134 6.12278608128 11.1000000000000
4.97239160025 6.12278608169 11.1000000000000
4.97239160016 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
}}
sageへの入力:
#pre{{
a_plt = plot(g, [x,-5, 5])
# 2αEw(m_N)をln(a)の関数として定義
def aEw(ln_a):
a = (e^ln_a).n()
A = a*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# m_Nいくつかの点が複素数になるためabs(a*m_N*m_N)とした
return abs(a*m_N*m_N)
Ew_plt = plot(aEw, [x,-5, 5], color='gray')
(a_plt + Ew_plt).show()
}}
&ref(sage0-9.png);
*** ガウス基底関数を使った残差と対数エビデンスのグラフ [#...
同様に、ガウス基底関数を使った残差と対数エビデンスのグラ...
残差の形状がやや違いますが、残差の最低値の付近に対数エビ...
部分はほぼ合っています。
sageへの入力:
#pre{{
b_phi_2 = beta*Phi_t*Phi
Phi_t_t = Phi_t* t
def ln_pt(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
E_mX = beta/2 * res*res + a/2 * m_N*m_N
return ((M/2 * ln_a + N/2 * ln(beta) - E_mX -1/2*ln(A...
# 残差の定義
def res2(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
return abs(res*res)
# 結果のプロット
ln_pt_plt = plot(ln_pt, [x,-5, 5], color='red')
res_2_plt = plot(res2, [x,-5, 5])
(ln_pt_plt + res_2_plt).show()
}}
&ref(sage0-10.png);
*** ガウス基底関数を回帰結果 [#rf54c993]
おまけに、ガウス基底関数の回帰結果と分散の分布、オリジナ...
プロットした図を示します。ガウス基底を使った回帰も良く合...
sageへの入力:
#pre{{
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t
f = lambda x : (sum((Wml[i]*_phi(x, i)) for i in range(0,...
# 分散
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)....
s_sqr = 1/beta + phi_x * S * phi_x
return sqrt(s_sqr)
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red')
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-11.png);
** コメント [#n7c3a48f]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/04からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/95/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第3章-エビデンス近似をSageで試す [#vecd04fc]
[[sage/PRML-線形回帰]]では、
与えられたモデルパラメータM, α、βに対して解いていましたが、
実際には自分でベストなモデルパラメータM, α、βを求める必要...
[[パターン認識と機械学習>http://www.amazon.co.jp/dp/44311...
の３章ではエビデンス関数を使って最適なパラメータを求める...
ここでは、Sageを使ってエビデンス関数の評価、最適なモデル...
** エビデンス関数の評価 [#xf7b5eab]
エビデンス関数の対数表現は、式(3.86)
$$
\ln p(t|\alpha, \beta) = \frac{M}{2} \ln \alpha + \frac{...
$$
で与えられます。
ここで、$E(m_N)$は、式(3.82)
$$
E(m_N) = \frac{\beta}{2}|| t = \Phi m_N ||^2 + \frac{\alp...
$$
$m_N$は、式(3.84)
$$
m_N = \beta A^{-1} \Phi^T t
$$
Aは、式(3.81)
$$
A = \alpha I + \beta \Phi^T \Phi
$$
目指すは、図3.14です。
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
*** 準備 [#ac2b67a1]
[[sage/PRML-線形回帰]]
で使ったデータと同じものを座標Xと目的値tにセットし、関数Φ...
sageへの入力:
#pre{{
# PRML fig.3.14の再現
# 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...
}}
sageへの入力:
#pre{{
# Φ関数定義
def _phi(x, j):
return x^j
}}
*** エビデンスの計算 [#p34968e0]
まずは、α、βを固定値
- $\alpha = 5 \times 10^{-3}$
- $\beta = 11.1$
として、多項式の次数をMを0から9まで変えた値を計算し、プロ...
一番良いエビデンスは、M=4あたりになります。
理由は分かりませんが、図3.14のようなM=3以降の急激な下降は...
sageへの入力:
#pre{{
# α=5*10^-3, β=11.1を固定して計算すると
ln_p_List = []
N = len(t)
alpha = 5*10^-3
beta = 11.1
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta*A.inverse() * Phi_t * t
# エビデンスの対数
res = (t - Phi*m_N)
E_mN = beta/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/...
ln_p_List += [ln_p_t]
list_plot(ln_p_List, ymin = -26, ymax = -5).show()
# 結果は、M=3以降でfig.3.14のような急激な降下は見られない
}}
&ref(sage0.png);
*** β_MLを使ってエビデンスを計算 [#rf539535]
そこで、βを平均値の分散$\beta_{ML}$を使って計算してみまし...
$$
\frac{1}{\beta_{ML}} = \frac{1}{N} \sum (y(x_n, W_ML) - t...
$$
期待に反し、M=0でも高い値となり、あまり良くありません。
sageへの入力:
#pre{{
# 1/β_ML = 1/N Σ{y(x_n,W_ML) - t_n}^2から計算したβ_MLでエ...
beta_ML_List = []
ln_p_List = []
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# エビデンスの対数
res = (t - Phi*m_N)
beta_ML = N / (res*res)
E_mN = beta_ML/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta_ML) - E_mN -...
ln_p_List += [ln_p_t]
beta_ML_List += [beta_ML]
list_plot(ln_p_List, ymin = -26, ymax = 0).show()
#print beta_ML_List
}}
&ref(sage0-1.png);
*** どうして図3.14と合わないのか？ [#q2979465]
どうして図3.14のような結果にならないのかと調べたのですが...
方がいらして、私と同じ傾向になったとの記述がありました。
[[基底関数を色々変えて、線形回帰のエビデンスを計算してみ...
私の推測では、ln |A|の計算をA.norm()と間違えたのではない...
sageへの入力:
#pre{{
# 当初、式3.85, 3.86のミスプリと考えたのですが、
# →http://d.hatena.ne.jp/n_shuyo/20090715/evidence#c
# n_shuyoさんのコメントで 正則行列の行列式は、逆行列の...
# 式3.85, 3.86 は正しい！
# 私の推測：図3.14の計算で ln |A|とするところをln A.norm(...
ln_p_List = []
for M in range(10):
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1...
Phi_t = Phi.transpose()
# m_N, Aを定義
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# エビデンスの対数
res = (t - Phi*m_N)
E_mN = beta/2* res * res + alpha/2 * m_N * m_N
ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/...
ln_p_List += [ln_p_t]
list_plot(ln_p_List)
}}
&ref(sage0-2.png);
*** 最適なα、βを求める [#tf19c14f]
α、βを固定にした場合に、エビデンスが最大はM=4だったので、...
最適なα、βを求めてみます。
エビデンスを最大にするαは、式(3.92)
$$
\alpha = \frac{\gamma}{m_N^T m_N}
$$
ここで\(\gamma\)は、式(3.87)
$$
(\beta \Phi^T \Phi) u_i = \lambda_i u_i
$$
の固有値から、式(3.91)
$$
\gamma = \sum_{i} \frac{\lambda_i}{\alpha + \lambda_i}
$$
で求まります。
βについても、式(3.95)
$$
\frac{1}{\beta} = \frac{1}{N - \gamma} \sum_{n=1}^{N} ( t...
$$
\(\gamma\)は、\(\alpha\)に依存し、\(m_N\)も\(\alpha\)に依...
できません。
+ 適当な \(\alpha, \beta\)を初期値とする
+ \(m_N, \gamma\)を求め
+ 式(3.95)から\(\beta\)を求める
+ 式(3.92)から\(\alpha\)を求め、2番目からを繰り返す
以下の例では、
- $\alpha = 5 \times 10^{-3}$
- $\beta = 11.1$
から20回の計算を繰り返し、\(\gamma, \alpha, \beta\)の値を...
sageへの入力:
#pre{{
# M=4の場合のα、βを決める
M=4
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# Φ^T Φは固定なので、先に固有値を計算
B = Phi_t*Phi
lambs = B.eigenvalues()
# 初期化
alpha = 5*10^(-3)
beta = 11.1
gamma = M
for i in range(20):
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
gamma = sum((lamb*beta) /(alpha + (lamb*beta)) for la...
alpha = gamma/(m_N*m_N)
res = (t - Phi*m_N)
beta = (N - gamma) / (res*res)
print gamma, alpha, beta
}}
sageの出力：
#pre{{
4.02295319874 0.0160029572678 15.5703330431
3.86957616939 0.0192233025727 14.6079344569
3.81804019954 0.0208658867531 14.0162355337
3.79079625631 0.021876053777 13.6574090928
3.77410819657 0.0225480023858 13.4248103487
3.76308615077 0.0230152846069 13.2666278994
3.75546956914 0.0233495381967 13.1554382026
3.75004882783 0.0235932327887 13.0754487343
3.74611222743 0.0237732924978 13.0169432556
3.74321234677 0.0239076177514 12.973632867
3.7410540217 0.0240085312871 12.9412853436
3.73943541844 0.0240847392383 12.9169657574
3.7382147331 0.0241425141272 12.8985911035
3.73729026269 0.0241864427411 12.8846563016
3.73658790674 0.0242199173426 12.8740587014
3.73605302024 0.0242454685376 12.8659818046
3.73564493065 0.0242649966702 12.8598160009
3.73533314922 0.0242799360403 12.8551032449
3.73509469602 0.0242913734086 12.8514976848
3.73491217789 0.0243001346625 12.848737196
}}
*** γと2αEw(mN)のグラフ [#qa302662]
M=4の多項式回帰に対して、図3.16(a)
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
のように\(\gammaと2\alpha E_w(m_N)\)のグラフを\(\ln \alph...
形が図3.16とは異なり、曲線の交差する点もα=0.024だとln αの...
交差しています。（何かちがう！）
sageへの入力:
#pre{{
# γをln(a)の関数として定義
def g(ln_a):
return abs(sum((lamb*beta) /(e^ln_a + (lamb*beta)) fo...
a_plt = plot(g, [x,-6, 6], color='red')
}}
sageへの入力:
#pre{{
# 2αEw(m_N)をln(a)の関数として定義
def aEw(ln_a):
a = (e^ln_a).n()
A = a*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# m_Nいくつかの点が複素数になるためabs(a*m_N*m_N)とした
return abs(a*m_N*m_N)
Ew_plt = plot(aEw, [x,-6, 6], color='blue')
(a_plt + Ew_plt).show()
}}
&ref(sage0-3.png);
*** 残差と対数エビデンスのグラフ [#z3abc719]
同様に、M=4の多項式回帰に対して、図3.16(b)
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
のように残差と対数エビデンスを\(\ln \alpha\)の関係を図化...
やはり形が、図3.16とは異なります。
sageへの入力:
#pre{{
b_phi_2 = beta*Phi_t*Phi
Phi_t_t = Phi_t* t
def ln_pt(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
E_mX = beta/2 * res*res + a/2 * m_N*m_N
# E_mXいくつかの点が複素数になるためabs(E_mX)とした
return (M/2 * ln_a + N/2 * ln(beta) - abs(E_mX) -1/2*...
#
def res2(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
return abs(res*res)
}}
sageへの入力:
#pre{{
ln_pt_plt = plot(ln_pt, [x,-6, 6], color='red')
res_2_plt = plot(res2, [x,-6, 6])
(ln_pt_plt + res_2_plt).show()
}}
&ref(sage0-4.png);
*** M=4に対する最適な解 [#bff6d673]
M=4の多項式回帰に対する最適なα、βを使って\(sin(2 \pi x)\)...
sageへの入力:
#pre{{
# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
data_plt = list_plot(zip(X, t))
}}
sageへの入力:
#pre{{
# 求まったα、βでベイズ的なフィッティングを再計算
alpha = 0.0243001346625
beta = 12.848737196
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t
f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1)))
# 分散
def s(x):
phi_x = vector([x^i for i in range(M+1)])
S = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi)....
s_sqr = 1/beta + phi_x * S * phi_x
return sqrt(s_sqr)
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red')
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-5.png);
M=9の最適解（図1.17）
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
と比べると若干フィッティングが良くありませんが、与えられ...
** どうしてエビデンス関数が合わないのか [#ic963cc3]
つぎにどうして、図3.16と合わないのか、調べてみます。
図3.16の説明をよく読むと、基底関数が多項式ではなく、ガウ...
*** 基底関数のチェック [#z831c5cd]
ガウス基底関数がどのようなものなのか、図3.1を再現しながら...
以下に、
- 多項式の基底関数
- ガウス基底関数
$$
\phi_j(x) = exp \left \{ - \frac{(x - \mu_j)^2}{2 s^2} \r...
$$
- シグモイド基底関数
$$
\phi_j(x) = \sigma \left( \frac{x - \mu_j}{s} \right)
$$
$$
\sigma(a) = \frac{1}{1 + exp(-a)}
$$
を表示します。
ガウス基底関数、シグモイド基底関数では、μの値として-1から...
とり、それを基底関数の数で等分した値をとるみたいです。
[[線形回帰モデルとか>http://d.hatena.ne.jp/smly/20100630/...
sageへの入力:
#pre{{
# 基底関数のチェック
from pylab import linspace
M = 11
mu = linspace(-1, 1, M)
s = 0.1
s_sq = (s)^2
# 多項式基底関数定義
def _phi_poly(x, j):
return x^j
# ガウス基底関数
def _phi_gauss(x, j):
return e^(-1*(x - mu[j])^2/(2* s_sq))
# ロジスティック基底関数
def _logi_sig(x):
return 1/(1 + e^(-1*x))
def _phi_sigmoid(x, j):
return _logi_sig((x - mu[j])/s)
}}
sageへの入力:
#pre{{
x = var('x')
poly_plt =Graphics()
for j in range(1, M+1):
f = lambda x : _phi_poly(x, j)
poly_plt += plot(f, [x, -1, 1])
poly_plt.show()
}}
&ref(sage0-6.png);
sageへの入力:
#pre{{
gauss_plt =Graphics()
for j in range(M):
f = lambda x : _phi_gauss(x, j)
gauss_plt += plot(f, [x, -1, 1])
gauss_plt.show()
}}
&ref(sage0-7.png);
sageへの入力:
#pre{{
sigmoid_plt =Graphics()
for j in range(M):
f = lambda x : _phi_sigmoid(x, j)
sigmoid_plt += plot(f, [x, -1, 1])
sigmoid_plt.show()
}}
&ref(sage0-8.png);
*** ガウス基底関数を使ったγと2αEw(mN)のグラフ [#p159acdf]
M=9個のガウス基底関数を使ったγと2αEw(mN)のグラフを以下に...
関数の形状と交差点の位置が図3.16と合います。
sageへの入力:
#pre{{
# ガウス基底関数で三角関数の例題を近似
# j=0に対応するため_phiを以下のように定義
def _phi(x, j):
if j == 0:
return 1
else:
return _phi_gauss(x, j-1)
# 定数をセット
M=9
mu = linspace(0, 1, M)
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] ...
Phi_t = Phi.transpose()
# Φ^T Φは固定なので、先に固有値を計算
B = Phi_t*Phi
lambs = B.eigenvalues()
# 初期化
alpha = 5*10^(-3)
beta = 11.1
gamma = M
for i in range(20):
A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
gamma = sum((lamb*beta) /(alpha + (lamb*beta)) for la...
alpha = gamma/(m_N*m_N)
#res = (t - Phi*m_N)
#beta = (N - gamma) / (res*res)
print gamma, alpha, beta
}}
#pre{{
8.97526712356 2.74314585233 11.1000000000000
5.9805253866 5.87614916325 11.1000000000000
5.02659378174 6.10302152117 11.1000000000000
4.97666250731 6.12118986054 11.1000000000000
4.97273606757 6.12265708879 11.1000000000000
4.97241943407 6.1227756571 11.1000000000000
4.97239384954 6.12278523925 11.1000000000000
4.97239178194 6.12278601364 11.1000000000000
4.97239161484 6.12278607622 11.1000000000000
4.97239160134 6.12278608128 11.1000000000000
4.97239160025 6.12278608169 11.1000000000000
4.97239160016 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
4.97239160015 6.12278608172 11.1000000000000
}}
sageへの入力:
#pre{{
a_plt = plot(g, [x,-5, 5])
# 2αEw(m_N)をln(a)の関数として定義
def aEw(ln_a):
a = (e^ln_a).n()
A = a*matrix((M+1),(M+1),1) + beta*Phi_t * Phi
m_N = beta* A.inverse() * Phi_t* t
# m_Nいくつかの点が複素数になるためabs(a*m_N*m_N)とした
return abs(a*m_N*m_N)
Ew_plt = plot(aEw, [x,-5, 5], color='gray')
(a_plt + Ew_plt).show()
}}
&ref(sage0-9.png);
*** ガウス基底関数を使った残差と対数エビデンスのグラフ [#...
同様に、ガウス基底関数を使った残差と対数エビデンスのグラ...
残差の形状がやや違いますが、残差の最低値の付近に対数エビ...
部分はほぼ合っています。
sageへの入力:
#pre{{
b_phi_2 = beta*Phi_t*Phi
Phi_t_t = Phi_t* t
def ln_pt(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
E_mX = beta/2 * res*res + a/2 * m_N*m_N
return ((M/2 * ln_a + N/2 * ln(beta) - E_mX -1/2*ln(A...
# 残差の定義
def res2(ln_a):
a = (e^ln_a).n()
A = a*matrix(RDF, (M+1),(M+1),1) + b_phi_2
m_N = beta* A.inverse() * Phi_t_t
res = (t - Phi*m_N)
return abs(res*res)
# 結果のプロット
ln_pt_plt = plot(ln_pt, [x,-5, 5], color='red')
res_2_plt = plot(res2, [x,-5, 5])
(ln_pt_plt + res_2_plt).show()
}}
&ref(sage0-10.png);
*** ガウス基底関数を回帰結果 [#rf54c993]
おまけに、ガウス基底関数の回帰結果と分散の分布、オリジナ...
プロットした図を示します。ガウス基底を使った回帰も良く合...
sageへの入力:
#pre{{
# ムーア＿ベンローズの疑似逆行列
Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi...
# 平均の重み
Wml = beta*Phi_dag * t
f = lambda x : (sum((Wml[i]*_phi(x, i)) for i in range(0,...
# 分散
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)....
s_sqr = 1/beta + phi_x * S * phi_x
return sqrt(s_sqr)
s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcol...
s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcol...
y_plt = plot(f, [x, 0, 1], rgbcolor='red')
(y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymi...
}}
&ref(sage0-11.png);
** コメント [#n7c3a48f]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
