sage/PRML - 混合密度ネットワーク の変更点

http://www.pwv.co.jp/~take/TakeWiki/index.php?sage%2FPRML%20-%20%E6%B7%B7%E5%90%88%E5%AF%86%E5%BA%A6%E3%83%8D%E3%83%83%E3%83%88%E3%83%AF%E3%83%BC%E3%82%AF

[ トップ ]   [ 編集 | 差分 | バックアップ | 添付 | リロード ]   [ 新規 | 一覧 | 単語検索 | 最終更新 | ヘルプ ]

• 追加された行はこの色です。
• 削除された行はこの色です。
• sage/PRML - 混合密度ネットワーク へ行く。
• sage/PRML - 混合密度ネットワーク の差分を削除

---

#freeze
[[FrontPage]]

#contents

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

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

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

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



* 第5章-混合密度ネットワークをSageで試す [#d7fb5911]

PRMLの第5章-混合密度ネットワークの例題、図5.19

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.19a.jpg,center,40%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.19b.jpg,center,40%);

のデータは、XとYを入れ替えただけで多峰性によりフィッティングが失敗する例です。
（x=0.5のところで、yの値として3つの選択しがあります）


** テストデータ [#z43bbaa8]

区間（0, 1）に一様分布する変数$x_n$をサンプリングし、目的値\(t_n\)
を\(x_n + 0.3 sin(2\pi x_n)\)に区間(-0.1, 0.1)の一様乱数を
加えたデータを200点生成します。


sageへの入力:
#pre{{
# PRML fig.5.19のデータ
#f = lambda x : x + 0.3*sin(2*pi*x) + random()*0.2 - 0.1
#X = [random() for i in range(200)]
#t = [f(x).n() for x in X]
}}


sageへの入力:
#pre{{
# データの保存
# save(X, DATA+'X')
# save(t, DATA+'t')
# リストアするには、以下のコメントを外す
X = load(DATA+'X')
t = load(DATA+'t')
print 'Done'
}}
sageからの出力:
#pre{{
Done
}}


*** グラフのプロット [#bb7ade87]

生成したデータとx,yを入れ替えたグラフをプロットします。

sageへの入力:
#pre{{
# データ生成
data_plt = list_plot(zip(X,t))
data_plt.show()
}}

&ref(sage0.png);

sageへの入力:
#pre{{
# X, tを入れ替えた図
data2_plt = list_plot(zip(t,X))
data2_plt.show()
}}

&ref(sage0-1.png);


** ニューラルネットワークでの回帰 [#wc37435e]

生成したデータをニューラルネットワークを使って、フィッティングさせたグラフを
以下に示します。（入力層1個、隠れ層6個、出力層1個で計算）

図5.19のbの様にx,tを入れ替えた場合には上手くフィッティングしません。

sageへの入力:
#pre{{
# 入力1点、隠れ層3個、出力1点
N_in = 1 
N_h = 6     
N_out = 1
# 隠れ層の活性化関数
h_1(x) = tanh(x)
# 出力層の活性化関数
h_2(x) = x
}}


sageへの入力:
#pre{{
# ニューラルネットワーク用のクラスのロード
attach "NeuralNetwork.sage"
# Scaled Conjugate Gradients（スケール共役勾配法）のロード
attach "SCG.sage"
}}


sageへの入力:
#pre{{
# 出力関数
def _f(x, net):
    y = net.feedForward(x)
    return y[0]
}}


sageへの入力:
#pre{{
# 定数設定
LEARN_COUNT = 500
RESET_COUNT = 50
beta = 0.5
# ニューラルネットワークの生成
net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# wの初期値を保存
saved_w = net.getW()
# 逐次勾配降下法による学習
_learn_csg(net, X, t, LEARN_COUNT, RESET_COUNT)
}}
sageからの出力:
#pre{{
WARNING: Output truncated!  
2 0.961554425907 3.77799819103
3 1.04117097273 2.41177356492
途中省略
499 1.14719236412 0.281449667896
500 0.626481765163 0.281421869851
}}


