# set sample size and number of samples set.seed(2345) n <- 2 reps <- 10000 # perform random sampling samples <- replicate(reps, rnorm(n,2.5,1.118)) # 2 x 1000 sample matrix # compute sample variances sample_var <- rep(NA, reps) for (i in 1:reps) { sample_var[i] <- var(samples[,i]) } # check that 'sample_var' is a vector is.vector(sample_var) mean(sample_var) var(sample_var) hist(sample_var, freq=F, col="grey", xlab="", xlim=c(-1, 6), breaks=100) par(new=T) plot(density(sample_var), axes=F, main="", xlim=c(-1, 6), lwd=2, col="blue")