　今回は下記書籍の混合ベルヌーイ分布による手書き文字分類のコードを Ruby で実装してみたいと思います。

ITエンジニアのための機械学習理論入門
https://www.amazon.co.jp/IT-ebook/dp/B016Q22IX2/

　サンプルコードはこちらで公開されています。

一様乱数の複数生成

　Python版では numpy.random.rand メソッドを使って、一様乱数を複数生成しています。引数で生成数を指定できます。戻り値の型は numpy.ndarray になります。

　下記サイトも参考にさせていただきました。

>>> from numpy.random import randint, rand

>>> type(rand(28*28*K))
<type 'numpy.ndarray'> 

>>> rand(28*28*K)                                               
array([ 0.5395075 ,  0.98248117,  0.02652428, ...,  0.60618989, 
        0.55422211,  0.97170692])                    

>>> rand(28*28*K)*0.5+0.25                                     
array([ 0.27982902,  0.52512579,  0.6166757 , ...,  0.28828048,
        0.43710937,  0.70152255])

　Ruby版では複数の一様乱数を生成するようなメソッドは見つからなかったので、times と map で Random.rand の結果を Numo::NArray のオブジェクトとして返すようにしました。

instance method Random#rand (Ruby 2.4.0)

irb(main):002:0* require 'numo/narray'
=> true
irb(main):003:0> 
irb(main):004:0* K = 3
=> 3
irb(main):005:0> 
irb(main):006:0* random = Random.new
=> #<Random:0x007f5d8e5ba958>
irb(main):007:0> 
irb(main):008:0* (Numo::NArray[*(28*28*K).times.map { rand }] * 0.5 + 0.25).inspect
=> "Numo::DFloat#shape=[2352]\n[0.57563, 0.42476, 0.271377, 0.37791, 0.639573, 0.469103, 0.381991, ...]"
irb(main):009:0>

配列の次元の操作

　Python版では numpy.reshape を使って、一次元配列を二次元配列に変換しています。

numpy.reshape — NumPy v1.12 Manual

　下記サイトも参考にさせていただきました。

配列の次元数や大きさの操作 — 機械学習の Python との出会い

>>> rand(28*28*K)*0.5+0.25                                     
array([ 0.27982902,  0.52512579,  0.6166757 , ...,  0.28828048,
        0.43710937,  0.70152255])                              

>>> (rand(28*28*K)*0.5+0.25).reshape(K, 28*28)                  
array([[ 0.45648701,  0.52251086,  0.64335816, ...,  0.58305645,
         0.6682408 ,  0.43845499],                              
       [ 0.56658712,  0.35430961,  0.43453176, ...,  0.47584553,
         0.28571063,  0.47738208],                              
       [ 0.45287756,  0.36832049,  0.66361054, ...,  0.57060036,
         0.32932144,  0.43710831]])

　Ruby版では Numo::NArray に reshape メソッドがあるのでこちらを使用しました。

Class: Numo::NArray — Documentation by YARD 0.8.7.6

MF / Rubyの行列計算ライブラリ Numo::NArray 覚え書き

irb(main):008:0* (Numo::NArray[*(28*28*K).times.map { rand }] * 0.5 + 0.25).inspect
=> "Numo::DFloat#shape=[2352]\n[0.57563, 0.42476, 0.271377, 0.37791, 0.639573, 0.469103, 0.381991, ...]"
irb(main):009:0> 

irb(main):012:0* puts (Numo::NArray[*(28*28*K).times.map { rand }] * 0.5 + 0.25).reshape(K, 28 * 28).inspect         
Numo::DFloat#shape=[3,784]
[[0.34944, 0.607955, 0.691216, 0.556252, 0.291292, 0.466035, 0.680318, ...], 
 [0.707424, 0.483005, 0.541919, 0.708824, 0.287847, 0.266209, 0.354161, ...], 
 [0.627625, 0.688389, 0.510972, 0.673134, 0.579064, 0.355312, 0.513933, ...]]
=> nil

二次元配列の要素指定と合計値の計算

　Python版では ndarray の一次元目のインデックスを指定すると該当する配列が取得できます。そして sum() メソッドで合計値が取得できます。

