sage/PRML - RVM（関連ベクトルマシン）
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/27からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/106/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第7章-RVM（関連ベクトルマシン）をSageで試す [#z65b795c]
最初に「パターン認識と機械学習」（PRML）を鈴木由宇さんか...
RVMが今の流行で、それをとてもスマートに実装している人がい...
se-kichiさんの
[[「きちめも」>http://d.hatena.ne.jp/se-kichi/20100519/12...
でした。
これを見て、Sageを使ってPRMLの例題を試してみようと思いま...
ようやく「きちめも」のRVMまでPRMLを読み進むことができまし...
ここで紹介するのは、PRMLの図7.9、
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
です。エビデンス近似の場合の最も良くフィットしたMが3個で...
RVMで求まった関連ベクトルが3個というのも意味深いものを感...
** テストデータ [#t117dc02]
1章のSin曲線のフィッティングで使用したデータを今回も使い...
[[sage/PRML-線形回帰]]
参照してください。
sageへの入力:
#pre{{
# 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)
# データのプロット
x = var('x');
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green');
data_plt = list_plot(zip(X, t)); data_plt
(data_plt + sin_plt).show()
}}
&ref(sage0.png);
** RVMの更新式 [#u95eea73]
重みベクトルに対する事後確率は、式(7.81)
$$
p(w|t, X, \alpha, \beta) = \mathcal{N}(w|m, \Sigma)
$$
平均と共分散は、式(7.82), (7.83)
$$
m = \beta \Sigma \Phi^T t
$$
$$
\Sigma = ( A + \beta \Phi^T \Phi)^{-1}
$$
で与えられると、αとβのエビデンス近似は、対数尤度の式(7.85...
$$
\ln p(t|X, \alpha, \beta) = \ln \mathcal{N}(t|0, C)
$$
$$
= - \frac{1}{2} \left \{ N \ln (2\pi) + \ln |C| t^TC^{-1}...
$$
$$
C = \beta^{-1} I + \Phi A^{-1}\Phi^T
$$
を最大化することで、式(7.87), (7.88)となります。
$$
\alpha_i^{new} = \frac{\gamma_i}{m_i^2}
$$
$$
(\beta^{new})^{-1} = \frac{|| t - \Phi m ||^2}{N - \Sigma...
$$
ここで\(\gamma_i\)は、式(7.89)で定義されます。
$$
\gamma_i = 1 - \alpha_i \Sigma_{ii}
$$
予測値の条件付き確率分布は、式(7.76)、
$$
p(t|x, w, \beta) = \mathcal{N}(t|y(x), \beta^{-1})
$$
平均値は、式(7.77)、
$$
y(x) = w^T \phi(x)
$$
で、予測値の平均は、重みベクトルwを事後平均mとしたものに...
$$
p(t|x, X, t, \alpha^*, \beta^*) = \mathcal{N}(t|m^T \phi(...
$$
\(\sigma^2(x)\)は、式(7.91)
$$
\sigma^2(x) = (\beta^*)^{-1} + \phi(x)^T \Sigma \phi(x)
$$
となります。
** RVMの計算 [#iecef450]
更新式(7.87), (7.88)を使って
[[「きちめも」>http://d.hatena.ne.jp/se-kichi/20100519/12...
を参考にSageに移植してみました。
sageへの入力:
#pre{{
# PRMLの7.2.1を「きちめも」
# http://d.hatena.ne.jp/se-kichi/20100519/1274275285
# からsageに移行
from scipy.linalg import norm as _norm
import scipy as sp
# ガウスカーネルの定義
def kernel(x, y):
return exp(-_norm(x-y)**2/0.05)
# 対角要素のみを持つマトリックスを返す
def diagonal(m):
n = Sigma.ncols()
return diagonal_matrix([Sigma[i][i] for i in range(n)])
}}
sageへの入力:
#pre{{
# 初期設定
tol = 0.01
alpha_max = 1e10
N = len(X)
alpha = vector(RDF, [1 for i in range(N+1)])
beta = 1.0
# 計画行列Φ
Phi = matrix(RDF, [[1]+[kernel(xi, xj) for xj in X] for x...
Phi_T = Phi.transpose()
}}
sageへの入力:
#pre{{
loop = 10000
diff = 1
for n in range(loop):
if diff < tol: break
A = diagonal_matrix(alpha)
Sigma = (A + beta*Phi_T*Phi).inverse()
m = beta*Sigma*Phi_T*t
gamma = vector(sp.ones(N+1)) - diagonal(Sigma)*alpha
alpha_new = vector(RDF, [gamma[i]/(m[i]*m[i]) for i i...
beta_new = (N - gamma.sum()) / (t - Phi*m).norm()^2
alpha_new = alpha_new.apply_map(lambda x : x if x < a...
d_alpha, d_beta = alpha_new - alpha, beta_new - beta
diff = d_alpha.norm()^2 + d_beta^2
alpha, beta = alpha_new, beta_new
print alpha
}}
sageからの出力:
#pre{{
(1.00000000000000e10, 1.00000000000000e10, 1.000000000000...
1.00000000000000e10, 1.00000000000000e10, 1.0000000000000...
1.95795402710279, 1.00000000000000e10, 10.9873539203255)
}}
sageへの入力:
#pre{{
# 関連ベクトルのインデックスを取得
rv_index = [i-1 for i in range(1, N+1) if alpha[i] < alph...
# 関連ベクトル以外の平均を０とする
_m = vector([m[i] if alpha[i] < alpha_max else 0 for i in...
print rv_index
}}
sageからの出力:
#pre{{
[2, 7, 9]
}}
sageへの入力:
#pre{{
# 関連ベクトルのプロット
rv_plt = Graphics()
for i in rv_index:
rv_plt += point((X[i], t[i]), rgbcolor='red', pointsi...
}}
sageへの入力:
#pre{{
# 予測値の平均
def mean(x):
phi = vector(RDF, [1.0] + [kernel(x, xj) for xj in X])
return _m * phi
# 予測値の分散
def variance(x):
phi = vector(RDF, [1.0] + [kernel(x, xj) for xj in X])
return 1/beta + phi*Sigma*phi
}}
** 結果のプロット [#y0fd55f7]
オリジナルのsin曲線を緑で、データ点を青で、計算したフィッ...
関連ベクトルを大きめの赤い点でプロットしました。
sageへの入力:
#pre{{
# 結果のプロット
var('x')
y_plt = plot(lambda x : mean(x), [x, 0, 1], rgbcolor='re...
s_u_plt = plot(lambda x : mean(x) + sqrt(variance(x)), [x...
s_d_plt = plot(lambda x : mean(x) - sqrt(variance(x)), [x...
(y_plt + s_u_plt + s_d_plt + data_plt + sin_plt + rv_plt)...
}}
&ref(sage0-1.png);
** コメント [#e67e2cf7]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/27からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/106/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第7章-RVM（関連ベクトルマシン）をSageで試す [#z65b795c]
最初に「パターン認識と機械学習」（PRML）を鈴木由宇さんか...
RVMが今の流行で、それをとてもスマートに実装している人がい...
se-kichiさんの
[[「きちめも」>http://d.hatena.ne.jp/se-kichi/20100519/12...
でした。
これを見て、Sageを使ってPRMLの例題を試してみようと思いま...
ようやく「きちめも」のRVMまでPRMLを読み進むことができまし...
ここで紹介するのは、PRMLの図7.9、
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
です。エビデンス近似の場合の最も良くフィットしたMが3個で...
RVMで求まった関連ベクトルが3個というのも意味深いものを感...
** テストデータ [#t117dc02]
1章のSin曲線のフィッティングで使用したデータを今回も使い...
[[sage/PRML-線形回帰]]
参照してください。
sageへの入力:
#pre{{
# 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)
# データのプロット
x = var('x');
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green');
data_plt = list_plot(zip(X, t)); data_plt
(data_plt + sin_plt).show()
}}
&ref(sage0.png);
** RVMの更新式 [#u95eea73]
重みベクトルに対する事後確率は、式(7.81)
$$
p(w|t, X, \alpha, \beta) = \mathcal{N}(w|m, \Sigma)
$$
平均と共分散は、式(7.82), (7.83)
$$
m = \beta \Sigma \Phi^T t
$$
$$
\Sigma = ( A + \beta \Phi^T \Phi)^{-1}
$$
で与えられると、αとβのエビデンス近似は、対数尤度の式(7.85...
$$
\ln p(t|X, \alpha, \beta) = \ln \mathcal{N}(t|0, C)
$$
$$
= - \frac{1}{2} \left \{ N \ln (2\pi) + \ln |C| t^TC^{-1}...
$$
$$
C = \beta^{-1} I + \Phi A^{-1}\Phi^T
$$
を最大化することで、式(7.87), (7.88)となります。
$$
\alpha_i^{new} = \frac{\gamma_i}{m_i^2}
$$
$$
(\beta^{new})^{-1} = \frac{|| t - \Phi m ||^2}{N - \Sigma...
$$
ここで\(\gamma_i\)は、式(7.89)で定義されます。
$$
\gamma_i = 1 - \alpha_i \Sigma_{ii}
$$
予測値の条件付き確率分布は、式(7.76)、
$$
p(t|x, w, \beta) = \mathcal{N}(t|y(x), \beta^{-1})
$$
平均値は、式(7.77)、
$$
y(x) = w^T \phi(x)
$$
で、予測値の平均は、重みベクトルwを事後平均mとしたものに...
$$
p(t|x, X, t, \alpha^*, \beta^*) = \mathcal{N}(t|m^T \phi(...
$$
\(\sigma^2(x)\)は、式(7.91)
$$
\sigma^2(x) = (\beta^*)^{-1} + \phi(x)^T \Sigma \phi(x)
$$
となります。
** RVMの計算 [#iecef450]
更新式(7.87), (7.88)を使って
[[「きちめも」>http://d.hatena.ne.jp/se-kichi/20100519/12...
を参考にSageに移植してみました。
sageへの入力:
#pre{{
# PRMLの7.2.1を「きちめも」
# http://d.hatena.ne.jp/se-kichi/20100519/1274275285
# からsageに移行
from scipy.linalg import norm as _norm
import scipy as sp
# ガウスカーネルの定義
def kernel(x, y):
return exp(-_norm(x-y)**2/0.05)
# 対角要素のみを持つマトリックスを返す
def diagonal(m):
n = Sigma.ncols()
return diagonal_matrix([Sigma[i][i] for i in range(n)])
}}
sageへの入力:
#pre{{
# 初期設定
tol = 0.01
alpha_max = 1e10
N = len(X)
alpha = vector(RDF, [1 for i in range(N+1)])
beta = 1.0
# 計画行列Φ
Phi = matrix(RDF, [[1]+[kernel(xi, xj) for xj in X] for x...
Phi_T = Phi.transpose()
}}
sageへの入力:
#pre{{
loop = 10000
diff = 1
for n in range(loop):
if diff < tol: break
A = diagonal_matrix(alpha)
Sigma = (A + beta*Phi_T*Phi).inverse()
m = beta*Sigma*Phi_T*t
gamma = vector(sp.ones(N+1)) - diagonal(Sigma)*alpha
alpha_new = vector(RDF, [gamma[i]/(m[i]*m[i]) for i i...
beta_new = (N - gamma.sum()) / (t - Phi*m).norm()^2
alpha_new = alpha_new.apply_map(lambda x : x if x < a...
d_alpha, d_beta = alpha_new - alpha, beta_new - beta
diff = d_alpha.norm()^2 + d_beta^2
alpha, beta = alpha_new, beta_new
print alpha
}}
sageからの出力:
#pre{{
(1.00000000000000e10, 1.00000000000000e10, 1.000000000000...
1.00000000000000e10, 1.00000000000000e10, 1.0000000000000...
1.95795402710279, 1.00000000000000e10, 10.9873539203255)
}}
sageへの入力:
#pre{{
# 関連ベクトルのインデックスを取得
rv_index = [i-1 for i in range(1, N+1) if alpha[i] < alph...
# 関連ベクトル以外の平均を０とする
_m = vector([m[i] if alpha[i] < alpha_max else 0 for i in...
print rv_index
}}
sageからの出力:
#pre{{
[2, 7, 9]
}}
sageへの入力:
#pre{{
# 関連ベクトルのプロット
rv_plt = Graphics()
for i in rv_index:
rv_plt += point((X[i], t[i]), rgbcolor='red', pointsi...
}}
sageへの入力:
#pre{{
# 予測値の平均
def mean(x):
phi = vector(RDF, [1.0] + [kernel(x, xj) for xj in X])
return _m * phi
# 予測値の分散
def variance(x):
phi = vector(RDF, [1.0] + [kernel(x, xj) for xj in X])
return 1/beta + phi*Sigma*phi
}}
** 結果のプロット [#y0fd55f7]
オリジナルのsin曲線を緑で、データ点を青で、計算したフィッ...
関連ベクトルを大きめの赤い点でプロットしました。
sageへの入力:
#pre{{
# 結果のプロット
var('x')
y_plt = plot(lambda x : mean(x), [x, 0, 1], rgbcolor='re...
s_u_plt = plot(lambda x : mean(x) + sqrt(variance(x)), [x...
s_d_plt = plot(lambda x : mean(x) - sqrt(variance(x)), [x...
(y_plt + s_u_plt + s_d_plt + data_plt + sin_plt + rv_plt)...
}}
&ref(sage0-1.png);
** コメント [#e67e2cf7]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