sageへの入力:
#pre{{
var('x')
g = lambda x : _f([x], net)
g_plt = plot(g, [x, 0, 1], rgbcolor='red')
(g_plt + data_plt).show()
}}

&ref(sage0-2.png);


sageへの入力:
#pre{{
# X, tを入れ替えた場合
# ニューラルネットワークの生成
net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# wの初期値を保存
saved_w = net.getW()
# 逐次勾配降下法による学習
_learn_csg(net, t, X, LEARN_COUNT, RESET_COUNT)
}}
sageからの出力:
#pre{{
WARNING: Output truncated!  
2 1.00671085821 5.493956126
3 0.983878474346 5.09315153282
途中省略
487 -19.5656631037 4.5870585875
488 -117.393981949 4.5870585875
489 -156.525310456 4.5870585875
}}


sageへの入力:
#pre{{
g2_plt = plot(g, [x, 0, 1], rgbcolor='red')
(g2_plt + data2_plt).show()
}}

&ref(sage0-3.png);


** 混合密度ネットワーク [#t943f5eb]

混合モデルの条件付き密度関数\(p(t|x)\)を式(5.148)
$$
p(t|x) = \sum_{k=1}^{K} \pi_k(x) \mathcal{N}(t|\mu_k(x), \sigma_k^2(x) I)
$$
で表し、モデルパラメータ

- 混合係数 \(\pi_k(x)\)

- 平均 \(\mu_k(x)\)

- 分散 \(\sigma^2_k(x)\)
をニューラルネットワークの出力とします。

混合計数は、式(5.149)
$$
\sum_{k=1}^{K} \pi_k(x) = 1,　0 \le \pi_k(x) \le 1
$$
の制約を満たすために、ソフトマックス関数、式(5.150)
$$
\pi_k(x) = \frac{exp(a_k^{\pi})}{\sum_{k=1}^{K}}
$$
とし、分散$\sigma_k^2(x) \ge 0$を満たすように、式(5.151)
$$
\sigma_k(x) = exp(a_k^{\sigma})
$$
平均は、ニューラルネットワークの出力をそのまま使って、式(5.152)
$$
\mu_{kj}(x) = a_{kj}^{\mu}
$$
とします。

誤差関数を式(5.153)
$$
E(w) = - \sum_{n=1}^{N} ln \left \{ \sum_{k=1}^K \pi_k(x_n, w) \mathcal{N}(t_n|\mu_k(x_n, w),  \sigma_k^2(x_n, w) I) \right \}
$$
とすると、事後分布式(5.154)
$$
\gamma_{nk}(t_n | x_n) = \frac{\pi_k \mathcal{N}_{nk}}{\sum_{l=1}^{K} \pi_l \mathcal{N}_{nl}}
$$
を使って各係数の微分を求めます。

混合係数の微分は、式(5.155)
$$
\frac{\partial E_n}{\partial a_k^{\pi}} = \pi_k - \gamma_{nk}
$$
平均の微分は、式(5.156)
$$
\frac{\partial E_n}{\partial a_{kl}^{\mu}} = \gamma_{nk} \left ( \frac{\mu_{kl} - t_{nl}}{\sigma_k^2} \right )
$$
分散の微分は、式(5.157)
$$
\frac{\partial E_n}{\partial a_{k}^{\sigma}} = \gamma_{nk} \left ( L - \frac{||t_{n} - \mu_k||^2}{\sigma_k^2} \right )
$$			

重みwへの部分は、
$$
\frac{\partial E_n}{\partial w_{ji}^{1}} = \delta_j x_i
$$
ここで、\(\delta_j\)は、
$$
\delta_j = (1 - z_j^2) \sum_{k=1}^2 \left ( w_k^{pi} \frac{\partial E_n}{\partial a_k^{\pi}} + w_k^{\mu} \frac{\partial E_n}{\partial w_{k}^{\sigma}} + \sum_{l=1}^{L} w_{kl}^{\mu} \frac{\partial E_n}{\partial w_{kl}^{\mu}} \right )
$$
混合係数の重み\(w_{kj}^{\pi}\)微分は、
$$
\frac{\partial E_n}{\partial w_{kj}^{\pi}} = \frac{\partial E_n}{\partial a_k^{\pi}} z_j
$$
平均の微分の重み\(w_{kj}^{\pi}\)微分は、
$$
\frac{\partial E_n}{\partial w_{kj}^{\mu}} = \sum_{l=1}^{L} \frac{\partial E_n}{\partial w_{kl}^{\mu}} z_j
$$
分散の微分の重み\(w_{kj}^{\sigma}\)微分は、
$$
\frac{\partial E_n}{\partial w_{kj}^{\sigma}} = \frac{\partial E_n}{\partial w_{k}^{\sigma}} z_j
$$		


*** MixtureDensityNetworkクラスの定義 [#c47c1e89]

MixtureDensityNetworkクラスを以下のように定義します。


sageへの入力:
#pre{{
class MixtureDensityNetwork:
    """ 混合密度ネットワーク用クラス（バイアス項を含む） """
    def __init__(self, N_in, N_h, K, L, h_1, dif_h_1, beta0, beta1, beta2, beta3):
        # N_in : 入力ユニットの数（バイアス項を除く）
        # N_h  : 隠れ層の数（バイアス項を除く）
        # K    : 混合密度の数
        # L    : 出力の次元
        # h_1  : 隠れ層の活性化関数
        # dif_h_1: 隠れ層の活性化関数の微分をz_jで表現した関数
        # beta0: w_1の初期分布β
        # beta1: w_pi_kの初期分布β
        # beta2: w_sig_kの初期分布β
        # beta3: w_mu_kの初期分布β 
        self.N_in = N_in+1; self.N_h = N_h+1 
        self.K = K; self.L = L; self.h_1 = h_1
        self.dif_h_1 = dif_h_1
        self.w_1 = matrix(RDF, [[gauss(0,beta0) for i in range(self.N_in)] for j in range(self.N_h-1)])
        self.w_pi_k = matrix(RDF, [[gauss(0,beta1) for j in range(self.N_h)] for k in range(self.K)]) 
        self.w_sig_k = matrix(RDF, [[gauss(0,beta2) for j in range(self.N_h)] for k in range(self.K)])
        self.w_mu_kl = [matrix(RDF, [[gauss(0,beta3) for j in range(self.N_h)] for k in range(self.K)]) for l in range(self.L)] 
        self.z_j = vector(RDF, self.N_h); self.z_j[0] = 1
        self.a_pi_k = vector(RDF, self.K)
        self.a_sig_k = vector(RDF, self.K)
        self.a_mu_kl = matrix(RDF, self.L, self.K, 0)
        self.d_k = vector(RDF, self.L)
        self.gam_nk = vector(RDF, self.K)
        self.Ident = identity_matrix(RDF,self.L)
        self.N_nk = vector(RDF, self.K)
    # ベクトルのガウス分布を定義
    def gauss_v(self, v, mu, sigma):
        d = len(v)
        val = -(v - mu)*(v - mu)/(2*sigma)
        a = n(2*pi*sigma)
        c = n(a**(-d*0.5))
        return n(c * (e**val))
    # 予測分布
    def getP(self, x, t):
        self.feedForward(x)
        self.backPropagate(t)
        return n(self.pi_k*self.N_nk)
    # フィードフォワード
    def feedForward(self, x): 
        # x : 入力値 
        self.x_i = vector(RDF, flatten([1, x])) 
        # フィードフォワード処理
        for j in range(1, self.N_h): 
            self.z_j[j] = self.h_1(self.w_1.row(j-1)*self.x_i).n()
        self.a_pi_k = self.w_pi_k*self.z_j
        self.a_sig_k = self.w_sig_k*self.z_j
        for l in range(self.L):
            self.a_mu_kl[l] = self.w_mu_kl[l]*self.z_j
        # 微分係数用の変数セット
        self.e_a_pi = vector(RDF, self.a_pi_k).apply_map(lambda x : n(e**x))
        self.e_a_sig_k = vector(RDF, self.a_sig_k).apply_map(lambda x : n(e**x))
        self.sig_k_sq = vector(RDF, self.e_a_sig_k).apply_map(lambda x : n(x*x))
        sum_e_a_pi_k = sum(self.e_a_pi.list())
        self.pi_k = vector(RDF, self.e_a_pi).apply_map(lambda x :x/sum_e_a_pi_k)
        ret = flatten([self.pi_k.list(), self.a_mu_kl.list(), self.a_sig_k.list()])
        return ret
    # バックプロパゲート処理 
    def backPropagate(self, t): 
        # t : 観測値
        v = vector(flatten([t]))
        self.t = vector(RDF, flatten([t]))
        mu_kl = matrix(RDF, self.a_mu_kl)
        sum_pi_N = 0
        for k in range(self.K):
            self.N_nk[k] = self.gauss_v(v, mu_kl.column(k), self.sig_k_sq[k])
            self.gam_nk[k] = self.pi_k[k]*self.N_nk[k]
            sum_pi_N += self.gam_nk[k]
        # avoid 0 devide.
        if sum_pi_N == 0: sum_pi_N = 1
        self.gam_nk = vector(RDF, self.gam_nk).apply_map(lambda x :x/sum_pi_N)
        self.dEn_da_pi_k = self.pi_k - self.gam_nk
        self.dEn_da_mu_k = matrix(RDF, self.K, self.L, 0)
        self.dEn_da_sig_k = vector(RDF, self.K)
        for k in range(self.K):
            d = mu_kl.column(k) - v
            self.dEn_da_mu_k[k] = self.gam_nk[k]*(d/self.sig_k_sq[k])
            self.dEn_da_sig_k[k] = self.gam_nk[k]*(self.L - d*d/self.sig_k_sq[k])
        # d_jを計算
        self.d_j = vector(RDF, self.N_h)
        for j in range(1, self.N_h):
            self.d_j[j] += self.dEn_da_pi_k*self.w_pi_k.column(j)
            self.d_j[j] += self.dEn_da_sig_k*self.w_sig_k.column(j)
            for l in range(self.L):
                self.d_j[j] += self.dEn_da_mu_k.column(l)*self.w_mu_kl[l].column(j)
            self.d_j[j] *= self.dif_h_1(self.z_j[j])
    # En(w)を返す
    def En(self): 
        return -ln(n(self.pi_k*self.N_nk))
    # ▽Enを返す
    def gradEn(self): 
        dE_1 = [[self.d_j[j]*self.x_i[i] for i in range(self.N_in)] for j in range(1, self.N_h)]
        dE_pi = [[self.dEn_da_pi_k[k]*self.z_j[j] for j in range(self.N_h)] for k in range(self.K)]
        dE_sig = [[self.dEn_da_sig_k[k]*self.z_j[j] for j in range(self.N_h)] for k in range(self.K)]  
        dE_mu = [[[self.dEn_da_mu_k[k][l]*self.z_j[j] for j in range(self.N_h)] for k in range(self.K)] for l in range(self.L)]
        return vector(RDF, flatten(dE_1 + dE_pi + dE_sig + dE_mu)) 
    # wを返す
    def getW(self): 
        return vector(RDF, flatten(self.w_1.list() + self.w_pi_k.list() + self.w_sig_k.list() 
        + [self.w_mu_kl[l].list() for l in range(self.L)]))
    # wをセットする
    def setW(self, w): 
        s_pi = self.N_in*(self.N_h-1)
        s_sig = self.N_in*(self.N_h-1)+self.N_h*self.K
        s_mu = self.N_in*(self.N_h-1)+2*self.N_h*self.K
        def index(l):
            return s_mu+(l*self.N_h*self.K)
        self.w_1 = matrix(RDF, self.N_h-1, self.N_in, w.list()[0:s_pi])   
        self.w_pi_k = matrix(RDF, self.K, self.N_h, w.list()[s_pi:s_sig])
        self.w_sig_k = matrix(RDF, self.K, self.N_h, w.list()[s_sig:s_mu])
        self.w_mu_kl = [matrix(RDF, self.K, self.N_h, w.list()[index(l):index(l+1)]) for l in range(self.L)]
}}


