５つの正規分布をバラまいたのに、データに合わせて最終的に３つになっているのがわかると思います。
本稿ではVBEMアルゴリズムの更新式とその実装について書きます。
ハイパーパラメータの初期化
このアルゴリズムには $\alpha_{0}$ $\beta_{0}$ $\boldsymbol{m}_{0}$ $\boldsymbol{W}_{0}$ $\nu_{0}$ というハイパーパラメータがあり、当然これらの値は結果に影響します。
これらは後に示すソースコードの set_hyper_parameters メソッド内で初期化を行っているので、各値はそちらを参照してください。ハイパーパラメータの調整は非常に職人的なところがあるようで、私も自信をもって説明できることがあまりないです。
このアルゴリズムの利点としては、対象のデータに混合正規分布を当てはめる際、混合要素数が未知であっても更新過程で適切な要素数になる（ある混合要素がどこかの点だけを極端に説明するようなことが起きづらい）ということが挙げられます。
その最終的な混合要素数を主に調整するのが $\alpha_{0}$ で、 $\alpha_{0}$ が小さいほど混合係数が0になる解が選ばれやすく、冒頭の例のように非零の混合係数が3つあるような解を得るためには $\alpha_{0} = 1$ とするといい感じになります。
VBEMアルゴリズムの更新式
以下に変分混合ガウス分布におけるVBEMアルゴリズムの更新式を書いていきます。導出過程はとても長いのでPRMLを参照してください。アルゴリズムとしては各ハイパーパラメータと負担率を初期化した後に、これらをひたすら更新し続けます。
N_k &=& \sum_{n}^{} r_{nk} \\
\boldsymbol{\overline{x}}_{k} &=& \frac{1}{N_{k}}\sum_{n}^{} r_{nk}\boldsymbol{x}_{n} \\
S_{kd} &=& \frac{1}{N_{k}}\sum_{n}^{} r_{nk}(\boldsymbol{x}_{n} - \boldsymbol{\overline{x}}_{k})(\boldsymbol{x}_{n} - \boldsymbol{\overline{x}}_{k})^{\mathrm{T}}
\end{eqnarray}
\alpha_{k} &=& \alpha_{0} + N_{k} \\
\beta_{k} &=& \beta_{0} + N_{k} \\
m_{k} &=& \frac{\beta_{0}m_{0} + N_{k}\boldsymbol{\overline{x}}_{k}}{\beta_{k}} \\
\boldsymbol{W}_{k}^{-1} &=& N_{k}\boldsymbol{S_{k}} + \frac{\beta_{0}N_{k}}{\beta_{0}+N_{k}}(\boldsymbol{\overline{x}}_{k} - \boldsymbol{m}_{0})(\boldsymbol{\overline{x}}_{k} - \boldsymbol{m}_{0})^{\mathrm{T}} \\
\nu_{k} &=& \nu_{0} + N_{k} \\
\end{eqnarray}
\mathbb{E}[(\boldsymbol{x}_{n} - \boldsymbol{\mu}_{k})^{\mathrm{T}}\boldsymbol{\Lambda}_{k}(\boldsymbol{x}_{n} - \boldsymbol{\mu}_{k})] &=& D\beta_{k}^{-1} + \nu_{k}(\boldsymbol{x}_{n} - \boldsymbol{m}_{k})^{\mathrm{T}}\boldsymbol{W}_{k}(\boldsymbol{x}_{n} - \boldsymbol{m}_{k}) \\
\mathbb{E}[\mathrm{ln}\,|\boldsymbol{\Lambda}_{k}|] &=& \sum_{i=1}^{D}\psi(\frac{\nu_{k}+1-i}{2}) + D\mathrm{ln}\,2 + \mathrm{ln}\,|\boldsymbol{W}_{k}| \\
\mathbb{E}[{\rm ln}\,\pi_{k}] &=& \psi(\alpha_{k}) - \psi(\sum_{j=1}^{K} \alpha_{j})
\end{eqnarray}
{\rm ln}\,\rho_{nk} &=& \mathbb{E}[{\rm ln}\,\pi_{k}] + \frac{1}{2}\mathbb{E}[\mathrm{ln}\,|\boldsymbol{\Lambda}_{k}|] - \frac{D}{2}\mathrm{ln}(2\pi) - \frac{1}{2}\mathbb{E}[(\boldsymbol{x}_{n} - \boldsymbol{\mu}_{k})^{\mathrm{T}}\boldsymbol{\Lambda}_{k}(\boldsymbol{x}_{n} - \boldsymbol{\mu}_{k})] \\
r_{nk} &=& \frac{\rho_{nk}}{\sum_{j=1}^{K} \rho_{nj}}
\end{eqnarray}
ソースコード
アルゴリズムの計算を行っているファイルのみ以下に示します。その他描画を行っているファイルなどは https://gist.github.com/seinosuke/4f565647f33c815e18123404a224a3e1 にあります。そこにある３つのファイルを同じディレクトリにおいて main.rb を実行すれば冒頭のようなGIFが出力されます。
Rubyは2.2.0ですが、そんなに無茶なことはしてないので多分どのバージョンでも動くと思います。
class VBEM
attr_reader :m, :W, :r, :n, :k
attr_reader :E_lambda, :E_pi
alias :E_mu :m
def initialize(options = {})
@dim = options[:X].first.size
@k = options[:k]
@X = options[:X].map { |x| Matrix[x.map(&:to_f)] }
@num = options[:X].size
set_hyper_parameters
rand_init
end
# 1ステップ進める
def update
vb_m_step
vb_e_step
p @E_pi
end
# ハイパーパラメータを初期化
def set_hyper_parameters
@alpha_0 = Matrix[ Array.new(@k) { 1e-0 } ]
@beta_0 = Matrix[ Array.new(@k) { 1e-3 } ]
@m_0 = Matrix[ Array.new(@dim) { 0.0 } ]
@W_0 = Matrix[ *Matrix.I(@dim).to_a.map { |row| row.map(&:to_f) } ]
@invW_0 = @W_0.inv
@nu_0 = Matrix[ Array.new(@k) { @dim } ]
end
# ランダムに負担率を初期化
def rand_init
@r = Array.new(@num) { Array.new(@k) { rand } }
@r = @r.map do |row_n|
sum = row_n.inject(:+)
row_n.map { |v| v / sum }
end
end
############################################################
# VB Mstep
############################################################
def vb_m_step
@n = Matrix[ @r.transpose.map { |row| row.inject(:+) } ]
@alpha = @alpha_0 + @n
@beta = @beta_0 + @n
@nu = @nu_0 + @n
# 後の計算に使う (K mat(1 D))
@x_ = @k.times.map do |k|
sum = Matrix[Array.new(@dim) { 0.0 }]
@num.times { |n| sum += @r[n][k] * @X[n] }
sum / @n[0, k]
end
# 後の計算に使う (K mat(D D))
@S = @k.times.map do |k|
sum = Matrix.zero(@dim)
@num.times do |n|
tmp = (@X[n] - @x_[k]).t * (@X[n] - @x_[k])
sum += @r[n][k] * tmp
end
sum / @n[0, k]
end
# 平均のパラメータを更新 (K mat(1 D))
@m = @k.times.map do |k|
ret = Matrix[Array.new(@dim) { 0.0 }]
ret += @m_0 * @beta_0[0, k]
ret += @n[0, k]*@x_[k]
ret / @beta[0, k]
end
# 精度のパラメータを更新 (K mat(D D))
@W = @k.times.map do |k|
ret = Matrix.zero(@dim)
ret += @invW_0 + @n[0, k]*@S[k]
tmp = (@beta_0[0, k] * @n[0, k]) / (@beta_0[0, k] + @n[0, k])
ret += tmp * (@x_[k] - @m_0).t * (@x_[k] - @m_0)
ret.inv
end
end
############################################################
# VB Estep
############################################################
def vb_e_step
# PIの期待値 (1 K)
@E_pi = @k.times.map do |k|
@alpha[0, k] / @alpha.to_a.flatten.inject(:+)
end
# lnPIの期待値 (1 K)
@E_lnpi = @k.times.map do |k|
Utils.digamma(@alpha[0, k]) - Utils.digamma(@alpha.to_a.flatten.inject(:+))
end
# Λの期待値 (K mat(D D))
@E_lambda = @k.times.map do |k|
@nu[0, k] * @W[k]
end
# lnΛの期待値 (1 K)
@E_lnlambda = @k.times.map do |k|
ret = @dim * Math.log(2.0) + Math.log(@W[k].det)
ret += (1..@dim).inject(0.0) { |sum, i| Utils.digamma((@nu[0, k]+1-i) / 2.0) }
ret
end
# 後の計算に使う mat(N K)
@E_x_mWx_m = @num.times.map do |n|
@k.times.map do |k|
ret = @dim / @beta[0, k]
d = @X[n] - @m[k]
ret += @nu[0, k] * (d * @W[k] * d.t)[0, 0]
ret
end
end
@E_x_mWx_m = Matrix[*@E_x_mWx_m]
# lnrhoの計算で使う mat(N K)
@E_lnN = @num.times.map do |n|
@k.times.map do |k|
ret = @E_lnlambda[k] - @dim*Math.log(2.0*PI) - @E_x_mWx_m[n, k]
ret / 2.0
end
end
@E_lnN = Matrix[*@E_lnN]
# 負担率を更新 mat(N K)
@lnrho = @E_lnN + Matrix[ *Array.new(@num) { @E_lnpi } ]
@r = @lnrho.to_a.map do |row|
logsumexp = Math.log( row.inject(0.0) { |sum, v| sum + Math.exp(v) } )
row.map { |v| Math.exp(v - logsumexp) }
end.map do |row|
sum = row.inject(:+)
row.map do |v|
v = Utils.max(v, 1e-10)
v / sum
end
end
end
endおわりに
gnuplotで等高線をxy平面図でみたいとき、splotしてnosurfaceしてview0,0で上から覗く以外にうまいことやる方法はないもんでしょうか。
参考
C.M.ビショップ(2012) 「パターン認識と機械学習 下」 丸善出版.
動く変分混合ガウス分布（実装編）- 動く PRML シリーズ（２）
http://web.science.mq.edu.au/~mjohnson/code/digamma.c
0 件のコメント:
コメントを投稿