ニューラルネットワークの数値微分による学習アルゴリズムの実装

 今回も引き続き「ゼロから作るDeepLearning」をベースに、数値微分による学習アルゴリズムをRubyで実装してみました。

www.oreilly.co.jp

 最初に書いておくと、今回の数値微分での実装は、実装はシンプルなもののその分処理に時間がかかり、手元の環境では繰り返し処理を終わらせることができませんでした。次回以降で書籍の次の章で解説されている、誤差逆伝播法での実装に変更してみたいと思います。

2層ニューラルネットワークのクラス

 書籍のPython実装に基づき、2層のニューラルネットワークをTwoLayerNetという一つのクラスとして実装します。まずはコード全体です。

require 'numo/narray'
require './softmax.rb'
require './sigmoid.rb'
require './cross_entropy_error.rb'
require './numerical_gradient.rb'

class TwoLayerNet
  def initialize(input_size, hidden_size, output_size, weight_init_std = 0.01)
    @params = {}
    @params[:w1] = weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm
    @params[:b1] = Numo::DFloat.zeros(hidden_size)
    @params[:w2] = weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm
    @params[:b2] = Numo::DFloat.zeros(output_size)
  end

  def params
    @params
  end

  def predict(x)
    w1 = @params[:w1]
    w2 = @params[:w2]
    b1 = @params[:b1]
    b2 = @params[:b2]

    a1 = x.dot(w1) + b1
    z1 = sigmoid(a1)
    a2 = z1.dot(w2) + b2
    softmax(a2)
  end

  def loss(x, t)
    y = predict(x)
    cross_entropy_error(y, t)
  end

  def accuracy(x, t)
    y = predict(x)
    y = y.max_index(1) % 10
    t = t.max_index(1) % 10

    y.eq(t).cast_to(Numo::UInt32).sum / x.shape[0].to_f
  end

  def numerical_gradients(x, t)
    loss_w = lambda {|w| loss(x, t) }

    grads = {}
    grads[:w1] = numerical_gradient(loss_w, @params[:w1])
    grads[:b1] = numerical_gradient(loss_w, @params[:b1])
    grads[:w2] = numerical_gradient(loss_w, @params[:w2])
    grads[:b2] = numerical_gradient(loss_w, @params[:b2])

    grads
  end
end

初期化

 まずはinitializeメソッドで重みとバイアスのパラメータを初期化し、インスタンス変数にHashとして保持します。合わせて外部からパラメータを参照するためのメソッドも用意しておきます。

def initialize(input_size, hidden_size, output_size, weight_init_std = 0.01)
  @params = {}
  @params[:w1] = weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm
  @params[:b1] = Numo::DFloat.zeros(hidden_size)
  @params[:w2] = weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm
  @params[:b2] = Numo::DFloat.zeros(output_size)
end

def params
  @params
end

 initializeメソッドの引数には、入力層のニューロン数、隠れ層のニューロン数、出力層のニューロン数を渡します。

 重みの初期値の生成は、まず Numo::DFloat.new で “入力層ニューロン数 x 隠れ層ニューロン数” もしくは “隠れ層ニューロン数 x 出力層ニューロン数” の行列を用意し、rand_normメソッドで標準正規分布に従う乱数を設定します。

Class: Numo::DFloat — Documentation by YARD 0.9.8

 バイアスの初期値は Numo::DFloat.zeros メソッドで全てゼロの配列を用意します。

Class: Numo::NArray — Documentation by YARD 0.9.8

推論処理

 画像データを引数にとり、推論処理を行います。

def predict(x)
  w1 = @params[:w1]
  w2 = @params[:w2]
  b1 = @params[:b1]
  b2 = @params[:b2]

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)
  a2 = z1.dot(w2) + b2
  softmax(a2)
end

 処理内容としては前回の記事で実装したものと同じで、隠れ層の活性化関数にシグモイド関数、出力層の活性化関数にソフトマックス関数を使用しています。

