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));