<file_sep>潜在変数モデルと学習法に関して
================================================================================

________________________________________________________________________________
はじめに
--------------------------------------------------------------------------------
### 自己紹介
はじめまして！いしはたです！

### 前置き
この記事は Machine Learning Advent Calendar の 12月8日 の記事として書かれています。  
思いの外長くなった割に内容が無いという悲惨な状況ですがご容赦ください。  

### 内容
今日の内容ですが、生成的確率モデリングに関して述べたいと思います。  
具体的には潜在変数（隠れ変数）を含む確率モデルの学習アルゴリズムについて簡単に解説します。  

### 想定する読者

対象は「__機械学習はじめてみよう！__」とか「__確率モデルで遊んでみたい！__」というビギナー向けのつもりです。  
難易度としては高校時代の同級生が「へー」って言ってくれるレベルを目指します。  
場合によってはあえて正確な表現を避け、簡単な言い回しを用いる予定です。  
なので玄人の方が読むと「__何を今更…__」とか「__それは違うだろ！__」となると思いますがご容赦ください。

________________________________________________________________________________
潜在変数モデル
--------------------------------------------------------------------------------
### 確率モデルと潜在変数
まず「__潜在変数モデル__」というものについて、簡単な実例を交えて簡単に解説します。  
（簡単を繰り替えしてごめんなさい。）  