>>> mu[0]                                                              
array([ 0.50779561,  0.40787245,  0.48671511,  0.65496467,  0.3637279 ,
        0.39219046,  0.53240974,  0.33735545,  0.63518755,  0.26163205,
        0.28325156,  0.4229063 ,  0.3944319 ,  0.66258652,  0.33099521,
        0.6460817 ,  0.58168761,  0.63767735,  0.59634922,  0.61818323,
〜〜〜 以下略 〜〜〜

>>> mu[0].sum()   
399.46258992135597

　Ruby版では Numo::NArray で同様に配列を取得するには一次元目だけでなく二次元目のインデックス範囲も指定する必要があります。そして合計値は sum メソッドで取得できます。

http://ruby-numo.github.io/narray/narray/Numo/DFloat.html#%5B%5D-instance_method

http://ruby-numo.github.io/narray/narray/Numo/DFloat.html#sum-instance_method

irb(main):030:0* mu[0,0..-1]
=> Numo::DFloat(view)#shape=[784]
[0.329509, 0.738836, 0.516475, 0.461981, 0.519925, 0.577925, 0.477174, ...]
irb(main):031:0> 
irb(main):032:0* mu[0,0..-1].sum
=> 392.4870691553733

画像データの表示

　Python版では配列データとして保持している画像データを画像として表示するために、 matplotlab の imshow を使用しています。一つの手書き文字データのピクセルデータの一次元配列を 28 x 28 の二次元配列に変換し、カラーマップにグレースケールを指定しています。また、元データはピクセルが黒だったら 1、白だったら 0 というデータなので、グレースケールとしては逆の指定になるため _r をつけた plt.cm.gray_r を指定しています。

subplot.imshow(mu[k].reshape(28,28), cmap=plt.cm.gray_r)

Image tutorial — Matplotlib 2.0.0 documentation

NumPyで画像処理 | mwSoft

pyplot — Matplotlib 2.0.0 documentation

　Ruby版では Nyaplot で簡単に画像を扱うことはできなそうだったので、 RMagick で画像データを生成し、 IRuby.display で画像を Jupyter Notebook 上に表示しました。imshow では渡したデータの上限値と下限値を自動的にスケールの上限と下限として扱ってくれるのですが、 import_pixel では自動的に変換はしてくれないので、最大値が 1.0 になるように map の中で変換して、 Float データであることを示すために、StorateType のパラメータには FloatPixel を指定しています。また、白黒を反転するために map の中では 1 から値を引いたものを使用するようにしています。

rate = 1.0 / mu[k, 0..-1].max
image = Magick::Image.new(28, 28)
image.import_pixels(0, 0, 28, 28, 'I', mu[k, 0..-1].map {|v| 1 - v * rate }.to_a, Magick::FloatPixel)
IRuby.display image

RMagick 2.12.0: class Image (instance methods e-o)

3階建ての2階角部屋２nd : ruby/tk + narray + rmagick でpixel操作(RGB)

ゼロ埋め配列データの生成

　Python版では numpy.zeros でパラメータで指定した要素数のゼロ埋めの配列を生成しています。

numpy.zeros — NumPy v1.12 Manual

mu = np.zeros((K, 28*28))

　Ruby 版では Numo::DFloat.zeros を使用しました。

Class: Numo::NArray — Documentation by YARD 0.8.7.6

mu = Numo::DFloat.zeros(K, 28 * 28)

配列データの最大値の要素のインデックス取得

　Python版では配列データの最大値を持つ要素のインデックスを取得するために、 numpy.argmax を使用しています。

for index, line in resp.iterrows():
    cls.append(np.argmax(line[0:]))

numpy.argmax — NumPy v1.12 Manual

NumPy配列の最大値/最小値に関する関数 | hydroculのメモ

　Ruby版では Numo::DFloat#max_index というメソッドがあるのでこちらを使用しました。

Class: Numo::DFloat — Documentation by YARD 0.8.7.6

cls = resp.map.with_index do |line, index|
  Numo::NArray[*line].max_index
end

スクリプト全体

　ここまでの内容を踏まえて、コード全体としては下記のように実装しました。

# 混合ベルヌーイ分布による手書き文字分類

require 'csv'
require 'numo/narray'
require 'rmagick'

K = 3  # 分類する文字数
N = 10 # 反復回数

# 分類結果の表示
def show_figure(df, mu, cls)
  K.times do |c|
    rate = 1.0 / mu[c, 0..-1].max
    image = Magick::Image.new(28, 28)
    image.import_pixels(0, 0, 28, 28, 'I', mu[c, 0..-1].map {|v| 1 - v * rate }.to_a, Magick::FloatPixel)
    puts "Master #{c}"
    IRuby.display image

    cnt = 0
    cls.each_with_index do |v, i|
      if v == c
        rate = 1.0 / df[i].max
        image = Magick::Image.new(28, 28)
        image.import_pixels(0, 0, 28, 28, 'I', df[i].map {|d| 1 - d * rate }, Magick::FloatPixel)
        IRuby.display image
        cnt += 1
      end
      
      break if cnt == 6
    end
  end
end

# ベルヌーイ分布
def bern(x, mu)
  r = 1.0
  x.zip(mu).each do |x_i, mu_i|
    if x_i == 1
      r *= mu_i
      next
    end

    r *= (1.0 - mu_i)
  end

  r
end

# トレーニングセットの読み込み
df = CSV.read('sample-images.txt').map {|line| line.map(&:to_i) }
data_num = df.size

# 初期パラメータの設定
mix = [1.0 / K] * K
random = Random.new
mu = (Numo::NArray[*(28 * 28 * K).times.map { rand }] * 0.5 + 0.25).reshape(K, 28 * 28)

K.times do |k|
  mu[k, 0..-1] /= mu[k, 0..-1].sum
end

puts "initial"
K.times do |k|
  rate = 1.0 / mu[k, 0..-1].max
  image = Magick::Image.new(28, 28)
  image.import_pixels(0, 0, 28, 28, 'I', mu[k, 0..-1].map {|v| 1 - v * rate }.to_a, Magick::FloatPixel)
  IRuby.display image
end

resp = nil
N.times do |iter_num|
  puts "iter_num #{iter_num}"

  # E Phase
  resp = []
  df.each_with_index do |line, index|
    tmp = []
    K.times do |k|
      a = mix[k] * bern(line, mu[k, 0..-1])
      if a == 0.0
        tmp << 0.0
      else
        s = 0.0
        K.times do |kk|
          s += mix[kk] * bern(line, mu[kk, 0..-1])
        end
        tmp << a / s
      end
    end
    resp << tmp
  end

  # M Phase
  mu = Numo::DFloat.zeros(K, 28 * 28)
  K.times do |k|
    nk = resp.transpose[k].inject(:+)
    mix[k] = nk / data_num
    df.each_with_index do |line, index|
      mu[k, 0..-1] = mu[k, 0..-1] + Numo::NArray[*line].cast_to(Numo::DFloat) * resp[index][k]
    end
    mu[k, 0..-1] /= nk

    rate = 1.0 / mu[k, 0..-1].max
    image = Magick::Image.new(28, 28)
    image.import_pixels(0, 0, 28, 28, 'I', mu[k, 0..-1].map {|v| 1 - v * rate }.to_a, Magick::FloatPixel)
    IRuby.display image
  end
end

# トレーニングセットの文字を分類
cls = resp.map.with_index do |line, index|
  Numo::NArray[*line].max_index
end

# 分類結果の表示
show_figure(df, mu, cls)

　これを Jupyter Notebook 上で実行すると、ランダムに作成された３つの画像生成器が表示され、それぞれがEMアルゴリズムで更新された画像が10回目まで順次表示されます。その後で、EMアルゴリズムで生成された画像生成器によってテストデータを分類した結果を６個ずつ表示します。Python 版で matplotlab でやっているように見やすく表形式で表示したかったのですが、 gnuplot 等調べてみたものの時間がかかりそうだったので、今回は分類自体はちゃんとできてそうだという確認までにしました。表示された画像を見やすく並べると下記のようになります。まずは画像生成器を更新していくところから。

f:id:akanuma-hiroaki:20170213224743j:plain