def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end
def softmax(a)
  if a.ndim == 2
    a = a.transpose
    a = a - a.max(0)
    y = Numo::DFloat::Math.exp(a) / Numo::DFloat::Math.exp(a).sum(0)
    return y.transpose
  end

  c = a.max
  exp_a = Numo::DFloat::Math.exp(a - c)
  sum_exp_a = exp_a.sum
  exp_a / sum_exp_a
end

損失関数

 入力データ(画像データ)と教師データ(正解ラベル)を引数にとり、損失関数の値を求めます。

def loss(x, t)
  y = predict(x)
  cross_entropy_error(y, t)
end

 画像データからpredictメソッドで推論を行なった結果と正解ラベルから交差エントロピー誤差を求めています。

def cross_entropy_error(y, t)
  if y.ndim == 1
    t = t.reshape(1, t.size)
    y = y.reshape(1, y.size)
  end

  batch_size = y.shape[0]
  -(t * (Numo::DFloat::Math.log(y))).sum / batch_size # one-hot表現用
end

認識精度の計算

 入力データ(画像データ)と教師データ(正解ラベル)を引数にとり、認識精度を計算します。

def accuracy(x, t)
  y = predict(x)
  y = y.max_index(1) % 10
  t = t.max_index(1) % 10

  y.eq(t).cast_to(Numo::UInt32).sum / x.shape[0].to_f
end

 画像データからpredictメソッドで推論を行なった結果は0〜9の数字である確度(確率)の配列として返されるため、どの数字である確率が最も高いかを max_index メソッドを用いて取得します。推論はバッチ処理で行うため、複数の画像データに対する結果として二次元配列として返ってくるので、二次元目を基準に最大値を求めるため、max_indexの引数には1を渡しています(0次元目、1次元目と考えるため)。ただし、max_index では多次元配列の場合、全ての要素数に対してのインデックス値(10 x 10の配列だったら0〜99の値)として戻ってくるため、10で割ったあまりを求めることで、0〜9のラベルを表すようにしています。

y = y.max_index(1) % 10

 これは正解ラベルについても同様で、 one-hot表現で渡されている複数データについての正解ラベルは二次元配列なので、10で割ることで0〜9のラベルを表すようにしています。

t = t.max_index(1) % 10

 そして最後に推論結果のラベル配列と正解ラベル配列をeqメソッドで比較し、正解数(配列の値が1になっている数)を入力データの数で割ることで、正解率を求めています。この辺りのeqメソッドの使い方等は前回記事と同様です。

重みパラメータに対する勾配の計算

 入力データ(画像データ)と教師データ(正解ラベル)を引数にとり、各パラメータに対する勾配を計算します。

def numerical_gradients(x, t)
  loss_w = lambda {|w| loss(x, t) }

  grads = {}
  grads[:w1] = numerical_gradient(loss_w, @params[:w1])
  grads[:b1] = numerical_gradient(loss_w, @params[:b1])
  grads[:w2] = numerical_gradient(loss_w, @params[:w2])
  grads[:b2] = numerical_gradient(loss_w, @params[:b2])

  grads
end

 損失関数を計算結果を取得するlambdaを用意し、勾配計算用のnumerical_gradientメソッドに各パラメータ共に渡し、結果をHashに格納して返します。

def numerical_gradient(f, x)
  h = 1e-4
  grad = Numo::DFloat.zeros(x.shape)

  x.size.times do |i|
    tmp_val = x[i]

    x[i] = tmp_val + h
    fxh1 = f.call(x)

    x[i] = tmp_val - h
    fxh2 = f.call(x)

    grad[i] = (fxh1 - fxh2) / (2 * h)
    x[i] = tmp_val
  end

  grad
end

 numerical_gradientメソッドでは中心差分での数値微分によって各パラメータの勾配を計算します。まず結果の格納用に入力パラメータと同じ形のゼロ配列を用意します。そして入力パラメータの各要素について (f(x + h) - f(x - h)) / 2h を計算して結果格納用の配列に格納し、最後にその配列を返しています。Numo::NArrayの多次元配列では単一のインデックスで配列の内容を参照する場合、全要素をフラットに並べた時のインデックス値を意味するので、x.size で全要素数を取得して、timesでその要素数分ループを回して、順次配列の内容を参照するようにしています。

