[R] Using R to illustrate the Central Limit Theorem

(Ted Harding) Ted.Harding at nessie.mcc.ac.uk
Thu Apr 21 20:48:05 CEST 2005


On 21-Apr-05 Paul Smith wrote:
> Dear All
> 
> I am totally new to R and I would like to know whether R is able and
> appropriate to illustrate to my students the Central Limit Theorem,
> using for instance 100 independent variables with uniform distribution
> and showing that their sum is a variable with an approximated normal
> distribution.
> 
> Thanks in advance,
> 
> Paul

Similar to Francisco's suggestion:

  m<-numeric(10000);
  for(k in (1:20)){
    for(i in(1:10000)){m[i]<-(mean(runif(k))-0.5)*sqrt(12*k)}
    hist(m,breaks=0.3*(-15:15),xlim=c(-4,4),main=sprintf("%d",k))
  }

(On my slowish laptop, this ticks over at a satidfactory rate,
about 1 plot per second. If your mahine is much faster, then
simply increase 10000 to a larger number.)

The real problem with demos like this, starting with the
uniform distribution, is that the result is, to the eye,
already approximately normal when k=3, and it's only out
in the tails that the improvement shows for larger values
of k.

This was in fact the way we used to simulate a normal
distribution in the old days: look up 3 numbers in
Kendall & Babington-Smith's "Tables of Random Sampling
Numbers", which are in effect pages full of integers
uniform on 00-99, and take their mean.

It's the one book I ever encountered which contained
absolutely no information -- at least, none that I ever
spotted.

A more dramatic illustration of the CLT effect might be
obtained if, instead of runif(k), you used rbinom(k,1,p)
for p > 0.5, say:

  m<-numeric(10000);
  p<-0.75; for(j in (1:50)){ k<-j*j
    for(i in(1:10000)){m[i]<-(mean(rbinom(k,1,p))-p)/sqrt(p*(1-p)/k)}
    hist(m,breaks=41,xlim=c(-4,4),main=sprintf("%d",k))
  }

Cheers,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 21-Apr-05                                       Time: 19:48:05
------------------------------ XFMail ------------------------------




More information about the R-help mailing list