ビュフォンの針をRのanimationで可視化モンテカルロ法

この記事は[R Advent Calendar 2015](http://qiita.com/advent-calendar/2015/r-rstudio)の3日目の記事です。

## はじめに

数値解析の分野では、Rを使用せずに違う言語で描画プログラムを書いた経験のある方が多数いらっしゃると思います。
Rには`animation`パッケージというものがあり、プロットを結合してアニメーションとして出力することができます。

今回はそれを使った数値解析・モンテカルロ法の話です。

## ビュフォンの針とは

[ビュフォンの針（Buffon's needle problem）](
https://ja.wikipedia.org/wiki/ビュフォンの針
)というのをご存知でしょうか。

長さ `a` の針を等間隔 `d` で引かれた平行線群のなかに落としたとき、針と平行線の交わる確率がどのようになるかという問題です。

<p align=center>
<img src=https://qiita-image-store.s3.amazonaws.com/0/3353/5b704aab-0693-2abb-2f21-ab81294f4d86.jpeg />
</p>

このとき

- 針の中心は平行線の存在する空間に一様に分布
- 針と平行線のなす角度の分布は区間 (0, π) 上の一様分布

を満たすことができれば

```math
\frac{2a}{\pi d}
```

となります。つまり、 `d=2a` とおけば


```math
\frac{2a}{\pi d} = \frac{1}{\pi}
```

となり、円周率（の逆数）を得ることができます。

※`a=l` と記述されるのが一般ですが、見難いため `a` にしています。


今回は、この問題をRで実装し、シミュレートして確認します。

## 実装

### Rコード

```r
buffon.needle = function(n = 100, sec = 0.1, l = 1.0, d = 1.0) {
  x <- y <- xd <- yd <- phi <- ctr <- numeric(0)
  sum <- 0; PI <- rep(NA, n)
  layout(matrix(c(1, 2), 2))
  for (j in 1:n) {
    # Draw parallel lines
    xlim <- c(-0.5 * l, 1.5 * l); ylim <- c(0, 2 * d)
    plot(1, xlim = xlim, ylim = ylim, type = 'n',
         xlab = '', ylab = '', axes = FALSE)
    abline(h = c(0.5 * d, 1.5 * d), lwd = 2, col = 'blue')
    
    # Simulate
    phi <- c(phi, runif(1, 0, pi))
    ctr <- c(ctr, runif(1, 0, 0.5 * d))
    x = c(x, runif(1, 0, l))
    y = c(y, sample(c(0.5 * d + ctr[j], 1.5 * d - ctr[j]), 1))
    xd = c(xd, 0.5 * l * cos(phi[j]))
    yd = c(yd, 0.5 * l * sin(phi[j]))
    
    # Draw Needle
    segments(x - xd, y - yd, x + xd, y + yd, col = 'gray')
    
    # Draw Mote Carlo
    xx = seq(0, pi, length = 200)
    plot(xx, 0.5 * l * sin(xx), ylim = c(0, 0.5 * d), type = 'l', bty = 'l',
         xlab = '', ylab = '', col = 'lightgray')
    idx = as.numeric(ctr > 0.5 * l * sin(phi)) + 1
    points(phi, ctr, col = c('red', 'blue')[idx])
    sum = sum + (ctr[j] <= 0.5 * l * sin(phi[j]))
    if (n > 0) {
      PI[j] = 2 * l * j/(d * sum)
    }
    ani.pause(sec)
  }
  
  return(PI)
}
```

`animation`の話だからといってやたらめったら変わったことはしていません。
今回は `ani.pause()` を使用して、段階的に表示することだけをしています。
実際の現場では `ani.start()`, `ani.stop()`, `ani.record()`, `ani.replay()` を使って録画することもあると思いますので、適したときに使用してください。

## 実行

animationパッケージを `library` か `require` でロードしたのち、ユーザ定義関数を呼び出します。

```r
> library('animation')
> buffon.needle()
```

`buffon.needle` 関数は引数を与えることもできます。

```r
> buffon.needle = function(n = 100, sec = 0.1, l = 1.0, d = 1.0)
```

## 結果

### アニメーション

針をプロットするのと、モンテカルロ法のアニメーションを見ることができます。

![output2.gif](https://qiita-image-store.s3.amazonaws.com/0/3353/cbc5bc89-fd73-db5b-a7cf-ae1b4c38e908.gif)

ビュフォンの針では下記の関数の積分値が解となるので、その積分＝面積をモンテカルロ法で近似解を得ている。

```math
\frac{L}{2}sin \theta
```

### 出力

`a=l=1`, `d=1`として関数を実行しているため、近似解は円周率に収束していくはずとなります。
その収束していく流れがわかるようにステップ毎の近似解を返却しています。

```r
> buffon.needle(30)
 [1] 2.000000 2.000000 2.000000 2.666667 2.500000 2.400000 2.800000 3.200000 3.000000 2.857143 3.142857
[12] 3.000000 2.888889 2.800000 2.727273 2.666667 2.833333 2.769231 2.714286 2.857143 2.800000 2.750000
[23] 2.705882 2.823529 2.941176 3.058824 3.176471 3.111111 3.052632 3.157895
```

## おわりに

Rというより数値解析の話が強かったかもしれませんが、数値解析方面をやっている方はどんどんRを使っていってほしいです。
僕も大学生のときにちょっと使っていた程度ですが、もっと使っていれば楽にできたことたくさんあったよなぁ。と今では思います。