RでA/Bテスト:ベイズ的アプローチ

# --- シナリオ設定 ---
# 「真の」CVR
# そもそも知り得ないパラメータ
# つまりB案勝利のシミュレーションであり、
# 「もしB案の真のCVRがA案より2パーセントポイント高かったとしたら、今回のテスト手法はそれを正しく検知できるか？」のシミュレーションになります。
true_cvr_A <- 0.10
true_cvr_B <- 0.12

# 事前分布 (Beta分布のパラメータ alpha, beta)
# 「CVRは10%くらいだろう」という弱い信念を表現 (平均 1/(1+9) = 0.1)
alpha_prior <- 1
beta_prior <- 9

# 1日あたりの訪問者数 (各グループ)
visitors_per_day <- 1500

# --- 意思決定ルール ---
# BがAを上回る確率の閾値
prob_threshold <- 0.95

seed <- 20250613
set.seed(seed)

# 毎日データを観測し、信念（事後分布）を更新していくプロセスをシミュレートします。

# 結果を記録するためのデータフレーム
results_log <- data.frame()

# 事後分布のパラメータを事前分布で初期化
alpha_A <- alpha_prior
beta_A <- beta_prior
alpha_B <- alpha_prior
beta_B <- beta_prior

cat("シミュレーション開始...\n")

# シミュレーションループ (最大30日間)
for (day in 1:30) {
  # 今日のデータを生成 (二項分布に従う)
  conversions_A_today <- rbinom(1, visitors_per_day, true_cvr_A)
  conversions_B_today <- rbinom(1, visitors_per_day, true_cvr_B)

  # 事後分布のパラメータを更新
  # 事後分布 ∝ 尤度 × 事前分布
  # Beta分布は二項分布の共役事前分布
  # 事前分布 Beta(alpha_prior, beta_prior)
  # 観測結果 n 回の試行のうち k 回成功
  # 事後分布 Beta(alpha_prior + k, beta_prior + (n - k))

  alpha_A <- alpha_A + conversions_A_today
  beta_A <- beta_A + (visitors_per_day - conversions_A_today)

  alpha_B <- alpha_B + conversions_B_today
  beta_B <- beta_B + (visitors_per_day - conversions_B_today)

  # --- 意思決定ルールの評価 ---
  # モンテカルロシミュレーションで事後分布から大量にサンプリング
  n_samples <- 100000
  posterior_samples_A <- rbeta(n_samples, alpha_A, beta_A)
  posterior_samples_B <- rbeta(n_samples, alpha_B, beta_B)

  # BがAを上回る確率 P(CVR_B > CVR_A) を計算
  prob_B_beats_A <- mean(posterior_samples_B > posterior_samples_A)

  # 期待損失を計算 (Bを選んだ場合の間違った意思決定による損失)
  loss_if_choose_B <- pmax(0, posterior_samples_A - posterior_samples_B)
  expected_loss_B <- mean(loss_if_choose_B)

  # ログに記録
  current_log <- data.frame(
    day = day,
    prob_B_beats_A = prob_B_beats_A,
    expected_loss_B = expected_loss_B
  )
  results_log <- rbind(results_log, current_log)

  # 終了条件のチェック
  if (prob_B_beats_A > prob_threshold) {
    cat(paste0("テスト終了 (", day, "日目): \n"))
    cat(paste0("  BがAを上回る確率: ", round(prob_B_beats_A, 3), " ( > 閾値 ", prob_threshold, ")\n"))
    break
  }

  # 最終日まで終了しなかった場合
  if (day == 30) {
    cat("30日間テストを継続しましたが、明確な結論は出ませんでした。\n")
  }
}

cat("\n--- 最終的な結論 ---\n")

# 最終日のデータを取得
final_day_stats <- results_log[nrow(results_log), ]
total_visitors <- visitors_per_day * final_day_stats$day

# AとBの事後分布の要約統計量
summary_A <- quantile(posterior_samples_A, probs = c(0.025, 0.5, 0.975))
summary_B <- quantile(posterior_samples_B, probs = c(0.025, 0.5, 0.975))

cat(paste0("テストは ", day, " 日間で終了しました。\n"))
cat(paste0("総訪問者数は各グループ ", total_visitors, " 人です。\n\n"))

cat(paste0("新しいデザインBのCVRが既存のAを上回る確率は ", round(final_day_stats$prob_B_beats_A * 100, 1), "% です。\n"))
cat("これは我々が設定した95%の信頼度閾値を超えています。\n\n")

cat("A（既存）のCVRの95%信用区間:", paste0(round(summary_A[1] * 100, 1), "% ~ ", round(summary_A[3] * 100, 1), "% (中央値: ", round(summary_A[2] * 100, 1), "%)\n"))
cat("B（新規）のCVRの95%信用区間:", paste0(round(summary_B[1] * 100, 1), "% ~ ", round(summary_B[3] * 100, 1), "% (中央値: ", round(summary_B[2] * 100, 1), "%)\n\n"))
cat(paste0("B案を採用した場合の期待損失は ", round(final_day_stats$expected_loss_B * 100, 4), " パーセントポイントと非常に小さく、リスクは低いと判断できます。\n"))
cat("ビジネス判断: 新しいデザインのボタンBを全面的に採用することを強く推奨します。\n")