\begin{Schunk}
\begin{Sinput}
 mint <- -4; maxt <- 4
 x <- seq(mint,maxt,length.out=1000)
 myt <- dt(x, df=50)
 plot(x, myt, type="l", xlab="t", ylab="probability density",  ylim=c(0,0.40), main="")
 abline(h=0, lwd=0.5)
 critValH <- qt(0.025, df=50, lower.tail=F)
 statAtCritH <- dt(critVal, df=50)
 lines(c(critValH,critValH), c(0, statAtCritH), lty=4)
 xvals <- seq(critValH, maxt, length.out=50)
 polygon( x=c(xvals, xvals[50], sort(xvals,decreasing=T), critValH), 
        y=c(dt(xvals,df=50), 0, rep(0,50), 0), col=gray(.90)) 
 ##Stupidly repeat same code for lower side
 critValL <- qt(0.025, df=50, lower.tail=T)
 statAtCritL <- dt(critVal, df=50)
 lines(c(critValL,critValL), c(0, statAtCritL), lty=4)
 xvals <- seq(mint, critValL,  length.out=50)
 polygon( x=c(xvals, xvals[50], mint), 
        y=c(dt(xvals,df=50), 0, 0), col=gray(.90)) 
\end{Sinput}
\end{Schunk}