*** ランダムな重みを使った計算 [#c05dba7b]

MixtureDensityNetworkの計算にランダムな重みを使った結果を以下に示します。

式(5.148)の分布が平均や分散に強く依存しているため、すべてランダムな重みで計算すると
なかなか収束しません。


sageへの入力:
#pre{{
# 定数設定 
beta = 0.5
L=1
K=3
N_h = 5 
LEARN_COUNT = 30
RESET_COUNT = 25
# ニューラルネットワークにバグあり、h'(x) = 1 -z^2を使えば、a_jを保持する必要なし！
# もう一つh_1, dif_h_1の定義はlambda式を使う
h_1 = lambda x : tanh(x)
dif_h_1 = lambda z : 1 - z*z
}}


sageへの入力:
#pre{{
# ニューラルネットワークの生成 
net = MixtureDensityNetwork(1, 5, 3, 1, h_1, dif_h_1, 0.5, 0.5, 0.5, 0.5) 
# wの初期値を保存 
saved_w = net.getW() 
# スケール共役勾配法による学習 
_learn_csg(net, t, X, LEARN_COUNT, RESET_COUNT)
}}
sageからの出力:
#pre{{
2 0.853606593211 220.691284013135
2 -1.02380954822 295.681734779221
4 49619.0430098 220.688261274285
途中省略
29 0.97008460868 -5.34286542507738
30 0.799581594259 -10.4327483758505
}}


