information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その4です。
TokyoR発表資料
その1 「"The R Foundation for Statistical Computing"を実感する」
その2 「Rの統計関連(dpqr)関数チートシート」
その3 「Rの統計関連(dpqr)関数の書き方」
誕生日のパラドックスとは
$n$人いるときに，その中に同じ誕生日である組が存在する確率です。ただし、閏年は考慮せず、1年は365日とします。何がパラドックスなのかは後ほど。
解答例(余事象から考える）
求める確率を$P(A)$とする。まず余事象$1 - P(A) = P(A^c)$を考える。
例えば2人がある部屋にいるとする。その2人が同じ誕生日でない確率は$\dfrac{364}{365}$
その後部屋に入ってきた3人目が前にいた2人と同じ誕生日でない確率は$\dfrac{363}{365}$
よって、部屋にいる3人の誕生日が違う確率は$\dfrac{364}{365}\times \dfrac{363}{365}$
これを汎用化して考えると、$n$人が部屋にいて同じ誕生日である二人組が存在しない確率は
P(A^c) = \dfrac{{}_{364} P_{n-1}}{365^{n-1}}
となる。ここで分子について計算すると、
{}_{364} P_{n-1} = \dfrac{364!}{\{364-(n-1)\}!} = \dfrac{364!}{(365-n)!}
よって、
P(A^c) = \dfrac{{}_{364} P_{n-1}}{365^{n-1}} = \dfrac{364!}{365^{n-1}(365-n)!}
$n$が大きい場合は、計算がとんでもないことになるのを私は知っている（←偉そう）。そこで対数を取って計算すると、
\ln P(A^c) = \ln \left\{ \dfrac{364!}{365^{n-1}(365-n)!} \right\}= \ln 364! - (n-1)\ln365 - \ln \{(365-n)! \}
$\ln(364!) = \ln \Gamma(365)$、$\ln {(365-n)!} = \ln \Gamma(366-n)$なので、
\ln P(A^c) = \ln \Gamma(365) - \ln \Gamma(366-n) + (1-n)\ln365
Rでのコーディング例
最後の式で求められた$\ln P(A^c)$を指数変換して$P(A^c)$（$n$人がいて同じ誕生日の組み合わせが１つもない確率）を求める。そして、真に求める確率は$P(A) = 1 - P(A^c)$なので、以下のようにコーディングできる。
library(conflicted)
library(tidyverse)
conflicts_prefer(dplyr::filter)
ln_gamma_365 <- lgamma(365)
ln_365 <- log(365)
df <- tibble(n = 2:365) |>
mutate(ln_p_c = ln_gamma_365 - lgamma(366 - n) + (1 - n) * ln_365) |>
mutate(p_c = exp(ln_p_c)) |>
mutate(p = 1 - p_c)
> head(df)
# A tibble: 6 × 4
n ln_p_c p_c p
<int> <dbl> <dbl> <dbl>
1 2 -0.00274 0.997 0.00274
2 3 -0.00824 0.992 0.00820
3 4 -0.0165 0.984 0.0164
4 5 -0.0275 0.973 0.0271
5 6 -0.0413 0.960 0.0405
6 7 -0.0579 0.944 0.0562
ln_p_c列が$\ln P(A^c)$の値、
p_c列は$P(A^c)$の値、
p列が求める値$P(A)$です。
誕生日のパラドックスの解を求めるRの関数
♪こんな計算しなくても、Rにはpbirthday()関数が用意されています。$n=2,n=3$のときの確率はそれぞれ下記のとおりで、私の自作した計算コードとも一致しています。
> pbirthday(2, classes = 365, coincident = 2)
[1] 0.002739726
> pbirthday(3, classes = 365, coincident = 2)
[1] 0.008204166
pbirthday()関数を使えば、先ほど私が作ったデータフレームも簡単に計算できます。せっかくデータフレームにしたので、グラフも作ってみます。
library(conflicted)
library(tidyverse)
conflicts_prefer(dplyr::filter)
df <- tibble(n = 2:365) |>
mutate(p = map_dbl(n, ~ pbirthday(.x, classes = 365, coincident = 2)))
g <- ggplot(data = filter(.data = df, p < 0.999), aes(x = n, y = p)) +
geom_line(linewidth = 1.5)
ggsave("birthday.png",
plot = g,
width = 2039,
height = 1447,
units = c("px"))
このグラフを見ると、50人を超えたところですでに確率が限りなく1に近くなっているのがわかります。
dpqrチートシートの記事を読まれた方はピンとくると思いますが、累積確率$p$のときの$n$を求める関数もちゃんと用意されていて、それがqbirthday()関数です。
> qbirthday(prob = 0.5, classes = 365, coincident = 2)
[1] 23
> qbirthday(prob = 0.99, classes = 365, coincident = 2)
[1] 57
確率が0.5を超えるのが23人のとき、0.99を超えるのはなんと57人のときです。直感的にはそんなように思えませんが、その人間の直感との乖離を「パラドックス」と呼んでいるのです。
こんなところまでカバーしているRの関数
Rの素晴らしいところは、baseのRに誕生日のパラドックスの計算関数まであるということです。しかも3人や4人の組がいる場合の計算もできます（coincident =を3や4に設定すれば良い）し、誕生日でなくても同様の計算もできます(classes=の値を365以外に設定できる)。やはり、
The R Foundation for Statistical Computing
ですね。