円周率を推定するのに円周率を使いたくない

• ビュフォンの針を使った円周率の計算機推定の話題(こちら)がある。
• 簡単に言うと、等長の針を、針の長さより長い幅の縞模様の上に投げて、針が１つの縞に納まるか、２つの縞に跨るかの割合から円周率が推定できるという話です。
• （１）針の頭の位置を一様乱数で発生し、針の向き(角度)を0-2πの一様乱数で発生させれば、針の位置は決まりますから、シミュレーションは成功です。ですが、針の向きの一様乱数発生に円周率を使っているので、これは「ずる」です。

# 針の長さをLとする
# n本の針
# 針の頭が[0,1] x [0,1]の正方形に一様に分布するように発生する
# 針の方向は[0,2pi) に一様に分布するように発生する
# x軸方向の座標の小数点以下切り捨ての値が同じか異なるかを確認する
# 切り捨て値で２分した割合を計算する

L <- 5/6
n <- 10^4
head.x <- runif(n)
head.y <- runif(n)
# 方向一様乱数を発生するのにpiを利用するのを許せば…
theta <- runif(n)*2*pi
# 三角関数を使うことも許せば…
tail.x <- head.x + L*cos(theta)

fl.head <- floor(head.x)
fl.tail <- floor(tail.x)
check <- fl.head==fl.tail
table(check)
# ２縞に跨った本数 n-sum(check) が占める割合と針の長さL（縞の幅に対する相対的長さ）とから次の式でpiの推定値が出る
2 * L * n / (n-sum(check))

• 別の方法。針の頭の位置を基準として、針先の相対位置が決まればよいのです。それは、単位円周上の一様分布が発生できればよいです。方向に関して一様な分布がランダムに作れれば、それを単位半径に縮めさえすればよいので、結局、方向一様な2次元分布を原点を中心に作れればよいことがわかります。x,yの両方に正規乱数を与えると2次元正規分布の乱点が作れます。

# 方向一様乱数を発生するのにpiを使うのが嫌なので
# 単位円周上の一様乱数をpiを陽に使わずに発生してみる
# 2次元正規分布は全方向対称な分布なので…
needle <- cbind(rnorm(n),rnorm(n))
plot(needle)
# 円周上にするには…
r <- sqrt(needle[,1]^2+needle[,2]^2)
needle.1 <- cbind(needle[,1]/r,needle[,2]/r)
plot(needle.1)

tail.x <- head.x + L * needle.1[,1]
fl.head <- floor(head.x)
fl.tail <- floor(tail.x)
check <- fl.head==fl.tail
table(check)
2 * L * n / (n-sum(check))

• （３）正規乱数を作る、とは言っても、正規分布の関数表現には、指数関数やpiが現れますから、これも「ずる」な感じです。方向一様な2次元分布を作ればよいので、[-1,1] x [-1,1]という正方形に一様乱数を発生させ、そのうち、半径1以下の点だけを取れば、単位円板上の一様乱点が得られますので、これを単位円周に移動させればよいですね。

# 正規乱数を発生するのも計算機的にずるだと思えば
# 正方形の乱数を発生させ、そのうち、単位円内に含まれる点だけを取り出し
# それを単位円周に移動すればよい

X <- runif(n,min=-1,max=1)
Y <- runif(n,min=-1,max=1)
plot(X,Y)
R <- sqrt(X^2+Y^2)
ok <- which(R<=1)
X.ok <- X[ok]
Y.ok <- Y[ok]
plot(X.ok,Y.ok)

# 後は、正規乱数の場合と同じ
needle <- cbind(X.ok,Y.ok)
r <- sqrt(needle[,1]^2+needle[,2]^2)
needle.1 <- cbind(needle[,1]/r,needle[,2]/r)
plot(needle.1)

tail.x <- head.x + L * needle.1[,1]
fl.head <- floor(head.x)
fl.tail <- floor(tail.x)
check <- fl.head==fl.tail
table(check)
2 * L * n / (n-sum(check))