40人の中に、同じ誕生日の人がいる確率を求めよう。 ただし、

1. R を使って計算する

(1) 1組も同じ誕生日がいない場合の数

誕生日が重複しないような選び方

  • 1人目は365日のどれを選んでもよい
  • 2人目は1人目が選んだ誕生日以外の364日から選ぶ
  • 3人目は2人目までが選んだ誕生日以外の363日から選ぶ
  • \(\cdots\)
  • 40人目は39人目までが選んだ誕生日以外の326日から選ぶ

よって、組み合わせは \(365 \times 364 \times \cdots \times 326\) 通り。

これをRで計算しよう。

まず、365から（from = 365)）1ずつ数が減る (by = -1)、長さ（要素の数）が40 (length.out = 40) の数列を作り、A という名前をつけて保存する（現時点ではコードの意味は理解する必要はない。第1クォータの終わりまでに理解すればよい）。

A <- seq(from = 365, by = -1, length.out = 40)

このAの中身を確認してみよう。

A
##  [1] 365 364 363 362 361 360 359 358 357 356 355 354 353 352 351 350 349 348 347
## [20] 346 345 344 343 342 341 340 339 338 337 336 335 334 333 332 331 330 329 328
## [39] 327 326

必要な要素が揃ったので、このAの中身をすべて掛け合わせればよい。 Rで次のコードを実行すると、かけあわせた結果がBとして保存される。

B <- prod(A)

Bの中身を確認してみよう。

B
## [1] 3.374555e+101

e+101 というのは、\(\times 10^{101}\)という意味である。

(2) すべての組み合わせ

40人それぞれについて、誕生日の選び方は365通りあるので、365を40回かければよい。つまり、365の40乗を求めればよい。Rで次のコマンドを実行すると、365の40乗がCとして保存される。

C <- 365 ^ 40

Cの中身を確認してみよう。

C
## [1] 3.10252e+102

これは、\(3.10252 \times 10^{102}\)という意味である。つまり、3.10252の小数点を102桁右に動かした数を表している。

(3) 40人の中に1組も同じ誕生日がいない確率

上で求めたB（1組も同じ誕生日がいない場合の数）をC（すべての組み合わせ）で割れば、1組も同じ誕生日のペアがいない確率を求められる。この結果を Dとして保存してみよう。

D <- B / C

Dの中身を確認してみよう。

D
## [1] 0.1087682

(4) 40人の中に少なくとも1組は同じ誕生日のペアがいる確率

40人の中に少なくとも1組は同じ誕生日のペアがいる確率は、すべての確率（1 = 100%）から、1組も同じ誕生日のペアがいない確率 (D) を引けば求められる。よって、Rでこの確率を求めると、

1 - D
## [1] 0.8912318

となる。

このように、Rを使いこなせれば、手計算で求めるのがとても面倒な問題でも、楽に答えを求めることができる。

以上の内容を1行で書くと、次のように書ける。

1 - prod(seq(from = 365, by = -1, length.out = 40) / 365)
## [1] 0.8912318

2. Rでもっと簡単に答えを求める

Rには、様々な関数 (functions) があらかじめ用意されており、目的にふさわしい関数を利用すると、答えを簡単に求めることができる。

誕生日問題を解くための関数として、pbirthday() という関数があるので、この関数を使って答えを求めてみよう。この関数では、全体の人数 (n) が何人いるかを指定すれば、その中に少なくとも1組は同じ誕生日のペアがいる確率を求めてくれる。私たちは40人の場合について知りたいので、n=40を指定する。以下のRコードを実行してみよう。

pbirthday(n = 40)
## [1] 0.8912318

先ほどと同じ答えが出た。

このように、Rの関数を使いこなせるようになると、さらに簡単に答えを求めることができる。

3. 可視化

最後に、全体の人数が変わると確率がどう変わるか可視化して確かめてみよう。 クラスの人数 \(n\) を1から50まで動かすと、同じ誕生日の人がいる確率 \(p\)はどのように変わるだろうか？

n <- 1:50
p <- sapply(n, pbirthday)
plot(n, p, type = "l", col = "royalblue", lwd = 2)



授業の内容に戻る