FrontPage
2011/06/27からのアクセス回数
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
また、Sageのサーバを公開しているサイト( 、 )にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
最初に「パターン認識と機械学習」(PRML)を鈴木由宇さんから紹介して頂いた時に、 RVMが今の流行で、それをとてもスマートに実装している人がいると教えて頂いたのが、 se-kichiさんの http://d.hatena.ne.jp/se-kichi/20100519/1274275285 でした。
これを見て、Sageを使ってPRMLの例題を試してみようと思いました。 ようやく「きちめも」のRVMまでPRMLを読み進むことができました。
ここで紹介するのは、PRMLの図7.9、
です。エビデンス近似の場合の最も良くフィットしたMが3個であり、 RVMで求まった関連ベクトルが3個というのも意味深いものを感じます。
1章のSin曲線のフィッティングで使用したデータを今回も使います。 sage/PRML-線形回帰 参照してください。
sageへの入力:
# PRMLのsin曲線のデータ
data = matrix([
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() }}
重みベクトルに対する事後確率は、式(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), (7.86) $$ \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}t \right \}. $$ $$ 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_i \gamma_i} $$ ここで\(\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としたものになります。式(7.90) $$ p(t|x, X, t, \alpha^*, \beta^*) = \mathcal{N}(t|m^T \phi(x), \sigma^2(x)) $$ \(\sigma^2(x)\)は、式(7.91) $$ \sigma^2(x) = (\beta^*)^{-1} + \phi(x)^T \Sigma \phi(x) $$ となります。
更新式(7.87), (7.88)を使って http://d.hatena.ne.jp/se-kichi/20100519/1274275285 を参考にSageに移植してみました。
sageへの入力:
# 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):
# 対角要素のみを持つマトリックスを返す
def diagonal(m):
}}
sageへの入力:
# 初期設定
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 xi in X]) Phi_T = Phi.transpose() }}
sageへの入力:
loop = 10000 diff = 1 for n in range(loop):
print alpha }} sageからの出力:
(1.00000000000000e10, 1.00000000000000e10, 1.00000000000000e10, 0.899348917265024, 1.00000000000000e10, 1.00000000000000e10, 1.00000000000000e10, 1.00000000000000e10, 1.95795402710279, 1.00000000000000e10, 10.9873539203255) }}
sageへの入力:
# 関連ベクトルのインデックスを取得
rv_index = [i-1 for i in range(1, N+1) if alpha[i] < alpha_max]
# 関連ベクトル以外の平均を0とする
_m = vector([m[i] if alpha[i] < alpha_max else 0 for i in range(N+1)]) print rv_index }} sageからの出力:
[2, 7, 9] }}
sageへの入力:
# 関連ベクトルのプロット
rv_plt = Graphics() for i in rv_index:
}}
sageへの入力:
# 予測値の平均
def mean(x):
# 予測値の分散
def variance(x):
}}
オリジナルのsin曲線を緑で、データ点を青で、計算したフィッティング曲線を緑で、予測値の標準偏差をグレイで、 関連ベクトルを大きめの赤い点でプロットしました。
sageへの入力:
# 結果のプロット
var('x') y_plt = plot(lambda x : mean(x), [x, 0, 1], rgbcolor='red') s_u_plt = plot(lambda x : mean(x) + sqrt(variance(x)), [x, 0, 1], rgbcolor='grey'); s_d_plt = plot(lambda x : mean(x) - sqrt(variance(x)), [x, 0, 1], rgbcolor='grey'); (y_plt + s_u_plt + s_d_plt + data_plt + sin_plt + rv_plt).show() }}
皆様のご意見、ご希望をお待ちしております。