Example R programs and commands Relationship between normal, chi-squared, and F distributions # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. # All other indented lines are the R program output. ########## Chi squared is the sum of standard normals M <- 1000; # large number of samples, for good histograms xisq11 <- vector(length=M); # to hold M samples of Chi-sq, df=11 for(i in 1:M) # M samples of ... xisq11[i] <- sum(rnorm(11)^2); # ...11 squared std.normals summed # Plot the histogram of the samples with narrow bins: tmax <- max(xisq11); # largest value among the samples tmax <- ceiling(tmax); # least integer greater than tmax hist(xisq11,freq=FALSE,breaks=seq(0,tmax,1)); # plot the histogram # Compare with plot of Chi-squared on 11 degrees of freedom: tmax <- max(xisq11); # largest value among the samples t <- seq(0,tmax,by=0.1); # plot range lines(t,dchisq(t,df=11)); # plot of chi-squared on 11 degrees of freedom ########## F density squared is the ratio of mean Chi-squareds M <- 1000 # large number of samples for good histograms # Generate the numerator samples: ndf <- 18; # numerator degrees of freedom ximsqnum <- vector(length=M); # to hold M samples of mean Chi-sq for(i in 1:M) # M samples of ... ximsqnum[i] <- mean(rnorm(ndf)^2); # ...ndf squared std.normals averaged # Generate the denominator samples: ddf <- 7; # denominator degrees of freedom ximsqden <- vector(length=M); # to hold M samples of mean Chi-sq for(i in 1:M) # M samples of ... ximsqden[i] <- mean(rnorm(ddf)^2); # ...ddf squared std.normals averaged # Plot the histogram of the ratios in narrow bins: tmax <- max(ximsqnum/ximsqden); # largest value among the samples tmax <- ceiling(tmax); # least integer greater than tmax hist(ximsqnum/ximsqden,freq=FALSE,breaks=seq(0,tmax,by=0.5)); # plot the histogram # Compare with plot of F density on ndf/ddf degrees of freedom: t <- seq(0,tmax,by=0.1); # plot range lines(t,df(t,df1=ndf,df2=ddf)); # plot of the F density on ndf over ddf d.o.f. ### EXTRA: Plot F(ndf,ddf) and F(ddf,ndf) side-by-side dev.off(); # clear the graphics par(mfrow=c(1,2)); # (1,2) means one row with two columns tmax <- max(ximsqnum/ximsqden); hist(ximsqnum/ximsqden,freq=FALSE,breaks=0:ceiling(tmax)); # F(ndf,ddf) lines(t,df(t,df1=ndf,df2=ddf)); tmax <- max(ximsqden/ximsqnum); hist(ximsqden/ximsqnum,freq=FALSE,breaks=0:ceiling(tmax)); # F(ddf,ndf) lines(t,df(t,df1=ddf,df2=ndf));