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) var(z2sum) 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) var(sampvar_scale) # 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) var(zmean01) ########## 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, (n-1)/(n-3)=9/7) tmean <- zmean/(zstd/sqrt(size.sample)) mean(tmean) var(tmean) # 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()