clt.illustration <- function(n=30, N=10000, from = NA){
if(is.na(from)) {
graphics.off() # clean graphical device
from <- 1 # we start with 1 random variable
}
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")[1]
x.max <- par("usr")[2]
# we generate 500 equidistant points on this range
x <- seq(x.min, x.max, length.out = 500)
range.between.2.cut.offs <- h$breaks[2]-h$breaks[1]
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
}
}