sage/PRML- エビデンス近似のバックアップ差分(No.2)

バックアップ一覧
現在との差分を表示
ソースを表示
バックアップを表示
sage/PRML- エビデンス近似へ行く。
- 1 (2011-06-04 (土) 18:25:29)
- 2 (2011-06-04 (土) 19:27:59)
追加された行はこの色です。
削除された行はこの色です。
[[FrontPage]]

#contents

2011/06/04からのアクセス回数 &counter;

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://sage.math.canterbury.ac.nz/home/pub/95/

また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。

* 第3章-エビデンス近似をSageで試す [#x3862501]

[[第1章-Sageを使って線形回帰を試してみる>http://sage.math.canterbury.ac.nz/home/pub/86/]]では、
* 第3章-エビデンス近似をSageで試す [#vecd04fc]

[[sage/PRML-線形回帰]]では、
与えられたモデルパラメータM, α、βに対して解いていましたが、
実際には自分でベストなモデルパラメータM, α、βを求める必要があります。

[[パターン認識と機械学習>http://www.amazon.co.jp/dp/443110013X/]]
の３章ではエビデンス関数を使って最適なパラメータを求める方法が説明されています。

ここでは、Sageを使ってエビデンス関数の評価、最適なモデルパラメータの推定を試してみます。


** エビデンス関数の評価 [#k3992ae4]
** エビデンス関数の評価 [#xf7b5eab]

エビデンス関数の対数表現は、式(3.86)
$$
\ln p(t|\alpha, \beta) =  \frac{M}{2} \ln \alpha + \frac{N}{2} \ln \beta -E(m_N) - \frac{1}{2} \ln | A | + \frac{N}{2} \ln (2 \pi)
$$
で与えられます。

ここで、$E(m_N)$は、式(3.82)
$$
E(m_N) = \frac{\beta}{2}|| t = \Phi m_N ||^2 + \frac{\alpha}{2} m_N^T m_N
$$
$m_N$は、式(3.84)
$$
m_N = \beta A^{-1} \Phi^T t
$$
Aは、式(3.81)
$$
A = \alpha I + \beta \Phi^T \Phi
$$

目指すは、図3.14です。		
<p align="center">
<img src="http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.14.jpg" width="50%"></img>
目指すは、図3.14です。

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.14.jpg,center,50%);

*** 準備 [#s196841f]
*** 準備 [#ac2b67a1]

<a href="http://sage.math.canterbury.ac.nz/home/pub/86/">第1章-Sageを使って線形回帰を試してみる</a>
[[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 X.list()]);
}}
#pre{{
}}


sageへの入力:
#pre{{
# Φ関数定義
def _phi(x, j):
    return x^j
}}
#pre{{
}}

*** エビデンスの計算 [#p34968e0]

*** エビデンスの計算 [#zbfb813f]

まずは、α、βを固定値

- $\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))] for x in X.list()])
    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/2 * ln(A.det()) - (N/2 * ln(2*pi))
    ln_p_List += [ln_p_t]
list_plot(ln_p_List, ymin = -26, ymax = -5).show()
# 結果は、M=3以降でfig.3.14のような急激な降下は見られない
}}
#pre{{
}}
&ref(sage0.png);

*** β_MLを使ってエビデンスを計算 [#rf539535]

*** β_MLを使ってエビデンスを計算 [#h8575640]

そこで、βを平均値の分散$\beta_{ML}$を使って計算してみました。
$$
\frac{1}{\beta_{ML}} = \frac{1}{N} \sum (y(x_n, W_ML) - t_n)^2 
$$

期待に反し、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))] for x in X.list()]) 
    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 - 1/2 * ln(A.det()) - N/2 * ln(2*pi)
    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
}}
#pre{{
}}

&ref(sage0-1.png);

