Python による同名の記事を書いた後，R ではもっと簡潔に書けるだろうと思ったので，やってみた。
確かに簡潔になるけど，読む人は苦労するだろうなあ。
n = 1000000 # 取りあえずの試行回数（実際はこれより少なくなる）
d = 10 # 平行線の間隔
l = d / 2 # 針の長さ
# 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sine.theta を生成する
xy = matrix(runif(2*n, 0, 1), 2) # 一辺 1 の正方形内の一様乱数による座標点
r = colSums(xy^2) # 原点からの距離
xy = xy[, r <= 1] # 半径 1 の円の内部にある点に限る
r = r[r <= 1] # 半径 1 の円の内部にある点に限る
sine.theta = xy[2,] /sqrt(r) # 三角関数の定義により
n = length(r) # sine.theta の個数
y = runif(n, 0, l) # 同じ個数の[0, l) の一様乱数（針の中心から，平行線までの距離）
2 * l / d / mean(y <= l / 2 * sine.theta) # 条件を満たすものの逆数が π の推定値
