Zenn
🐈

「ゼロからできるMCMC」本をJuliaで書いてみると，Cの代わりにJuliaも良いなという気持ちになりました

2020/09/24に公開
Julia
tech

まとめ

  • PythonもいいけどJuliaも良いよ！

まえがき

私はJuliaを0.2とか0.3ぐらいで知った後に，だいたい0.5の頃からお仕事プログラムの一部で利用しています．0.7/1.0に切り替わったときにそこそこの破壊的な更新があったと思いますが，1.0以降は比較的安定して使えるという印象を持っています．最近Twitterを見ていると，Juliaは1.0以降はそこそこ汎用的な言語で高速であり，チューニングを細かくしたい場合でなければ，あえてCythonやC++を使わなくてもいいのでは？という意見がある様子です．

以前たまたま読んでいた本がC/C++ (とはいってもほとんどCです) でMCMCを1から理解できるように詳しく説明された良い本があり，手元に本が戻ってきたので (出社したから持ち帰りました)，Juliaで書いてみようと思います．こちら題材の本です．

githubのレポジトリはこちらです．

この記事は，Juliaが優れているよ！といったことを宣伝するためではなく，「Pythonは使ったことあるけどC++はそこまで使えないから，なんとか高速化できないかな？」みたいな人がJuliaを選択肢に入れられると良いなと思って書きました ．

（こんなJuliaの書き方は下手くそだぞ！！！みたいなのは，はい，申し訳ないです，ごめんなさい)

πの近似値を求める (C/Julia/特に工夫していない脳死Python)

よくあるやつですね．Cの実装を見てみます．これは本の実装をほぼ○コピしたものです．

#include <stdio>
#include <stdlib.h>
#include <time.h>

const int n_iter = 1000000;
const int print_step = 10000;

int main() {
  srand((unsigned)time(NULL));
  int n_in = 0;
  for (int i = 1; i <= n_iter; i++) {
    double x = (double)rand() / RAND_MAX;
    double y = (double)rand() / RAND_MAX;
    if (x * x + y * y < 1.0f)
      n_in++;

    if (i % print_step == 0)
      printf("%d   %.10f\n", i, (double)n_in/n_iter);
  }
}

次はJuliaでの実装です．Pythonだとループを直接書くのにビクビクしますが (大げさです)，ここは特に脳死した状態でやってみます．LinearAlgebra: normは2乗ノルムの計算に使っていますがなくても良いです．

using LinearAlgebra: norm

function compute(;n_iter=1_000_000, print_step=1_000)
    n_in = 0

    for i in 1:n_iter
        xy = rand(2)
        (norm(xy) < 1) && (n_in += 1)
        (i % print_step == 0) && println("$i $(n_in / i)")
    end

    println("true: $(pi/4)")
end

compute()

いずれの実装も0.785ぐらいになっていましたね．ついでにPythonでも工夫せずにfor文で書くとこういう雰囲気です．

import random

n_in = 0
n_iter = 1000000
print_step = 1000

for i in range(1, n_iter + 1):
    x = random.random()
    y = random.random()
    if x * x + y * y <= 1:
        n_in += 1

    if i % print_step == 0:
        print("{} {}".format(i, n_in / i))

print(n_in / n_iter)

Pythonがスッと頭に入っている人はJuliaもスッと入ってきそうですね！（要検討）

補足

繰り返しですが，別にC++/Python/Juliaで，どれを選ぶべき！Pythonは遅い！みたな話をするつもりは特にありません．ただ，Python使えてC++勉強する（他にもCythonを書く，とかもそうですが），というようなステップを踏まなくても，同じJupyter notebook上ですぐに使えるJuliaを使うのも選択肢なのかな，という話がしたかったです．ちなみにJuliaからはPythonのパッケージも簡単に呼べるので，JuliaでつらいところはPythonで書けば良いです．

私はPythonを使っているときはmatplotlibを使いますが，JuliaのときはPlots.jl/grバックエンドを使っていることもあります．

例えばガウス分布からのサンプリングをヘイスティング法で計算して描画する部分のコード (本の4章です) はこのような感じになっています．ほとんど好みぐらいの差でしょうか．好きなものを使えばよいと思います．

function plot_values_in_hist(values; n_bins=50)
    f = plot(size=(400, 300))
    histogram!(f, values, bins=n_bins, normalize=true)
    savefig(f, "chapter4/histogram1_n$(length(values)).png")
end

オマケ

ところで実装するときですが，VS Codeで編集しています．ターミナルを右側に出してOhMyREPL.jlとRevise.jlを読み込みます (OhMyREPL.jlはいつも)．ここで編集したコードをinclude("")で呼び出しています．このときエラーはいた場合に再読込が反映されないのがあるのですが，これって何が問題なんですかね？includetじゃないからダメとかそういう感じでしょうか？VS Codeでincremental実装するときの最高のやり方を識者の方教えてください…．

Discussion