sageへの入力:
#pre{{
# 結果のプロット 30回目
var('x') 
b_plt = plot(lambda x : net.feedForward(x)[0], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[1], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[2], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
b_plt = plot(lambda x : net.feedForward(x)[3], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[4], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[5], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
}}

&ref(sage0-4.png);

&ref(sage1.png);

*** k-meansクラスター分析の結果を重みの初期値とする [#gf6f8c67]

収束をよくするために、重みのバイアス項の初期値をk-meansクラスター分析の結果を
使うことにします。

データをクラスタ分析した結果を以下に示します。


sageへの入力:
#pre{{
def euclid(x, y):
    return sqrt(sum([(x[i]-y[i])**2 for i in range(len(x))]))
# k-meansでデータを分類し、平均の重みを計算する
# 集合知から引用
def kcluster(rows,distance,k=3): 
    # Determine the minimum and maximum values for each point 
    ranges=[(min([row[i] for row in rows]),max([row[i] for row in rows])) 
    for i in range(len(rows[0]))]
    
    # Create k randomly placed centroids 
    clusters=[[random()*(ranges[i][1]-ranges[i][0])+ranges[i][0] 
    for i in range(len(rows[0]))] for j in range(k)]
    sig = [[] for i in range(k)]
    lastmatches=None 
    for t in range(100):
        print 'Iteration %d' % t 
        bestmatches=[[] for i in range(k)]
        
        # Find which centroid is the closest for each row 
        for j in range(len(rows)):
            row=rows[j] 
            bestmatch=0 
            for i in range(k):
                d=distance(clusters[i],row)
                if d<distance(clusters[bestmatch],row): 
                    bestmatch=i 
            bestmatches[bestmatch].append(j)
        # If the results are the same as last time, this is complete 
        if bestmatches==lastmatches: break 
        lastmatches=bestmatches
        # Move the centroids to the average of their members 
        for i in range(k):
            v = [0.0]*len(rows[0])
            avgs=[0.0]*len(rows[0]) 
            n = len(bestmatches[i])
            if len(bestmatches[i])>0:
                sumX = [0.0]*len(rows[0])
                sumXX = [0.0]*len(rows[0])
                for rowid in bestmatches[i]:
                    for m in range(len(rows[rowid])):
                        avgs[m]+=rows[rowid][m] 
                        sumX[m] += rows[rowid][m]
                        sumXX[m] += rows[rowid][m]*rows[rowid][m]
                    
                for j in range(len(avgs)):
                    avgs[j]/= n
                    v[j] = sumXX[j]/n - (sumX[j]/n)*(sumX[j]/n)
                sig[i] = v
                clusters[i]=avgs
    print clusters
    print sig
    return bestmatches
data = zip(t,X)  
group = kcluster(data, distance=euclid, k=3)
}}
sageからの出力:
#pre{{
Iteration 0
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
[[0.755995147678366, 0.906635821415391], [0.489412268723978, 0.535817638206947], 
[0.325834642861869, 0.133408609458438]]
[[0.0184016725675581, 0.00227560546930994], [0.00822710110832403, 0.0257078799298951], 
[0.0280504201990768, 0.00434725503090522]]
}}


