読者です 読者をやめる 読者になる 読者になる

30代無職のプログラミング入門

暇つぶしにプログラミングを独学してみる

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

コードはRです。解説はできないので端折りますが、『高校数学+α』のpp.518-9の説明が分かりやすかったです。著者のサイトで読めます→ http://tad311.xsrv.jp/hsmath/

library(tidyverse)

N = 10000

df <- tibble(
  x = runif(N, -1, 1),
  y = runif(N, -1, 1),
  inside =  (x**2 + y**2) <= 1
)

pi_monte = sum(inside) * 4/N
error = abs((pi_monte - pi)/pi) * 100

結果:

> print(pi_monte)
[1] 3.1296
> print(error)
[1] 0.381738

3.1296なので、円周率っぽい数字にはなってる。誤差は0.381738%。

ggplot2で作画。

ggplot(data = df, mapping = aes(x = x, y = y)) +
  geom_point(mapping = aes(color=inside)) +
  coord_fixed() # 散布図のサイズを正方形に

f:id:unEmployed:20170418212314p:plain