MNISTデータセットをNArray配列として扱う

　RubyからMNISTのデータセットを扱えるように、書籍「ゼロから作るDeepLearning」のコードをベースにMNISTのデータセットをNArray配列としてロードするコードを書いてみました。

www.oreilly.co.jp

コード全体

　まずはコード全体を掲載します。

require 'open-uri'
require 'zlib'
require 'fileutils'
require 'numo/narray'

URL_BASE = 'http://yann.lecun.com/exdb/mnist/'
KEY_FILES = {
  train_img:   'train-images-idx3-ubyte.gz',
  train_label: 'train-labels-idx1-ubyte.gz',
  test_img:    't10k-images-idx3-ubyte.gz',
  test_label:   't10k-labels-idx1-ubyte.gz'
}

DATASET_DIR = "#{File.dirname(__FILE__)}/dataset"
SAVE_FILE = "#{DATASET_DIR}/mnist.dump"

IMG_SIZE = 784

def download(file_name)
  puts "Downloading #{file_name} ..."
  open(URL_BASE + file_name) do |file|
    open("#{DATASET_DIR}/#{file_name}", 'w+b') do |out|
      out.write(file.read)
    end
  end
  puts "Done."
end

def download_mnist
  FileUtils.mkdir_p(DATASET_DIR)
  KEY_FILES.each do |k, file|
    download(file)
  end
end