sageへの入力:
#pre{{
c_plt = Graphics()
r_grp = [ data[i] for i in group[0]]
b_grp = [ data[i] for i in group[1]]
g_grp = [ data[i] for i in group[2]]
for pt in b_grp:
    c_plt += point(pt, rgbcolor="blue")
for pt in g_grp:
    c_plt += point(pt, rgbcolor="green")
for pt in r_grp:
    c_plt += point(pt, rgbcolor="red")
c_plt.show()
}}

&ref(sage0-5.png);


*** 初期値を変えて回帰分析 [#t7ffd506]

初期値を変えて回帰分析をし直します。

今回は、バイアス項を

- \(\pi_1 = 0.4, \pi_2 = 0.2, \pi_3 = 0.4\)の対数値

- \(\sigma_1 = 0.025, \sigma_2 = 0.125, \sigma_3 = 0.025\)の対数値

- \(\mu_1 = 0.13, \mu_2 = 0.53, \mu_3 = 0.91\)
とし、重みの分散も\(w_1 = 1, w_{\pi} = 0.25, w_{\sigma} = 0.1, w_{\mu} = 0.1\)と
しました。


sageへの入力:
#pre{{
# ニューラルネットワークの生成 
net = MixtureDensityNetwork(1, 5, 3, 1, h_1, dif_h_1, 1, 0.25, 0.1, 0.1)  
# wの初期値を保存 
saved_w = net.getW() 
tmp = saved_w
tmp[10] = -0.9  # pi_1 ln(0.4)
tmp[16] = -1.6  # pi_2 ln(0.2)
tmp[22] = -0.9  # pi_3 ln(0.4)
tmp[28] = -3.69 # sig_1 ln(0.025) Fig.5.21から4σ = 0.5と読み取る
tmp[34] = -2.08 # sig_2 ln(0.125) Fig.5.21から4σ = 0.25と読み取る
tmp[40] = -3.69 # sig_3 ln(0.025) Fig.5.21から4σ = 0.5と読み取る
tmp[46] = 0.13  # mu_1
tmp[52] = 0.53  # mu_2
tmp[58] = 0.91  # mu_3
# スケール共役勾配法による学習 
net.setW(tmp)
_learn_csg(net, t, X, LEARN_COUNT, RESET_COUNT)
}}
sageからの出力:
#pre{{
2 1.00140531678 162.746607216641
3 0.978118009568 94.6561189652474
4 1.00006922311 74.3749728302636
途中省略
28 0.989452644705 -76.7534334638356
29 1.32320485458 -81.6847672263680
30 0.876928985373 -86.8153622565841
}}