ミニバッチ学習と評価の実装

 TwoLayerNetクラスを使ってMNISTデータセットを用いた学習と評価を行います。学習の実装は、訓練データから無作為に一部のデータを取り出して入力データとする、ミニバッチ学習で行います。まずはコード全体です。

require 'numo/narray'
require 'numo/gnuplot'
require './mnist.rb'
require './two_layer_net.rb'

# データの読み込み
x_train, t_train, x_test, t_test = load_mnist(true, true, true)

network = TwoLayerNet.new(784, 50, 10)

iters_num = 10_000 # 繰り返し回数
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

train_loss_list = []
train_acc_list = []
test_acc_list = []

iter_per_epoch = [train_size / batch_size, 1].max

iters_num.times do |i|
  batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
  x_batch = x_train[batch_mask, true]
  t_batch = t_train[batch_mask, true]

  # 勾配の計算
  grad = network.numerical_gradients(x_batch, t_batch)

  # パラメータの更新
  %i(w1 b1 w2 b2).each do |key|
    network.params[key] -= learning_rate * grad[key]
  end

  loss = network.loss(x_batch, t_batch)
  train_loss_list << loss

  next if i % iter_per_epoch != 0

  train_acc = network.accuracy(x_train, t_train)
  test_acc = network.accuracy(x_test, t_test)
  train_acc_list << train_acc
  test_acc_list << test_acc
  puts "train acc, test acc | #{train_acc}, #{test_acc}"
end

# グラフの描画
x = (0..(train_acc_list.size - 1)).to_a
Numo.gnuplot do
  plot x, train_acc_list, { w: :lines, t: 'train acc', lc_rgb: 'blue' },
       x, test_acc_list, { w: :lines, t: 'test acc', lc_rgb: 'green' }
  set xlabel: 'epochs'
  set ylabel: 'accuracy'
  set yrange: 0..1
end

MNISTデータの読み込みとパラメータ初期化

 まずは以前実装したMNISTデータのロード処理を使ってMNISTデータを読み込んだ後、TwoLayerNetの初期化と、繰り返し数等の設定を行います。

x_train, t_train, x_test, t_test = load_mnist(true, true, true)

network = TwoLayerNet.new(784, 50, 10)

iters_num = 10_000 # 繰り返し回数
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

 ニューロン数の構成は、入力層は28 x 28の画像データなので 784、出力層は0-9の数字を表すので 10、隠れ層は50としています。また、確率勾配降下法によるパラメータ更新の繰り返し数は10,000回で、ミニバッチのサイズは100としています。

ミニバッチデータの取得

 繰り返し処理の中では、まず60,000件の画像データの中からミニバッチデータを取得しています。

batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
x_batch = x_train[batch_mask, true]
t_batch = t_train[batch_mask, true]

 Numo::Int32.new でバッチサイズ分の要素数の配列を用意し、randメソッドで0以上60,000未満のランダム値を生成しています。これによって、60,000件の画像データのどのインデックス値のデータを取得するかを決めています。そしてそのインデックス値に該当する画像データと正解ラベルを取得していますが、それぞれ 60,000 x 784、60,000 x 10 の二次元配列になっているため、ランダムに生成したインデックス値の配列が一次元目を表すことを意味するように、二次元目はtrueを指定し、該当する二次元目のデータは全て取得するようにしています。

勾配の計算とパラメータの更新

 ミニバッチデータをTwoLayerNetの勾配計算メソッドに渡して勾配データを取得し、学習率をかけたものを各パラメータから引くことで、各パラメータを更新します。

# 勾配の計算
grad = network.numerical_gradients(x_batch, t_batch)

# パラメータの更新
%i(w1 b1 w2 b2).each do |key|
  network.params[key] -= learning_rate * grad[key]
