今回は下記書籍のベイズ推定による回帰分析のコードを Ruby で実装してみたいと思います。今回で下記書籍のコードの置き換えは終了です。
サンプルコードはこちらで公開されています。
単位行列の生成
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 上で実行すると下記のようなグラフが表示されます。
コードは下記にも公開してあります。
