Rで円周率計算

Rでモンテカルロ法（大量に試行して精度の高い結果を得ようとする方法）によって、円周率を計算することを行った。

set.seed(1)
points.num <- 1000
coords <- matrix(runif(2*points.num),ncol=2)
dist.from.origin <- apply(coords,1,function(x){sqrt(sum(x^2))})
myPI <- sum(dist.from.origin<1)/points.num*4

コードを説明する。set.seed(1)は乱数のシードを固定するためのものである。毎回同じ結果を得るためだけに利用している。points.numはシミュレーションに用いる点の数である。coordsはpoints.num*2の行列であり、1行が1組の(x,y)座標を表している．座標は[0,1]×[0,1]の領域で一様乱数を発生させて得られた点になる。dist.from.originは、原点からの距離を計算した結果である。定義どおりに距離を計算しているが、最終的に原点からの距離が1未満の点の数を求めるだけなので、実はsqrtは不要である。0以上の実数aに対し、a*a<1とsqrt(a*a)<1は同値だからである。myPIが円周率の近似解になる。円周率は原点からの距離が1未満である点の割合を4倍したものになる。4倍するのは、一様乱数を発生させた領域が第1象限のみに限られているからである。

上記のコードを走らせるとmyPIは3.148という値が得られる。実際の円周率3.141592...に近い値が得られていることが確認できる。精度を上げようとしてpoints.numを大きくすると、急激に時間がかかるようになる。そこで、次回はsnowパッケージを用いて、並列処理をさせてみることにしよう。

追記(2011/06/08) 並列処理のコードを実装しました。こちらよりご覧ください。

追記(2011/11/09) 上のコードをいじった高速版のコードを以下に挙げておきます。@a_bickyさんのツイートと、とあるブログからの指摘がありました。applyとrowSumsの速度比較結果についてはまた後日。

set.seed(1)
points.num <- 1000
coords <- matrix(runif(2*points.num),ncol=2)
dist.from.origin <- rowSums(coords^2) # 上で説明した通りsqrtは不要
myPI <- sum(dist.from.origin<1)/points.num*4

追記(2011/12/11)
applyとrowSumsの速度比較を行ってみました。詳細はこちら。