第5章-混合密度ネットワークをSageで試す
PRMLの第5章-混合密度ネットワークの例題、図5.19
テストデータ
区間（0, 1）に一様分布する変数$x_n$をサンプリングし、目的値$t_n$ を$x_n + 0.3 sin(2\pi x_n)$に区間(-0.1, 0.1)の一様乱数を 加えたデータを200点生成します。
# 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]
|
|
# データの保存
# save(X, DATA+'X')
# save(t, DATA+'t')
# リストアするには、以下のコメントを外す
X = load(DATA+'X')
t = load(DATA+'t')
print 'Done'
Done Done |
グラフのプロット
生成したデータとx,yを入れ替えたグラフをプロットします。
# データ生成
data_plt = list_plot(zip(X,t))
data_plt.show()
|
|
# X, tを入れ替えた図
data2_plt = list_plot(zip(t,X))
data2_plt.show()
|
|
ニューラルネットワークでの回帰
生成したデータをニューラルネットワークを使って、フィッティングさせたグラフを 以下に示します。（入力層1個、隠れ層6個、出力層1個で計算）
図5.19のbの様にx,tを入れ替えた場合には上手くフィッティングしません。
# 入力1点、隠れ層3個、出力1点
N_in = 1
N_h = 6
N_out = 1
# 隠れ層の活性化関数
h_1(x) = tanh(x)
# 出力層の活性化関数
h_2(x) = x
|
|
# ニューラルネットワーク用のクラスのロード
attach "NeuralNetwork.sage"
# Scaled Conjugate Gradients（スケール共役勾配法）のロード
attach "SCG.sage"
|
|
# 出力関数
def _f(x, net):
y = net.feedForward(x)
return y[0]
|
|
# 定数設定
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)
WARNING: Output truncated! full_output.txt 2 0.961554425907 3.77799819103 3 1.04117097273 2.41177356492 途中省略 499 1.14719236412 0.281449667896 500 0.626481765163 0.281421869851 WARNING: Output truncated! full_output.txt 2 0.961554425907 3.77799819103 3 1.04117097273 2.41177356492 途中省略 499 1.14719236412 0.281449667896 500 0.626481765163 0.281421869851 |
var('x')
g = lambda x : _f([x], net)
g_plt = plot(g, [x, 0, 1], rgbcolor='red')
(g_plt + data_plt).show()
|
|
# 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)
WARNING: Output truncated! full_output.txt 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 WARNING: Output truncated! full_output.txt 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 |
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) $$ で表し、モデルパラメータ
- 混合係数 $\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クラスの定義
MixtureDensityNetworkクラスを以下のように定義します。
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)の分布が平均や分散に強く依存しているため、すべてランダムな重みで計算すると なかなか収束しません。
# 定数設定
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
|
|
# ニューラルネットワークの生成
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)
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 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 |
# 結果のプロット 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クラスター分析の結果を重みの初期値とする
収束をよくするために、重みのバイアス項の初期値をk-meansクラスター分析の結果を 使うことにします。
データをクラスタ分析した結果を以下に示します。
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)
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]] 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]] |
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()
|
|
初期値を変えて回帰分析
初期値を変えて回帰分析をし直します。
今回は、バイアス項を
- $\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$
# ニューラルネットワークの生成
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)
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 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 |
# 結果のプロット 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()
|
|
tmp_w = net.getW()
save(tmp_w, DATA+'tmp_w')
|
|
_learn_csg(net, t, X, 100, RESET_COUNT)
2 0.999254054027 -87.3044338521603 3 0.982374480725 -91.6889548386661 途中省略 99 1.01213525834 -188.881647537428 100 1.03434125073 -189.734542888085 2 0.999254054027 -87.3044338521603 3 0.982374480725 -91.6889548386661 途中省略 99 1.01213525834 -188.881647537428 100 1.03434125073 -189.734542888085 |
# 結果のプロット 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()
|
|
tmp_w = net.getW()
save(tmp_w, DATA+'tmp_w')
_learn_csg(net, t, X, 100, RESET_COUNT)
2 1.00010153793 -189.851273562036 3 1.00013157633 -190.041029348136 途中省略 99 0.962650552011 -206.426993641372 100 1.03840031165 -206.497416181699 2 1.00010153793 -189.851273562036 3 1.00013157633 -190.041029348136 途中省略 99 0.962650552011 -206.426993641372 100 1.03840031165 -206.497416181699 |
# 結果のプロット 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()
|
|
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1])
|
|
w_230 = net.getW()
save(w_230, DATA+'w_230')
_learn_csg(net, t, X, 100, RESET_COUNT)
2 1.00024395346 -206.557574296743 3 1.00161969024 -206.711515824961 途中省略 99 1.03276294743 -212.624652019776 100 0.971359659233 -212.656180208402 2 1.00024395346 -206.557574296743 3 1.00161969024 -206.711515824961 途中省略 99 1.03276294743 -212.624652019776 100 0.971359659233 -212.656180208402 |
最終結果
最終結果は、PRMLの図5.21
# 結果のプロット 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()
|
|
var('x y')
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=True, contours=20)
|
|
contour_plot(lambda x, y : net.getP(x, y), [x, 0, 1], [y, 0, 1], fill=False, cmap='hsv', contours=20)
|
|
w_330 = net.getW()
save(w_330, DATA+'w_330')
|
|
予測値
混合密度の平均として、PRMLでは条件付き密度の近似的な条件付きモード
# π_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]
|
|
# 最終的な解をプロット
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} $$ を使って重みを近似計算してみました。これと手計算の値を比較しながら、プログラムをデバッグしました。
# ナブラを数値計算で求めて確認する
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
|
|
# 入力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)
|
|
ニューラルネットワークの場合
最初に単純なニューラルネットワークの場合について近似計算とネットワークの計算を比較しました。 非常に良い一致を得ました。
mlp.feedForward(0)
mlp.backPropagate(1)
mlp.gradEn()
(0.0668611517486, 0.0, 0.0207058923852, 0.0, -1.05811442086, 0.258067961374) (0.0668611517486, 0.0, 0.0207058923852, 0.0, -1.05811442086, 0.258067961374) |
grad_w(mlp, 0, 1)
(0.0, 0.0, 0.0207058923896, 0.0, -1.05811442086, 0.258067961278) (0.0, 0.0, 0.0207058923896, 0.0, -1.05811442086, 0.258067961278) |
混合密度ネットワークの場合
次に混合密度ネットワークの場合についても近似計算とネットワークの計算を比較しました。 複雑な処理をしているにも関わらず、かなり一致することに驚きました。
# 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)
|
|
mdn.feedForward(0)
mdn.backPropagate(1)
mdn.gradEn()
(-0.024679341566, -0.0, 0.0, -0.0, -0.199641572236, 0.162818802495, -1.20402804974, 0.981951820123) (-0.024679341566, -0.0, 0.0, -0.0, -0.199641572236, 0.162818802495, -1.20402804974, 0.981951820123) |
grad_w(mdn, 0, 1)
(-0.0246793415704, 0.0, 0.0, 0.0, -0.199641572274, 0.162818802441, -1.20402804982, 0.981951820078) (-0.0246793415704, 0.0, 0.0, 0.0, -0.199641572274, 0.162818802441, -1.20402804982, 0.981951820078) |
|
|