# How to compute simulated differences of means repetitively?

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$.

-

Let's go for the one-line solution:

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

-
Yes, your solution is clearly more aesthetic than mine (+1). I just wanted to keep it simple ;-) –  ocram Nov 10 '11 at 11:00
sorry chi... not thinking in R... feel free to roll that back. –  John Nov 10 '11 at 15:18
Much clean R code, thanks @John. –  chl Nov 10 '11 at 16:12
yes, it looks better but now '-' is replicated instead of vectorized... I suppose the performance shouldn't matter much... you're welcome? :) –  John Nov 10 '11 at 16:14
What about rnorm(1000,mean=69.5-63.9,sd=sqrt((2.9)^2/100+(2.7)^2/100)) ? –  Xi'an Nov 11 '11 at 6:52
show 1 more comment

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).

-
(+1) Those are good comments, indeed. (I was thinking in terms of bootstrap when replying, even if this is not the question nor did I resample the same dataset.) –  chl Nov 10 '11 at 13:46

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 :-)

-
i would like to substract y from x –  Dbr Nov 10 '11 at 10:53
Oh, the minus sign was not clear (before @chl's edit). Well, this is a simple adaptation now.... –  ocram Nov 10 '11 at 10:59
(+1) Keeping it simple has several merits, especially from a pedagogical point of view. –  chl Nov 10 '11 at 12:05