【Ruby】 EMアルゴリズムでクラスタリング

今回はEMアルゴリズムでクラスタリングをやってみました。

ソースコード

以下に示すアルゴリズム部分の ema.rb と結果を描画する main.rb があります。 main.rb の方はいつものように (https://gist.github.com/seinosuke/7438fad6c92f25e2a8f4) においてあります。

クラス名悩みましたが、EMとするとRubyだとEventMachineっぽいなって思ったのでEMAとしました。

class EMA
  attr_reader :all_patterns, :k

  def initialize(options = {})
    @dimension = options[:dimension]
    @k = options[:k]
    orig_patterns = options[:patterns].map do |pattern|
      Matrix[pattern.map(&:to_f)]
    end
    @means = orig_patterns.sample(@k)

    @all_patterns = Array.new(@k) { [] }
    orig_patterns.each do |pattern|
      c = @means.map.with_index do |m, i|
        d = (m - pattern).inject(0.0) { |sum, v| sum + v*v }
        {:i => i, :d => d}
      end.min_by { |h| h[:d] }[:i]
      @all_patterns[c] << pattern
    end

    @means = @all_patterns.map do |patterns|
      patterns.inject(Matrix[Array.new(@dimension) { 0.0 }], :+)
      .map { |v| v / patterns.size.to_f }
    end

    @sigmas = @all_patterns.map.with_index do |patterns, i|
      patterns.inject(Matrix.zero(@dimension)) do |sum, x|
        sum + (x - @means[i]).t * (x - @means[i])
      end / patterns.size.to_f
    end
  end

  def update
    e_step
    m_step
  end

  # 発生している確率が最も高いクラスタへ振り分け直す
  def e_step
    patterns = @all_patterns.flatten
    @all_patterns = Array.new(@k) { [] }
    patterns.each do |pattern|
      c = @k.times.map do |i|
        {:i => i, :prob => gauss(i, pattern)}
      end.max_by { |h| h[:prob] }[:i]
      @all_patterns[c] << pattern
    end
  end

  # 更新されたクラスタ内で平均と共分散行列を計算し直す
  def m_step
    @means = @all_patterns.map.with_index do |patterns, i|
      patterns.inject(Matrix[Array.new(@dimension) { 0.0 }]) do |sum, x|
        sum + (x * prob(i, x))
      end / patterns.inject(0.0) { |sum, x| sum + prob(i, x) }
    end

    @sigmas = @all_patterns.map.with_index do |patterns, i|
      patterns.inject(Matrix.zero(@dimension)) do |sum, x|
        sum + (x - @means[i]).t * (x - @means[i])
      end / patterns.inject(0.0) { |sum, x| sum + prob(i, x) }
    end
  end

  # i番目の正規分布についてxにおける値を返す
  def gauss(i, x)
    Math.exp(-0.5 * ((x - @means[i]) * @sigmas[i].inv * (x - @means[i]).t)[0, 0]) /
      ( (Math.sqrt(2.0*Math::PI)**@dimension) * Math.sqrt(@sigmas[i].det) )
  end

  # xがi番目のクラスタから発生している確率 P(Wi|x)
  def prob(i, x)
    gauss(i, x) / @k.times.inject(0.0) { |sum, j| sum + gauss(j, x) }
  end

  # gnuplot用に2次元正規分布の文字列を返す
  # 特にアルゴリズムには関係ない
  def gauss_str(i)
    s1 = Math.sqrt(@sigmas[i][0, 0])
    s2 = Math.sqrt(@sigmas[i][1, 1])
    m1 = @means[i][0, 0]
    m2 = @means[i][0, 1]
    rho = @sigmas[i][0, 1] / (s1 * s2)
    "exp( -( ((x-#{m1})/#{s1})**2 - #{2.0*rho}*((x-#{m1})/#{s1})*((y-#{m2})/#{s2}) + ((y-#{m2})/#{s2})**2 )" <<
      " / " <<
      " #{2.0*(1.0 - rho**2)} )" <<
      " / " <<
      "#{( (Math.sqrt(2.0*Math::PI)**@dimension) * Math.sqrt(@sigmas[i].det) )}"
  end
end

処理の流れ

k-means法でセントロイドに最も近いパターンをそのクラスタに属させるという処理をしますが、そこを各クラスタが正規分布で近似できるとし、あるパターンについて各分布での値を計算して最も大きな値となる分布のクラスタに属しているとしています。

いくつのクラスタに分けるかという数字は初期化の際に与えます。この後紹介する結果では3つのクラスタに分ける例を2次元のデータセットに対して行います。

行列で書いたので一応何次元のデータセットに対してでもクラスタリングはできますが、処理の途中経過をgifで出力する部分は2次元のものにしか対応してません。

実行結果 その1

• クラスタリング前のデータ
• 実行中のgif
• クラスタリング後のデータ

という順でふたつのデータセットに対する実行結果を以下に示します。ちなみに、gnuplotのバージョンは4.6です。

 

図1 クラスタリング前のデータ

 

図2 実行中 

 

図3 クラスタリング後のデータ

実行結果 その2

 

 図4 クラスタリング前のデータ

 

 図5 実行中

 

図6 クラスタリング後のデータ
 

おわりに

うまく書けてるか自信が無かったのですが、なんとかクラスタリングできてるっぽいのでよかったです。何かおかしな点などあったらご指摘ください。

2016年初投稿でした。
今年もよろしくお願いします。

※追記　2016/09/08
改良版です。こちらもよろしくお願いします。
【Ruby】 VBEMアルゴリズム（PRML10章）