Rで安定した男系男子継承に関する分析(ゴルトン･ワトソン過程)

library(dplyr)
library(knitr)

# ゴルトン・ワトソン過程 によるシミュレーションの関数
Galton_Watson_process <- function(NUM_CHILDREN) {
  # --- 1. パラメータ設定 ---
  PROB_MALE <- 105 / 205 # 男子が生まれる確率 (p)
  MAX_GENERATION <- 20 # 計算する最大世代数
  # 引数 NUM_CHILDREN は 1夫婦あたりの子供の数 (c)

  # --- 2. 確率計算の準備 ---
  # 結果を格納するベクトルを初期化
  # survival_prob[n] は第n世代で存続している確率を格納します
  survival_prob <- numeric(MAX_GENERATION)

  # 第1世代は定義により存続しているので、存続確率 = 1.0
  survival_prob[1] <- 1.0

  # --- 3. 漸化式による計算 (ループ処理) ---
  # 第2世代からMAX_GENERATIONまで、存続確率を順次計算
  for (n in 2:MAX_GENERATION) {
    # 前世代(n-1)の存続確率 S_{n-1} を取得
    s_prev <- survival_prob[n - 1]

    # 存続確率の漸化式 S_n = 1 - (1 - p * S_{n-1})^c を計算
    s_current <- 1 - (1 - PROB_MALE * s_prev)^NUM_CHILDREN

    # 計算結果を格納
    survival_prob[n] <- s_current
  }

  # --- 4. 結果の表示 ---
  # 世代と計算された存続確率をデータフレームにまとめる
  results_df <- data.frame(
    Generation = 1:MAX_GENERATION,
    Survival_Probability = survival_prob
  )

  # 存続確率をパーセンテージ形式に変換して追加
  results_df$Survival_Percentage <- sprintf("%.1f %%", results_df$Survival_Probability * 100)
  colnames(results_df)[c(1, 3)] <- c("世代", paste0("子供の数:", NUM_CHILDREN))

  return(results_df[c(1, 3)])
}

results_df <- lapply(seq(6), function(x) Galton_Watson_process(x)) %>% Reduce(function(x, y) plyr::join(x, y), .)

# kable()関数を使って、結果をテーブル形式で出力
kable(
  results_df,
  align = "r"
)