# 11-19-2016 figs = function(pdf=F,w=6,h=4) { # N(0,1) if(pdf) { pdf("fig1.pdf",F,width=w,height=h) } x = seq(-4,4,.01); density = dnorm(x) plot(x,density,type="l",lwd=3,cex.lab=1.5,main="N(0,1)", xlab=bquote(italic(x))); abline(h=0,lty=2) if(pdf) { dev.off() } # T_r if(pdf) { pdf("fig2.pdf",F,width=w,height=h) } x = seq(-4,4,.01); density = dnorm(x) plot(x,density,type="l",lwd=2,cex.lab=1.5, main=bquote(italic(T)[italic(r)]), col=1, xlab=bquote(italic(x))); abline(h=0,lty=2) df=c(1,2,3,5); for( k in 1:4) { lines(x,dt(x,df[k]),lwd=2,col=k) } legend(2.5,.35,rev(c("1","2","3","5","inf")),lwd=2,col=5:1) if(pdf) { dev.off() } # chi-squared if(pdf) { pdf("fig3.pdf",F,width=w,height=h) } x = seq(0,15,.02); density = dchisq(x,2) plot(x,density,type="l",lwd=2,cex.lab=1.5, main=bquote(italic(chi)^2), col=1, xlab=bquote(italic(x))); abline(h=0,lty=2) for( k in 1:4) { lines(x,dchisq(x,2+2*k),lwd=2,col=k+1) } legend( 6.5 , 0.47,c("2","4","6","8","10"),lwd=2,col=1:5) if(pdf) { dev.off() } # F if(pdf) { pdf("fig4.pdf",F,width=w,height=h) } x = seq(0,3,.01); density = df(x,50,50) plot(x,density,type="l",lwd=2,cex.lab=1.5, main=bquote(italic(F)[italic(rs)]), col=1, xlab=bquote(italic(x))); abline(h=0,v=1,lty=2) for( k in 1:4) { df=50-10*k; lines(x,df(x,df,df),lwd=2,col=k+1) } legend( 2 , 1.3,rev(c("10,10","20,20","30,30","40,40","50,50")), lwd=2,col=1:5) if(pdf) { dev.off() } # F 20,20 simulation n=21 if(pdf) { pdf("fig5.pdf",F,width=w,height=h) } x = seq(0,3,.01); density = df(x,20,20) plot(x,density,type="l",lwd=2,cex.lab=1.5, main=bquote(italic(F)["20,20"]), ylim=c(0,1), col=1, xlab=bquote(italic(x))); abline(h=0,v=1,lty=3) set.seed(129) nsim=1e5; X = scale( matrix( rnorm(nsim*21), 21, nsim), T, F ) Y = scale( matrix( rnorm(nsim*21), 21, nsim), T, F ) # xs2 = apply(X,2,"var"); ys2 = apply(Y,2,"var") # fs = xs2/ys2; hist( fs, seq(0,10,.05), prob=T,add=T,col=rgb(.8,0,0,.5) ) aa = c(rep(1,21) %*% X^2)/20; bb = c(rep(1,21) %*% Y^2)/20; fs = aa/bb; print(range(fs)) hist( fs, seq(0,20,.05), prob=T,add=T,col=rgb(.8,0,0,.5) ) text( 2.5, 0.90, bquote(paste(10^5," repetitions")), cex=1.5) text( 2.5, 0.70, bquote(n[y]==21), cex=1.5 ) text( 2.5, 0.50, bquote(n[y]==21), cex=1.5 ) if(pdf) { dev.off() } } figs()