RでA/Bテスト:多腕バンディット問題

# 各バナー（腕）の真のクリック率を設定
# Arm 3が最も高い(0.15)が、アルゴリズムはこの事実を知らない
true_probs <- c(0.05, 0.08, 0.15, 0.11, 0.03)
arms <- length(true_probs)
trials <- 2000

# ε-Greedy 法の実装
epsilon_greedy <- function(true_probs, trials, epsilon = 0.1) {
  arms <- length(true_probs)
  counts <- rep(0, arms) # 各腕を引いた回数
  values <- rep(0, arms) # 各腕の報酬の合計

  total_reward <- 0
  rewards_history <- numeric(trials) # 試行ごとの累積報酬を記録

  for (t in 1:trials) {
    if (runif(n = 1, min = 0, max = 1) < epsilon) {
      # 探索(Exploration): ランダムに腕を選ぶ
      chosen_arm <- sample(1:arms, 1)
    } else {
      # 活用(Exploitation): 現時点で平均報酬が最大の腕を選ぶ
      # 平均が0の場合や同率一位がある場合も考慮し、ランダムに一つ選ぶ
      avg_values <- ifelse(counts > 0, values / counts, 0)
      best_arms <- which(avg_values == max(avg_values))
      chosen_arm <- sample(best_arms, 1)
    }

    # 選んだ腕で報酬を得る (ベルヌーイ試行)
    reward <- rbinom(1, 1, prob = true_probs[chosen_arm])

    # 結果を更新
    counts[chosen_arm] <- counts[chosen_arm] + 1
    values[chosen_arm] <- values[chosen_arm] + reward
    total_reward <- total_reward + reward
    rewards_history[t] <- total_reward
  }

  return(list(total_reward = total_reward, rewards_history = rewards_history, counts = counts))
}

# UCB 法の実装
ucb <- function(true_probs, trials, c_param = 2) {
  arms <- length(true_probs)
  counts <- rep(0, arms)
  values <- rep(0, arms)

  total_reward <- 0
  rewards_history <- numeric(trials)

  # 最初に各腕を1回ずつ試す
  for (i in 1:arms) {
    reward <- rbinom(1, 1, prob = true_probs[i])
    counts[i] <- 1
    values[i] <- reward
    total_reward <- total_reward + reward
    rewards_history[i] <- total_reward
  }

  for (t in (arms + 1):trials) {
    avg_values <- values / counts
    # UCBスコアを計算
    ucb_scores <- avg_values + sqrt(c_param * log(t) / counts)

    # UCBスコアが最大の腕を選ぶ
    best_arms <- which(ucb_scores == max(ucb_scores))
    chosen_arm <- sample(best_arms, 1)

    # 報酬を得て更新
    reward <- rbinom(1, 1, prob = true_probs[chosen_arm])
    counts[chosen_arm] <- counts[chosen_arm] + 1
    values[chosen_arm] <- values[chosen_arm] + reward
    total_reward <- total_reward + reward
    rewards_history[t] <- total_reward
  }

  return(list(total_reward = total_reward, rewards_history = rewards_history, counts = counts))
}

# トンプソンサンプリングの実装
thompson_sampling <- function(true_probs, trials) {
  arms <- length(true_probs)
  # ベータ分布のパラメータ (alpha:成功+1, beta:失敗+1)
  alpha <- rep(1, arms)
  beta <- rep(1, arms)

  total_reward <- 0
  rewards_history <- numeric(trials)

  for (t in 1:trials) {
    # 各腕の現在のベータ分布から値をサンプリング
    sampled_probs <- rbeta(arms, alpha, beta)

    # サンプル値が最大の腕を選ぶ
    best_arms <- which(sampled_probs == max(sampled_probs))
    chosen_arm <- sample(best_arms, 1)

    # 報酬を得る
    reward <- rbinom(1, 1, prob = true_probs[chosen_arm])

    # ベータ分布のパラメータを更新
    if (reward == 1) {
      alpha[chosen_arm] <- alpha[chosen_arm] + 1
    } else {
      beta[chosen_arm] <- beta[chosen_arm] + 1
    }

    total_reward <- total_reward + reward
    rewards_history[t] <- total_reward
  }

  return(list(total_reward = total_reward, rewards_history = rewards_history, counts = alpha + beta - 2))
}

# シミュレーション設定
num_simulations <- 100

# シナリオ1: 差が明確な場合
true_probs_easy <- c(0.05, 0.08, 0.15, 0.11, 0.03)

# シナリオ2: 差が僅かな場合
true_probs_hard <- c(0.10, 0.12, 0.15, 0.14, 0.09)

# シミュレーション実行関数
run_simulations <- function(probs) {
  # replicateを使って各アルゴリズムを100回実行
  results_greedy <- replicate(num_simulations, epsilon_greedy(probs, trials)$total_reward, simplify = FALSE)
  results_ucb <- replicate(num_simulations, ucb(probs, trials)$total_reward, simplify = FALSE)
  results_ts <- replicate(num_simulations, thompson_sampling(probs, trials)$total_reward, simplify = FALSE)

  # 結果をデータフレームにまとめる
  results_df <- data.frame(
    Simulation = 1:num_simulations,
    EpsilonGreedy = unlist(results_greedy),
    UCB = unlist(results_ucb),
    ThompsonSampling = unlist(results_ts)
  )
  return(results_df)
}

seed <- 20250615
set.seed(seed)

# 各シナリオでシミュレーションを実行
results_easy_df <- run_simulations(true_probs_easy)
results_hard_df <- run_simulations(true_probs_hard)

library(dplyr)

# 平均総報酬を計算
avg_rewards_easy <- results_easy_df %>%
  summarise(across(c(EpsilonGreedy, UCB, ThompsonSampling), mean))

avg_rewards_hard <- results_hard_df %>%
  summarise(across(c(EpsilonGreedy, UCB, ThompsonSampling), mean))

cat("\n--- シナリオ1 (差が明確な場合) の平均総報酬 ---\n")
print(avg_rewards_easy)

cat("\n--- シナリオ2 (差が僅かな場合) の平均総報酬 ---\n")
print(avg_rewards_hard)