I got my hands on the 1978 JASA paper The Randomized Play-the-Winner Rule in Medical Trials and saw that the authors analyzed an "inverse" stopping rule:
"We continue taking observations until either $S_A(n) + F_B(n)$ or $S_B(n) + F_A(n)$ is equal to some given value $r$."
Above, $S_i(n)$ and $F_i(n)$ are the number of sucesses and failures observed for treatment $i$ after trial $n$, for both $i=A, B$.
They talk about how you can view the sequential experiment on a 2d lattice, with state and transition probabilities, and they prove that the probability of "correct selection" (chosing the better treatment) tends towards 1 as r tends towards infinity. Still, I think they had to run a simulation with different values of $r$ and the sucess probabilities. (This was the late 70s. I'd be really curious to learn whether this was analytical or a Monte Carlo experiment using a computer.)
In the table they present below, $p_A$ and $p_B$ are the success probabilities of treatment A and B, respectively. $Pr(CS)$ is the probability of correct stopping (note A is always superior) and $ASN$ is the average sample number at which $r$ was reached. Ignore the "PW" columns, and note that when I said "Play the Winner" in my original post, I meant "Randomized Play the Winner" (RPW).
Below is code that implements a Monte Carlo experiment to replicate these numbers. Again, I'm curious if there's a more elegant way to get these. Before posting the code, I'll conclude with the procedure:
- Choose r based on probability of correct selection and average sample number, given hypothetical $p_A$ and $p_B$.
- Run the Randomized Play the Winner strategy until one treatment has $r$ balls in the urn.
simulate_trial <- function(p_a, p_b, stop_value) {
# Initial counts
r_a <- 0
r_b <- 0
# Keep track of the history
history_a <- c()
history_b <- c()
while (r_a < stop_value && r_b < stop_value) {
# Calculate the current probability to draw treatment A
prob_a <- r_a / (r_a + r_b)
# If r_a and r_b are both 0, we need to randomly select a treatment since prob_a will be NaN
if (is.nan(prob_a)) {
prob_a <- 0.5
}
# Randomly draw a treatment
treatment_drawn <- ifelse(runif(1) < prob_a, "A", "B")
if (treatment_drawn == "A") {
# Check for success or failure
if (runif(1) < p_a) {
# Success for treatment A
r_a <- r_a + 1
} else {
# Failure for treatment A
r_b <- r_b + 1
}
} else {
# Check for success or failure
if (runif(1) < p_b) {
# Success for treatment B
r_b <- r_b + 1
} else {
# Failure for treatment B
r_a <- r_a + 1
}
}
# Append to history
history_a <- c(history_a, r_a)
history_b <- c(history_b, r_b)
}
return(list(history_a = history_a, history_b = history_b))
}
# Simulation parameters
grid <- merge(
data.frame(r=c(5, 10, 15, 30, 50)),
data.frame(p_a=c(.3, .7, .9, .5, .8, .7), p_b = c(.1, .5, .7, .1, .4, .1)),
all=TRUE
)[, c(2, 3, 1)]
M <- 100000 # Monte Carlo reps
for (j in 1:nrow(grid)) {
p_a <- grid[j, "p_a"]
p_b <- grid[j, "p_b"]
r <- grid[j, "r"]
cat("running with p_a =", p_a, ", p_b =", p_b, ", r =", r, "for", M, "reps\n")
correct_stopping <- numeric(M)
sample_size <- numeric(M)
for (rep in 1:M) {
result <- simulate_trial(p_a, p_b, r)
correct_stopping[rep] <- max(result$history_a) == r
sample_size[rep] <- length(result$history_a)
}
grid[j, "pr_correct_stopping"] <- mean(correct_stopping)
grid[j, "average_sample_number"] <- mean(sample_size)
}
print(grid, digits=3)
