ニューラルネットワークの誤差逆伝播法による学習アルゴリズムの実装

　今回も引き続き「ゼロから作るDeepLearning」をベースに、前回数値微分で実装した学習アルゴリズムの誤差逆伝播法版をRubyで実装してみました。計算の内容等は書籍を参照いただくとして、Rubyで実装した際のポイントを説明していきたいと思います。

www.oreilly.co.jp

活性化関数レイヤ実装

　今回はニューラルネットワークの各レイヤをそれぞれ一つのクラスとして実装しています。活性化関数レイヤはReLUレイヤとSigmoidレイヤです。

ReLUレイヤ

　ReLUレイヤは下記のように実装しました。

class Relu
  def initialize
    @mask = nil
  end

  def forward(x:)
    @mask = (x <= 0)
    out = x.copy
    out[@mask] = 0
    out
  end

  def backward(dout:)
    dout[@mask] = 0
    dout
  end
end

　順伝播時のforwardメソッドの引数としてはNArray配列を期待しています。NArray配列 x に対して (x <= 0) をすると0以下の要素が1、それ以外が0のBit配列を返しますので、それをマスクに使い、入力に対して0以下の要素を0に置き換えた結果を返しています。また、マスクは逆伝播時にも使うのでインスタンス変数に保持します。

　逆伝播時のbackwardメソッドでは順伝播時に保持したマスクを元に、上流からの入力に対してマスクの要素が1になっている要素に0を設定しています。

Sigmoidレイヤ

　Sigmoidレイヤの実装は下記のように行いました。

require './sigmoid.rb'

class Sigmoid
  def initialize
    @out = nil
  end

  def forward(x:)
    @out = sigmoid(x)
    @out
  end

  def backword(dout:)
    dout * (1.0 - @out) * @out
  end
end

　順伝播時は以前実装したsigmoidメソッドを呼んでいるだけで、その結果をインスタンス変数に保持しておきます。

　逆伝播時は順伝播時の出力を元に計算を行なった結果を返します。

AffineレイヤとSoftmaxレイヤの実装

　ニューラルネットワークの順伝播において重みの計算とバイアスの加算を行なっていた層をAffineレイヤとして実装し、出力層ではソフトマックス関数を用いて出力を正規化するSoftmaxレイヤを実装し、それぞれの順伝播、逆伝播の処理を実装します。

Affineレイヤ

　Affineレイヤの実装は下記の通りです。

class Affine
  def initialize(w:, b:)
    @w = w
    @b = b

    @x = nil
    @original_x_shape = nil

    # 重み・バイアスパラメータの微分
    @dw = nil
    @db = nil
  end

  def dw
    @dw
  end

  def db
    @db
  end

  def forward(x:)
    # テンソル対応
    @original_x_shape = x.shape
    x = x.reshape(x.shape[0], nil)
    @x = x

    @x.dot(@w.first) + @b.first
  end

  def backward(dout:)
    dx = dout.dot(@w.first.transpose)
    @dw = @x.transpose.dot(dout)
    @db = dout.sum(0)

    dx.reshape(*@original_x_shape)
  end
end

　initializeメソッドでは重み、バイアスパラメータをインスタンス変数に格納し、それ以外にも処理に必要になるインスタンス変数を定義しています。

　順伝播時は入力の行列の形と入力行列を保持し、入力行列と重みパラメータの内積にバイアスを加算した結果を返します。

　逆伝播時はまず内積の逆伝播の計算として重みパラメータと順伝播時の入力値の転置行列を使った計算を行なっています。NArray行列の転置行列はtransposeメソッドで取得できます。

Numo::NArray#transpose
http://ruby-numo.github.io/narray/narray/Numo/NArray.html#transpose-instance_method

　バイアスの逆伝播の計算では行列の0番目の軸（データ単位）に対しての合計を求めるため、sumメソッドのパラメータに0を指定しています。

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

　最後に入力値の逆伝播の計算結果を入力値と同じ行列の形にreshapeして返しています。入力値の形状は変数に配列データとして格納していますが、reshapeメソッドのパラメータは配列ではないので、* で配列を展開して渡しています。

Softmaxレイヤ

　Softmaxレイヤは損失関数である交差エントロピー誤差の計算も含めて、SoftmaxWithLossクラスとして下記のように実装しました。

require './softmax.rb'
require './cross_entropy_error.rb'

class SoftmaxWithLoss
  def initialize
    @loss = nil
    @y = nil # softmaxの出力
    @t = nil # 教師データ
  end

  def forward(x:, t:)
    @t = t
    @y = softmax(x)
    @loss = cross_entropy_error(@y, @t)

    @loss
  end

  def backward(dout: 1)
    batch_size = @t.shape[0]
    if @t.size == @y.size # 教師データがon-hot-vectorの場合
      return (@y - @t) / batch_size
    end

    dx = @y.copy

    (0..(batch_size - 1)).to_a.zip(@t).each do |index_array|
      dx[*index_array] -= 1
    end

    dx / batch_size
  end
end

