# How to run and plot simulation in R of sum of 20 random variables?

I would like to run a simulation where the sum of 20 random variables are generated 1000 times and plotted in a histogram.

I generate the random numbers using runif(20,0,1) then sum these random numbers using the sum function, but I was wondering if there is an easier method to do this 1000 then to write a for loop.

-
Try help(replicate). –  chl Jul 22 '11 at 15:12
thanks! it worked –  DBR Jul 22 '11 at 15:17

X = matrix(runif(20000), 1000, 20)
S1 = apply(X, 1, sum)
S2 = rowSums(X)
hist(S1) ## same as hist(S2)


Are two ways to do this without a for loop. X contains 1000 rows, each containing 20 uniform(0,1) variables. S1 and S2 will both contain the same data; 1000 independent realizations of the sum of 20 independent uniform(0,1) random variables.

-
Thanks a lot! a further related question: I would like to run 100 simulations of x-y, where x has a mean of 60.9 and sd of 2.9 and y has mean63.7 and sd2.7. do I need to use the qnorm function? –  DBR Jul 22 '11 at 15:40
What distribution do x,y have? If it's normal then you could do X = rnorm(100, mean=60.9, sd=2.9); Y = rnorm(100, mean=63.7, sd=2.7); Z=X-Y; hist(Z). qnorm() tells you the quantiles (percentiles) of the normal distribution so I don't think that is what you want. –  Macro Jul 22 '11 at 15:42
@Marco good answer. Worth mentioning that by generating 20000 uniform random variables you incur overhead of a single call to runif(), instead of repeatedly calling that function to produce 20 random variables at a time. The same numbers are generated but doing it your way is recommended and far more efficient. Also, rowSums() is preferred over the apply() version as it makes use of C code optimised for that particular task whereas apply() needs to operate more generically. –  Gavin Simpson Jul 22 '11 at 15:53
@Macro But I would like to run 1000 simulations of x-y where 100 samples are re-drawn every time from the two distributions specified. –  DBR Jul 22 '11 at 16:22
@DBR, well, if your variables are normally distributed, by additivity of normal distributions, $X-Y$ is normal with mean $60.9-63.7$ and standard deviation $\sqrt{ 2.9^{2} + 2.7^{2} }$ so you can just use the same trick I specified in my answer except using the normal instead of uniform and the mean/sd I specified here. –  Macro Jul 22 '11 at 16:29
y <- replicate(1000,sum(runif(20,0,1)))