#############################################
# Regressin example: Problem 8.1. from GCSR
#############################################

# 1. Read in the data.
#    X = matrix with columns x1,x2,x3
# 2. Compute the m.l.e. beta-hat and s2. Remember
#    p(beta|sig2,y) = N(beta-hat, sig2*V-beta)
#    p(sig2|y)      = Inv-X(n-k,s2)
# 3. Compute V-beta. See the explainations 1.-3.on the top of 
#    page 238 for an explaination. 
#    Summary: we find a matrix L, such that V-beta=LL'.
# 4. Now simulate values for sig2 ~ Inv-X(n-k,s2) and
#                            beta ~ N(beta-hat, sig2*V-beta)
#    Note: to get x~N(m,V), V=LL'
#          a) generate z~N(0,I)
#          b) x = L*z+m (compare with the "Lemmas" at the beginning 
#             of the homework 9 solutions.
# 5. Now consider various summaries and plots of the posterior
#    simulations for beta and sig2 to answer the questions.

#############################################
# 1. read in the data
#
data <- scan("radon.data",what=list(1,1,1,1))
				       # read in data from "radon.data"

y <- log(data[[4]])                       # y is the last data column
X <- cbind(data[[1]],data[[2]],data[[3]]) # design matrix
   # note: x0 = intercept (not in design matrix - Splus will
   #            add it automatically!)n
   #       x1 = indicator for Blue Earth
   #       x2 = indicator for Clay
   #       x3 = indicator for first floor
k <- ncol(X)+1    # number pars (+1 for intercept)
n_ length(y)      # number of observations

#############################################
# 2. get the posterior momemts, i.e. beta-hat, V-beta and s2
#
lf <- lsfit(X,y)  # get the m.l.e. beta-hat (remember, this is at the 
                  # same time the posterior mean for (beta|sigma,y)
print(lf)         # see what you have

bhat <- lf$coef   # beta hat

#############################################
# 3. Get V-beta
# see the comments 1.-3. on the top of page 238 for an explaination of the 
# following 2 lines.
# Summary: L is such that V-beta = LL'
R <- qr.R(lf$qr)  # QR decomp of X
L <- solve(R)     # chol decomp of V_beta

s2 <- sum(lf$resid^2)/(n-k)   # s2 in the expression for the marginal 
                              # posterior p(sigma^2|y) = Inv-X(n-k,s2)

#############################################
# 4. now generate (sig2,beta) ~ p(sig2,beta|y) 
#                          = p(sig2|y) * p(beta|sig2,y)
#

# generate sig2 ~ Inv-Chi(n-k,s2)
sig2 <- (n-k)*s2/rchisq(1000,n-k)
# and beta ~ N(bhat, sig2*LL')
# (i)  z[j] ~ N(0,sig2[j]*I), columns j=1...1000
z <- matrix(rnorm(4000,sd=sqrt(sig2)),byrow=T,ncol=1000)
# (ii) b[j] = bhat+L*z[j],    columns j=1...1000
b <- bhat + L %*% z

#############################################
# 5. various summaries and plots
#
options(digits=2)
apply(b,1,summary)

boxplot(b[1,],b[2,],b[3,],b[4,],
	names=c("CONST","BLUE E","CLAY","1ST FL"))
abline(h=0)

boxplot(b[1,]+b[2,],b[1,]+b[3,],b[1,],names=c("BLUE","CLAY","GOODH"))

cts <- cbind(b[1,]+b[2,], b[1,]+b[3,], b[1,])
sum(b[2,]>0)/1000 # = 0.56
sum(b[3,]>0)/1000 # = 0.4
sum(b[2,]-b[3,]>0)/1000 # 0.68