　順伝播時は以前実装したsoftmaxメソッドとcross_entropy_errorメソッドを使用しています。入力値をsoftmaxメソッドで正規化し、その結果と教師データをcross_entropy_errorメソッドに渡して交差エントロピー誤差を計算しています。

　逆伝播時はデータ一個あたりの誤差を伝播するために誤差をデータ数で割っています。

　行列の要素の参照方法については、Pythonのndarrayの場合は配列を二つ渡すと、それぞれを行・列のインデックスとして要素を参照してくれますが、NArray行列で同じような指定をすると、一つ目の配列で指定した全ての行で、二つ目の配列に指定した全ての要素が取得されてしまいます。

>>> array
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])
>>> array[[0, 2], [1, 3]]
array([ 2, 12])

irb(main):021:0> array                
=> Numo::UInt32#shape=[3,4]           
[[1, 2, 3, 4],                        
 [5, 6, 7, 8],                        
 [9, 10, 11, 12]]                     
irb(main):022:0> array[[0, 2], [1, 3]]
=> Numo::UInt32(view)#shape=[2,2]     
[[2, 4],                              
 [10, 12]]                            

　そこでzipメソッドを使ってそれぞれの配列から一つずつデータを取得し、直接一つの要素を参照するようにしました。

irb(main):024:0* [0, 2].zip([1, 3]).each do |array_index| 
irb(main):025:1*   puts array[*array_index]               
irb(main):026:1> end                                      
2                                                         
12                                                        
=> [[0, 1], [2, 3]]                                       

誤差逆伝播法でのニューラルネットワーク実装

　上記で実装したレイヤを組み合わせて、誤差逆伝播法でのニューラルネットワークを実装します。下記のように一つのクラスとして実装します。

require 'numo/narray'
require './numerical_gradient.rb'
require './layers.rb'

class TwoLayerNet
  def initialize(input_size:, hidden_size:, output_size:, weight_init_std: 0.01)
    # 重みの初期化
    Numo::NArray.srand
    @params = {
      w1: [weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm],
      b1: [Numo::DFloat.zeros(hidden_size)],
      w2: [weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm],
      b2: [Numo::DFloat.zeros(output_size)]
    }

    # レイヤの生成
    @layers = {
      affine1: Affine.new(w: @params[:w1], b: @params[:b1]),
      relu1:   Relu.new,
      affine2: Affine.new(w: @params[:w2], b: @params[:b2])
    }
    @last_layer = SoftmaxWithLoss.new
  end

  def params
    @params
  end

  def predict(x:)
    @layers.values.inject(x) do |x, layer|
      x = layer.forward(x: x)
    end
  end

  # x: 入力データ, t: 教師データ
  def loss(x:, t:)
    y = predict(x: x)
    @last_layer.forward(x: y, t: t)
  end

  def accuracy(x:, t:)
    y = predict(x: x)
    y = y.max_index(1) % 10
    if t.ndim != 1
      t = t.max_index(1) % 10
    end

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

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

    {
      w1: numerical_gradient(loss_w, @params[:w1].first),
      b1: numerical_gradient(loss_w, @params[:b1].first),
      w2: numerical_gradient(loss_w, @params[:w2].first),
      b2: numerical_gradient(loss_w, @params[:b2].first)
    }
  end

  def gradient(x:, t:)
    # forward
    loss(x: x, t: t)

    # backward
    dout = 1
    dout = @last_layer.backward(dout: dout)

    layers = @layers.values.reverse
    layers.inject(dout) do |dout, layer|
      dout = layer.backward(dout: dout)
    end

    {
      w1: @layers[:affine1].dw,
      b1: @layers[:affine1].db,
      w2: @layers[:affine2].dw,
      b2: @layers[:affine2].db
    }
  end
end

　initializeでは重みとバイアスパラメータの初期化と、各レイヤのインスタンスの生成を行なっています。重みとバイアスのパラメータを配列として保持しているのは、学習結果をTwoLayerNetクラス内のパラメータに反映したものを各レイヤで参照できるようにするためです。Affineクラスのインスタンス生成時にパラメータを w: @params[:w1] という形で渡していますが、Rubyでは参照の値渡しになるため、TwoLayerNet側で @params[:w1] に計算結果を代入しても、 Affineインスタンス側では初期化時に渡されたパラメータを参照し続けます。そこで参照先を配列にして、計算結果はその配列の中身を更新する形にしています。書籍のPythonコードでは参照渡しになるため、配列として保持しなくてもTwoLayerNet側での変更がAffineインスタンス側で参照されています。

　それと、パラメータをランダムに生成する前に、Numo::NArray.srandメソッドでseedが変わるようにしています。これをしないとスクリプト実行時に毎回同じ内容のパラメータが作成されてしまいます。

　また、ニューラルネットワークでは各レイヤが実行される順番が重要なので、Pythonの場合は通常のDictionaryではなくOrderedDictを使っていますが、RubyのHashでは順番が保持されるため、通常のHashをそのまま使っています。

