πを明示的には使わない，ビュフォンの針のシミュレーション

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) # 条件を満たすものの逆数が π の推定値

言語により演算子の優先順位には違いがある

R では，

> 123 * 345 %% 10
[1] 615
> (123 * 345) %% 10
[1] 5
> 123 * (345 %% 10)
[1] 615

Python では，

>>> 123 * 345 % 10
5
>>> (123 * 345) % 10
5
>>> 123 * (345 % 10)
615

ということ

次の 13 日の金曜日は？

ツェラーの公式を使う。ただ，Python のときと違いがあるので注意。

Zeller = function(y, m, d) {
 w = c("Sat", "Sun", "Mon", "Tue", "Wed", "Thu", "Fri")
 if (m == 1 || m == 2) {
  m = m + 12
  y = y - 1
 }
 C = y%/%100
 Y = y%%100
 D = (d + (26 * (m + 1))%/%10 + Y + Y%/%4 - 2 * C + C%/%4)%%7 + 1
 return(w[D])
}

for (y in 2020:2030) {
 for (m in 1:12) {
  if (Zeller(y, m, 13) == "Fri") {
   print(paste0(y, "年", m, "月"))
  }
 }
}

[1] "2020年3月"
[1] "2020年11月"
[1] "2021年8月"
[1] "2022年5月"
[1] "2023年1月"
[1] "2023年10月"
[1] "2024年9月"
[1] "2024年12月"
[1] "2025年6月"
[1] "2026年2月"
[1] "2026年3月"
[1] "2026年11月"
[1] "2027年8月"
[1] "2028年10月"
[1] "2029年4月"
[1] "2029年7月"
[1] "2030年9月"
[1] "2030年12月"

Python の場合と違うところ。大事！！

(26 * (m + 1))%/%10 を 26 * (m + 1) %/% 10 としてしまうと間違い

 

アベの餌食にならなくて本当によかったね！！

本日放送 2020/03/05

ニンゲン観察！モニタリング 　TBS

最後の卒業式

合唱　ゆず　も　一緒に

早く卒業式をやって，ほんとうによかった

日本全国で，想い出に残る卒業式を企画していたたくさんの卒業生の気持ちを思うと，怒りの念がふつふつと沸き起こる