誕生日のパラドックスをRの関数で求める方法

:::note info
information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その4です。
TokyoR発表資料
[その1 「"The R Foundation for Statistical Computing"を実感する」](https://qiita.com/hige_megane/private/e49ef15222bba6be14d7)
[その2 「Rの統計関連(dpqr)関数チートシート」](https://qiita.com/hige_megane/private/8383b8926cb97f5040cf)
[その3 「Rの統計関連(dpqr)関数の書き方」](https://qiita.com/hige_megane/private/416c0d6b5742dfcf4050)
:::

## 誕生日のパラドックスとは

$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$人が部屋にいて同じ誕生日である二人組が存在しない確率は

```math
P(A^c) = \dfrac{{}_{364} P_{n-1}}{365^{n-1}}
```
となる。ここで分子について計算すると、

```math
{}_{364} P_{n-1} = \dfrac{364!}{\{364-(n-1)\}!} = \dfrac{364!}{(365-n)!}
```

よって、
```math
P(A^c) = \dfrac{{}_{364} P_{n-1}}{365^{n-1}} = \dfrac{364!}{365^{n-1}(365-n)!}
```

$n$が大きい場合は、計算がとんでもないことになるのを私は知っている（←偉そう）。そこで対数を取って計算すると、

```math
\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)$なので、

```math
\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)$なので、以下のようにコーディングできる。

```birthday_paradox_cul.r
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.r
> pbirthday(2, classes = 365, coincident = 2)
[1] 0.002739726
> pbirthday(3, classes = 365, coincident = 2)
[1] 0.008204166
```

`pbirthday()`関数を使えば、先ほど私が作ったデータフレームも簡単に計算できます。せっかくデータフレームにしたので、グラフも作ってみます。
```birthday.r
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"))
```
![birthday.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1015152/7b647493-a219-4e44-a8f2-ba683fdb1f4f.png)

このグラフを見ると、50人を超えたところですでに確率が限りなく1に近くなっているのがわかります。
dpqrチートシートの記事を読まれた方はピンとくると思いますが、累積確率$p$のときの$n$を求める関数もちゃんと用意されていて、それが`qbirthday()`関数です。

```qbirthday.r
> 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

ですね。

# ENJOY!