sageへの入力:
#pre{{
# 結果のプロット 30回目
var('x') 
b_plt = plot(lambda x : net.feedForward(x)[0], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[1], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[2], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show(ymax=1, ymin=-1)
b_plt = plot(lambda x : net.feedForward(x)[3], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[4], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[5], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
}}

&ref(sage0-6.png);

&ref(sage1-1.png);


sageへの入力:
#pre{{
tmp_w = net.getW()
save(tmp_w, DATA+'tmp_w')
}}


sageへの入力:
#pre{{
_learn_csg(net, t, X, 100, RESET_COUNT)
}}
sageからの出力:
#pre{{
2 0.999254054027 -87.3044338521603
3 0.982374480725 -91.6889548386661
途中省略
99 1.01213525834 -188.881647537428
100 1.03434125073 -189.734542888085
}}


sageへの入力:
#pre{{
# 結果のプロット 130回目
var('x') 
b_plt = plot(lambda x : net.feedForward(x)[0], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[1], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[2], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show(ymax=1, ymin=-1)
b_plt = plot(lambda x : net.feedForward(x)[3], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[4], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[5], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
}}

&ref(sage0-7.png);

&ref(sage1-2.png);



sageへの入力:
#pre{{
tmp_w = net.getW()
save(tmp_w, DATA+'tmp_w')
_learn_csg(net, t, X, 100, RESET_COUNT)
}}
sageからの出力:
#pre{{
2 1.00010153793 -189.851273562036
3 1.00013157633 -190.041029348136
途中省略
99 0.962650552011 -206.426993641372
100 1.03840031165 -206.497416181699
}}


