R version 2.5.0 (2007-04-23)
Copyright (C) 2007 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pdf("clt.pdf")
> # central limit theorem simulation
> # population: uniform(a, b) with mean=(a+b)/2; var=(b-a)^2/12
> # draw a random sample of n=36 
> a=0; b=9;  # uniform(a=0, b=9) so mu=4.5 and sigma^2=81/12=6.75, sigma=2.60
> n=36;  # sample size n=36
> x=runif(n, min=a, max=b)
> stem(x)  # stem plot

  The decimal point is at the |

  0 | 247
  1 | 233458
  2 | 137
  3 | 399
  4 | 26789
  5 | 0125899
  6 | 25
  7 | 037
  8 | 0136

> mean(x)  # sample mean
[1] 4.354623
> sd(x) # sample sd
[1] 2.495644
> 
> # do it again, noting that the sample mean and sd change
> x=runif(n, min=a, max=b)
> stem(x)  # stem plot

  The decimal point is at the |

  0 | 114599
  1 | 557
  2 | 117
  3 | 6778
  4 | 1259
  5 | 234
  6 | 567
  7 | 04699
  8 | 14445

> mean(x)  # sample mean
[1] 4.505213
> sd(x) # sample sd
[1] 2.815882
> 
> # do it again, noting that the sample mean and sd change
> x=runif(n, min=a, max=b)
> stem(x)  # stem plot

  The decimal point is at the |

  0 | 2
  1 | 00458
  2 | 222579
  3 | 09
  4 | 036
  5 | 026779
  6 | 78
  7 | 123566
  8 | 1123
  9 | 0

> mean(x)  # sample mean
[1] 4.830192
> sd(x) # sample sd
[1] 2.58289
> 
> # repeat the procedure K times, save the K sample means in xbar
> K=1000
> xbar = rep(0, K)
> for(i in 1:K) { x=runif(n, min=a, max=b); xbar[i] = mean(x) }
> # now examine the distribution of K sample means (xbar). 
> stem(xbar, scale=4)

  The decimal point is 1 digit(s) to the left of the |

  31 | 2
  31 | 
  32 | 
  32 | 
  33 | 1
  33 | 69
  34 | 0
  34 | 899
  35 | 124
  35 | 67789
  36 | 013
  36 | 557888
  37 | 0000001123344
  37 | 5556666677788888999
  38 | 0111112222333333344
  38 | 55566666677889
  39 | 00000011111112333333334444
  39 | 556678888999
  40 | 00000000111122223333444444
  40 | 555555555666667777788889999999
  41 | 0000001122222223333344
  41 | 55556677777777777788888888888899999
  42 | 000000011111111222333333333444444444
  42 | 55555555555556666666777777788888888899999999
  43 | 00000000000000000001111111222222222223333333444444444444
  43 | 555666666777777777888888888999999
  44 | 00000000000000111122222222222222223333333333344444444
  44 | 5555555556666666677777888888888889999999
  45 | 0000000000111111111122222222223333333333344444444
  45 | 555555555555555666666677888888999999
  46 | 000000000011111222222222222333333334444444444444
  46 | 55555555566666667777777888888888899999999999
  47 | 00000000111112222333333344444444
  47 | 5555556666666667777777788889999999
  48 | 000000000011111122222233334444444
  48 | 5555555566666666666777777777788888888999999999
  49 | 000000011122222223333334444444
  49 | 555555555556666667778888888999999
  50 | 0000001111122222233344
  50 | 555556666666677788999999
  51 | 0000001122333444
  51 | 55556678999
  52 | 13333
  52 | 55555666789
  53 | 034
  53 | 667788
  54 | 00223
  54 | 59
  55 | 01
  55 | 9
  56 | 124
  56 | 
  57 | 0

> hist(xbar)  #  approximately normal
> mean(xbar)
[1] 4.487563
> sd(xbar)
[1] 0.4221001
> 
> # look at the Normal Probability Plot
> qqnorm(xbar)  # note the vertical axis is sample x
> 
> # what is the expected value of sd(xbar)?
> sqrt(6.75/n)
[1] 0.4330127
> 
> # first create a function to make the plot as in the textbook
> qqnorm2 = function(x)
+ {
+ temp=qqnorm(x, plot=F)
+ plot(temp$y, temp$x, xlab="Sample  x(j)", ylab="Normal quantile zj", main="Normal Probability Plot")
+ }
> # the normal probability plot  as in the textbook
> qqnorm2(xbar)   # close to a straight line
> 
> 
> dev.off()
null device 
          1 
>