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:
- I get two full pills. As on the first morning, I split one and place the unused half back into the bottle.
- I get one full pill and one half-pill, both of which I consume.
- I get two half-pills. In this case, I take out another pill at random. If it’s a half-pill, then I consume all three halves. But if it’s a full pill, I split it and place the unused half back in the bottle.
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)
}
}
}
}
}
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.
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.