p        = .06
n_events = 40
t_max    = 1000
n_pop    = 1e3
t_death  = matrix(nrow = n_pop, ncol = 1)
for(i in 1:n_pop){
  dat        = replicate(n_events, sample(0:1, t_max, replace = T, prob = c(1-p, p)))
  t_death[i] = which(apply(apply(dat, 2, function(x) cumsum(x) > 0), 1, all))[1]
}

hist(t_death, breaks = seq(0, max(t_death) + 5, by = 5),
     freq = F,  col = "Grey", xlab = "Age",
     main = c("Age at Death", paste0("Mean = ", mean(t_death))))

lines(2:max(t_death), diff((1 - (1 - p)^(1:max(t_death)))^n_events), lwd = 3, col = "Red")