確率モデルとはざっくり言えば、__同時分布__ です。  
__生成的__ 確率モデルとはざっくり言えば、
__データを生成できる（サンプリングできる）__ 確率モデルです。  
潜在変数モデルとはざっくり言えば、生成的確率モデルのうち、  
観測可能である __観測変数 x__ と観測不能である __潜在変数 z__ を持つものです。  
ここでモデルの同時分布を具体的に決めるために、モデルは __パラメータ θ__ を持つとします。  
そしてその同時分布を p(x, z | θ）と書くことにします。

生成的確率モデルの素敵なところは、データの生成過程を想像できるので、解釈しやすいところです。  
潜在変数モデルの素敵なところは、複雑なデータが大量にあっても、  
それらを潜在変数という不思議なものでなんとなく説明してくれるところです。  

潜在変数モデルの使い道としてもっともポピュラーなものがクラスタリングです。  
データがいっぱいあるけど、実はそれらは潜在的にグループに分かれているという解釈です。  
例えば、テレビの視聴履歴を集めたとします。  
すると潜在グループとして、__オタク__ or __リア充__ などが考えられます。  
データはその潜在グループに依存して生成されると考えます。  
例えばオタク野郎はアニメばっかり観てて、リア充様はドラマとかスポーツとか観てるんじゃねえの？(激怒)  

この潜在変数を推定するといろいろ便利なんです。  
例えば Amazon の購入履歴を考えます。  
クラスタリングの結果「こいつらは潜在的に萌豚えだな」と分かれば、
アニメBD売りたい放題とかなるわけです。  

では具体的にはどうやってその潜在変数を推定するのでしょうか。

### 学習と推論
潜在変数モデル p(x, z | θ) では x は既知で、 z と θ は未知とします。  
よって x から z, θ を推定する必要があります。

ここで推定ってなんやねん！となることがあります。  
ここでは単純に z と θ の具体的な値を知りたいという意味で「推定」と言っています。  
しかし個人的な感覚だと z を推定することを推論、  
θ を推定する事を学習と呼ぶ気もします。  
なぜ同じ未知の変数なのに扱いが違うのか。それはそれらが別のタイプの変数だからです。  
潜在変数 z は確率変数であるのに対して、θはパラメータであり、単なる変数なのです。  
前置きはまぁ置いといて推定してみましょう。

#### 準備
一気に両方推定するのは難しそうなので１つずつ推定しましょう。  

まず x, z が既知のとき、θを推定することを考えます。  
これはモデル中の全確率変数が観測可能であるときに、  そのパラメータを推定することに対応します。  
このようなデータを __完全データ__ と呼びます。  
よってもはや x,z と区別することなく、p(x | θ) と書いても良いのです。  
完全データからのパラメータの推定は実は中学生くらいでもできます。  
例えばコインを用意し、１００回振ります。そしたら表 70 回、裏 30 回でした。  
するとこのコインのパラメータは 表が出る確率 70 % , 裏が出る確率 30 % としたくなります。  
表が55回で裏が45回なら表55%,裏45%としたくなります。  
これは __最尤推定__ と呼ばれるれっきとした推定法なのです。  
具体的には p(x | θ) を最大化するように θ を決めるのです。  
要するに今起きたことは、『最も尤もらしい出来事』だったと考えるのです。  
最尤推定はまぁ素敵な性質がいっぱいあるのですがここでは割愛します。  

次に x, θ が既知のとき、z を推定することを考えます。  
これは同時分布 p(x, z | θ) が既知であるときに z を推定することに対応します。  
先のコインの例で言えば、コインの表裏の出る確率がわかっている時に、  
今から振るコインの面を推定することに対応します。  
例えば表がでる確率が 99 % で裏の確率が 1 % だったとします。  
次にどっちが出ると思う？と聞かれたらまぁ「表」と答えるでしょう。  
では表55%, 裏45%ではどうでしょう？まぁ一つ選べと言われれば「表」でしょう。  
つまり、p(x, z | θ) を最大化するように z を決めるのです。  
さっきの最尤推定と同じアイディアですね。

しかし、多くの人は 表55%,裏45% の状況で「表！」と叫ぶことを躊躇すると思います。  
それに対し、θを推定するときに表55%,裏45%とすることはあまり抵抗がなかったと思います。  
この違いは z が確率変数であることを知っているからなのです。  
これがさっき z と θ の推定を推論と学習と言いわけたい気持ちの答えであり、  
同時に今日のメインテーマでもあります。  

とりあえずこの問題は棚上げし、z,θ を同時に推定する方法を考えましょう。  
今、x,z から θ を、そして x,θ からzを推定する方法を知っています。  
では x から z,θ を推定するにはどうすればいいか。  
答えは簡単。交互に推定するのです！  


#### MM algorithm (Viterbi Training, K-means etc)

まず適当に θ を決めます。  
すると x,θ がわかるので z を p(x,z | θ) を最大化するように推定します。  
すると x,z がわかるので θ を p(x,z | θ) を最大化するように推定します。  
これを繰り返すといつか z,θ は変化しなくなり、その z,θ が推定結果となります。  

「え？これでいいの？」って感じですがいいのです。  
これは Viterbi Training と呼ばれたりする方法ですが、ここでは MM algorithm と呼びます。  
なぜなら z を p(x,z | θ) を最大化(Maximize)するように決め、  
次に θ を p(x,z | θ) を最大化(Maximize)するように決めるからです。  
k-means というクラスタリング手法はこの一例です。  

初心者の人でも「これならできそうだ！」と思ったでしょう。  
できます。暇な時にやってください。  

この手法では z を推定するとき思い切りを持って p(x, z | θ) を最大化するように決めました。  
しかし、表55%,裏45%のときに「表」と答えるのを躊躇する方は納得出来ないでしょう。  
そこで z をもっとソフトに推定する方法があります。  
それがかの有名な __EM algorithm__ なのです！  

#### EM algorithm

表55%,裏45%のときになんて答えればしっくりくるのか。  
きっとその回答の一つが「表の期待値0.55(どやぁ」です。  
麻雀でもパチンコでも競馬でも皆さん期待値で考えますよね？  
「期待値が1を超えないギャンブルをする奴はクズ」とかいう人もいますが彼らは夢を買っているのです（多分。  

つまり、z を推定するときに p(x, z | θ) を最大化(Maximize)するのでなく、  
期待値(Expectation)を取ればいいんじゃないの？  
それでいいんです。  

まず適当に θ を決めます。  
すると p(z |x, θ) がわかるので z の期待値を出します。  
すると x と z の期待値がわかるので、それを用いて θ を更新します。  
これを繰り返すと p(x | θ) がどんどん増加していきます。  
これが EM algorithm なのです。  

「なんでp(x | θ)が増えるんや！」という声が聞こえますが、Markdown で数式書きたくないので割愛します。  
コインの例で言えば、  
表回数 55, 裏回数 45のとき、表55%,裏45%と計算したのと同様、  
表期待値 0.55, 裏期待値 0.45 のとき、表55%,裏45%と計算できますよね？  

これで「表」と叫ぶことが不安な方々も学習できるようになりました。  
しかし次の問題も生まれます。  
「z の期待値考えるなら、θも期待値考えたくない？」  
MM, EM algorithm がああるなら ME, EE algorithm もあるよね？  
貪欲ですね〜。  
でもθは確率変数じゃないんです。  
ではどうするか。  
θも確率変数としちゃおうよ！  
これがかの有名な「ベイズ」の始まりである…。

#### ME algorithm と EE algorithm

__ベイズ的な__ 確率モデルとはモデルパラメータ θ に __事前分布 p(θ)__ を導入することで、  
θ も確率変数として扱うモデルのことです。  
これにより潜在変数モデルは p(x, z,  θ) = p(x, z | θ)p(θ) はとなります。  
(ベイズは奥が深いのでここでは詳細は述べません。)  
この拡張により θ も p(x, z | θ) を最大化(Maximize)するのではなく、  
p(θ | x) で期待値 (Expectation) を取れるようになります。  
これによって MM, EM algorithm と同様に ME, EE algorithm が構成できますね。  

しかし実際は、ベイズ的なモデルで p(z | x), p(θ | x) を計算することは難しく、  
それらの近似分布 q(z), q(θ) を変分法で求める __変分ベイズ法__ が使われます。  
これはとても説明がめんどくせいので興味のある人は自分で調べてください。  

### まとめ

登場した４つの学習法をまとめると、  

1. MM algorithm : z を最大化, θ も最大化
2. EM algorithm : z を期待値, θ は最大化
3. ME algorithm : z を最大化, θ は期待値
4. EE algorithm : z は期待値, θ も期待値

となります。  
ではどれがいいの？  
難しい質問です。  
答えは「場合による」です。  

一般的に、以下の特徴が挙げられます。  

* M* algorithm は収束が早い。
* M* algorithm は初期値依存性が強い。
* E*, *E algorithm は汎化性能が高い。
* E* algorithm は遅い。
* *E algorithm は近似誤差が含まれる。

上の特徴を踏まえ、データやモデルから適切に学習法を選ぶ必要があります。  
これらを適切に使い分けられるのが今流行りの __データサイエンティスト__ なんじゃないかなぁ（しらね。  

________________________________________________________________________________
実験
================================================================================
ここでは実際に紹介した４手法を利用してクラスタリングを行ってみます。

### データ
今回は本をクラスタリングしてみます。  
データは私が趣味で作っている[本推薦 AI bot][book_rec_ai] のデータを利用します。  
まず bot がフォローしているユーザの本棚中の出現頻度 [Top 200][titles] の本を持ってきます。  
次にユーザのうち、それらの本のうち大体半分くらいはもっているユーザ 65 人を取ってきます。  
最後にユーザと本の関係を行列で表現します。  
縦軸を本 i、横軸をユーザ j、行列の値 xij を  

- 0:本を持ってない
- 1:本を持ってる
- 2:本を持ってる、かつ、高評価

の３値で表現します。  
[この行列][relation]の各列をデータと思って本をクラスタリングしてみます。  
クラスタ数は 10 として各学習法を 10 回ずつ初期値を変えて実行し、  
もっとも尤度の高い結果をクラスタリング結果としました。  

### モデル
今回はクラスタリングに naive Bayes model を用います。  
本 i に関するデータは xi = {xi1,...,xi65} によって表現されます。  
この本 i のクラスを zi とすれば、naive bayes model は  
p(xi | zi, θ) = Πj p(xij | zi, θ) と分解することを許します。  
つまりクラスが与えられたとき、各属性は条件付き独立と仮定するモデルです。  
naive Baeys model に対する MM, EM, ME, EE の[実装][nbc]も一応載っけておきます。  
この実装は「遅い、長い、汚い」の三拍子そろっていますが、唯一の見どころは、  
z, θ の更新ルーチンに渡す関数を変えるだけで MM, EM, ME, EE が実現できるところでしょうか。  
ようするにこれら４つの手法は紙一重ということが伝えたかったのです（フヒヒ、サーセンw  

### 結果

各学習法のクラスタリング結果です。

1. [MM][mm]
2. [EM][em]
3. [ME][me]
4. [EE][ee]

ぶっちゃけおおきな違いはないですね。  
よくよく見ると、漫画の扱いで個性が出ています。  
MM が聖☆おにいさんを３月のライオンやよつばと！とくっつけているのに対して、  
EE は聖☆おにいさんとテルマエロマエをくっつけてますね。  
漫画好きの私としては EE の分け方に賛成です。  

まぁこのように微妙に結果変わるので適切にアルゴリズム選んでねって話です。

________________________________________________________________________________
最後に
================================================================================
長々と書きましたが、後半だれてきてるのが目に見えますね。すいません。  

この記事では潜在変数モデルの学習法がまぁいろいろあって、それらがどういう関係にあるかを述べました。  
みなさんいろんなモデルを作るのが好きなようですが、MMは実装しやすく、早いのでお勧めです。  

もし最後まで読んでくださった方がいれば、大変お疲れ様でした＆ありがとうございました。  
またどこかでお会いしましょう。  

[book_rec_ai]: https://twitter.com/book_rec_ai "book_rec_ai"
[relation]: https://github.com/masakazu-ishihata/advent2013/blob/master/relation.dat "relation.dat"
[titles]: https://github.com/masakazu-ishihata/advent2013/blob/master/titles.dat "titles.dat"
[nbc]: https://github.com/masakazu-ishihata/advent2013/blob/master/nbc.rb "nbc.rb"
[mm]: https://github.com/masakazu-ishihata/advent2013/blob/master/mm.txt "mm.txt"
[em]: https://github.com/masakazu-ishihata/advent2013/blob/master/em.txt "em.txt"
[me]: https://github.com/masakazu-ishihata/advent2013/blob/master/me.txt "me.txt"
[ee]: https://github.com/masakazu-ishihata/advent2013/blob/master/ee.txt "ee.txt"

<file_sep>#!/usr/bin/env ruby
# -*- coding: utf-8 -*-

require "optparse"
require './nbc.rb'

################################################################################
# default
################################################################################
@ifile = "relation.dat"  # 入力ファイル
@ofile = "result.dat"    # 出力ファイル
@m = 0                   # 学習法
@a = 1.1
@b = 1.1
@k = 5
@r = 1
@@mm, @@em, @@me, @@ee = [0, 1, 2, 3]

################################################################################
# Arguments
################################################################################
OptionParser.new { |opts|
  # options
  opts.on("-h","--help","Show this message") {
    puts opts
    exit
  }
  opts.on("-i [INPUT]"){ |f|
    @ifile = f
  }
  opts.on("-o [OUTPUT]"){ |f|
    @ofile = f
  }
  opts.on("-k [# clasters]"){ |f|
    @k = f.to_i
  }
  opts.on("-r [# restart]"){ |f|
    @r = f.to_i
  }
  opts.on("-a [alpha]"){ |f|
    @a = f.to_f
  }
  opts.on("-b [beta]"){ |f|
    @b = f.to_f
  }
  opts.on("-1", "--mm"){
    @m = @@mm
  }
  opts.on("-2", "--em"){
    @m = @@em
  }
  opts.on("-3", "--me"){
    @m = @@me
  }
  opts.on("-4", "--ee"){
    @m = @@ee
  }
  # parse
  opts.parse!(ARGV)
}

################################################################################
# load data
################################################################################
class MyData
  #### new ####
  attr_reader :d    # データの次元
  attr_reader :n    # n[d] = d次元の属性値の種類
  attr_reader :data # データ
  def initialize(_file)
    # データの読み取り
    # データ形式 : i a_i,1 ... a_i,d
    @raw_data = []
    open(_file).read.split("\n").each do |line|
      ary = line.split(" ")
      name = ary.shift # 名前
      attr = ary       # 属性
      @raw_data.push([name, attr])
    end

    # データの加工
    @N = @raw_data.size                     # データ数
    @d = @raw_data[0][1].size               # データの次元
    @ns = Array.new(@N)                     # 各データ点の名前
    @vs = Array.new(@d){|d| Hash.new(nil)}  # 各次元の値

    @data = []
    for i in 0..@N-1
      name, attr = @raw_data[i]
      # 名前
      @ns[i] = name

      # 属性
      a = Array.new(@d)
      for d in 0..@d-1
        v = attr[d]
        @vs[d][v] = @vs[d].size if @vs[d][v] == nil
        a[d] = @vs[d][v]
      end
      @data.push(a)
    end
    @n = Array.new(@d){|d| @vs[d].size}
  end

  #### clastering ####
  def clastering(_file, opt = {})
    opt = {:m => 0, :a => 1.1, :b => 1.1, :k => 5, :r => 1}.merge(opt)
    p opt
    # init
    m = NaiveBayesModel.new(opt[:k], @n)
    m.a = opt[:a]
    m.b = opt[:b]
    m.restart = opt[:r]
    @data.each{|a| m.add(a) }

    # learn
    m.rand_p
    case opt[:m]
    when @@mm then m.mm_learn
    when @@em then m.em_learn
    when @@me then m.me_learn
    when @@ee then m.ee_learn
    end

    # result
    f = open(_file, "w")
    m.result.sort{|a,b| b.size <=> a.size}.each do |c|
      f.puts "#{c.map{|i| @ns[i]}.join(" ")}"
    end
    f.close
  end
end

################################################################################
# main
################################################################################
d = MyData.new(@ifile)
d.clastering(@ofile, {:m => @m, :a => @a, :b => @b, :k => @k, :r => @r})
<file_sep>#!/usr/bin/env ruby
# -*- coding: utf-8 -*-

require './samplers.rb'

################################################################################
# 便利なやつら
################################################################################
#### Array まわり ####
class Array
  def product; inject(:*); end
  def sum; inject(:+); end
  def normalize
    s = self.sum
    s == 0.0 || s == 1.0 ? self.clone : map{|i| i / s.to_f}
  end
end

#### 各種関数 ####
def exp(x); Math.exp(x); end
def log(x); Math.log(x); end
def delta(x, y); x == y ? 1 : 0; end

# digamma 関数 (テイラー近似)
# copied from https://github.com/csauper/content-attitude
def dig(x)
  x = x + 6
  p = 1.0 / (x * x)
  return (((0.004166666666667*p-0.003968253986254)*p+
           0.008333333333333)*p-0.083333333333333)*p+
    Math.log(x) - 0.5 / x - 1.0 / (x-1) - 1.0 / (x-2) -
    1.0 / (x-3) - 1.0 / (x-4) - 1.0 / (x-5) - 1.0 / (x-6)
end

#### 収束判定用 : 平均二乗誤差 ####
def mse(_a, _b)
  e = 0
  for i in 0.._a.size-1
    e += (_a[i] - _b[i]) ** 2
  end
  e
end

################################################################################
# Naive Bayes Model
################################################################################
class NaiveBayesModel
  #### for each instacne / class / dimension / value ####
  def each_i;     for i in 0..@x.size-1; yield i; end; end
  def each_k;     for k in 0..@k-1;      yield k; end; end
  def each_d;     for d in 0..@d-1;      yield d; end; end
  def each_v(_d); for v in 0..@n[_d]-1;  yield v; end; end

  ########################################
  # 初期化
  ########################################
  #### new ####
  attr_reader   :d, :n, :k     # モデル定数
  attr_accessor :a, :b         # ハイパーパラメータ
  attr_reader   :pk, :pvgk     # パラメータ
  attr_reader   :x, :z         # 確率変数
  attr_accessor :restart, :iter, :radius # 学習設定
  def initialize(_k, _n)
    # モデル定数
    @d = _n.size  # 観測ベクトルの次元数
    @n = _n.clone # 各次元の取りうる値
    @k = _k       # クラスタ数

    # ハイパーパラメータ
    @a = 1.0
    @b = 1.0

    # パラメータ
    # @pk[k] = p(k)
    # @pvgk[k][d][v] = p(a_d = v | k)
    @pk   = Array.new(@k){|k| 1 / @k.to_f}
    @pvgk = Array.new(@k){|k| Array.new(@d){|d| Array.new(@n[d]){|v| 1/@n[d].to_f}}}

    # 確率変数
    @x = [] # 観測変数集合 [ a_i = [a_i,1,...,a_i,@d],... ]
    @z = [] # 潜在変数集合 [ k_i,...]

    # 学習用期待値 (もしくはカウント)
    @ek   = Array.new(@k){|k| @a-1}
    @evgk = Array.new(@k){|k| Array.new(@d){|d| Array.new(@n[d]){|v| @b-1}}}

    # 学習設定
    @restart = 1        # ランダムリスタート
    @iter    = 100      # 繰り返し数
    @radius  = 1e-5     # 収束半径
  end

  #### パラメータの初期化 ####
  # rand_p ランダムに初期化
  # max_p  最頻値に初期化
  def rand_p
    @pk = Dirichlet.new_simple(@k, @a).sample
    each_k{|k| each_d{|d| @pvgk[k][d] = Dirichlet.new_simple(@n[d], @b).sample } }
  end
  def max_p
    @pk = @ek.normalize
    each_k{|k| each_d{|d| @pvgk[k][d] = @evgk[k][d].normalize } }
  end

  #### 潜在変数の初期化 ####
  # rand_z ランダムに初期化
  # max_z  最頻値に初期化
  def rand_z; each_i{|i| @z[i] = Categorical.new(@pk).sample}; end
  def max_z;  each_i{|i| @z[i] = max_k(@x[i])}; end
  def max_k(_a)
    mk, mp = [nil, nil]
    each_k{|k| mk, mp = [k, pkga(k, _a)] if mk == nil || mp < pkga(k, _a) }
    mk
  end

  #### 期待値の初期化 ####
  # ハイパーパラメータを考慮した擬似カウント
  def init_e
    each_k{|k| @ek[k] = @a-1}
    each_k{|k| each_d{|d| each_v(d){|v| @evgk[k][d][v] = @b-1}}}
  end

  ########################################
  # データ生成、入力
  ########################################
  #### sample ####
  def sample; sample_with_class[1]; end
  def sample_with_class
    k = Categorical.new(@pk).sample
    a = Array.new(@d){|d| Categorical.new(@pvgk[k][d]).sample }
    [k, a]
  end
  #### add ####
  def add(_a); add_with_class(-1, _a); end
  def add_with_class(_k, _a); @z.push(_k); @x.push(_a); end

  ########################################
  # 確率値 (場合によっては期待値)
  ########################################
  # pka   : p(_k, _a)
  # pkga  : p(_k | _a)
  def pka(_k, _a);  @pk[_k] * pagk(_a, _k); end
  def pkga(_k, _a); pka(_k, _a) / pa(_a);  end

  # pa    : p(_a)
  # pagk  : p(_a | _k)
  def pa(_a);       Array.new(@k){|k| pka(k, _a) }.sum;              end
  def pagk(_a, _k); Array.new(@d){|d| @pvgk[_k][d][_a[d]] }.product; end

  # export/import parameters
  def export_par; [@pk, @pvgk].flatten; end
  def import_par(_ary)
    # load pk
    each_k do |k|
      @pk[k] = _ary.shift
    end
    # load pvgk
    each_k do |k|
      each_d do |d|
        each_v(d) do |v|
          @pvgk[k][d][v] = _ary.shift
        end
      end
    end
  end

  ########################################
  # 学習 : z, p をそれぞれ e, m で推定
  ########################################
  #### z の推定 ####
  # zm_step : z を最頻値で推定
  # ze_step : z を期待値で推定
  def zm_step; max_z; z_step( lambda{|i, k| delta(k, @z[i])} ); end
  def ze_step;        z_step( lambda{|i, k| pkga(k, @x[i])}  ); end
  def z_step(_f)
    init_e
    each_i do |i|
      each_k do |k|
        # ek
        @ek[k] += _f.call(i, k)
        # evgk
        each_d do |d|
          @evgk[k][d][ @x[i][d] ] += _f.call(i, k)
        end
      end
    end
  end

  #### p の推定 ####
  # pm_step : p を最頻値で推定
  # pe_step : p を期待値で推定
  def pm_step; p_step( lambda{|x| log(x)} ); end
  def pe_step; p_step( lambda{|x| dig(x)} ); end
  def p_step(_f)
    sum_ek = @ek.sum
    each_k do |k|
      # pk
      @pk[k] = exp( _f.call(@ek[k]) - _f.call(sum_ek) )
      # pvgk
      each_d do |d|
        sum_evgk = @evgk[k][d].sum
        each_v(d) do |v|
          @pvgk[k][d][v] = exp( _f.call(@evgk[k][d][v]) - _f.call(sum_evgk) )
        end
      end
    end
  end

  #### 学習本体 ####
  def mm_learn; puts "MM"; learn{zm_step; pm_step}; end
  def em_learn; puts "EM"; learn{ze_step; pm_step}; end
  def me_learn; puts "ME"; learn{zm_step; pe_step}; end
  def ee_learn; puts "EE"; learn{ze_step; pe_step}; end
  def learn
    best_p, best_l = [nil, nil]
    for r in 1..@restart
      puts "Try #{r}"

      # パラメータの初期化
      rand_p
      p0 = export_par

      # _n 回 | 収束するまで繰り返す
      for i in 1..@iter
        # learn
        t1 = Time.now
        yield
        t2 = Time.now

        # 収束判定
        p1 = export_par
        e = mse(p0, p1)
        printf "%5d %.5e %.5e\n", i, e, t2-t1
        break if e < @radius
        p0 = p1
      end

      # prediction
      max_p
      max_z

      # loglikelihood
      l = loglikelihood
      printf "LL = %e\n", l

      best_p, best_l = [p1, l] if best_l == nil || best_l < l
    end

    # best result
    import_par(best_p)
    max_z
  end

  #### log likelihood ####
  def loglikelihood
    ll = 0
    each_i{|i| ll += log( pa( @x[i] ) )}
    ll
  end
  def loglikelihood_with_class
    ll = 0
    each_i{|i| ll += log( pagk(@x[i], @z[i]) ) }
    ll
  end

  #### clastering result ####
  def result
    cs = Array.new(@k){|k| []}
    each_i{|i| cs[ @z[i] ].push(i) }
    cs
  end

  ########################################
  # show
  ########################################
  def show
    show_cons
    show_hpar
    show_par
    show_var
    show_exp
  end
  #### constants ####
  def show_cons
    puts "Model Constants"
    puts "D = #{@d}, N = #{@n}, K = #{@k}"
  end
  #### hyper parameters ####
  def show_hpar
    puts "Hyper Parameters"
    puts "a = #{@a}"
    each_d{|d| puts "b = #{@b}"}
  end
  #### parameters ####
  def show_par
    puts "Parameters"
    puts "p(k) = #{@pk}"
    each_k do |k|
      each_d do |d|
        puts "p(a_#{d} | k = #{k}) = #{@pvgk[k][d]}"
      end
    end
  end
  #### variables ####
  def show_var
    puts "Variables"
    each_i do |i|
      puts "#{@z[i]}, #{@x[i]}"
    end
  end
  #### expectations ####
  def show_exp
    puts "Expectations"
    puts "ek = #{@ek}"
    each_k do |k|
      each_d do |d|
        puts "ea[#{k}][#{d}] = #{@evgk[k][d]}"
      end
    end
  end
end
