Processing math: 100%
　

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

最終更新：2020/11/7

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

Julia Version 1.5.0

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

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

Julia 
  1.  
  2. using Distributions
  3. using Random
  4. Random.seed!(12345)
  5.  
  6. niter = 10000;
  7.  
  8. xx = rand(Uniform(), niter);
  9. yy = rand(Uniform(), niter);
  10.  
  11. 4*sum(@. (xx^2 + yy^2) < 1.0) / niter

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

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


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

Julia 
  1.  
  2. using Distributions
  3. using Random: seed!
  4.  
  5. seed!(12345)
  6.  
  7. niter = 1000;
  8. xx = rand(Uniform(), niter);
  9. res =  sum(@. sqrt(1 - xx^2))/niter;
  10. println("$(4*res)")

でした.

期待値と積分 ( p.30 )

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

Julia 
  1. using Random: seed!
  2. using  BenchmarkTools
  3. seed!(12345)
  4.  
  5. function main(niter)
  6.     n_in = zero(Int64)
  7.     z    = 0.0
  8.     for _ in 1:niter
  9.         x, y = rand(), rand()
  10.         if x^2 + y^2 ≤ 1
  11.             n_in += 1
  12.             z += sqrt(1.0-x*x-y*y)
  13.         end
  14.     end
  15.     println(2π*z/n_in)
  16. end
  17.  
  18. # result
  19. julia> 	@btime main(10^5)
  20. 	922.678 μs (0 allocations: 0 bytes)
  21. 	4.180887719852872

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

ボックスミュラー法を用いた⟨x2⟩ の計算 ( p.33 )

ボックスミュラー法を用いて⟨x2⟩ を計算してみる. 1になるはずである.

Julia 
  1. using Random: seed!
  2. using  BenchmarkTools
  3. seed!(12345)
  4.  
  5. function main(niter)
  6.     z = 0.0
  7.     for _ in 1:niter
  8.         p, q = rand(), rand()
  9.         x = sqrt(-2log(p))*cos(2π*q)
  10.         z += x*x
  11.     end
  12.  
  13.     z/niter
  14. end
  15.  
  16. # result
  17. julia> 	@btime main(10^5)
  18.   	2.394 ms (0 allocations: 0 bytes)
  19. 	0.9918424337202467

だいたい1になった.

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

この記事のTOP    BACK    TOP