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


###################################################
### chunk number 2: lm12
###################################################
lm1 <- lm(optden ~ 1 + carb, Formaldehyde)
set.seed(1234)                      # allow for reproducible "random" numbers
badDat <- within(data.frame(x1=1:20, x2=rnorm(20,mean=6,sd=0.2),
                            x4=rexp(20,rate=0.02),
                            y=runif(20,min=18,max=24)),
                 x3 <- x1 + 2*x2)   # create linear combination
lm2 <- lm(y ~ x1 + x2 + x3 + x4, badDat)
lm3 <- lm(count ~ spray, InsectSprays)
lmlst <- list(lm1=lm1, lm2=lm2, lm3=lm3)
mmlst <- lapply(lmlst, model.matrix)


###################################################
### chunk number 3: ranks
###################################################
sapply(lmlst, function(fm) c(rank=fm[["rank"]], p=length(coef(fm))))


###################################################
### chunk number 4: rdiag
###################################################
lapply(lmlst, function(fm) diag(fm[["qr"]][["qr"]]))
lapply(mmlst, function(mm) svd(mm, nu=0, nv=0)[["d"]])
sapply(lmlst, kappa, exact=TRUE)


###################################################
### chunk number 5: Moore-Penrose
###################################################
X <- mmlst[[3]]
lm3qr <- lmlst[[3]]$qr
Q1 <- qr.Q(lm3qr)
R <- qr.R(lm3qr)
Xpinv <- backsolve(R, t(Q1))
zapsmall(Xpinv %*% X)
all.equal(X %*% Xpinv %*% X, X, check.attr=FALSE)
all.equal(Xpinv %*% X %*% Xpinv, Xpinv, check.attr=FALSE)


###################################################
### chunk number 6: rank-deficient-Moore-Penrose
###################################################
X <- mmlst[[2]]
lm2qr <- lmlst[[2]]$qr
SVD <- svd(X)
(rr <- lm2qr$rank)                      # rank
(rrind <- seq_len(rr))                  # safer than 1:rr
(dpinv <- c(1/SVD$d[rrind], rep(0, ncol(X) - rr))) 
str(Xpinv1 <- SVD$v %*% (dpinv * t(SVD$u)))
zapsmall(Xpinv1 %*% X)
all.equal(X %*% Xpinv1 %*% X, X, check.attr=FALSE)
all.equal(Xpinv1 %*% X %*% Xpinv1, Xpinv1, check.attr=FALSE)

## An alternative construction is to reduce the SVD components to the first 4 columns
str(SVDred <- list(d=SVD$d[rrind], u=SVD$u[,rrind], v=SVD$v[,rrind]))
str(Xpinv2 <- with(SVDred, v %*% (1/d * t(u))))
all.equal(Xpinv2, Xpinv1)

## Finally, we can use a similar construction on the QR decomposition
## taking into account the rearrangement of the columns of X
Xpiv <- X[, lm2qr$pivot]
str(Xpinv3 <- rbind(backsolve(qr.R(lm2qr)[rrind, rrind], t(qr.Q(lm2qr)[, rrind])), 0))
all.equal(Xpiv %*% Xpinv3 %*% Xpiv, Xpiv, check.attr=FALSE)
all.equal(Xpinv3 %*% Xpiv %*% Xpinv3, Xpinv3, check.attr=FALSE)


###################################################
### chunk number 7: simulate
###################################################
str(Ymat <- data.matrix(unname(simulate(lm1, 10000))))
fits <- lm(Ymat ~ carb, Formaldehyde)
str(coefs <- coef(fits))


###################################################
### chunk number 8: coefs
###################################################
str(coefs <- data.frame(t(coef(fits)), check.names=FALSE))


###################################################
### chunk number 9: truecoef
###################################################
printCoefmat(coef(summary(lm1)))


###################################################
### chunk number 10: coefMeans
###################################################
sapply(coefs, mean)


###################################################
### chunk number 11: coefsds
###################################################
sapply(coefs, sd)


###################################################
### chunk number 12: coefCorr
###################################################
summary(lm1, corr=TRUE)$correlation
cor(coefs)


###################################################
### chunk number 13: coefVcov
###################################################
vcov(lm1)
var(coefs)


###################################################
### chunk number 14: coefDens
###################################################
print(qplot(`(Intercept)`, data=coefs, geom="density"),
      TRUE, vp=viewport(0.25,0.5,0.5,1))
print(qplot(carb, data=coefs, geom="density"),
      FALSE, vp=viewport(0.75,0.5,0.5,1))


###################################################
### chunk number 15: coefQQ
###################################################
print(qplot(sample=quantile(`(Intercept)`,ppoints(200)), data=coefs),
      TRUE, vp=viewport(0.25,0.5,0.5,1))
print(qplot(sample=quantile(carb, ppoints(200)), data=coefs),
      FALSE, vp=viewport(0.75,0.5,0.5,1))


###################################################
### chunk number 16: coef2ddens
###################################################
print(ggplot(coefs, aes(x=carb, y=`(Intercept)`)) + stat_binhex() + stat_density2d())