Advertisement
Guest User

Untitled

a guest
Jul 27th, 2023
83
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 1.05 KB | None | 0 0
  1. # A simulation of asteroid impacts.  Asteroids are modeled as a normal
  2. # distribution N(m, 1) with an unknown mean.  Impacts above a certain threshold are
  3. # deadly, everything else is harmless.
  4. # Even though observations above the threshold are impossible, the survivors
  5. # are still able to get an unbiased estimate of mean asteroid size.
  6. sim <- function(n_earths, n_impacts, extinction_thresh=1) {
  7.   means <- rnorm(n_earths)
  8.   impacts <- matrix(
  9.     rnorm(
  10.       n_earths * n_impacts,
  11.       rep(means, each = n_impacts)),
  12.     byrow = TRUE,
  13.     nrow = n_earths)
  14.  
  15.   alive <- apply(impacts, 1, max) < extinction_thresh
  16.   sample_means <- rowMeans(impacts)
  17.   # The survivors use a N(0, 1) prior.
  18.   post_means <- sample_means / (1 / n_impacts + 1)
  19.   biases <- (post_means - means)[alive]
  20.  
  21.   plot(density(biases), main = "Bias density")
  22.   x <- seq(min(biases), max(biases), length.out = 100)
  23.   y <- dnorm(x, 0, sd(biases))
  24.   lines(x, y, lty = "dashed", col = "red")
  25.   abline(v = 0, lty = "dotted")
  26. }
  27.  
  28. sim(1000, 3) # Should be roughly normal, centered at 0.
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement