How to compute simulated differences of means repetitively?

Question

I would like to run a 1000 simulations of $X-Y$ where $X$ is the average of 100 random observations with mean of 69.5 and SD 2.9 and $Y$ is the average of 100 random observations with mean of 63.9 and SD 2.7.

I did the following:

x <- rnorm(100, 69.5, 2.9)
y <- rnorm(100, 63.9, 2.7)

And now I would like to somehow run a 1000 simulations of $X-Y$.

Answer 1

Let's go for the one-line solution:

replicate(1000, mean(rnorm(100, 69.5, 2.9)) - mean(rnorm(100, 63.9, 2.7)))

Answer 2

To say something about the math of it rather than the programming: it would be useful to think about whether you can actually avoid the repeated simulations altogether in this case by using elementary properties of the Gaussian normal distribution. In this toy example it doesn't matter, but in general it's vastly more efficient to exploit an analytic solution if it exists. In this case here, the code for the solution becomes much shorter than even chl's example.

Hint 1: What is the distribution of the average of 100 normal random variables with a mean of 69.5 and SD 2.9?

Hint 2: Apply hint 1 to the average of 1000 instantiations of random variable X - Y, for which the closed form is also elementary (if X and Y are normal).

Answer 3

Do you mean something like that?

x <- numeric(1000)
y <- numeric(1000)
for(i in 1:1000)
{
  x[i] <- mean(rnorm(n=100, mean=69.5, sd=2.9))
  y[i] <- mean(rnorm(n=100, mean=63.9, sd=2.7))
}

My trial gives

Of note, the normal distribution was used... but it was an arbitrary choice... Following this assumption, I've superimposed the density of $\bar{X}$ and that of $\bar{Y}$ on the histograms. Can you find these densities? In particular, they show that the simulations work pretty well :-)