sageへの入力:
#pre{{
# 結果のプロット 230回目
var('x') 
b_plt = plot(lambda x : net.feedForward(x)[0], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[1], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[2], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
b_plt = plot(lambda x : net.feedForward(x)[3], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[4], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[5], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
}}

&ref(sage0-8.png);

&ref(sage1-3.png);



sageへの入力:
#pre{{
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1])
}}

&ref(sage0-9.png);


sageへの入力:
#pre{{
w_230 = net.getW()
save(w_230, DATA+'w_230')
_learn_csg(net, t, X, 100, RESET_COUNT)
}}
sageからの出力:
#pre{{
2 1.00024395346 -206.557574296743
3 1.00161969024 -206.711515824961
途中省略
99 1.03276294743 -212.624652019776
100 0.971359659233 -212.656180208402
}}


*** 最終結果 [#vec210c8]

最終結果は、PRMLの図5.21

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.21a.jpg,center,40%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.21b.jpg,center,40%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.21c.jpg,center,40%);

の結果と良く似たものになりました。


sageへの入力:
#pre{{
# 結果のプロット 330回目
var('x') 
b_plt = plot(lambda x : net.feedForward(x)[0], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[1], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[2], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
b_plt = plot(lambda x : net.feedForward(x)[3], [x, 0, 1], color='blue') 
g_plt = plot(lambda x : net.feedForward(x)[4], [x, 0, 1], color='green') 
r_plt = plot(lambda x : net.feedForward(x)[5], [x, 0, 1], color='red')
(b_plt+g_plt+r_plt).show()
}}

&ref(sage0-10.png);

&ref(sage1-4.png);



sageへの入力:
#pre{{
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=True, contours=20)
}}

&ref(sage0-11.png);

sageへの入力:
#pre{{
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=False, cmap='hsv', contours=20)
}}

&ref(sage0-12.png);


sageへの入力:
#pre{{
w_330 = net.getW()
save(w_330, DATA+'w_330')
}}


*** 予測値 [#tb0a7dba]

混合密度の平均として、PRMLでは条件付き密度の近似的な条件付きモード

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.21d.jpg,center,40%);

をプロットしていますが、ここでは本文の説明にある混合係数が最大の要素の平均をとりました。


sageへの入力:
#pre{{
# π_kが最大の分布の平均
def probAve(net, x):
    net.feedForward(x)
    max_pi_index = [ k for k in range(net.K) if max(net.pi_k) == net.pi_k[k]][0]
    return net.a_mu_kl[0][max_pi_index]
}}



sageへの入力:
#pre{{
# 最終的な解をプロット
N = 100
probAve_plt = list_plot([[x/N, probAve(net, x/N)] for x in range(100)], rgbcolor='red')
(probAve_plt + data2_plt).show()
}}

