Processing math: 0%

2016年9月7日水曜日

【Ruby】 VBEMアルゴリズム（PRML10章）

今回はPRML10章中の「変分混合ガウス分布」の内容をRubyで実装してみました。
５つの正規分布をバラまいたのに、データに合わせて最終的に３つになっているのがわかると思います。
本稿では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を参照してください。アルゴリズムとしては各ハイパーパラメータと負担率を初期化した後に、これらをひたすら更新し続けます。

負担率から求められる統計量を更新
\begin{eqnarray} 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}

変分Mステップ
\begin{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}

変分Eステップ
$\psi(\cdot)$はディガンマ関数である.
\begin{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}

負担率を更新
\begin{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ですが、そんなに無茶なことはしてないので多分どのバージョンでも動くと思います。

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
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 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)