def load_img(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = gz.each_byte.to_a[16..-1].each_slice(IMG_SIZE).to_a
    data = Numo::UInt8[*data]
  end
  puts "Done"
  data
end

def load_label(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = Numo::UInt8[*gz.each_byte.to_a[8..-1]]
  end
  puts "Done"
  data
end

def convert_narray
  dataset = {}
  dataset[:train_img] = load_img(KEY_FILES[:train_img])
  dataset[:train_label] = load_label(KEY_FILES[:train_label])
  dataset[:test_img] = load_img(KEY_FILES[:test_img])
  dataset[:test_label] = load_label(KEY_FILES[:test_label])
  dataset
end

def init_mnist
  download_mnist
  dataset = convert_narray
  puts "Creating dump file ..."
  File.write(SAVE_FILE, Marshal.dump(dataset))
  puts "Done!"
end

def change_one_hot_label(x)
  one_hot_arrays = x.to_a.map do |v|
    one_hot_array = Array.new(10, 0)
    one_hot_array[v] = 1
    one_hot_array
  end
  Numo::UInt8[*one_hot_arrays]
end

def load_mnist(normalize = true, flatten = true, one_hot_label = false)
  unless File.exist?(SAVE_FILE)
    init_mnist
  end

  dataset = Marshal.load(File.read(SAVE_FILE))

  if normalize
    %i(train_img test_img).each do |key|
      dataset[key] = dataset[key].cast_to(Numo::DFloat)
      dataset[key] /= 255.0
    end
  end

  if one_hot_label
    dataset[:train_label] = change_one_hot_label(dataset[:train_label])
    dataset[:test_label] = change_one_hot_label(dataset[:test_label])
  end

  unless flatten
    %i(train_img test_img).each do |key|
      dataset[key] = dataset[key].reshape(dataset[key].shape[0], 28, 28)
    end
  end

  [dataset[:train_img], dataset[:train_label], dataset[:test_img], dataset[:test_label]]
end

データのダウンロード

　初回実行時はMNISTデータがホスティングされている Yann LeCun’s MNIST からデータをダウンロードします。対象ファイルは下記のように定義しています。

URL_BASE = 'http://yann.lecun.com/exdb/mnist/'
KEY_FILES = {
  train_img:   'train-images-idx3-ubyte.gz',
  train_label: 'train-labels-idx1-ubyte.gz',
  test_img:    't10k-images-idx3-ubyte.gz',
  test_label:   't10k-labels-idx1-ubyte.gz'
}

　これをデータセット用のディレクトリにダウンロードします。

def download(file_name)
  puts "Downloading #{file_name} ..."
  open(URL_BASE + file_name) do |file|
    open("#{DATASET_DIR}/#{file_name}", 'w+b') do |out|
      out.write(file.read)
    end
  end
  puts "Done."
end

def download_mnist
  FileUtils.mkdir_p(DATASET_DIR)
  KEY_FILES.each do |k, file|
    download(file)
  end
end

画像データのロード

　ダウンロードしたMNISTデータセットの画像データをロードします。データ形式は一般的な画像フォーマットではなく、Yann LeCun’s MNISTに下記のように説明があります。

TRAINING SET IMAGE FILE (train-images-idx3-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000803(2051) magic number 
0004     32 bit integer  60000            number of images 
0008     32 bit integer  28               number of rows 
0012     32 bit integer  28               number of columns 
0016     unsigned byte   ??               pixel 
0017     unsigned byte   ??               pixel 
........ 
xxxx     unsigned byte   ??               pixel
Pixels are organized row-wise. Pixel values are 0 to 255. 0 means background (white), 255 means foreground (black).

TEST SET IMAGE FILE (t10k-images-idx3-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000803(2051) magic number 
0004     32 bit integer  10000            number of images 
0008     32 bit integer  28               number of rows 
0012     32 bit integer  28               number of columns 
0016     unsigned byte   ??               pixel 
0017     unsigned byte   ??               pixel 
........ 
xxxx     unsigned byte   ??               pixel

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

TensorFlow : MNIST データ・ダウンロード (コード解説) – TensorFlow

　上記内容と書籍のコードを元に、下記のように実装しました。

def load_img(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = gz.each_byte.to_a[16..-1].each_slice(IMG_SIZE).to_a
    data = Numo::UInt8[*data]
  end
  puts "Done"
  data
end

　ダウンロードしたgzファイルをopenし、byte単位で読み込みます。先頭16byteは画像以外の情報なので、16byte以降から末尾までを読み込み、画像サイズ（28 x 28 = 784）ごとの配列で分割して二次元配列にして、それをNArrayの配列データに変換します。

ラベルデータのロード

　ラベルデータのフォーマットは下記の通りです。

TRAINING SET LABEL FILE (train-labels-idx1-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000801(2049) magic number (MSB first) 
0004     32 bit integer  60000            number of items 
0008     unsigned byte   ??               label 
0009     unsigned byte   ??               label 
........ 
xxxx     unsigned byte   ??               label
The labels values are 0 to 9.

TEST SET LABEL FILE (t10k-labels-idx1-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000801(2049) magic number (MSB first) 
0004     32 bit integer  10000            number of items 
0008     unsigned byte   ??               label 
0009     unsigned byte   ??               label 
........ 
xxxx     unsigned byte   ??               label
The labels values are 0 to 9.

　実装は下記の通りです。基本的には画像データの場合と同じですが、画像データと違って正解ラベルの一次配列データなので、画像サイズで分割するような処理はありません。

def load_label(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = Numo::UInt8[*gz.each_byte.to_a[8..-1]]
  end
  puts "Done"
  data
end

MNISTデータのセーブとロード

　毎回MNISTデータをサイトからダウンロードしたらgzファイルからロードするのは時間がかかるので、初回実行時にデータをローカルに保持するようにします。今回はひとまず汎用性はあまり考慮せず、このスクリプトからのみ扱う前提で、 Marshal.dump したものをファイルに保存します。

File.write(SAVE_FILE, Marshal.dump(dataset))

　次回実行時はdumpファイルが存在すればダウンロードや変換処理はスキップして、dumpファイルをロードします。

dataset = Marshal.load(File.read(SAVE_FILE))

オプション

　基本的な処理はここまでですが、書籍の内容に沿って、3つのオプションを用意します。

　まず normalize オプションは画像データの各ピクセルデータを正規化するかどうかの指定で、元のデータは 0 - 255 の範囲の数値ですが、true を指定するとこれを 0.0 - 1.0 の値に正規化します。

if normalize
  %i(train_img test_img).each do |key|
    dataset[key] = dataset[key].cast_to(Numo::DFloat)
    dataset[key] /= 255.0
  end
end

　次に one_hot_label オプションは、正解ラベルデータを one_hot 表現の配列データとして返すかどうかの指定です。元データは画像データがどの数字なのかを表す一次元配列で、例えば画像データが 3, 5, 2 であれば、[3, 5, 2] という形ですが、 one_hot 表現だと、それぞれについて要素数10の配列を用意し、正解のインデックスのみ1で、それ以外は0の二次元配列になります。

[
 [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], # 3が正解なのでインデックス3の要素だけ1
 [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], # 5が正解なのでインデックス5の要素だけ1
 [0, 0, 1, 0, 0, 0, 0, 0, 0, 0]  # 2が正解なのでインデックス2の要素だけ1
]

def change_one_hot_label(x)
  one_hot_arrays = x.to_a.map do |v|
    one_hot_array = Array.new(10, 0)
    one_hot_array[v] = 1
    one_hot_array
  end
  Numo::UInt8[*one_hot_arrays]
end

if one_hot_label
  dataset[:train_label] = change_one_hot_label(dataset[:train_label])
  dataset[:test_label] = change_one_hot_label(dataset[:test_label])
end

　最後に flatten オプションは各画像のデータを 28 x 28 の二次元配列として返すか、784要素の一次元配列として返すかの指定になります。元データは784要素の一次元配列なので、flatten オプションに false が指定されれば 28 x 28 の二次元配列に変換します。

unless flatten
  %i(train_img test_img).each do |key|
    dataset[key] = dataset[key].reshape(dataset[key].shape[0], 28, 28)
  end
end

ロードしたMNISTデータを参照するサンプル

　上記コードによって読み込んだMNISTデータを参照してみます。上記のロード用コードは mnist.rb として保存しています。

require 'rmagick'
require './mnist.rb'

x_train, t_train, x_test, t_test = load_mnist

img = x_train[0, true]
label = t_train[0]
puts img.shape
puts img.max
puts img.min
puts label

image = Magick::Image.new(28, 28)
image.import_pixels(0, 0, 28, 28, 'I', img, Magick::FloatPixel)
image.display

　画像データとラベルデータからそれぞれ一つ目を取り出し、ラベルはターミナルに表示しています。画像データは ImageMagick(RMagick)で 28 x 28 の画像オブジェクトを作成した上で、 import_pixelsメソッドでインポートして、displayメソッドで別ウィンドウに表示します。先頭データは5なので、5の画像が表示されます。

　今回実装したコードは下記にも公開しています。

github.com

3層ニューラルネットワーク実装

　書籍「ゼロから作るDeepLearning」で紹介されている3層ニューラルネットワークをRubyで実装してみました。

www.oreilly.co.jp

サンプルコード全体

　まずは実装したコード全体を掲載します。

# 3層ニューラルネットワーク実装
# 各層のニューロン数は下記の通り
#  第0層（入力層）：2
#  第1層（隠れ層）：3
#  第2層（隠れ層）：2
#  第3層（出力層）：2

require 'numo/narray'

# 活性化関数としてシグモイド関数を使う
def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end

# 出力層の活性化関数として恒等関数を使う
def identity_function(x)
  x
end

# 重みとバイアスの初期化
def init_network
  network = {}
  network['w1'] = Numo::DFloat[[0.1, 0.3, 0.5], [0.2, 0.4, 0.6]]
  network['b1'] = Numo::DFloat[0.1, 0.2, 0.3]
  network['w2'] = Numo::DFloat[[0.1, 0.4], [0.2, 0.5], [0.3, 0.6]]
  network['b2'] = Numo::DFloat[0.1, 0.2]
  network['w3'] = Numo::DFloat[[0.1, 0.3], [0.2, 0.4]]
  network['b3'] = Numo::DFloat[0.1, 0.2]
  network
end

# 入力から出力までの処理
def forward(network, x)
  w1 = network['w1']
  w2 = network['w2']
  w3 = network['w3']
  b1 = network['b1']
  b2 = network['b2']
  b3 = network['b3']

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)
  a2 = z1.dot(w2) + b2
  z2 = sigmoid(a2)
  a3 = z2.dot(w3) + b3
  identity_function(a3)
end

network = init_network     # 重みとバイアスのデータを用意
x = Numo::DFloat[1.0, 0.5] # 入力層として2つのニューロンを用意
y = forward(network, x)    # 入力層、重み、バイアスのデータを渡して出力層のデータを取得
puts y.inspect             # 出力層のデータを表示

重みとバイアスのデータの用意

　まずは各層の重みとバイアスのデータを用意します。今回のニューロン数の構成としては下記のようになります。

• 第0層（入力層）：2
• 第1層（隠れ層1）：3
• 第2層（隠れ層2）：2
• 第3層（出力層）：2

　init_newtorkメソッドで上記構成の第1層〜第3層の各層への処理に使うデータを作成しています。w1, w2, w3 はそれぞれ第1層、第2層、第3層への重みの値で、b1, b2, b3 はバイアスの値になります。Numo::DFloatの配列として用意します。

def init_network
  network = {}
  network['w1'] = Numo::DFloat[[0.1, 0.3, 0.5], [0.2, 0.4, 0.6]]
  network['b1'] = Numo::DFloat[0.1, 0.2, 0.3]
  network['w2'] = Numo::DFloat[[0.1, 0.4], [0.2, 0.5], [0.3, 0.6]]
  network['b2'] = Numo::DFloat[0.1, 0.2]
  network['w3'] = Numo::DFloat[[0.1, 0.3], [0.2, 0.4]]
  network['b3'] = Numo::DFloat[0.1, 0.2]
  network
end

入力層のデータ

　入力層（第0層）のニューロン数は2なので、Numo::DFloatの配列を下記のように用意します。

x = Numo::DFloat[1.0, 0.5]

活性化関数

　単純パーセプトロンでは活性化関数としてステップ関数（ある閾値を界に出力が切り替わる関数）が使われますが、ニューラルネットワーク（多層パーセプトロン）では活性化関数としては入力に対して連続的に出力が変化する（グラフにすると滑らかな曲線になる）関数が使われるということで、今回はシグモイド関数を用いています。

def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end

　Ruby標準のMathモジュールのexpメソッドは引数に配列を渡すことができませんが、Numo::DFloatのMathモジュールのexpメソッドは配列に対応しているので、そちらを使用しています。

Module: Numo::DFloat::Math — Documentation by YARD 0.9.8

恒等関数

　出力層の活性化関数は解く問題の性質によって選択するようですが、今回はひとまずのサンプル実装ということで恒等関数を使用しています。入力値をそのまま返しています。

def identity_function(x)
  x
end

各層の伝達処理

　各層の伝達処理では、前の層のニューロンからの入力値に重みを掛け合わせ、バイアスを加算したものをシグモイド関数で処理します。例えば第0層（入力層）から第1層（隠れ層1）への伝達処理は下記のようになっています。

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)

　第0層からの入力値 x と、第1層へ伝達する際の重みの値 w1 はいずれもNumo::DFloatの配列として作成されているので、dotメソッドで内積を計算し、第1層へ伝達する際のバイアス値 b1 を加算しています。そしてそれをシグモイド関数で処理し、第1層の値としてz1を得ています。これを後続の層でも繰り返しますが、出力層の場合だけシグモイド関数の代わりに恒等関数を用います。

Class: Numo::NArray — Documentation by YARD 0.9.8

  a2 = z1.dot(w2) + b2
  z2 = sigmoid(a2)
  a3 = z2.dot(w3) + b3
  identity_function(a3)

スクリプト実行結果

　上記のスクリプトを実行すると下記のような出力になります。出力層のニューロンの値として、2つの要素を持つNumo::DFloatの配列が出力されます。

$ ruby three_layer_newral_network.rb 
Numo::DFloat#shape=[2]
[0.316827, 0.696279]

　コードは下記にも公開しました。

github.com

NumPyとNumo::NArray, Matrixでの演算の比較

　最近「ゼロから作るDeepLearning」を読み始めました。

www.oreilly.co.jp

　この本ではプログラミング言語としてはPythonを使用していて、配列や行列の演算にはNumPyが使われています。第1章ではNumPyでの基本的な演算について説明されているのですが、その内容をRubyのNumo::NArrayと、Ruby標準クラスのMatrixでの演算と比較してみました。NumPyでの内容は書籍記載の内容と同様に実行しています。Rubyのバージョンは 2.3.1p112 です。

使用準備

　NumPyは外部ライブラリなのでimportして使用します。

>>> import numpy as np

　Numo::NArrayを使用するにはまずgemをインストールします。

$ git clone git://github.com/ruby-numo/narray
$ cd narray
$ gem build numo-narray.gemspec
$ gem install numo-narray-0.9.0.3.gem

　そしてrequireでロードします。

irb(main):001:0> require 'numo/narray'
=> true

　Matrixは標準クラスですが、requireでロードする必要があります。

irb(main):001:0> require 'matrix'
=> true

配列の生成

　NumPyで配列を生成するには np.array() メソッドで引数に配列を渡すと numpy.ndarray クラスの配列が生成されます。

>>> x = np.array([1.0, 2.0, 3.0])
>>> x
array([ 1.,  2.,  3.])
>>> type(x)
<class 'numpy.ndarray'>

　Numo::NArrayでは下記のようにして配列を生成できます。

irb(main):002:0> x = Numo::NArray[1.0, 2.0, 3.0]
=> Numo::DFloat#shape=[3]
[1, 2, 3]

　データの型は自動的に判別してくれるようですが、下記のようにすることで明示的に型を指定して生成することも可能です。

irb(main):024:0* x = Numo::DFloat[1, 2, 3]
=> Numo::DFloat#shape=[3]
[1, 2, 3]

　Matrixの場合は下記のようになります。

irb(main):002:0> x = Matrix[[1.0, 2.0, 3.0]]
=> Matrix[[1.0, 2.0, 3.0]]

算術演算

　NumPyでの算術演算の例は下記の通りです。

>>> x = np.array([1.0, 2.0, 3.0])
>>> y = np.array([2.0, 4.0, 6.0])
>>> x + y
array([ 3.,  6.,  9.])
>>> x - y
array([-1., -2., -3.])
>>> x * y
array([  2.,   8.,  18.])
>>> x / y
array([ 0.5,  0.5,  0.5])
>>> x / 2.0
array([ 0.5,  1. ,  1.5])

　NumPyでは要素数が同じ配列での演算は各要素に対して行われ、配列と単一の値での計算の時にはブロードキャストの機能によって、単一の値と配列の各値との演算が行われます。

　次に Numo::NArray での例です。

irb(main):026:0* x = Numo::NArray[1.0, 2.0, 3.0]
=> Numo::DFloat#shape=[3]
[1, 2, 3]
irb(main):028:0> y = Numo::NArray[2.0, 4.0, 6.0]                                                                                                                                                                                              
=> Numo::DFloat#shape=[3]
[2, 4, 6]
irb(main):029:0> x + y
=> Numo::DFloat#shape=[3]
[3, 6, 9]
irb(main):030:0> x - y
=> Numo::DFloat#shape=[3]
[-1, -2, -3]
irb(main):031:0> x * y
=> Numo::DFloat#shape=[3]
[2, 8, 18]
irb(main):032:0> x / y
=> Numo::DFloat#shape=[3]
[0.5, 0.5, 0.5]
irb(main):033:0> x / 2.0
=> Numo::DFloat#shape=[3]
[0.5, 1, 1.5]

　こちらも NumPy と同様で、ブロードキャストも同じように行われています。

　Matrixでは下記のようになります。

irb(main):004:0* x = Matrix[[1.0, 2.0, 3.0]]
=> Matrix[[1.0, 2.0, 3.0]]
irb(main):005:0> y = Matrix[[2.0, 4.0, 6.0]]
=> Matrix[[2.0, 4.0, 6.0]]
irb(main):006:0> x + y
=> Matrix[[3.0, 6.0, 9.0]]
irb(main):007:0> x - y
=> Matrix[[-1.0, -2.0, -3.0]]
irb(main):008:0> x * y
ExceptionForMatrix::ErrDimensionMismatch: Matrix dimension mismatch
        from /usr/local/rbenv/versions/2.3.1/lib/ruby/2.3.0/matrix.rb:966:in `*'
        from (irb):8
        from /usr/local/rbenv/versions/2.3.1/bin/irb:11:in `<main>'
irb(main):009:0> x / y
ExceptionForMatrix::ErrDimensionMismatch: Matrix dimension mismatch
        from /usr/local/rbenv/versions/2.3.1/lib/ruby/2.3.0/matrix.rb:1062:in `inverse'
        from /usr/local/rbenv/versions/2.3.1/lib/ruby/2.3.0/matrix.rb:1049:in `/'
        from (irb):9
        from /usr/local/rbenv/versions/2.3.1/bin/irb:11:in `<main>'
irb(main):010:0> x / 2.0
=> Matrix[[0.5, 1.0, 1.5]]

　配列同士の乗算と除算では、ErrDimensionMismatch が発生してしまいます。Matrixの場合、行列同士の乗算は単純にそれぞれの要素同士を計算するのではなく、行列の積の演算になるようなので、次元数が合わずにエラーになっています。

instance method Matrix#* (Ruby 2.4.0)

　除算については逆行列に対しての乗算になるようで、逆行列を求めようとしているところで同様にエラーになっています。

instance method Matrix#/ (Ruby 2.4.0)

N次元配列

　NumPyでのN次元配列の生成と演算は下記のようになります。

>>> A = np.array([[1, 2], [3, 4]])
>>> A
array([[1, 2],
       [3, 4]])
>>> A.shape
(2, 2)
>>> A.dtype
dtype('int64')
>>> 
>>> B = np.array([[3, 0], [0, 6]])
>>> A + B
array([[ 4,  2],
       [ 3, 10]])
>>> A * B
array([[ 3,  0],
       [ 0, 24]])
>>> A * 10
array([[10, 20],
       [30, 40]])

　一次元配列の時と同様に、行列同士の演算、もしくは行列と単一の値での演算が行われます。行列の形状は shape メソッドで参照しています。

　次に Numo::NArray での例です。

irb(main):034:0> A = Numo::NArray[[1, 2], [3, 4]]
=> Numo::Int32#shape=[2,2]
[[1, 2], 
 [3, 4]]
irb(main):035:0> A.shape
=> [2, 2]
irb(main):036:0> 
irb(main):037:0* B = Numo::NArray[[3, 0], [0, 6]]
=> Numo::Int32#shape=[2,2]
[[3, 0], 
 [0, 6]]
irb(main):038:0> A + B
=> Numo::Int32#shape=[2,2]
[[4, 2], 
 [3, 10]]
irb(main):039:0> A * B
=> Numo::Int32#shape=[2,2]
[[3, 0], 
 [0, 24]]
irb(main):008:0* A * 10
=> Numo::Int32#shape=[2,2]
[[10, 20], 
 [30, 40]]

　NumPyと同様の演算が行われます。shapeメソッドも同様に用意されています。

　Matrixでは下記のようになります。

irb(main):016:0* A = Matrix[[1, 2], [3, 4]]
=> Matrix[[1, 2], [3, 4]]
irb(main):021:0> A.row_size
=> 2
irb(main):022:0> B = Matrix[[3, 0], [0, 6]]
=> Matrix[[3, 0], [0, 6]]
irb(main):023:0> A + B
=> Matrix[[4, 2], [3, 10]]
irb(main):024:0> A * B
=> Matrix[[3, 12], [9, 24]]
irb(main):025:0> A * 10
=> Matrix[[10, 20], [30, 40]]
irb(main):026:0> A.column_size
=> 2
irb(main):027:0> A.row_size
=> 2

　前述の通り、Matrixの乗算は行列の積になるので、NumPyとNumo::NArrayとは結果が異なっています。配列の形状を参照するためのメソッドは用意されていないようですが、 row_size, column_size メソッドで行と列の値をそれぞれ取ることはできるようになっています。

要素へのアクセス

　NumPyでの行列の要素へのアクセスの例は下記のようになります。

>>> X = np.array([[51, 55], [14, 19], [0, 4]])
>>> X
array([[51, 55],
       [14, 19],
       [ 0,  4]])
>>> X[0]
array([51, 55])
>>> X[0][1]
55
>>> for row in X:
...   print(row)
... 
[51 55]
[14 19]
[0 4]
>>> X = X.flatten()
>>> X
array([51, 55, 14, 19,  0,  4])
>>> X[np.array([0, 2, 4])]
array([51, 14,  0])
>>> X > 15
array([ True,  True, False,  True, False, False], dtype=bool)
>>> X[X > 15]
array([51, 55, 19])

　通常の多次元配列と同じようにインデックスを指定してアクセスすることができます。さらにインデックスとして配列を渡すことによって、複数の要素を指定することもできます。NumPyの配列に対して不等号などでの演算を行うと、Booleanの配列が生成されますので、それをインデックスとして渡すことで、Trueに対応する要素のみ抜き出すこともできます。

　次に Numo::NArray での例です。

irb(main):042:0* X = Numo::NArray[[51, 55], [14, 19], [0, 4]]
=> Numo::Int32#shape=[3,2]
[[51, 55], 
 [14, 19], 
 [0, 4]]
irb(main):043:0> X[0]
=> 51
irb(main):052:0> X[0..X.shape[1]-1]
=> Numo::Int32(view)#shape=[2]
[51, 55]
irb(main):055:0> X.to_a[0]
=> [51, 55]
irb(main):059:0> X.slice(0, 1)
=> Numo::Int32(view)#shape=[1,1]
[[55]]
irb(main):074:0* X.each do |row|
irb(main):075:1*   puts row.inspect
irb(main):076:1> end
51
55
14
19
0
4
=> Numo::Int32#shape=[3,2]
[[51, 55], 
 [14, 19], 
 [0, 4]]
irb(main):077:0> 
irb(main):078:0* X.to_a.each do |row|
irb(main):079:1*   puts row.inspect
irb(main):080:1> end
[51, 55]
[14, 19]
[0, 4]
=> [[51, 55], [14, 19], [0, 4]]
irb(main):088:0> X = X.flatten
(irb):88: warning: already initialized constant X
(irb):42: warning: previous definition of X was here
=> Numo::Int32(view)#shape=[6]
[51, 55, 14, 19, 0, 4]
irb(main):089:0> X
=> Numo::Int32(view)#shape=[6]
[51, 55, 14, 19, 0, 4]
irb(main):090:0> X[[0, 2, 4]]
=> Numo::Int32(view)#shape=[3]
[51, 14, 0]
irb(main):091:0> X > 15
=> Numo::Bit#shape=[6]
[1, 1, 0, 1, 0, 0]
irb(main):092:0> X[X > 15]
=> Numo::Int32(view)#shape=[3]
[51, 55, 19]

　Numo::NArrayの配列は、NumPyの配列とは違ってインデックスを指定すると全要素に対してのインデックス指定となります。なので、ある行の要素を指定したい時には、明確にどこからどこまでの要素かを指定するか、通常の配列に変換した上で行数を指定する必要があります。ある特定の要素を指定するには slice メソッドが利用できます。eachでの参照でも同様に各要素が参照されるので、行ごとに参照したい時には通常の配列に変換後にeachで参照するなどの工夫が必要そうです。不等号等での演算については、結果がBooleanではなく1か0での配列として返されますが、NumPyと同様に扱うことができます。

　最後にMatrixでの例です。

irb(main):052:0> X = Matrix[[51, 55], [14, 19], [0, 4]]                                                                                                                                                                                       
=> Matrix[[51, 55], [14, 19], [0, 4]]
irb(main):060:0> X.row(0)
=> Vector[51, 55]
irb(main):061:0> X.row(0)[1]
=> 55
irb(main):078:0> X.each_slice(X.column_size) do |row|                                                                                                                                                                                         
irb(main):079:1*   puts row.inspect
irb(main):080:1> end
[51, 55]
[14, 19]
[0, 4]
=> nil

　Matrixの行列に対して行を指定して参照するには row メソッドを使用します。単一の要素を指定する場合も、通常の多次元配列のような参照はできないので、row で行をVectorとして取得してからインデックスを指定します。eachでの参照については Numo::NArray と同様なので、通常の配列に変換してから each で参照するか、 each_slice で要素数を明示して参照するなどの工夫が必要そうです。また、flattenメソッドや、不等号での演算には対応していないようです。

　ここまで試した限りでは、Rubyで行列の演算をするにはひとまず Matrix よりも Numo::NArray の方が便利そうなので、書籍を読み進めつつ、 Numo::NArray 中心にRubyでの実装を試してみたいと思います。

ベイズ推定による回帰分析のコード

　今回は下記書籍のベイズ推定による回帰分析のコードを Ruby で実装してみたいと思います。今回で下記書籍のコードの置き換えは終了です。

www.amazon.co.jp

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

github.com

単位行列の生成

　Python版では numpy.identity メソッドで単位行列を生成しています。引数で指定した行数の正方行列になります。

>>> import numpy as np
>>> 
>>> np.identity(3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> 

numpy.identity — NumPy v1.12 Manual

こちらも参考にさせていただきました。
Pythonの数値計算ライブラリ NumPy入門 « Rest Term

　Ruby版では Numo::DFloat.eye メソッドを使用しました。

irb(main):039:0* require 'numo/narray'
=> true                               
irb(main):040:0>                      
irb(main):041:0* Numo::DFloat.eye(3)  
=> Numo::DFloat#shape=[3,3]           
[[1, 0, 0],                           
 [0, 1, 0],                           
 [0, 0, 1]]                           

Class: Numo::NArray — Documentation by YARD 0.9.8

対角線要素の取得

　Python版では numpy.diagonal メソッドで行列の対角線要素を配列として取得しています。

>>> np.array([[1,2,3],[4,5,6],[7,8,9]]).diagonal()
array([1, 5, 9])
>>> 

numpy.diagonal — NumPy v1.12 Manual

　Ruby版では Numo::NArray#diagonal メソッドを使用しました。

irb(main):079:0* Numo::NArray[[1,2,3],[4,5,6],[7,8,9]].diagonal
=> Numo::Int32(view)#shape=[3]                                 
[1, 5, 9]                                                      

Class: Numo::NArray — Documentation by YARD 0.9.8

行列の積

　Python版では numpy.dot メソッドで二次元配列の積を求めています。

>>> DataFrame([1,2,3])                               
   0                                                 
0  1                                                 
1  2                                                 
2  3                                                 
>>>                                                  
>>> DataFrame([1,2,3]).T                             
   0  1  2                                           
0  1  2  3                                           
>>>                                                  
>>> np.dot(DataFrame([1,2,3]), DataFrame([1,2,3]).T) 
array([[1, 2, 3],                                    
       [2, 4, 6],                                    
       [3, 6, 9]])                                   
>>>                                                  

　Ruby版ではMatrixを*でかけることで積を求めました。

irb(main):084:0* Daru::DataFrame.new(x: [1,2,3])
=> #<Daru::DataFrame(3x1)>
       x
   0   1
   1   2
   2   3
irb(main):085:0> 
irb(main):086:0* Daru::DataFrame.new(x: [1,2,3]).transpose
=> #<Daru::DataFrame(1x3)>
       0   1   2
   x   1   2   3
irb(main):087:0> 
irb(main):088:0* Daru::DataFrame.new(x: [1,2,3]).to_matrix * Daru::DataFrame.new(x: [1,2,3]).transpose.to_matrix
=> Matrix[[1, 2, 3], [2, 4, 6], [3, 6, 9]]

負の値の平方根

　Python版の下記コードにおいて、 numpy.sqrt メソッドで平方根を求めていますが、対象が負の値になることがあります。

deviation = np.sqrt(1.0/beta + np.dot(np.dot(phi_x0.T, s), phi_x0))

　numpy.sqrt メソッドは負の値に対して実行してもRuntimeWarningだけでエラーにはならず、nanを返します。

>>> np.sqrt(-1)                                              
__main__:1: RuntimeWarning: invalid value encountered in sqrt
nan                                                          
>>>                                                          

　Ruby版では Math.sqrt メソッドは負の値に対して実行するとエラーになってしまいます。

irb(main):002:0* Math.sqrt(-1)
Math::DomainError: Numerical argument is out of domain - "sqrt"
        from (irb):2:in `sqrt'
        from (irb):2
        from /usr/local/rbenv/versions/2.3.1/bin/irb:11:in `<main>'
irb(main):003:0> 

　なので対象が負の値の場合は明示的に Float::NAN を返すようにしました。

deviation = (1.0 / BETA + phi_x0.transpose.dot(Numo::NArray[*s.to_a]).dot(phi_x0)).map {|v| v < 0 ? Float::NAN : Math.sqrt(v) }

constant Float::NAN (Ruby 2.0.0)

多次元正規分布乱数の取得

　Python版では多項式のサンプル生成に numpy.random.multivariate_normal メソッドを使用しています。

ws_samples = DataFrame(multivariate_normal(mean,sigma,4))

　multivariate_normal は平均の配列と共分散行列を引数に取り、多次元正規分布乱数を生成してくれます。

　以前パーセプトロンの二項分類のコードをRuby版で実装した時にも出てきたメソッドで、その時は共分散を考慮しなくてもなんとかなるケースだったのですが、今回はそうもいかず、Rubyでの代替ができなかったので、Rubyのコードの中からPythonのコードを実行して結果を利用するようにしました。この書籍のコード置き換え最後にして（自分的には）禁じ手を使ってしまいました。。

  command = "python -c 'import numpy; print numpy.random.multivariate_normal(#{mean.to_a.inspect}, #{sigma.to_a.inspect}, 4).tolist()'"
  output, std_error, status = Open3.capture3(command)
  ws_samples = Daru::DataFrame.rows(eval(output))

　Open3.capture3 メソッドでは標準出力を返すので、 print で multivariate_normal の結果を出力するようにしています。また、戻り値は文字列なので、Rubyで配列として扱われるように eval しています。

　そして配列から Daru::DataFrame を生成するには Daru::DataFrame.rows メソッドを使用しました。

Method: Daru::DataFrame.rows — Documentation for daru (0.1.5)

スクリプト全体

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

# ベイズ推定による回帰分析

require 'daru'
require 'numo/narray'
require 'open3'
require 'nyaplot'

DATASET_NUMS = [4, 5, 10, 100]

BETA = 1.0 / (0.3) ** 2 # 真の分布の分散
ALPHA = 1.0 / 100 ** 2  # 事前分布の分散
ORDER = 9               # 多項式の字数

def normal_rand(mu = 0, sigma = 1.0)
  random = Random.new
  (Math.sqrt(-2 * Math.log(random.rand)) * Math.sin(2 * Math::PI * random.rand) * sigma) + mu
end

# データセット {x_n,y_n} (n=1...N) を用意
def create_dataset(num)
  dataset = Daru::DataFrame.new({'x': [], 'y': []})  

  num.times do |i|
    x = i.to_f / (num - 1).to_f
    y = Math.sin(2 * Math::PI * x) + normal_rand(0, 0.3)
    dataset.add_row(Daru::Vector.new([x, y], index: [:x, :y]))
  end
  
  return dataset
end

# 事後分布に基づく推定曲線、および、事後分布の平均と分散を計算
def resolve(dataset, m)
  t = dataset.y

  columns = {}
  (m+1).times do |i|
    columns["x**#{i}"] = dataset.x ** i
  end
  phis = Daru::DataFrame.new(columns)

  phiphi = nil
  phis.each_row_with_index do |line, index|
    phi = Daru::DataFrame.new(x: line)
    if index == 0
      phiphi = phi.to_matrix * phi.transpose.to_matrix
    else
      phiphi += phi.to_matrix * phi.transpose.to_matrix
    end
  end

  s_inv = Matrix[*(ALPHA * Numo::DFloat.eye(m + 1))] + BETA * phiphi
  s = s_inv.inv # 事後分布の共分散行列

  # 平均 m(x)
  mean_fun = lambda {|x0|
    phi_x0 = Numo::NArray[*((m + 1).times.map {|i| (x0 ** i).to_a })]
    tmp = 0
    phis.each_row_with_index do |line, index|
      if index == 0
        tmp = t[index] * Numo::NArray[*line.to_a]
        next
      end

      tmp += t[index] * Numo::NArray[*line.to_a]
    end
    BETA * phi_x0.transpose.dot(Numo::NArray[*s.to_a]).dot(tmp)
  }

  # 標準偏差 s(x)
  deviation_fun = lambda {|x0|
    phi_x0 = Numo::NArray[*((m + 1).times.map {|i| (x0 ** i).to_a })]
    deviation = (1.0 / BETA + phi_x0.transpose.dot(Numo::NArray[*s.to_a]).dot(phi_x0)).map {|v| v < 0 ? Float::NAN : Math.sqrt(v) }
    deviation.diagonal
  }

  tmp = nil
  phis.each_row_with_index do |line, index|
    if index == 0
      tmp = t[index] * Numo::NArray[*line.to_a]
      next
    end

    tmp += t[index] * Numo::NArray[*line.to_a]
  end
  mean = BETA * Numo::NArray[*s.to_a].dot(tmp).flatten

  return mean_fun, deviation_fun, mean, s
end

fig1 = Nyaplot::Frame.new
fig2 = Nyaplot::Frame.new

DATASET_NUMS.each do |num|
  train_set = create_dataset(num)
  mean_fun, deviation_fun, mean, sigma = resolve(train_set, ORDER)
  command = "python -c 'import numpy; print numpy.random.multivariate_normal(#{mean.to_a.inspect}, #{sigma.to_a.inspect}, 4).tolist()'"
  output, std_error, status = Open3.capture3(command)
  ws_samples = Daru::DataFrame.rows(eval(output))
  
  # トレーニングセットを表示
  plot1 = Nyaplot::Plot.new
  scatter1 = plot1.add(:scatter, train_set.x.to_a, train_set.y.to_a)
  scatter1.color('blue')
  scatter1.title('train_set')
  
  plot1.configure do
    x_label("N=#{num}")
    y_label('')
    xrange([-0.05, 1.05])
    yrange([-2, 2])
    legend(true)
    height(300)
    width(490)
  end
  
  plot2 = Nyaplot::Plot.new
  scatter2 = plot2.add(:scatter, train_set.x.to_a, train_set.y.to_a)
  scatter2.color('blue')
  scatter2.title('train_set')

  plot2.configure do
    x_label("N=#{num}")
    y_label('')
    xrange([-0.05, 1.05])
    yrange([-2, 2])
    legend(true)
    height(300)
    width(490)
  end
  
  linex = Numo::NArray[*(0..1).step(0.01).to_a]

  # 真の曲線を表示
  liney = (2 * Math::PI * linex).map {|v| Math.sin(v) }
  collect_line = plot1.add(:line, linex, liney)
  collect_line.color('green')
  collect_line.title('collect')

  # 平均と標準偏差の曲線を表示
  m = mean_fun.call(linex)
  d = deviation_fun.call(linex)
  mean_line = plot1.add(:line, linex, m)
  mean_line.color('red')
  mean_line.title('mean')
  lower_std_line = plot1.add(:line, linex, m - d)
  lower_std_line.color('black')
  lower_std_line.title('')
  upper_std_line = plot1.add(:line, linex, m + d)
  upper_std_line.color('black')
  upper_std_line.title('')

  # 多項式のサンプルを表示
  liney = m
  mean_line = plot2.add(:line, linex, liney)
  mean_line.color('red')
  mean_line.title('mean')
  
  f = lambda {|x, ws|
    y = 0
    ws.each_with_index do |w, i|
      y += w * (x ** i.to_i)
    end
    y
  }

  ws_samples.each_row do |ws|
    liney = f.call(linex, ws)
    sample_line = plot2.add(:line, linex, liney)
    sample_line.color('pink')
    sample_line.title('sample')
  end
  
  fig1.add(plot1)
  fig2.add(plot2)
end

fig1.show
fig2.show

　これを Jupyter Notebook 上で実行すると下記のようなグラフが表示されます。

　コードは下記にも公開してあります。

github.com

ベイズ推定による正規分布の推定のコード

　今回は下記書籍のベイズ推定による正規分布の推定のコードを Ruby で実装してみたいと思います。

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

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

github.com

スクリプト全体

　今回は使っているメソッド等は特に新しい要素はなかったので、スクリプト全体だけ掲載します。

# ベイズ推定による正規分布の推定

require 'rubystats/normal_distribution'
require 'nyaplot'

class Array
  def mean
    self.inject(:+) / self.size.to_f
  end
end

def normal_rand(mu = 0, sigma = 1.0)
  random = Random.new
  (Math.sqrt(-2 * Math.log(random.rand)) * Math.sin(2 * Math::PI * random.rand) * sigma) + mu
end

DATA_COUNTS = [2, 4, 10, 100]

# 真の分布
MU_TRUE = 2.0
BETA_TRUE = 1.0

# 事前分布
MU_0 = -2.0
BETA_0 = 1.0

ds = 100.times.map do
  normal_rand(MU_TRUE, 1.0 / BETA_TRUE)
end

fig1 = Nyaplot::Frame.new
fig2 = Nyaplot::Frame.new

DATA_COUNTS.each_with_index do |n, c|
  train_set = ds[0..n-1]
  mu_ML = train_set.mean
  mu_N = (BETA_TRUE * mu_ML + BETA_0 * MU_0 / n) / (BETA_TRUE + BETA_0 / n)
  beta_N = BETA_0 + n * BETA_TRUE

# 平均μの推定結果を表示
  plot = Nyaplot::Plot.new
  linex = (-10..10).step(0.01).to_a

  # 平均μの確率分布
  sigma = 1.0 / beta_N
  mu_est = Rubystats::NormalDistribution.new(mu_N, Math.sqrt(sigma))
  label = "mu_N=%.2f var=%.2f" % [mu_N, sigma]
  liney = linex.map {|x| mu_est.pdf(x) }
  line = plot.add(:line, linex, liney)
  line.title(label)
  line.color('red')

  # トレーニングセットを表示
  scatter = plot.add(:scatter, train_set, [0.2] * n)
  scatter.title('train_set')
  scatter.color('blue')

  plot.configure do
    x_label("N=#{n}")
    y_label('')
    xrange([-5, 5])
    yrange([0, liney.max * 1.1])
    legend(true)
    height(300)
    width(490)
  end

  fig1.add(plot)

# 次に得られるデータの推定分布を表示
  plot = Nyaplot::Plot.new
  
  # 真の分布を表示
  orig = Rubystats::NormalDistribution.new(MU_TRUE, Math.sqrt(1.0 / BETA_TRUE))
  liney = linex.map {|x| orig.pdf(x)}
  y_max = liney.max * 1.1
  line = plot.add(:line, linex, liney)
  line.title('original')
  line.color('green')
  
  # 推定分布を表示
  sigma = 1.0 / BETA_TRUE + 1.0 / beta_N
  mu_est = Rubystats::NormalDistribution.new(mu_N, Math.sqrt(sigma))
  label = "mu_N=%.2f var=%.2f" % [mu_N, sigma]
  liney = linex.map {|x| mu_est.pdf(x) }
  line = plot.add(:line, linex, liney)
  line.title(label)
  line.color('red')
  
  # トレーニングセットを表示
  scatter = plot.add(:scatter, train_set, train_set.map {|x| orig.pdf(x) })
  scatter.title('train_set')
  scatter.color('blue')

  plot.configure do
    x_label("N=#{n}")
    y_label('')
    xrange([-5, 5])
    yrange([0, y_max])
    legend(true)
    height(300)
    width(490)
  end

  fig2.add(plot)
end

fig1.show
fig2.show

　これを Jupyter Notebook 上で実行すると下記のようなグラフが表示されます。

　コードは下記にも公開してあります。

github.com

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

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

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

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

github.com

一様乱数の複数生成

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

numpy.random.rand — NumPy v1.12 Manual

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

qiita.com

>>> 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 等調べてみたものの時間がかかりそうだったので、今回は分類自体はちゃんとできてそうだという確認までにしました。表示された画像を見やすく並べると下記のようになります。まずは画像生成器を更新していくところから。

　次に画像生成器によって分類された結果です。一番左の画像は画像生成器のデータで、分類されたデータがその右に並んでいます。

　実行する度に結果は異なり、誤判定されるデータもありますが、元のPython版と比較して大体同じ程度の精度にはなっているようなので、ひとまずは問題なさそうかと思います。

　コードは下記に公開しました。

github.com

　ちなみに上記スクリプトを実行する前にトレーニングセットを絞り込むための準備用のスクリプトがあるのですが、こちらも一応Rubyで実装して下記に置いてあります。

ml4se_study/07-prep_data.rb at master · h-akanuma/ml4se_study · GitHub

k平均法による画像の減色処理のコード

　今回は下記書籍のk平均法による画像の減色処理のコードを Ruby で実装してみたいと思います。

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

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

github.com

画像ファイルの読み込み

　Python版ではPILライブラリのImageクラスを使って画像ファイルを読み込み、RGBデータとして取得して配列にしています。

Python Imaging Library (PIL)

The Python Imaging Library Handbook

from PIL import Image

im = Image.open("photo.jpg")
pixels = list(im.convert('RGB').getdata())

　Ruby版では ImageMagick をインストールし、RMagick を使って画像を扱いました。

github.com

RMagick 2.12.0 User's Guide and Reference

https://rmagick.github.io/image2.html#export_pixels

　下記サイトを参考にさせていただき、 export_pixels メソッドで取得した 16bit のRGBデータを 8bit に変換しました。

d.hatena.ne.jp

irb(main):001:0> require 'rmagick'
=> true
irb(main):002:0> img = Magick::ImageList.new('./photo.jpg')
=> [./photo.jpg JPEG 520x346 520x346+0+0 DirectClass 8-bit 46kb]
scene=0
irb(main):003:0> img.export_pixels.map {|p| p / 257 }.each_slice(3).to_a

乱数の取得

　Python版ではランダムな整数を取得するために numpy.random.randint を使っています。

from numpy.random import randint

randint(256)

numpy.random.randint — NumPy v1.12 Manual

　このメソッドは引数に渡した整数未満のランダムな整数を返すと言うことで、引数もその範囲に含む random.randint とは動きが異なるようです。

qiita.com

　Ruby版では Random#rand メソッドが同様の挙動になります。

random = Random.new

random.rand(256)

class Random (Ruby 2.4.0)

Rubyでランダムな数値を得る方法 -- ぺけみさお

ゼロ徐算によるエラーの回避

　Python版ではk平均法によって新しい代表点を求める際に、下記のような計算を行なっています。

for i in range(k):
    center_new[i] = center_new[i] / num_points[i]

　このとき、特にクラスタ数が多いケースだと、あるクラスタに属する点がないときに num_points[i] が 0 になり、ゼロ徐算が発生します。 center_new[i] は ndarray のオブジェクトですが、通常の数値計算のゼロ徐算だとエラーになるのに対し、ndarrayに対するゼロ徐算はWarningで、エラーにはならないようです。

>>> np.array([100, 100, 100]) / 0                               
__main__:1: RuntimeWarning: divide by zero encountered in divide
array([0, 0, 0])                                                

　Ruby では Numo::NArray に対するゼロ徐算はエラーになります。

irb(main):085:0* Numo::NArray[100,100,100] / 0
ZeroDivisionError: error in NArray operation
        from (irb):85:in `/'
        from (irb):85
        from /usr/local/rbenv/versions/2.3.1/bin/irb:11:in `<main>'

　このエラーを回避するために、 num_points[i] が 0 の場合は、 Numo::NArray[0, 0, 0] を返すようにしました。

k.times do |i|
  if num_points[i] == 0
    center_new[i] = Numo::NArray[0, 0, 0]
    next
  end

  center_new[i] = center_new[i] / num_points[i]
end

画像ファイルの変更、保存

　Python版では読み込み時と同様にPILを使用し、元の画像にputdataで計算結果を反映し、ビットマップ形式の別ファイルとして保存しています。

im.putdata(result) # Update image
im.save("output%02d.bmp" % k, "BMP")

　Ruby版では Magick::Image#import_pixels で計算結果を反映して、 Magick::Image#write で別ファイルとして保存しました。読み込みの時に export_pixels で取得した配列データを 8bit に変換してRGB区切りの配列データにしたので、それとは逆に、一次元配列に戻して16bitに変換したものを渡しています。

img.import_pixels(0, 0, img.columns, img.rows, "RGB", result.flatten.map {|p| p * 257 })
img.write("bmp:output%02d.bmp" % k)

https://rmagick.github.io/image2.html#import_pixels

https://rmagick.github.io/image3.html#write

スクリプト全体

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

# k平均法による画像の減色処理

require 'rmagick'
require 'numo/narray'

Colors = [2, 3, 5, 16]  # 減色後の色数

# k平均法による減色処理
def run_kmeans(pixels, k)
  cls = Array.new(pixels.size).fill(0)

  # 代表色の初期値をランダムに設定
  random = Random.new
  center = k.times.map { Numo::NArray[random.rand(256), random.rand(256), random.rand(256)] }
  puts "Initial centers: #{center.map(&:to_a)}"
  puts "========================"
  distortion = 0.0

  # 最大50回のIterationを実施
  50.times do |iter_num|
    center_new = Array.new(k).fill(Numo::NArray[0, 0, 0])
    num_points = Array.new(k).fill(0)
    distortion_new = 0.0

    # E Phase: 各データが属するグループ（代表色）を計算
    pixels.each_with_index do |point, pix|
      min_dist = 256 * 256 * 3
      point = Numo::NArray[*point]
      k.times do |i|
        d = (point - center[i]).cast_to(Numo::Int64).to_a.inject(0) {|sum, x| sum + (x ** 2) }
        if d < min_dist
          min_dist = d
          cls[pix] = i
        end

        center_new[cls[pix]] += point
        num_points[cls[pix]] += 1
        distortion_new += min_dist
      end
    end

    # M Phase: 新しい代表色を計算
    k.times do |i|
      if num_points[i] == 0
        center_new[i] = Numo::NArray[0, 0, 0]
        next
      end

      center_new[i] = center_new[i] / num_points[i]
    end
    center = center_new
    puts center.map(&:to_a).inspect
    puts "Distortion: J=#{distortion_new}"

    # Distortion(J)の変化が0.1%未満になったら終了
    if iter_num > 0 && distortion - distortion_new < distortion * 0.001
      break
    end
    distortion = distortion_new
  end

  # 画像データの各ピクセルを代表色で置き換え
  pixels.each_with_index do |point, pix|
    pixels[pix] = center[cls[pix]].to_a
  end

  pixels
end

Colors.each do |k|
  puts ""
  puts "========================"
  puts "Number of clusters: K=#{k}"

  # 画像ファイルの読み込み
  img = Magick::ImageList.new('./photo.jpg')
  pixels = img.export_pixels.map {|p| p / 257 }.each_slice(3).to_a

  # k平均法による減色処理
  result = run_kmeans(pixels, k)
  img.import_pixels(0, 0, img.columns, img.rows, "RGB", result.flatten.map {|p| p * 257 })
  img.write("bmp:output%02d.bmp" % k)
end

　これを実行するとターミナル上では下記のように、初期のランダムな代表点と、計算ごとに変化していく代表点と二乗歪み（J）の値が出力されます。

$ ruby 06-k_means.rb 

========================
Number of clusters: K=2
Initial centers: [[73, 209, 26], [101, 193, 89]]
========================
[[147, 136, 114], [177, 155, 141]]
Distortion: J=11686169170.0
[[140, 126, 102], [231, 210, 215]]
Distortion: J=5777512344.0
[[143, 126, 103], [237, 225, 228]]
Distortion: J=4710119998.0
[[143, 126, 104], [238, 228, 231]]
Distortion: J=4670715192.0
[[143, 126, 104], [238, 228, 231]]
Distortion: J=4668708398.0

========================
Number of clusters: K=3
Initial centers: [[214, 192, 131], [132, 203, 96], [118, 69, 71]]
========================
[[193, 169, 164], [99, 112, 65], [104, 83, 52]]
Distortion: J=8555090758.0
[[192, 173, 165], [102, 106, 61], [105, 55, 41]]
Distortion: J=6458497492.0
[[192, 173, 165], [101, 104, 60], [110, 46, 40]]
Distortion: J=6391712615.0
[[192, 173, 165], [100, 103, 59], [120, 41, 42]]
Distortion: J=6384285256.0
[[192, 173, 165], [96, 101, 56], [153, 33, 52]]
Distortion: J=6369689981.0
[[191, 173, 165], [94, 100, 55], [182, 35, 64]]
Distortion: J=6295123108.0
[[191, 174, 165], [93, 99, 54], [185, 38, 68]]
Distortion: J=6267370834.0
[[191, 174, 165], [93, 99, 54], [186, 40, 70]]
Distortion: J=6265803130.0

========================
Number of clusters: K=5
Initial centers: [[194, 128, 53], [20, 175, 19], [143, 69, 162], [10, 6, 13], [232, 145, 213]]
========================
[[164, 135, 117], [67, 100, 43], [219, 196, 202], [49, 63, 25], [235, 218, 223]]
Distortion: J=17277755636.0
[[174, 141, 131], [84, 99, 48], [235, 220, 224], [71, 50, 25], [245, 243, 244]]
Distortion: J=7740485201.0
[[179, 148, 140], [91, 101, 52], [237, 227, 230], [89, 45, 29], [250, 249, 249]]
Distortion: J=7006558605.0
[[182, 151, 144], [94, 102, 54], [238, 230, 233], [103, 42, 34], [251, 251, 251]]
Distortion: J=6813616639.0
[[183, 153, 146], [94, 102, 54], [239, 232, 234], [118, 38, 38], [252, 252, 252]]
Distortion: J=6752640626.0
[[183, 154, 147], [90, 99, 52], [239, 232, 235], [164, 29, 51], [252, 252, 252]]
Distortion: J=6699438034.0
[[182, 155, 147], [89, 98, 52], [239, 233, 235], [181, 32, 61], [252, 252, 252]]
Distortion: J=6536049096.0
[[182, 156, 147], [89, 98, 51], [239, 233, 235], [183, 35, 64], [252, 252, 252]]
Distortion: J=6513101512.0
[[182, 156, 147], [89, 98, 51], [239, 233, 235], [184, 36, 65], [252, 252, 252]]
Distortion: J=6510842397.0

========================
Number of clusters: K=16
Initial centers: [[240, 182, 207], [150, 87, 15], [121, 236, 188], [140, 204, 22], [66, 131, 79], [157, 120, 116], [17, 196, 63], [164, 6, 24], [34, 40, 245], [45, 113, 4], [229, 65, 229], [117, 158, 111], [114, 36, 146], [51, 62, 230], [
48, 35, 252], [14, 78, 55]]
========================
[[225, 210, 211], [118, 82, 56], [151, 171, 141], [124, 151, 86], [82, 111, 55], [176, 118, 115], [0, 0, 0], [176, 22, 51], [0, 0, 0], [49, 76, 26], [234, 96, 155], [125, 150, 99], [131, 57, 81], [0, 0, 0], [0, 0, 0], [38, 63, 27]]
Distortion: J=17994009460.0
[[226, 215, 214], [117, 78, 55], [169, 157, 139], [116, 141, 83], [78, 107, 48], [201, 99, 122], [24, 49, 15], [180, 25, 54], [0, 0, 0], [46, 72, 26], [224, 103, 146], [121, 147, 105], [153, 62, 72], [0, 0, 0], [0, 0, 0], [32, 57, 20]]
Distortion: J=12774225380.0
[[226, 216, 215], [115, 77, 53], [176, 153, 141], [112, 137, 80], [75, 105, 46], [210, 84, 119], [31, 57, 19], [179, 23, 52], [9, 21, 8], [49, 74, 28], [226, 116, 157], [125, 150, 107], [163, 67, 73], [0, 0, 0], [0, 0, 0], [33, 59, 20]]
Distortion: J=12102951219.0
[[226, 217, 216], [114, 77, 53], [181, 150, 142], [110, 136, 78], [74, 103, 44], [211, 74, 114], [34, 59, 20], [177, 22, 49], [14, 35, 10], [51, 76, 29], [226, 129, 164], [129, 152, 106], [168, 70, 77], [0, 0, 0], [0, 0, 0], [36, 62, 21]]
Distortion: J=11926470826.0
[[226, 217, 216], [113, 78, 53], [184, 148, 143], [109, 135, 77], [73, 102, 43], [210, 68, 109], [35, 61, 21], [174, 21, 47], [18, 41, 12], [52, 78, 30], [227, 135, 168], [131, 152, 106], [170, 73, 79], [0, 0, 0], [0, 0, 0], [39, 65, 21]]
Distortion: J=11849083461.0
[[226, 217, 216], [112, 79, 53], [186, 146, 144], [108, 135, 77], [72, 101, 43], [209, 64, 105], [37, 62, 21], [172, 20, 45], [20, 45, 13], [53, 79, 31], [227, 137, 170], [132, 152, 107], [169, 76, 80], [9, 18, 5], [0, 0, 0], [42, 67, 21]
]
Distortion: J=11816838760.0
[[226, 217, 216], [111, 80, 53], [188, 145, 144], [108, 135, 77], [72, 101, 42], [208, 60, 102], [38, 63, 22], [171, 19, 44], [22, 47, 14], [54, 79, 32], [227, 138, 171], [133, 152, 107], [168, 80, 81], [10, 27, 9], [0, 0, 0], [44, 69, 21
]]
Distortion: J=11800409417.0
[[226, 217, 216], [110, 81, 53], [189, 143, 144], [108, 135, 77], [71, 101, 42], [207, 58, 99], [39, 64, 22], [170, 19, 43], [23, 48, 14], [54, 79, 32], [228, 138, 171], [134, 153, 108], [167, 84, 82], [12, 33, 9], [0, 0, 0], [46, 71, 20]
]
Distortion: J=11789721844.0

　そして下記の元画像に対して、代表色を2色、3色、5色、16色に減色した画像が生成されます。元画像はPython版のサンプルコードと一緒に公開されているものをそのまま使わせていただきましたが、生成後の画像もPython版と同じようなものが生成できたので、問題なさそうです。

オリジナル画像

2色

3色

5色

16色

　コードは下記にも公開しました。

github.com