FrontPage

2011/06/05からのアクセス回数 3800

ここで紹介した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で試す

sage/PRML-逐次ベイズ学習では、 重みwの分布について解いていましたが、 実際には新しいxの値たいするtを予測することが目的となります。 ここではPRML3.3.2の予測分布、図3.8

Figure3.8a.jpg Figure3.8b.jpg

Figure3.8c.jpg Figure3.8d.jpg

をSageを使って試してみます。

準備

sage/PRML-線形回帰 で使ったデータと同じものを座標Xと目的値tにセットし、関数Φを定義します。

sageへの入力:

# 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への入力:

# データのプロット
x = var('x')
sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')
data_plt = list_plot(zip(X, t))

ガウス基底関数の定義

sage/PRML- エビデンス近似 のようにガウス基底関数を定義しますが、今回は0から1の範囲で定義します。

近似に使う場合には、_phiのようにj=0が1となるような項を追加します。

sageへの入力:

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への入力:

# ガウス基底関数で三角関数の例題を近似
# j=0に対応するため_phiを以下のように定義
def _phi(x, j):
    if j == 0:
        return 1
    else:
        return _phi_gauss(x, j-1)
# 初期化
alpha = 2
beta = 25

予測分布

予測分布は、式(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への入力:

# 予測分布をプロットする
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)

予測分布のプロット

以下にM=1, 2, 4, 10の予測分布を示します。パラメータは、

  • \(\alpha = 2\)
  • \(\beta = 25\)

としました。

sageへの入力:

# 1点の場合
X_ = vector([X[4]])
t_ = vector([t[4]])
_plot_predict(X_, t_)

sage0.png

sageへの入力:

# 2点の場合
X_ = vector([X[4], X[6]])
t_ = vector([t[4], t[6]])
_plot_predict(X_, t_)

sage0-1.png

sageへの入力:

# 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_)

sage0-2.png

sageへの入力:

# 10点の場合
_plot_predict(X, t)

sage0-3.png

コメント

選択肢 投票
おもしろかった 3  
そうでもない 0  
わかりずらい 1  

皆様のご意見、ご希望をお待ちしております。

  • 自分としては、PRMLの図3.8をある程度再現できた思っていたのですが、分かりづらいとのコメントはちょっときつかったです。 -- 竹本 浩? 2013-06-25 (火) 20:49:17

(Input image string)


添付ファイル: filesage0-3.png 599件 [詳細] filesage0-2.png 615件 [詳細] filesage0-1.png 577件 [詳細] filesage0.png 581件 [詳細]

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2017-06-11 (日) 17:13:15 (491d)
SmartDoc