---
title: "Montecarlo Method"
author: "Aya Naito"
date: "2015年10月28日"
output: html_document
---

Montecarlo法を用いて円周率$\pi$の近似値を求める。

　

半径が$r$の円と、その円に外接する正方形の面積を考える。
円の面積は$\pi r^{2}$であり、正方形の面積は$4r^{2}$である。

これらの面積比を$a$として計算すると、$a=\dfrac{\pi r^{2}}{4r^{2}}=\dfrac{\pi}{4}$となる。

つまり、$a$の値を推定することで円周率を推定できる。

ここで、この正方形内にランダムに$n$個の点を打ち、そのうち円内部に入った点の個数を$m$とする。

このとき、円内部に点が打たれる確率は$\dfrac{m}{n}$である。

よって、$\dfrac{\pi}{4}=\dfrac{m}{n}$すなわち$\pi=\dfrac{4m}{n}$となる。

これによって円周率を推定することができる。

以下のプログラムでは簡単のため、円の半径を$1$とし、第一象限のみで考えるものとする。

```{r}
montecarlo <- function(n){
  m <- 0 #円内部にある点の個数; 始めは0個とする
  cols <- NULL #plotする点に色をつけるためのもの
  range <- c(0,1) #範囲を第一象限とする
  x <- runif(n,0,1) #一様乱数でx座標を定める; runif(試行回数,開始値,終了値)
  y <- runif(n,0,1) #一様乱数でy座標を定める
  r <- x^2+y^2 #円の範囲
  for(i in r){ 
    if(i<=1){
      m <- m +1 #円内部であれば、円内部の点の個数を1つ増やす
      cols <- c(cols, "red") #更に赤色でプロットする
    } else{
      cols <- c(cols,"blue") #円外部であれば青色でプロットする
    }
  }
# forループは以下のようにも書ける
#  inner <- (r<=1)
#  m <- sum(inner)
#  cols <- ifelse(inner,"red","blue")
  return_data <- 4*m/n #円周率の計算
  plot(x,y,col=cols,xlim=range,ylim=range,asp=1)
  #点をプロットする,xとyの範囲は上で定めたもの,x軸とy軸の比率を1にする
  return(return_data)
}
```

$n=1000$として結果を出力すると、図と出力結果は以下である。

```{r, echo=FALSE}
montecarlo(1000)
```


また、$n=10000$として結果を出力すると結果は以下である。

```{r, echo=FALSE}
montecarlo(10000)
```



