###################################################
### chunk number 1: preliminaries
###################################################
options(width=85, show.signif.stars = FALSE,
        lattice.theme = function() canonical.theme("pdf", color = FALSE),
        str = strOptions(strict.width = "cut"))
.f3 <- "%.3f"


###################################################
### chunk number 2: stdnorm
###################################################
set.seed(1234)                          # reproducible "random" values
X <- rnorm(100000)                      # standard normal values
zapsmall(summary(V <- X^2))             # a very skew distribution
var(V)


###################################################
### chunk number 3: stdnormdens
###################################################
library(ggplot2)
qmax <- qchisq(1-1/400, df=1)
quants <- seq(1/400, qmax, len=200)
probs <- pchisq(quants, df=1)
print(qplot(x=quants, y=quantile(V, probs, names=FALSE), geom="line",
            xlab=expression(chi[1]^2 * " Quantiles"), ylab="Sample quantiles")+
      geom_segment(aes(x=0,y=0,xend=qmax,yend=qmax, linetype=2)))


###################################################
### chunk number 4: noncentral
###################################################
V1 <- rnorm(100000, mean=2)^2
zapsmall(summary(V1))
var(V1)


###################################################
### chunk number 5: noncentral
###################################################
pp <- ppoints(200)
print(qplot(qchisq(pp,df=1,ncp=4), quantile(V1,pp,names=FALSE),
            xlab=expression(chi[1]^2 * "(4) Quantiles"), ylab="Sample quantiles")
      + geom_abline(intercept=0,slope=1,linetype=2),
      TRUE, vp=viewport(0.25,0.5,0.5,1))
print(qplot(qchisq(pp,df=1), quantile(V1,pp,names=FALSE),
            xlab=expression(chi[1]^2 * " Quantiles"), ylab="Sample quantiles")
      + geom_abline(intercept=0,slope=1,linetype=2),
      FALSE, vp=viewport(0.75,0.5,0.5,1))


###################################################
### chunk number 6: RSSsim
###################################################
lm1 <- lm(optden ~ carb, Formaldehyde)
str(Ymat <- data.matrix(unname(simulate(lm1, 10000))))
str(RSS <- deviance(fits <- lm(Ymat ~ carb, Formaldehyde)))
fits[["df.residual"]]


###################################################
### chunk number 7: sigsq
###################################################
(sigsq <- deviance(lm1)/4)
summary(RSS)


###################################################
### chunk number 8: sigsqscale
###################################################
summary(RSSsc <- RSS/sigsq)
var(RSSsc)


###################################################
### chunk number 9: RSSqq
###################################################
print(qplot(qchisq(pp,df=4), quantile(RSSsc,pp,names=FALSE),
            xlab=expression(chi[4]^2 * " Quantiles"),
            ylab="Sample quantiles"),
      TRUE, vp=viewport(0.25,0.5,0.5,1))
print(qplot(pp, pchisq(quantile(RSSsc, pp, names=FALSE), df=4),
            xlab="Cumulative probability",
            ylab="Theoretical cdf. evaluated at observed quantiles"),
      FALSE, vp=viewport(0.75,0.5,0.5,1))


###################################################
### chunk number 10: RSSdens
###################################################
xv <- seq(min(RSSsc), max(RSSsc), len=200)
print(qplot(RSSsc, geom="density",
            xlab="Scaled residual sums of squares") +
      geom_line(aes(x=xv,y=dchisq(xv,df=4),linetype=2)),
      TRUE, vp=viewport(0.5,0.5,0.8,1))