Solution to 538’s aug25/2017 Riddler Express

timothy leffel (@lefft)

---

here is the problem (original post from 538 here):

You take half of a vitamin every morning. The vitamins are sold in a bottle of 100 (whole) tablets, so at first you have to cut the tablets in half. Every day you randomly pull one thing from the bottle — if it’s a whole tablet, you cut it in half and put the leftover half back in the bottle. If it’s a half-tablet, you take the vitamin. You just bought a fresh bottle. How many days, on average, will it be before you pull a half-tablet out of the bottle?

Extra credit: What if the halves are less likely to come up than the full tablets? They are smaller, after all.

---

note: extra credit is the case where weighted=TRUE

1. draw a pill from our bottle, with or without weights.

• weighted=TRUE means that full pills are inherently twice as likely to be drawn as half pills
• weighted=FALSE means that full pills and half pills have the same inherent likelihood of being drawn

draw_pill <- function(bottle, weighted){
  if (!weighted){
    weight_full <- sum(bottle=="full_pill")
    weight_half <- sum(bottle=="half_pill") 
    weights <- c(weight_full, weight_half) / (weight_full + weight_half)
  } else {
    weight_full <- sum(bottle=="full_pill") * 2 
    weight_half <- sum(bottle=="half_pill") 
    weights <- c(weight_full, weight_half) / (weight_full + weight_half)
  }
  return(sample(c("full_pill", "half_pill"), size=1, prob=weights))
}

2. after drawing a pill, update the bottle

• if you draw a full, put a half back in the bottle
• if you draw a half, just remove the half

update_bottle <- function(bottle, draw, day_num){
  full_pills_left <- sum(bottle=="full_pill")
  half_pills_left <- sum(bottle=="half_pill")
  
  if (draw=="full_pill"){
    bottle <- c(rep("full_pill", times=full_pills_left-1), 
                rep("half_pill", times=half_pills_left+1))
    return(bottle)
  }
  if (draw=="half_pill"){
    bottle <- c(rep("full_pill", times=full_pills_left), 
                rep("half_pill", times=half_pills_left-1))
    return(bottle)
  }
}

3. simulate pill-taking

• until you draw a half, then stop.
• return the number of days it took.

simulate_pill_taking <- function(weighted){
  my_bottle <- rep("full_pill", times=100)
  day_num <- 1
  draw <- draw_pill(my_bottle, weighted=weighted)
  while (draw != "half_pill"){
    draw <- draw_pill(my_bottle, weighted=weighted)
    my_bottle <- update_bottle(bottle=my_bottle, draw=draw, day_num=day_num)
    day_num <- day_num + 1
  }
  # uncomment if you want to see how many fulls are left after day `day_num`
  # message(paste0(sum(my_bottle=="full_pill"), " full pills left"))
  return(day_num)
}

4. run the simulation for weighted and unweighted, show results

• 10000 times apiece
• display the average number of days it takes
• then plot the results

num_sims <- 10000

# 10000 simulations, without weights
days_till_half_unweighted <- 
  replicate(n=num_sims, simulate_pill_taking(weighted=FALSE))

# 10000 simulations, with weights
days_till_half_weighted <- 
  replicate(n=num_sims, simulate_pill_taking(weighted=TRUE))

# collect the results
sim_data <- data.frame(
  num_days = c(days_till_half_weighted, days_till_half_unweighted), 
  sim_type = rep(c("weighted", "unweighted"), each=num_sims)
)

here’s the averages:

message(paste0(
  "if we're equally likely to draw a full or a half, it takes an average of ",
  round(mean(days_till_half_unweighted), 1), " days to draw a half pill"
))

## if we're equally likely to draw a full or a half, it takes an average of 14.1 days to draw a half pill

message(paste0(
  "if we're twice as likely to draw a full as a half, it takes an average of ",
  round(mean(days_till_half_weighted), 1), " days to draw a half pill"
))

## if we're twice as likely to draw a full as a half, it takes an average of 18.8 days to draw a half pill

here’s a density plot showing the distribution of days it takes to get a half pill, with and without weights.

ggplot(sim_data, aes(color=sim_type, x=num_days)) + geom_density()

aaaand a boxplot too, for the types who appreciate straining their eyes/minds.

ggplot(sim_data, aes(x=sim_type, y=num_days)) + geom_boxplot()