![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2011/06/20からのアクセス回数 6546
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://sage.math.canterbury.ac.nz/home/pub/104/
また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
PRMLの第5章-混合密度ネットワークの例題、図5.19
のデータは、XとYを入れ替えただけで多峰性によりフィッティングが失敗する例です。 （x=0.5のところで、yの値として3つの選択しがあります）
区間（0, 1）に一様分布する変数$x_n$をサンプリングし、目的値\(t_n\) を\(x_n + 0.3 sin(2\pi x_n)\)に区間(-0.1, 0.1)の一様乱数を 加えたデータを200点生成します。
sageへの入力:
# 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への入力:
# データの保存 # save(X, DATA+'X') # save(t, DATA+'t') # リストアするには、以下のコメントを外す X = load(DATA+'X') t = load(DATA+'t') print 'Done'
sageからの出力:
Done
生成したデータとx,yを入れ替えたグラフをプロットします。
sageへの入力:
# データ生成 data_plt = list_plot(zip(X,t)) data_plt.show()
sageへの入力:
# X, tを入れ替えた図 data2_plt = list_plot(zip(t,X)) data2_plt.show()
生成したデータをニューラルネットワークを使って、フィッティングさせたグラフを 以下に示します。（入力層1個、隠れ層6個、出力層1個で計算）
図5.19のbの様にx,tを入れ替えた場合には上手くフィッティングしません。
sageへの入力:
# 入力1点、隠れ層3個、出力1点 N_in = 1 N_h = 6 N_out = 1 # 隠れ層の活性化関数 h_1(x) = tanh(x) # 出力層の活性化関数 h_2(x) = x
sageへの入力:
# ニューラルネットワーク用のクラスのロード attach "NeuralNetwork.sage" # Scaled Conjugate Gradients（スケール共役勾配法）のロード attach "SCG.sage"
sageへの入力:
# 出力関数
def _f(x, net):
y = net.feedForward(x)
return y[0]
sageへの入力:
# 定数設定 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からの出力:
WARNING: Output truncated! 2 0.961554425907 3.77799819103 3 1.04117097273 2.41177356492 途中省略 499 1.14719236412 0.281449667896 500 0.626481765163 0.281421869851
sageへの入力:
var('x')
g = lambda x : _f([x], net)
g_plt = plot(g, [x, 0, 1], rgbcolor='red')
(g_plt + data_plt).show()
sageへの入力:
# 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からの出力:
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への入力:
g2_plt = plot(g, [x, 0, 1], rgbcolor='red') (g2_plt + data2_plt).show()
混合モデルの条件付き密度関数\(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) $$ で表し、モデルパラメータ
混合計数は、式(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クラスを以下のように定義します。
sageへの入力:
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)]
MixtureDensityNetworkの計算にランダムな重みを使った結果を以下に示します。
式(5.148)の分布が平均や分散に強く依存しているため、すべてランダムな重みで計算すると なかなか収束しません。
sageへの入力:
# 定数設定 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への入力:
# ニューラルネットワークの生成 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からの出力:
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への入力:
# 結果のプロット 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()
収束をよくするために、重みのバイアス項の初期値をk-meansクラスター分析の結果を 使うことにします。
データをクラスタ分析した結果を以下に示します。
sageへの入力:
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からの出力:
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への入力:
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()
初期値を変えて回帰分析をし直します。
今回は、バイアス項を
sageへの入力:
# ニューラルネットワークの生成 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からの出力:
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への入力:
# 結果のプロット 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()
sageへの入力:
tmp_w = net.getW() save(tmp_w, DATA+'tmp_w')
sageへの入力:
_learn_csg(net, t, X, 100, RESET_COUNT)
sageからの出力:
2 0.999254054027 -87.3044338521603 3 0.982374480725 -91.6889548386661 途中省略 99 1.01213525834 -188.881647537428 100 1.03434125073 -189.734542888085
sageへの入力:
# 結果のプロット 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()
sageへの入力:
tmp_w = net.getW() save(tmp_w, DATA+'tmp_w') _learn_csg(net, t, X, 100, RESET_COUNT)
sageからの出力:
2 1.00010153793 -189.851273562036 3 1.00013157633 -190.041029348136 途中省略 99 0.962650552011 -206.426993641372 100 1.03840031165 -206.497416181699
sageへの入力:
# 結果のプロット 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()
sageへの入力:
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1])
sageへの入力:
w_230 = net.getW() save(w_230, DATA+'w_230') _learn_csg(net, t, X, 100, RESET_COUNT)
sageからの出力:
2 1.00024395346 -206.557574296743 3 1.00161969024 -206.711515824961 途中省略 99 1.03276294743 -212.624652019776 100 0.971359659233 -212.656180208402
最終結果は、PRMLの図5.21
の結果と良く似たものになりました。
sageへの入力:
# 結果のプロット 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()
sageへの入力:
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=True, contours=20)
sageへの入力:
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=False, cmap='hsv', contours=20)
sageへの入力:
w_330 = net.getW() save(w_330, DATA+'w_330')
混合密度の平均として、PRMLでは条件付き密度の近似的な条件付きモード
をプロットしていますが、ここでは本文の説明にある混合係数が最大の要素の平均をとりました。
sageへの入力:
# π_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への入力:
# 最終的な解をプロット N = 100 probAve_plt = list_plot([[x/N, probAve(net, x/N)] for x in range(100)], rgbcolor='red') (probAve_plt + data2_plt).show()
今回、プログラムが収束しない現象続き、とても悩みました！
そこで、5.3.3逆伝搬の効率で紹介されている $$ \frac{\partial E_n}{\partial w_{ji}} = \frac{En(w_{ji} + \epsilon) - En(w_{ji} - \epsilon )}{2 \epsilon} $$ を使って重みを近似計算してみました。これと手計算の値を比較しながら、プログラムをデバッグしました。
sageへの入力:
# ナブラを数値計算で求めて確認する
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への入力:
# 入力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)
最初に単純なニューラルネットワークの場合について近似計算とネットワークの計算を比較しました。 非常に良い一致を得ました。
sageへの入力:
mlp.feedForward(0) mlp.backPropagate(1) mlp.gradEn()
sageからの出力:
(0.0668611517486, 0.0, 0.0207058923852, 0.0, -1.05811442086, 0.258067961374)
sageへの入力:
grad_w(mlp, 0, 1)
sageからの出力:
(0.0, 0.0, 0.0207058923896, 0.0, -1.05811442086, 0.258067961278)
次に混合密度ネットワークの場合についても近似計算とネットワークの計算を比較しました。 複雑な処理をしているにも関わらず、かなり一致することに驚きました。
sageへの入力:
# 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への入力:
mdn.feedForward(0) mdn.backPropagate(1) mdn.gradEn()
sageからの出力:
(-0.024679341566, -0.0, 0.0, -0.0, -0.199641572236, 0.162818802495, -1.20402804974, 0.981951820123)
sageへの入力:
grad_w(mdn, 0, 1)
sageからの出力:
(-0.0246793415704, 0.0, 0.0, 0.0, -0.199641572274, 0.162818802441, -1.20402804982, 0.981951820078)
皆様のご意見、ご希望をお待ちしております。
