# GdA 20/1/2016

# distribution of averages 
mu=0             # or rnorm(1) if you like     
m=100            # nr of averages  (set m=3 to understand the varius steps,
                 #                  printing the variables after each step)
N=1000           # size of a sample 
medie=rep(0,m)   # empty vector
for(i in 1:m) medie[i] = mean(rnorm(N, mu))      # unitary sigma 

m.medie = mean(medie)    # mean of the means 
s.medie = sd(medie)      # standard deviation of the means

cat( sprintf("sigma medie = %f;   1/sqrt(N) = %f\n", s.medie, 1/sqrt(N)) )

# differences of means
m.d <- outer(medie, medie, '-')         # matrix of differences
ut <- upper.tri(m.d)                    # matrix of logical variables
utr <- as.vector(upper.tri(m.d) * m.d)  # vector in which only the upper triangle
                                        # elements are different from zero
diff <- utr[ut]                         # only upper triangular elements are kept
sigma.diff <- sd(diff)                  # standard deviation of all differences
cat( sprintf("sigma diff medie = %f;   1/sqrt(N) * sqrt(2) = %f\n",
             sigma.diff, 1/sqrt(N)*sqrt(2)) )

# you can also check the mean of the differences -> \approx 0
# given any mean, all others are around it with sigma = 1/sqrt(N) * sqrt(2) !