&ref(sage0-13.png);


** プログラムのチェック [#y11969c3]

今回、プログラムが収束しない現象続き、とても悩みました！

そこで、5.3.3逆伝搬の効率で紹介されている
$$
\frac{\partial E_n}{\partial w_{ji}} = \frac{En(w_{ji} + \epsilon) - En(w_{ji} - \epsilon )}{2 \epsilon}
$$
を使って重みを近似計算してみました。これと手計算の値を比較しながら、プログラムをデバッグしました。


sageへの入力:
#pre{{
# ナブラを数値計算で求めて確認する
def grad_w(net, x, t):
    eps = 1e-6
    # 現在のwを保存
    saved_w = net.getW()
    grad = vector(RDF, len(saved_w))
    work_w = vector(RDF, saved_w)
    for i in range(len(work_w)):
        w = work_w[i]
        work_w[i] = w - eps
        net.setW(work_w)
        net.feedForward(x)
        net.backPropagate(t)
        E1 = net.En()
        work_w[i] = w + eps
        net.setW(work_w)
        net.feedForward(x)
        net.backPropagate(t)
        E2 = net.En()
        work_w[i] = w
        grad[i] = (E2 - E1)/(2*eps)
    net.setW(saved_w)
    return grad
}}


sageへの入力:
#pre{{
# 入力1点、隠れ層3個、出力1点
N_in = 1 
N_h = 6     
N_out = 1
# 隠れ層の活性化関数
h_1(x) = tanh(x)
# 出力層の活性化関数
h_2(x) = x
dif_h_1 = lambda z : 1 - z*z
# mlpでチェック
mlp = NeuralNetwork(1, 1, 1, h_1, h_2, 0.1)
}}


*** ニューラルネットワークの場合 [#ua7f6763]

最初に単純なニューラルネットワークの場合について近似計算とネットワークの計算を比較しました。
非常に良い一致を得ました。


sageへの入力:
#pre{{
mlp.feedForward(0)
mlp.backPropagate(1)
mlp.gradEn()
}}
sageからの出力:
#pre{{
(0.0668611517486, 0.0, 0.0207058923852, 0.0, -1.05811442086, 0.258067961374)
}}


sageへの入力:
#pre{{
grad_w(mlp, 0, 1)
}}
sageからの出力:
#pre{{
(0.0, 0.0, 0.0207058923896, 0.0, -1.05811442086, 0.258067961278)
}}


*** 混合密度ネットワークの場合 [#bd405dc8]

次に混合密度ネットワークの場合についても近似計算とネットワークの計算を比較しました。
複雑な処理をしているにも関わらず、かなり一致することに驚きました。


sageへの入力:
#pre{{
# mdnでチェック
mdn = MixtureDensityNetwork(1, 1, 1, 1, h_1, dif_h_1, 0.5, 0.2, 0.05, 0.05)
# テスト用の重み
w =vector(RDF, [-1.14339938022, 0.626829670558, 0.147625024973, 0.269208141326, -0.126939642185, 
-0.0395784292044, 0.0589154048316, 0.0677724926114])
mdn.setW(w)
}}


sageへの入力:
#pre{{
mdn.feedForward(0)
mdn.backPropagate(1)
mdn.gradEn()
}}
sageからの出力:
#pre{{
(-0.024679341566, -0.0, 0.0, -0.0, -0.199641572236, 0.162818802495, -1.20402804974, 0.981951820123)
}}


sageへの入力:
#pre{{
grad_w(mdn, 0, 1)
}}
sageからの出力:
#pre{{
(-0.0246793415704, 0.0, 0.0, 0.0, -0.199641572274, 0.162818802441, -1.20402804982, 0.981951820078)
}}


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

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

       

Site admin: anonymous

PukiWiki 1.4.7 Copyright © 2001-2006 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji. Powered by PHP 5.1.6. HTML convert time: 0.007 sec.