MCMCへの道（１）

 そういえば，機械学習をまともに勉強していなかった。いや，そもそも統計学も確率論もあれもこれもである。先日の水素原子波動関数の可視化についても進んでいないのだけれど，どうやら，メトロポリス・ヘイスティング法を使うのが望ましいらしい（by tsujimotter）。そういえば，モンテカルロ法もまともに勉強してこなかった。

そこで，反省してMCMC（マルコフ連鎖モンテカルロ法）の全体像をマスターすべく，いろいろ探したが，おこちゃま向けの適当なテキストがない。こうなるとYouTubeだのみだ。機械学習基礎理論独習というサイトがあったので，このPythonコードをJuliaで再現しながら学ぶことにした。

その１回目はMCMC法#1モンテカルロ法なので，これを再現してみた。

using BenchmarkTools
using Random
using Plots

function findpi(n)
  rng = MersenneTwister(0)
  count = 0
  for i in 1:n
    count += ifelse(rand(rng)^2+rand(rng)^2<=1.0,1,0)
  end
  return 4*count/n
end

#@btime findpi(10^6)
#@benchmark findpi(10^8)

function circle(n)
  x=rand(n)
  y=rand(n)
  c=zeros(n)
  for i in 1:n
    c[i]=ifelse(x[i]^2+y[i]^2<=1.0,0.0,1.0)
  end
  scatter(x,y,
  marker_z = c,
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.2,
  markercolor = :grey,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white,
  markerstrokestyle = :dot,
  aspect_ratio = 1.0,
  xlim=(0.0,1.0),
  ylim=(0.0,1.0),
  legend=:topright
  )
end

circle(10000)

jupyterlabを使っているが，多数のデータの描画で表示が重くなってかたまる現象がみられる。黒木さんは400msで10^8点のデータが計算できるとしていたが，こちらではメルセンヌツイスターアルゴリズムを用いて300msくらいかかってしまう。

図：単純なモンテカルロ法による円周率の計算