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 proceedure:

1. Choose r based on probability of correct selection and average sample number, given hypothetical $p_A$ and $p_B$.
2. 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)