Rの関数:pbirthday {stats}

Rの関数から pbirthday {stats} を確認します。

関数 pbirthday とは

pbirthday は、確率論においてよく知られる 「誕生日のパラドックス（Birthday Paradox）」 およびその一般化された問題の確率を計算する関数です。

誕生日のパラドックスとは、「部屋にいる集団の中に、誕生日が同じである2人の組が存在する確率は、直感よりも高い」という現象を指します。例えば、23人の集団であれば、確率は約50%を超えます。

この関数は、単に「ペア（2人）の一致」だけでなく、「3人が同じ誕生日である確率」や、「1年が365日ではない場合（カテゴリ数が異なる場合）」など、問題を一般化して計算することができます。

関数 pbirthday の引数

args(pbirthday)
function (n, classes = 365, coincident = 2) 
NULL
  1. n

    • 集団のサイズ（人数、サンプルの数）。
    • 整数または数値ベクトル。
    • 誕生日問題における「部屋にいる人数」に相当します。この人数が増えるほど、一致する確率は上昇します。
  2. classes

    • カテゴリの総数。
    • デフォルト: 365。
    • 誕生日問題では「1年の日数」に相当します。誕生日以外にも、例えば「サイコロをn回振って同じ目がk回出る確率」を計算したい場合は 6 を指定します。ハッシュ衝突の確率計算などにも応用可能です。
  3. coincident

    • 一致が必要な最小人数（\(k\)）。
    • デフォルト: 2。
    • 2: 「少なくとも1組のペア（2人）」が同じ誕生日である確率を計算します（標準的な誕生日問題）。
    • 3: 「少なくとも1組のトリオ（3人）」が同じ誕生日である確率を計算します。
    • 一般に、集団内に「同じカテゴリに属する個体が \(k\) 人以上いるグループ」が少なくとも1つ存在する確率を求めます。

シミュレーションコード

以下に、pbirthday の理論値が実際の現象と一致しているかを検証するためのシミュレーションコードを示します。

本シミュレーションでは、以下の2つのシナリオについて、人数（\(n\)）を増やしながら確率がどのように上昇していくかを確認します。

  1. 標準的な誕生日問題 (\(k=2\)):

    • 同じ誕生日の2人組がいる確率。
  2. 3人一致の問題 (\(k=3\)):

    • 同じ誕生日の3人組がいる確率。

乱数を用いたモンテカルロ法による結果と、pbirthday 関数による計算値を比較します。

シミュレーションの設定と実行

# パッケージの読み込み
library(ggplot2)
library(tidyr)

# 乱数シードの固定
seed <- 20260302
set.seed(seed)

# 誕生日のパラドックス：シミュレーションと理論値の比較

# --- シミュレーション設定 ---
# 人数 n を 5人から 100人まで変化させる
n_vec <- seq(5, 100, by = 5)
n_sim <- 2000 # 各人数における試行回数
days_in_year <- 365

# 結果格納用のデータフレーム
results <- data.frame()

# --- シミュレーション実行 ---
for (n in n_vec) {
  # 1. モンテカルロ・シミュレーション
  # n人の誕生日をランダムに生成し、重複をチェックする試行を n_sim 回繰り返す

  # k=2 (ペア) のチェック
  count_pair <- 0
  # k=3 (トリオ) のチェック
  count_trio <- 0

  for (i in 1:n_sim) {
    # 1〜365の乱数をn個生成 (復元抽出)
    birthdays <- sample(1:days_in_year, n, replace = TRUE)

    # 度数分布を計算
    tab <- table(birthdays)

    # 最大重複数を取得 (max_overlap)
    max_overlap <- max(tab)

    if (max_overlap >= 2) count_pair <- count_pair + 1
    if (max_overlap >= 3) count_trio <- count_trio + 1
  }

  prob_sim_k2 <- count_pair / n_sim
  prob_sim_k3 <- count_trio / n_sim

  # 2. pbirthday関数による理論値計算
  prob_theo_k2 <- pbirthday(n, classes = days_in_year, coincident = 2)
  prob_theo_k3 <- pbirthday(n, classes = days_in_year, coincident = 3)

  # 結果の保存
  results <- rbind(results, data.frame(
    N = n,
    Type = "k=2 (ペア)",
    Simulation = prob_sim_k2,
    Theoretical = prob_theo_k2
  ))

  results <- rbind(results, data.frame(
    N = n,
    Type = "k=3 (トリオ)",
    Simulation = prob_sim_k3,
    Theoretical = prob_theo_k3
  ))
}