end

損失関数の計算と記録

 ミニバッチデータから損失関数を計算し、繰り返しごとの結果を格納し、経過を記録します。

loss = network.loss(x_batch, t_batch)
train_loss_list << loss

認識精度の計算

 トレーニングデータとテストデータに対する認識精度を計算します。

next if i % iter_per_epoch != 0

train_acc = network.accuracy(x_train, t_train)
test_acc = network.accuracy(x_test, t_test)
train_acc_list << train_acc
test_acc_list << test_acc
puts "train acc, test acc | #{train_acc}, #{test_acc}"

 繰り返しごとに毎回認識精度を計算すると時間がかかってしまうため、1エポック(学習において全ての訓練データを使い切った時の回数。今回は60,000件の訓練データに対してミニバッチサイズが100なので、 60,000 / 100 = 600回)ごとに認識精度を計算してその値を結果配列に格納し、コンソールにも表示します。

認識精度推移のグラフ描画

 最後に認識精度の推移をグラフに描画します。

x = (0..(train_acc_list.size - 1)).to_a
Numo.gnuplot do
  plot x, train_acc_list, { w: :lines, t: 'train acc', lc_rgb: 'blue' },
       x, test_acc_list, { w: :lines, t: 'test acc', lc_rgb: 'green' }
  set xlabel: 'epochs'
  set ylabel: 'accuracy'
  set yrange: 0..1
end

 グラフの描画にはgnuplotを使い、Rubyでgnuplotを扱うためのgemとして、Numo::Gnuplot を使わせていただきました。

github.com

 こちらも参考にさせていただきました。

MF / 【Ruby】numo-gnuplotで遊ぶ

 以前に Jupyter Notebook 上で Nyaplot を使わせてもらっていたことはあるのですが、ターミナル上で動かして手軽にグラフを表示したいというときは Numo::Gnuplot が手軽に使えて良さそうです。

次は誤差逆伝播法

 冒頭でも書きましたが、今回実装はしたものの、繰り返しの一回の処理でもかなり時間がかかり、10,000回の繰り返しは現実的な時間では終わりそうもありませんでした。書籍でも言われていますが、数値微分による実装は実装はシンプルなものの処理にはかなり時間がかかるので、次は誤差逆伝播法での実装で試してみたいと思います。

 今回のコードは下記リポジトリで公開してあります。

github.com

ニューラルネットワークの推論処理

 前回に引き続き書籍「ゼロから作るDeepLearning」をベースに、前回NArray配列として扱うようにしたMNISTデータに対して推論処理を行うニューラルネットワークを実装してみます。

www.oreilly.co.jp

 ニューロンの構成は下記の通りです。

  • 入力層: 784 # 28 x 28 の画像データのピクセル数
  • 隠れ層1: 50
  • 隠れ層2: 100
  • 出力層: 10 # 結果として10クラス(0から9の数字)に分類する

サンプルコード全体

 まずはサンプルコードの全体を掲載しておきます。

require 'numo/narray'
require 'json'
require './sigmoid.rb'
require './softmax.rb'
require './mnist.rb'

def get_data
  x_train, t_train, x_test, t_test = load_mnist(true, true, false)
  [x_test, t_test]
end

def init_network
  nw = JSON.load(File.read('sample_weight.json'))
  network = {}
  nw.each do |k, v|
    network[k.to_sym] = Numo::DFloat[*v]
  end
  network
end

