> pdf("F:\\Rmisc\\normal_stuff.pdf")
> 
> 
> ################# Functions of Normal Random Variables  ##########################
> 
> 
> set.seed(12345)              # set the seed for reproducible results
> num.sample <- 10000          # Number of samples
> size.sample <- 10            # Sample size for each sample
> z2sum <- rep(0,num.sample)   # Create vector to save sum of squared Zs
> zstd <- rep(0,num.sample)    # Create vector to save SD of Zs
> zmean <- rep(0,num.sample)   # Create vector to save mean of Zs
> 
> for (i in 1:num.sample) {
+ z <- rnorm(size.sample,0,1)       # Generate size.sample N(0,1) RVs
+ z2sum[i]  <- sum(z^2)             # Compute the sum of the squares of size.sample N(0,1) RVs
+ zstd[i]  <- sd(z)                 # Compute the SD of size.sample N(0,1) RVs
+ zmean[i] <- mean(z)               # Compute the mean of size.sample N(0,1) RVs
+ }
> 
> # Compute mean and variance of z2sum values (Theory => size.sample, 2*size.sample)
> mean(z2sum)    
[1] 10.00628
> var(z2sum)     
[1] 20.56895
> 
> par(mfrow=c(2,2))
> # Histogram of z2sum values with overlayed Chi-square(size.sample) distribution
> hist(z2sum,breaks=0:ceiling(max(z2sum)),main="",xlab="sum(Z^2)")
> xv <- seq(0,ceiling(max(z2sum)),0.01)
> yv <- dchisq(xv,size.sample)*num.sample
> lines(xv,yv)
> 
> # Compute mean and variance of (n-1)s^2/sigma^2 values (Theory => size.sample-1, 2*(size.sample-1))
> sampvar_scale <- (size.sample - 1)*(zstd^2)/1
> mean(sampvar_scale)
[1] 9.006212
> var(sampvar_scale)
[1] 18.50489
> 
> # Histogram of (n-1)s^2/sigma^2 values with overlayed Chi-square(size.sample-1) distribution
> hist(sampvar_scale,breaks=0:ceiling(max(sampvar_scale)),main="",xlab="(n-1)s^2/sigma^2")
> xv <- seq(0,ceiling(max(sampvar_scale)),0.01)
> yv <- dchisq(xv,size.sample-1)*num.sample
> lines(xv,yv)
> 
> # Compute mean and variance of (zmean-mu)/(sigma/sqrt(n))  (Theory => 0, 1) 
> zmean01 <- (zmean-0)/(1/sqrt(size.sample))
> mean(zmean01)
[1] 0.002582605
> var(zmean01)
[1] 1.000166
> 
> ########## Note that the histograms for z and t will use densities, not frequencies
> ###### due to scaling 
> 
> # Histogram of (zmean-mu)/(sigma/sqrt(n)) values with overlayed N(0,1) distribution 
> rangez <- ceiling(max(abs(zmean01)))
> hist(zmean01,breaks=seq(-rangez,rangez,0.1),main="",xlab="zmean01",freq=F)
> xv <- seq((-rangez),rangez,0.01)
> yv <- dnorm(xv,0,1)
> lines(xv,yv)
> 
> # Compute mean and variance of zmean01/sqrt(sampvar_scale/(n-1))   (Theory => 0, 9/7)
> tmean <- zmean/(zstd/sqrt(size.sample))
> mean(tmean)
[1] 0.005366197
> var(tmean)
[1] 1.29794
> 
> # Histogram of zmean01/sqrt(sampvar_scale/(n-1)) values with overlayed t(size.sample-1) distribution
> ranget <- ceiling(max(abs(tmean)))
> hist(tmean,breaks=seq(-ranget,ranget,0.1),main="",xlab="tmean",freq=F)
> xv <- seq((-ranget),ranget,0.01)
> yv <- dt(xv,size.sample-1)
> lines(xv,yv)
> 
> 
> dev.off()
null device 
          1 
>