cat("--- 計算結果の一部を確認 ---\n\n")
cat("--- ペアの場合 ---\n")
print(results[results$N %in% c(20, 25, 30, 40, 50, 80, 90) & results$Type == "k=2 (ペア)", ])

cat("\n--- トリオの場合 ---\n")
print(results[results$N %in% c(20, 25, 30, 40, 50, 80, 90) & results$Type == "k=3 (トリオ)", ])
--- 計算結果の一部を確認 ---

--- ペアの場合 ---
    N       Type Simulation Theoretical
7  20 k=2 (ペア)     0.4115   0.4114384
9  25 k=2 (ペア)     0.5610   0.5686997
11 30 k=2 (ペア)     0.7150   0.7063162
15 40 k=2 (ペア)     0.8825   0.8912318
19 50 k=2 (ペア)     0.9715   0.9703736
31 80 k=2 (ペア)     1.0000   0.9999143
35 90 k=2 (ペア)     1.0000   0.9999938

--- トリオの場合 ---
    N         Type Simulation Theoretical
8  20 k=3 (トリオ)     0.0055 0.009560094
10 25 k=3 (トリオ)     0.0155 0.018399727
12 30 k=3 (トリオ)     0.0285 0.031265743
16 40 k=3 (トリオ)     0.0715 0.071120021
20 50 k=3 (トリオ)     0.1360 0.131678634
32 80 k=3 (トリオ)     0.4230 0.419738066
36 90 k=3 (トリオ)     0.5365 0.532106840

ペアの場合、N=25の時点で、トリオの場合、N=90の時点で確率は50%を超えています。

データの整形と可視化

# --- データの整形と可視化 ---
# ggplot用にデータを縦長に変換
df_plot <- pivot_longer(results,
  cols = c("Simulation", "Theoretical"),
  names_to = "Method", values_to = "Probability"
)

p1 <- ggplot(df_plot, aes(x = N, y = Probability)) +
  # 理論値のライン
  geom_line(
    data = subset(df_plot, Method == "Theoretical"),
    aes(color = Type, linetype = "理論値 (pbirthday)"), linewidth = 1
  ) +
  # シミュレーション値のポイント
  geom_point(
    data = subset(df_plot, Method == "Simulation"),
    aes(color = Type, shape = "シミュレーション値"), size = 2.5, alpha = 0.7
  ) +

  # 確率0.5のライン
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +

  # スケール設定
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +

  # 色と凡例の設定
  scale_color_manual(values = c("k=2 (ペア)" = "blue", "k=3 (トリオ)" = "red")) +
  scale_linetype_manual(name = "", values = c("理論値 (pbirthday)" = "solid")) +
  scale_shape_manual(name = "", values = c("シミュレーション値" = 16)) +
  labs(
    title = "誕生日のパラドックス (理論値とシミュレーション)",
    x = "人数 (n)",
    y = "確率"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom", legend.box = "vertical")

print(p1)
Figure 1

Figure 1 から人数が増えるにつれて、一致する確率が上昇していく過程を確認できます。

k=2 (青線) の場合、23人程度で確率が0.5を超え、50人を超えるとほぼ確実(>0.97)になります。

k=3 (赤線) の場合、3人が同じ誕生日になる確率は、ペアの場合よりも緩やかな上昇になりますが、それでも88人程度で確率が0.5を超えます。

なお、関数 qbirthday を利用して、確率から集団のサイズ(n)を求めることが可能です。

cat("--- 一致する確率が50%になる集団サイズ ---\n\n")
cat("--- ペア(k=2)の場合 ---\n")
qbirthday(prob = 0.5, classes = days_in_year, coincident = 2)
cat("\n--- トリオ(k=3)の場合 ---\n")
qbirthday(prob = 0.5, classes = days_in_year, coincident = 3)
--- 一致する確率が50%になる集団サイズ ---

--- ペア(k=2)の場合 ---
[1] 23

--- トリオ(k=3)の場合 ---
[1] 88

以上です。