Solution to fivethirtyeight Riddler’s puzzle, Can You Find Your Pills? (August 19, 2022).

I’ve been prescribed to take 1.5 pills of a certain medication every day for 10 days, so I have a bottle with 15 pills. Each morning, I take two pills out of the bottle at random.

On the first morning, these are guaranteed to be two full pills. I consume one of them, split the other in half using a precision blade, consume half of that second pill, and place the remaining half back into the bottle.

On subsequent mornings when I take out two pills, there are three possibilities:

Assume that each pill — whether it is a full pill or a half-pill — is equally likely to be taken out of the bottle.

On the 10th day, I again take out two pills and consume them. In a rush, I immediately throw the bottle in the trash before bothering to check whether I had just consumed full pills or half-pills. What’s the probability that I took the full dosage, meaning I don’t have to dig through the trash for a remaining half-pill?

Since rules were laid out, I will brute-force the simulation to return the final pill composition, the remaining sum.

sim_medications <- function(target = 1.5, bottle_size = 15) {
  bottle <- rep(1, bottle_size)
  end_date <- bottle_size / target

  for (day in 1:end_date) {
    if (day == 1) {
      bottle <- sample(bottle, bottle_size - 2, replace = FALSE)
      bottle <- append(bottle, 0.5)
    } else {
      two_pills <- sample(bottle, 2, replace = FALSE)
      if (day == end_date) {
        return(sum(two_pills))
      }

      for (i in seq_along(two_pills)) {
        bottle <- bottle[-match(two_pills[i], bottle)]
      }

      diff <- sum(two_pills) - target

      if (diff == 0) {
        next
      } else if (diff > 0) {
        bottle <- append(bottle, 0.5)
      } else {
        extra_pill <- sample(bottle, 1, replace = FALSE)
        three_pills <- c(two_pills, extra_pill)
        bottle <- bottle[-match(extra_pill, bottle)]

        diff <- sum(three_pills) - target

        if (diff > 0) {
          bottle <- append(bottle, 0.5)
        }
      }
    }
  }
}

Answer

sims <- replicate(1e5, sim_medications())
table(sims) / 1e5
## sims
##       1     1.5 
## 0.21067 0.78933

The answer to the Riddler’s question is 0.79.

Additional Notes

The above probability will approach to 0.5 as the bottle size increases.

bottle_sizes <- seq(15, 450, by = 15)

sims <- sapply(bottle_sizes, function(x) {
  replicate(
    1000,
    sim_medications(bottle_size = x)
  ) |> table() / 1000
})

df_prop_sims <- data.frame(bottle_size = bottle_sizes, t(sims))

df_prop_sims_long <- pivot_longer(df_prop_sims, cols = -bottle_size, names_to = "two_pills", values_to = "proportion") |>
  mutate(two_pills = parse_number(two_pills) |> factor())

ggplot(df_prop_sims_long, aes(x = bottle_size, y = proportion, col = two_pills)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Probability of Success by Bottle Size and Sum of Two Pills",
    x = "Bottle Size",
    y = "Probability",
    col = "Sum of Two Pills"
  )

This is because the limiting distribution reflects the accumulation of half pills as the bottle size increases. This also suggests that as the distribution of possible sums becomes more spread out, the probability of getting the target sum of two pills becomes more uncertain.