今回は下記書籍のロジスティック回帰とパーセプトロンの比較のコードを Ruby で実装してみたいと思います。
ITエンジニアのための機械学習理論入門 https://www.amazon.co.jp/IT-ebook/dp/B016Q22IX2/
サンプルコードはこちらで公開されています。
ロジスティック回帰のコード
データセットの用意やパーセプトロンのアルゴリズムの実装については前回までと基本的な内容は同じです。ロジスティック回帰のメソッドは下記のように実装しました。
def run_logistic(train_set, plot) w = Numo::NArray[[0], [0.1], [0.1]] phi = train_set[:x, :y] phi[:bias] = Array.new(phi.size).fill(1) t = train_set[:type] t = t.to_matrix # 最大30回のIterationを実施 30.times do # IRLS法によるパラメータの修正 y = [] phi.each_row do |line| a = Vector[*line.to_a].dot(w) y << (1.0 / (1.0 + Math.exp(-a))) end r = Matrix.diagonal(*(Numo::NArray[*y] * (1 - Numo::NArray[*y]).to_a)) y = Numo::NArray[y].transpose tmp1 = Matrix[*Numo::NArray[*phi.to_matrix.transpose.to_a].dot(Numo::NArray[*r.to_a]).dot(Numo::NArray[*phi.to_matrix.to_a]).to_a].inverse tmp2 = Numo::NArray[*phi.to_matrix.transpose.to_a].dot(y - Numo::NArray[*t.transpose.to_a]) w_new = w - Numo::NArray[*tmp1.to_a].dot(tmp2) # パラメータの変化が 0.1% 未満になったら終了 if (w_new - w).transpose.dot(w_new - w).flatten[0] < (0.001 * (w.transpose.dot(w))).flatten[0] w = w_new break end end # 分類誤差の計算 w0 = w[0] w1 = w[1] w2 = w[2] err = 0 train_set.each_row do |point| x = point.x y = point.y type = point[:type] type = type * 2 - 1 if type * (w0 + w1*x + w2*y) < 0 err += 1 end end err_rate = err * 100 / train_set.size # 結果を表示 xmin = train_set.x.min - 5 xmax = train_set.x.max + 10 linex = Numo::NArray[*(xmin.to_i-5..xmax.to_i+4).to_a] liney = -linex * w1 / w2 - w0 / w2 label = "ERR %.2f%%" % err_rate line_err = plot.add(:line, linex.cast_to(Numo::Int64).to_a, liney.cast_to(Numo::Int64).to_a) line_err.title(label) line_err.color('blue') end
内積の計算
python では内積の計算を numpy.dot メソッドで行なっています。
numpy.dot — NumPy v1.12 Manual
高速数値計算ライブラリ「Numpy」覚書き - Pashango’s Blog
ruby では Daru::DataFrame にはそのまま内積を計算する処理はなさそうだったので、Vectorに変換した上で dot メソッドを使いました。
instance method Vector#dot (Ruby 2.4.0)
Ruby2.2 ではアレが死ぬほど使いやすくなるの! - Qiita
a = Vector[*line.to_a].dot(w)
入力配列を対角要素とする配列
python では numpy.diag メソッドで、入力配列を対角要素とする配列を取得しています。
numpy.diag — NumPy v1.10 Manual
r = np.diag(y*(1-y))
ruby では Matrix.diagonal メソッドを使用しました。
singleton method Matrix.diagonal (Ruby 2.4.0)
r = Matrix.diagonal(*(Numo::NArray[*y] * (1 - Numo::NArray[*y]).to_a))
パラメータの計算
python での実装は下記のようになっています。
tmp1 = np.linalg.inv(np.dot(np.dot(phi.T, r),phi)) tmp2 = np.dot(phi.T, (y-t)) w_new = w - np.dot(tmp1, tmp2)
内積の計算用の dot メソッドや、転置行列を求める T メソッド、逆行列を求める numpy.linalog.inv メッソドなどが使われています。
ruby で同様の処理を実装しようとしたところ、 Daru::DataFrame、 Numo::NArray、 Matrix それぞれで持っているメソッドを使うための変換処理でだいぶ煩雑になってしまいました。我ながらなんともイケてないコードです。前回までの実装もそうですが、Daru::DataFrame、Numo::NArray、Ruby標準のMatrix/Vector のどれか一つで基本的に処理することができればもっとすっきりかけるしやりようはあるのだろうと思っているのですが、まだその辺の使い方がよくわかっていません。元のpythonのコードの流れをそのままトレースするのではなく、そもそもどんな処理なのかをちゃんと理解して書く必要があるというのを改めて感じました。
tmp1 = Matrix[*Numo::NArray[*phi.to_matrix.transpose.to_a].dot(Numo::NArray[*r.to_a]).dot(Numo::NArray[*phi.to_matrix.to_a]).to_a].inverse tmp2 = Numo::NArray[*phi.to_matrix.transpose.to_a].dot(y - Numo::NArray[*t.transpose.to_a]) w_new = w - Numo::NArray[*tmp1.to_a].dot(tmp2)
それ以降の処理についてはパーセプトロンの場合とあまり変わらないので割愛します。
スクリプト全体
コード全体では下記のように実装しました。
# ロジスティック回帰とパーセプトロンの比較 require 'daru' require 'nyaplot' require 'numo/narray' Variances = [5, 10, 30, 50] # 両クラス共通の分散(4種類の分散で計算を実施) 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 # データセット {x_n,y_n,type_n} を用意 def prepare_dataset(variance) n1 = 10 n2 = 20 mu1 = [7, 7] mu2 = [-3, -3] sigma = Math.sqrt(variance) df1 = n1.times.map do [normal_rand(mu1[0], sigma), normal_rand(mu1[1], sigma)] end df1 = df1.transpose df1 = Daru::DataFrame.new(x: df1[0], y: df1[1], type: Array.new(n1).fill(1)) df2 = n2.times.map do [normal_rand(mu2[0], sigma), normal_rand(mu2[1], sigma)] end df2 = df2.transpose df2 = Daru::DataFrame.new(x: df2[0], y: df2[1], type: Array.new(n2).fill(0)) df = df1.concat(df2) df = df.reindex(Daru::Index.new(df.index.to_a.shuffle)) df[:index] = (n1 + n2).times.to_a df.set_index(:index) end # ロジスティック回帰 def run_logistic(train_set, plot) w = Numo::NArray[[0], [0.1], [0.1]] phi = train_set[:x, :y] phi[:bias] = Array.new(phi.size).fill(1) t = train_set[:type] t = t.to_matrix # 最大30回のIterationを実施 30.times do # IRLS法によるパラメータの修正 y = [] phi.each_row do |line| a = Vector[*line.to_a].dot(w) y << (1.0 / (1.0 + Math.exp(-a))) end r = Matrix.diagonal(*(Numo::NArray[*y] * (1 - Numo::NArray[*y]).to_a)) y = Numo::NArray[y].transpose tmp1 = Matrix[*Numo::NArray[*phi.to_matrix.transpose.to_a].dot(Numo::NArray[*r.to_a]).dot(Numo::NArray[*phi.to_matrix.to_a]).to_a].inverse tmp2 = Numo::NArray[*phi.to_matrix.transpose.to_a].dot(y - Numo::NArray[*t.transpose.to_a]) w_new = w - Numo::NArray[*tmp1.to_a].dot(tmp2) # パラメータの変化が 0.1% 未満になったら終了 if (w_new - w).transpose.dot(w_new - w).flatten[0] < (0.001 * (w.transpose.dot(w))).flatten[0] w = w_new break end end # 分類誤差の計算 w0 = w[0] w1 = w[1] w2 = w[2] err = 0 train_set.each_row do |point| x = point.x y = point.y type = point[:type] type = type * 2 - 1 if type * (w0 + w1*x + w2*y) < 0 err += 1 end end err_rate = err * 100 / train_set.size # 結果を表示 xmin = train_set.x.min - 5 xmax = train_set.x.max + 10 linex = Numo::NArray[*(xmin.to_i-5..xmax.to_i+4).to_a] liney = -linex * w1 / w2 - w0 / w2 label = "ERR %.2f%%" % err_rate line_err = plot.add(:line, linex.cast_to(Numo::Int64).to_a, liney.cast_to(Numo::Int64).to_a) line_err.title(label) line_err.color('blue') end # パーセプトロン def run_perceptron(train_set, plot) w0 = 0.0 w1 = 0.0 w2 = 0.0 bias = 0.5 * (train_set.x.mean + train_set.y.mean) # Iterationを30回実施 30.times do # 確率的勾配降下法によるパラメータの修正 train_set.each_row do |point| x = point.x y = point.y type = point[:type] type = type * 2 - 1 if type * (w0*bias + w1*x + w2*y) <= 0 w0 += type * bias w1 += type * x w2 += type * y end end end # 分類誤差の計算 err = 0 train_set.each_row do |point| x = point.x y = point.y type = point[:type] type = type * 2 - 1 if type * (w0*bias + w1*x + w2*y) <= 0 err += 1 end end err_rate = err * 100 / train_set.size # 結果を表示 xmin = train_set.x.min - 5 xmax = train_set.x.max + 10 linex = Numo::NArray[*(xmin.to_i-5..xmax.to_i+4).to_a] liney = -linex * w1 / w2 - bias * w0 / w2 label = "ERR %.2f%%" % err_rate line_err = plot.add(:line, linex.cast_to(Numo::Int64).to_a, liney.cast_to(Numo::Int64).to_a) line_err.title(label) line_err.color('red') end def run_simulation(variance, plot) train_set = prepare_dataset(variance) train_set1 = train_set.filter_rows {|row| row[:type] == 1 } train_set2 = train_set.filter_rows {|row| row[:type] == 0 } ymin = train_set.y.min - 5 xmin = train_set.x.min - 5 ymax = train_set.y.max + 10 xmax = train_set.x.max + 10 plot.configure do x_label("Variance: #{variance}") y_label('') xrange([xmin, xmax]) yrange([ymin, ymax]) legend(true) height(300) width(490) end scatter_true = plot.add_with_df(train_set1.to_nyaplotdf, :scatter, :x, :y) scatter_true.color('green') scatter_true.title('1') scatter_false = plot.add_with_df(train_set2.to_nyaplotdf, :scatter, :x, :y) scatter_false.color('orange') scatter_false.title('0') run_logistic(train_set, plot) run_perceptron(train_set, plot) end fig = Nyaplot::Frame.new Variances.each do |variance| plot = Nyaplot::Plot.new run_simulation(variance, plot) fig.add(plot) end fig.show
これを Jupyter Notebook 上で実行すると下記のようなグラフが表示されます。
python版と比べるとロジスティック回帰の優位性が見えにくい結果になっている気がするので、どこか実装に間違いがある可能性もあるかなと思っていますが、ひとまず近い形のグラフの表示までは行けた形です。もし、ここが違っている、とお気付きの方はコメントいただけると嬉しいです。
コードは下記にも公開しました。