### Preliminars

The runif function generates draws from a uniform distribution:
runif(5, min = 2, max = 3)

##  2.686 2.684 2.099 2.413 2.590

The rnorm function generates draws from a normal distribution:
rnorm(5, mean = 0, sd = 1)

##  -2.13356  0.20484 -1.59384 -0.07009  1.03166

The density function computes kernel density estimates:
# estimation of the density of the uniform distribution
plot(density(runif(10000, min = 2, max = 3))) # estimation of the density of the uniform distribution
plot(density(rnorm(10000, mean = 0, sd = 1))) ### Detailed code

The following function performs an animated plot of the mean of $k$ uniform random variables, with $k$ varying from $1$ to $n$, and using $N$ individuals. The code was initially written by Bill Venables, updated by Ted Harding, and then by Emmanuel Rousseaux.
clt.illustration <- function(n=30, N=10000, from = NA){
if(is.na(from)) {
graphics.off() # clean graphical device
}
par(mfrow = c(1,2), pty = "s") # split the screen in 2, for plotting the density
# and the QQ-plot
for(k in from:n) {
# we store our k r.v. in a matrix with N rows and k columns
k.rv.unif <- matrix(runif(N*k), N, k)
# we compute the mean for each row to have the estimation of the mean r. v.
rv.mean <- rowMeans(k.rv.unif)
# centering and reduction
rv.mean.centered.reduced <- (rv.mean - 0.5)*sqrt(12*k) # Bienaymé formula
# we compute an histogram
h <- hist(
rv.mean.centered.reduced,
breaks = "FD", # we use Freedman-Diaconis algorithm to estimate cut-offs (a
# range-fixed cutting)
plot=FALSE
)
# we plot the histogram
if(k == 1) {s <- ''} else {s <- 's'}
plot(h,xlim = c(-4,4),main = paste0(k, ' random variable', s),
xlab = "Centered and reduced mean R.V. ",
ylab = "Frequency"
)
# then we want to plot a centered reduced normal law
# we first need to get min and max bounds on x abscissa
x.min <- par("usr")
x.max <- par("usr")
# we generate 500 equidistant points on this range
x <- seq(x.min, x.max, length.out = 500)
range.between.2.cut.offs <- h$breaks-h$breaks
y <- N*dnorm(x)*(range.between.2.cut.offs) #dnorm compute the value of the
# normal density function with mean = 0 and sd = 1 by default
lines(
x, # abscissa
y, # ordinate
col = "red"
)
# we plot the QQ-plot
qqnorm(
rv.mean.centered.reduced,
ylim = c(-4,4),
xlim = c(-4,4),
pch = ".",
col = "blue"
)
abline(0, 1, col = "red")
Sys.sleep(1) # 1s sleep to have time to watch
}
}


### Illustration

Here are snapshots for some values of $k$.   