*** どうして図3.14と合わないのか？ [#w40aef53]

*** どうして図3.14と合わないのか？ [#q2979465]
どうして図3.14のような結果にならないのかと調べたのですが、他にも同様の計算をした
方がいらして、私と同じ傾向になったとの記述がありました。
<a href="http://d.hatena.ne.jp/n_shuyo/20090715/evidence">基底関数を色々変えて、線形回帰のエビデンスを計算してみた</a>
[[基底関数を色々変えて、線形回帰のエビデンスを計算してみた>http://d.hatena.ne.jp/n_shuyo/20090715/evidence]]

私の推測では、ln |A|の計算をA.norm()と間違えたのではないかと思われます。以下にA.norm()としたときの図を示します。


sageへの入力:
#pre{{
# 当初、式3.85, 3.86のミスプリと考えたのですが、
# →http://d.hatena.ne.jp/n_shuyo/20090715/evidence#c
#    n_shuyoさんのコメントで 正則行列の行列式は、逆行列の行列式の逆数：|A| = 1/|A^-1|
#    式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))] for x in X.list()]) 
    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/2 * ln(A.norm()) - N/2 * ln(2*pi)
    ln_p_List += [ln_p_t]
list_plot(ln_p_List)
}}
#pre{{
}}

&ref(sage0-2.png);

*** 最適なα、βを求める [#na288652]

*** 最適なα、βを求める [#tf19c14f]

α、βを固定にした場合に、エビデンスが最大はM=4だったので、M=4に対する
最適なα、βを求めてみます。

エビデンスを最大にするαは、式(3.92)
$$
\alpha = \frac{\gamma}{m_N^T m_N}
$$
ここで$\gamma$は、式(3.87)
ここで\(\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_n - m_N^T \phi(x_n))^2
$$

$\gamma$は、$\alpha$に依存し、$m_N$も$\alpha$に依存するため、すぐに最適な値を得ることが
\(\gamma\)は、\(\alpha\)に依存し、\(m_N\)も\(\alpha\)に依存するため、すぐに最適な値を得ることが
できません。


+ 適当な $\alpha, \beta$を初期値とする
+ 適当な \(\alpha, \beta\)を初期値とする
+ \(m_N, \gamma\)を求め
+ 式(3.95)から\(\beta\)を求める
+ 式(3.92)から\(\alpha\)を求め、2番目からを繰り返す

+ $m_N$, $\gamma$を求め

+ 式(3.95)から$\beta$を求める

+ 式(3.92)から$\alpha$を求め、2番目からを繰り返す

以下の例では、

- $\alpha = 5 \times 10^{-3}$

- $\beta = 11.1$
から20回の計算を繰り返し、$\gamma, \alpha, \beta$の値をプリントアウトしてみました。
から20回の計算を繰り返し、\(\gamma, \alpha, \beta\)の値をプリントアウトしてみました。



sageへの入力:
#pre{{
# M=4の場合のα、βを決める
M=4
Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()]) 
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 lamb in lambs)
    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]


*** γと2αEw(mN)のグラフ [#e961b2ff]

M=4の多項式回帰に対して、図3.16(a)
<p align="center">
<img src="http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.16a.jpg" width="50%"></img>

のように$\gammaと2\alpha E_w(m_N)$のグラフを$\ln \alpha$の関係を図化してみます。
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.16a.jpg,center,50%);

のように\(\gammaと2\alpha E_w(m_N)\)のグラフを\(\ln \alpha\)の関係を図化してみます。

形が図3.16とは異なり、曲線の交差する点もα=0.024だとln αの値は負になるのに、図3.16は正の部分で
交差しています。（何かちがう！）


sageへの入力:
#pre{{
# γをln(a)の関数として定義
def g(ln_a):
    return abs(sum((lamb*beta) /(e^ln_a + (lamb*beta)) for lamb in lambs))
a_plt = plot(g, [x,-6, 6], color='red')
}}
#pre{{
}}


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()
}}
#pre{{
}}

&ref(sage0-3.png);

*** 残差と対数エビデンスのグラフ [#n0d1b0bd]

*** 残差と対数エビデンスのグラフ [#z3abc719]

同様に、M=4の多項式回帰に対して、図3.16(b)
<p align="center">
<img src="http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.16b.jpg" width="50%"></img>

のように残差と対数エビデンスを$\ln \alpha$の関係を図化してみます。
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure3.16b.jpg,center,50%);

のように残差と対数エビデンスを\(\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*ln(A.det()) - N/2*ln(2*pi)).n()
#
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)
}}
#pre{{
}}


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()
}}
#pre{{
}}

