function (k=3,n=5,Nsim=10000,a.recip=3,Fmin.observed=NULL) 
{
# k is the number of sample variances to be compared
# n is either the common sample size behind each variance, 
# n-1 degrees of freedom each, or it is a vector of length k
# giving the sample sizes behind each sample variance. 
# Nsim is the number of simulations to run (assuming normality 
# and equal variances in the k samples of size n).
# a=1/a.recip is the suggested cut-off for any given rule of thumb.
# When Fmin.observed is provided, the plotted histogram will indicate
# its position by an orange line and give the proportion of simulated
# values <= Fmin.observed. Otherwise only the location of a is indicated 
# by a blue line with the proportion of simulated Fmin values <= a.
#-------------------------------------------------------------------------- 
kn=length(n)
if(kn>1){if(kn!=k)return(paste("the n-vector (",paste(n,collapse=", "),
 ") does not have length k =",k))}
ax=1/a.recip
s2mat=NULL
for(i in 1:k){
if(length(n)==1){
s2=rchisq(Nsim,n-1)/(n-1)}else{
s2=rchisq(Nsim,n[i]-1)/(n[i]-1)
}
s2mat=cbind(s2mat,s2)
}
if(kn>1){nlab=paste("n = (",paste(n,collapse=", "), ")")}else{
nlab=paste("n = ",n)
}
if(length(Fmin.observed)>0){Flab=paste(", Fmin.observed =",format(signif(Fmin.observed,5)))}else{
Flab=""}
maxs2=apply(s2mat,1,max)
mins2=apply(s2mat,1,min)
Frat=mins2/maxs2
out=hist(Frat,
xlab=expression("min"~group("(",list(s[1]^2, ... ,s[k]^2),")")/
"max"~group("(",list(s[1]^2, ... ,s[k]^2),")")),
main=paste("k =",k,", ", nlab,",   Nsim = ",
  Nsim,",   a = 1/",a.recip,Flab),
nclass=100,cex.lab=1.3,xlim=c(0,1))
abline(v=ax,col="blue",lwd=2)
mh=max(out$counts)
text(ax-.01,.7*mh,format(signif(mean(Frat<=ax),3)),
adj=1,cex=1.3,col="blue",srt=90)
if(length(Fmin.observed)>0) {
text(Fmin.observed-.01,mh,
format(signif(mean(Frat<=Fmin.observed),3)),
adj=1,cex=1.3,col="orange",srt=90)
}
abline(v=Fmin.observed,lwd=2,col="orange")
}