# A simulation of asteroid impacts.  Asteroids are modeled as a normal
# distribution N(m, 1) with an unknown mean.  Impacts above a certain threshold are
# deadly, everything else is harmless.
# Even though observations above the threshold are impossible, the survivors
# are still able to get an unbiased estimate of mean asteroid size.
sim <- function(n_earths, n_impacts, extinction_thresh=1) {
  means <- rnorm(n_earths)
  impacts <- matrix(
    rnorm(
      n_earths * n_impacts,
      rep(means, each = n_impacts)),
    byrow = TRUE,
    nrow = n_earths)

  alive <- apply(impacts, 1, max) < extinction_thresh
  sample_means <- rowMeans(impacts)
  # The survivors use a N(0, 1) prior.
  post_means <- sample_means / (1 / n_impacts + 1)
  biases <- (post_means - means)[alive]

  plot(density(biases), main = "Bias density")
  x <- seq(min(biases), max(biases), length.out = 100)
  y <- dnorm(x, 0, sd(biases))
  lines(x, y, lty = "dashed", col = "red")
  abline(v = 0, lty = "dotted")
}

sim(1000, 3) # Should be roughly normal, centered at 0.