&ref(sage0-4.png);

*** M=4に対する最適な解 [#efabd56e]

M=4の多項式回帰に対する最適なα、βを使って$sin(2 \pi x)$、観測データ、フィッティング曲線、分散をグラフに表示してみます。
*** 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))
}}
#pre{{
}}


sageへの入力:
#pre{{
# 求まったα、βでベイズ的なフィッティングを再計算
alpha = 0.0243001346625 
beta = 12.848737196

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]*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).inverse()
    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],  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.25, ymax=1.25)
}}
#pre{{
}}
&ref(sage0-5.png);


M=9の最適解（図1.17）
<p align="center">
<img src="http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure1.17.jpg" width="50%"></img>

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure1.17.jpg,center,50%);

と比べると若干フィッティングが良くありませんが、与えられた点をスムーズに補完しているのがわかります。


** どうしてエビデンス関数が合わないのか [#g0cdb7ad]
** どうしてエビデンス関数が合わないのか [#ic963cc3]

つぎにどうして、図3.16と合わないのか、調べてみます。

図3.16の説明をよく読むと、基底関数が多項式ではなく、ガウス基底関数になっていました。

*** 基底関数のチェック [#z831c5cd]

*** 基底関数のチェック [#waf261c2]

ガウス基底関数がどのようなものなのか、図3.1を再現しながら、見てみましょう。

以下に、

- 多項式の基底関数

- ガウス基底関数
$$
\phi_j(x) = exp \left \{ - \frac{(x - \mu_j)^2}{2 s^2} \right \}
$$

- シグモイド基底関数
$$
\phi_j(x) = \sigma \left( \frac{x - \mu_j}{s} \right)
$$
$$
\sigma(a) = \frac{1}{1 + exp(-a)}
$$
を表示します。

ガウス基底関数、シグモイド基底関数では、μの値として-1から1（問題によって0から1の場合もある）に
とり、それを基底関数の数で等分した値をとるみたいです。
<a href="http://d.hatena.ne.jp/smly/20100630/1277904761">線形回帰モデルとか</a>を参考にしました。
[[線形回帰モデルとか>http://d.hatena.ne.jp/smly/20100630/1277904761]]を参考にしました。


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)
}}
#pre{{
}}


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()
}}
#pre{{
}}
&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()
}}
#pre{{
}}
&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()
}}
#pre{{
}}
&ref(sage0-8.png);


*** ガウス基底関数を使ったγと2αEw(mN)のグラフ [#s4fcb9f7]
*** ガウス基底関数を使ったγと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))] for x in X.list()]) 
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 lamb in lambs)
    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()
}}
#pre{{
}}

&ref(sage0-9.png);

*** ガウス基底関数を使った残差と対数エビデンスのグラフ [#x1175a39]

*** ガウス基底関数を使った残差と対数エビデンスのグラフ [#gdf2d0a8]

同様に、ガウス基底関数を使った残差と対数エビデンスのグラフを示します。

残差の形状がやや違いますが、残差の最低値の付近に対数エビデンスの最大値が位置する
部分はほぼ合っています。


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.det()) - N/2*ln(2*pi))).n()
# 残差の定義
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()
}}
#pre{{
}}
&ref(sage0-10.png);


*** ガウス基底関数を回帰結果 [#v73db688]
*** ガウス基底関数を回帰結果 [#rf54c993]

おまけに、ガウス基底関数の回帰結果と分散の分布、オリジナルのsin曲線、データ点を
プロットした図を示します。ガウス基底を使った回帰も良く合っていることが分かります。


sageへの入力:
#pre{{
# ムーア＿ベンローズの疑似逆行列
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)
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.25, ymax=1.25)
}}
#pre{{
}}
&ref(sage0-11.png);



** コメント [#n7c3a48f]
#vote(おもしろかった,そうでもない,わかりずらい)

皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
sage/PRML- エビデンス近似 のバックアップ差分(No.2)

sage/PRML- エビデンス近似のバックアップ差分(No.2)
