『ゼロからできるMCMC』 Chapter 2をJuliaで

最終更新：2020/11/7

『ゼロからできるMCMC』 Chapter 2のコードを, 練習を兼ねてJuliaで書いてみました. 詳細なことを書くと著作権的にだめな気もするので, 解説などはありません.

Julia Version 1.5.0

一様乱数を用いた円周率の計算 ( p.20 )

モンテカルロ法で円周率を求める.

Julia 

using Distributions
using Random
Random.seed!(12345)

niter = 10000;

xx = rand(Uniform(), niter);
yy = rand(Uniform(), niter);

4*sum(@. (xx^2 + yy^2) < 1.0) / niter

でした. 思ったより近づかない印象.

一様乱数を用いた積分 ( p.23 )


モンテカルロ法で円周率を求める.

Julia 

using Distributions
using Random: seed!

seed!(12345)

niter = 1000;
xx = rand(Uniform(), niter);
res =  sum(@. sqrt(1 - xx^2))/niter;
println("$(4*res)")

でした.

期待値と積分 ( p.30 )

モンテカルロ法で三次元球 ( 半径１ ) の体積を求める.

Julia 
using Random: seed!
using  BenchmarkTools
seed!(12345)

function main(niter)
    n_in = zero(Int64)
    z    = 0.0
    for _ in 1:niter
        x, y = rand(), rand()
        if x^2 + y^2 ≤ 1
            n_in += 1
            z += sqrt(1.0-x*x-y*y)
        end
    end
    println(2π*z/n_in)
end

# result
julia> 	@btime main(10^5)
	922.678 μs (0 allocations: 0 bytes)
	4.180887719852872

函数にした方がはやいらしいので. 早いのかは今ひとつわからないけど, 0 allocationsなのはいいことやと思う.

ボックスミュラー法を用いた\( \langle x^2 \rangle \) の計算 ( p.33 )

ボックスミュラー法を用いて\( \langle x^2 \rangle \) を計算してみる. 1になるはずである.

Julia 
using Random: seed!
using  BenchmarkTools
seed!(12345)

function main(niter)
    z = 0.0
    for _ in 1:niter
        p, q = rand(), rand()
        x = sqrt(-2log(p))*cos(2π*q)
        z += x*x
    end

    z/niter
end

# result
julia> 	@btime main(10^5)
  	2.394 ms (0 allocations: 0 bytes)
	0.9918424337202467

だいたい1になった.

函数から外すとどれくらい遅くなるのかと思い計算してみた. niter=\( 10^8 \) で, 函数化すると2.5秒, 函数から外すと11秒であった.

この記事のTOP    BACK    TOP