def predict(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
  softmax(a3)
end

x, t = get_data
network = init_network

batch_size = 100
accuracy_cnt = 0
x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  y_batch = predict(network, Numo::DFloat[*x_batch])
  p = y_batch.max_index(1) % 10
  accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum
end

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

MNISTデータの取得

 前回実装したMNISTデータのロード処理(mnist.rb)を使ってMNISTデータのテストデータを取得します。

def get_data
  x_train, t_train, x_test, t_test = load_mnist(true, true, false)
  [x_test, t_test]
end

重みパラメータのロード

 学習済みの重みとバイアスのパラメータは、書籍のサンプルコードで pickle ファイルとして提供されているものをあらかじめJSONに変換してファイルに保存しておき(sample_weight.json)、それを読み込んでいます。

def init_network
  nw = JSON.load(File.read('sample_weight.json'))
  network = {}
  nw.each do |k, v|
    network[k.to_sym] = Numo::DFloat[*v]
  end
  network
end

推論処理

 上記でロードしたテストデータと重みパラメータに対して、隠れ層での活性化関数にはシグモイド関数、出力層での活性化関数にはソフトマックス関数を使って推論処理を行います。シグモイド関数は sigmoid.rb、ソフトマックス関数は softmax.rb として保存して読み込んでおきます。

def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end
def softmax(a)
  c = a.max
  exp_a = Numo::DFloat::Math.exp(a - c)
  sum_exp_a = exp_a.sum
  exp_a / sum_exp_a
end
def predict(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
  softmax(a3)
end

処理の実行

 上記メソッドを使って処理を実行します。

x, t = get_data
network = init_network

batch_size = 100
accuracy_cnt = 0
x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  y_batch = predict(network, Numo::DFloat[*x_batch])
  p = y_batch.max_index(1) % 10
  accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum
end

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

 get_data で取得した画像データを100件ずつバッチ処理します。NArrayの二次元配列データはそのままループすると一次元配列として並べて各要素が参照されてしまうので、to_a で通常の配列データに変換した上で、 each_slice で100件ずつのまとまりにし、with_index でインデックスを取得します。

x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  ...
end

 ループの中では100件分の画像データを再度NArray配列に変換してpredictメソッドに渡し、推論処理を実行した結果を y_batch として受け取っています。predict では100件分の各画像データについて、0-9の各数字に対しての確度が確率として返されるので、100 x 10 の二次元配列になります。

y_batch = predict(network, Numo::DFloat[*x_batch])

 その y_batch に対して max_index メソッドで二次元目の配列データについて最も数値が大きい要素のインデックスを取得しています。ただし取得されるインデックスは 100 x 10 を一次元配列として並べた場合のインデックスになるので、10で割った余りを取得することで、0-9の分類に変換しています。

Class: Numo::Int32 — Documentation by YARD 0.9.8

p = y_batch.max_index(1) % 10

 処理している画像データに対応するラベルデータを取得します。each_sliceとwith_indexを使った場合、インデックスとしてはeach_sliceのまとまりごとに0から順番に振られるので、インデックスにバッチサイズ(100)をかけることで配列データの始点を特定し、それにさらにバッチサイズを加算したものから1を引くことで、終点を特定して、ラベルデータ配列 t から該当するラベルデータを取得します。

t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]

 それを分類結果データ p と eq メソッドで比較します。

Class: Numo::Int32 — Documentation by YARD 0.9.8

 eqメソッドでは配列の各要素を比較し、一致する場合は1を、一致しない場合は0を配列として返します。

irb(main):043:0* a
=> Numo::Int32#shape=[5]
[1, 3, 5, 7, 9]
irb(main):044:0> b
=> Numo::Int32#shape=[5]
[1, 2, 5, 7, 8]
irb(main):045:0> a.eq(b)
=> Numo::Bit#shape=[5]
[1, 0, 1, 1, 0]

 返される配列のデータ型はNumo::BitなのでNumo::UInt8に変換し、合計値を取得することで推論が正しかった要素の数が取得できます。

accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum

 そして最後に正解数を要素数で割ることで、正解率を計算しています。

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

 今回実装したコードはこちらにも公開しました。

github.com

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 上で実行すると下記のようなグラフが表示されます。

f:id:akanuma-hiroaki:20170301001605p:plain

f:id:akanuma-hiroaki:20170301001624p:plain

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

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 上で実行すると下記のようなグラフが表示されます。

f:id:akanuma-hiroaki:20170217073027p:plain

f:id:akanuma-hiroaki:20170217073042p:plain

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

github.com