　勾配の計算処理のgradientメソッドでは、まず順伝播の処理を行うためにlossメソッドを実行します。lossメソッドではpredictメソッドを呼び出し、injectメソッドで各レイヤのforward処理を実行し、順伝播の処理を行なっています。そして最後に last_layer変数に保持しているSoftmaxWithLossインスタンスの順伝播処理で交差エントロピー誤差を計算しています。

　続いて逆伝播処理ではまずSoftmaxWithLossインスタンスの逆伝播処理を行なったあと、各レイヤを逆順に逆伝播処理を行ない、計算結果を返しています。

　精度確認用のaccuracyメソッドでは、まず各レイヤの順伝播処理を行い、その結果一番確度の高い要素のインデックスと教師データを付き合わせて正解率を計算しています。max_indexメソッドでは各データの最大値を判定したいので、引数で1軸目を指定しています（0だと全ての要素の中からの最大値を判定する）。結果は行列全ての要素の中でのインデックス値を返すので、10で割ることで各データのインデックス値に変換しています。

　eqメソッドでは一致する要素は1、一致しない要素は0のNumo::Bit配列を返します。正解数としてBitが1の要素の合計を計算するため、Numo::Bit配列をcast_toメソッドでInt配列に変換しています。この時、各ビット値は1か0なので、Int8でも正しく変換できますが、合計を計算するときにInt8の範囲を超えると正しく合計値が取得できません。今回の入力データの件数は60,000件あるので、Int8だと範囲を超えてしまうため、UInt16に変換した上で合計値を取得しています。

誤差逆伝播法の勾配確認

　誤差逆伝播法で求めた勾配が正しいかどうかを確認するため、数値微分で求めた勾配と比較して確認します。

require './mnist.rb'
require './two_layer_net.rb'

x_train, t_train, x_test, t_test = load_mnist(normalize: true, one_hot_label: true)

network = TwoLayerNet.new(input_size: 784, hidden_size: 50, output_size: 10)

x_batch = x_train[0..2, true]
t_batch = t_train[0..2, true]

grad_numerical = network.numerical_gradients(x: x_batch, t: t_batch)
grad_backprop = network.gradient(x: x_batch, t: t_batch)

grad_numerical.keys.each do |key|
  diff = (grad_backprop[key] - grad_numerical[key]).abs.mean
  puts "#{key}: #{diff}"
end

　MNISTデータの先頭3件を使い、数値微分での計算と誤差逆伝播法での計算を一度行い、それぞれの結果の差を計算しています。実行結果は下記のようになり、ほぼ差がないことが確認できます。

[vagrant@localhost vagrant]$ ruby gradient_check.rb
w1: 2.616843054226442e-13                          
b1: 8.347379796928845e-13                          
w2: 1.0232621862468414e-12                         
b2: 1.2012612987666316e-10                         

誤差逆伝播法を使った学習

　前回の記事で実装したミニバッチ学習と評価を、誤差逆伝播法を使うように変更します。前回との違いはnumerical_gradientsメソッドではなくgradientメソッドを使うようにした点です。

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

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

network = TwoLayerNet.new(input_size: 784, hidden_size: 50, output_size: 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|
  Numo::NArray.srand
  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)
  grad = network.gradient(x: x_batch, t: t_batch)

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

  loss = network.loss(x: x_batch, t: t_batch)
  train_loss_list << loss

  next if i % iter_per_epoch != 0

  train_acc = network.accuracy(x: x_train, t: t_train)
  test_acc = network.accuracy(x: x_test, t: 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

　実行結果は下記のようになります。シェルからスクリプトファイルを実行してもグラフの描画が行われなかったので、irbからロードすることで実行し、グラフが表示されるようにしています。

irb(main):001:0> load './train_neuralnet.rb'
train acc, test acc | 0.12578333333333333, 0.1225
train acc, test acc | 0.9026833333333333, 0.9071
train acc, test acc | 0.92225, 0.923
train acc, test acc | 0.9348833333333333, 0.9331
train acc, test acc | 0.9436, 0.9403
train acc, test acc | 0.9513666666666667, 0.9503
train acc, test acc | 0.9568833333333333, 0.9565
train acc, test acc | 0.9614666666666667, 0.9594
train acc, test acc | 0.96415, 0.9611
train acc, test acc | 0.9659, 0.9616
train acc, test acc | 0.9677833333333333, 0.9622
train acc, test acc | 0.97175, 0.9655
train acc, test acc | 0.97275, 0.9666
train acc, test acc | 0.9745666666666667, 0.968
train acc, test acc | 0.9764166666666667, 0.9682
train acc, test acc | 0.9780666666666666, 0.9696
train acc, test acc | 0.9788166666666667, 0.9704
=> true

実行速度

　手元のVagrant環境で実行した結果、Ruby版とPython版のそれぞれの実行速度は下記のようになりました。

Ruby: 4分41秒 Python: 39秒

　Python版の方がかなり早いです。今回はPython版をベースにRuby版を実装したので、パフォーマンスチューニングが可能かやってみたいところです。

　コードは下記リポジトリでも公開しています。

github.com

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

　今回も引き続き「ゼロから作る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 上で実行すると下記のようなグラフが表示されます。

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

github.com