[R]Monte Carlo Method による円周率計算

お仕事で Monte Carlo 法でシミュレーションをやってほしいという話があったので(その後、流れたけど)、wikipedia の最初の数十行だけお勉強。

モンテカルロ法の概要

  1. Define a domain of possible inputs.
  2. Generate inputs randomly from the domain using a certain specified probability distribution.
  3. Perform a deterministic computation using the inputs.
  4. Aggregate the results of the individual computations into the final result.

このアルゴリズムを利用し、円周率を求めるには次のようにする。

  1. 1辺2の正方形を描く
  2. その中に単位円を描く
  3. 4分割し、右上のブロックにランダムに点を打っていく。
  4. 打った点の数と単位円のなかにある点の数から円周率を求める。

# イメージ図

円周率をもとめるのは Monte Carlo 法の Hello World 的な例らしい。
で、実際にシミュレートした結果得られたのが下のグラフ。
赤線が円周率で横軸が打った点の数、縦軸がシミュレートして求めた円周率。
打った点の数は上から順に100個、1000個、10000個。

シミュレートに利用したのは次のコード。勉強のために R で書いてみた。

within_circle <- function(x, y) {
  if (x ^2 + y ^ 2 < 1) return (1)
  else return (0)
}

sim_pi_generator <- function(N) {
  counter <- 0
  circle_area <- 0
  xrand <- runif(N)
  yrand <- runif(N)
  return (function() {
    counter <<- counter + 1
    x <- xrand[counter]
    y <- yrand[counter]
    if (within_circle(x, y)) circle_area <<- circle_area + 1
    return (4 * circle_area / counter)
  }
 )
}

N <- 1000 # number of trials
sim_pi <- sim_pi_generator(N)
plot(replicate(N, sim_pi()), type='l', xlab='number of trials', ylab='estimation of pi', main='Monte Carlo estimate of pi using R')
abline(h=pi, col="red")

参考リンク

Advertisements
Tagged with: , ,
Posted in algorithm, R
2 comments on “[R]Monte Carlo Method による円周率計算
  1. foo says:

    within_circle <- function(x, y) {
    if (x ^2 + y ^ 2 < 1) return (1)
    else return (0)
    }

    within_circle <- function(x, y) {
    return (x ^2 + y ^ 2 < 1)
    }
    で良いでしょう。戻り値は TRUE/FALSE ですが,引用されるのが if ( ) の条件式ですから。1/0 にする必要がない)

  2. foo says:

    またお邪魔します。可能な限り短く。
    sim <- function(n=1000)
    {
    d <- matrix(runif(2*n), n)^2
    pi.est <- 4*cumsum(rowSums(d) < 1)/1:n
    plot(1:n, pi.est, type="l")
    abline(h=pi, col="red")
    }
    sim()
    sim(100000)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 4 months ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 4 months ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 4 months ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 10 months ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 1 year ago
%d bloggers like this: