モンテカルロ法で π を求める

pi を，試行回数 n = 10, 100, 1000, ..., 1000000 として，モンテカルロ法で求めよという単純な問題。

ずっと前にここにも書いたが，R では基本的に
mean(1 >= colSums(matrix(runif(2*n)^2, 2)))*4
の 1 行で書ける（複数の関数を使うが）

junk = sapply(10^(1:6), function(n) {
  pi0 = mean(1 >= colSums(matrix(runif(2*n)^2, 2)))*4
  cat(sprintf("n = %7i pi = %.9f error = %10.7f\n", n, pi0, pi-pi0))
}

n =      10 pi = 3.200000000 error = -0.0584073
n =     100 pi = 3.120000000 error =  0.0215927
n =    1000 pi = 3.252000000 error = -0.1104073
n =   10000 pi = 3.150000000 error = -0.0084073
n =  100000 pi = 3.142480000 error = -0.0008873
n = 1000000 pi = 3.144192000 error = -0.0025993

なお，
「ランダムな点を増やしていくと、段々と Πに近づいていくはずです」
というのは，一見正しそうであるが，一時的に遠ざかることもあるということで，実行例のような n が小さいときの値の方が誤差が少ないということもある（あたりまえなんだ）。

また，
「精度は単精度で構いません」
とは，時代遅れ。今時，単精度でプログラムすることなんかないだろう。