library(MASS)
degreePlus <- 4 # degree of polynomials used plus one
#coefficients of cubic
# beta <- matrix(runif(degreePlus,min=-2,max=2),nrow=degreePlus)
beta <- c(1.3965849,1.9407724,-0.1215529,0.1535441)
numSamples <- 40
right <- 2
left <- -right
s <- seq(from=left,to=right,length.out=numSamples)
s <- matrix(s,nrow=numSamples)
X <- cbind(s^0,s,s^2,s^3) 
y <- X %*% beta + matrix(rnorm(numSamples),nrow=numSamples) 
M <- t(X) %*% X
Minv <- ginv(M) # MASS must be loaded
accurate.y <- X %*% beta
stdDevs <- sqrt(diag(X %*% Minv %*% t(X)))
qn <- qnorm(0.975)
barrier <- qn*stdDevs
maxy <- max(accurate.y+barrier)
miny <- min(accurate.y-barrier)
plot(s,accurate.y,type='l',ylim=c(miny-1,maxy+1))

# We now want to choose various values for $\hat\beta$ within the 95%
# confidence interval, and, for each draw the corresponding cubic in
# green.
rad <- sqrt(qchisq(0.95,df=degreePlus)) #the radius we need when
	# dealing with the standard normal distribution.

#numCurves <- 2
#betaSamples1 <- matrix(runif(numCurves*degreePlus,min=left,max=right),
#	ncol=degreePlus)
#fac <- sqrt(apply(betaSamples1,1,function(x){t(x) %*% M %*% x}))
#betaBase <- matrix(rep(beta,numCurves),ncol=degreePlus,byrow=T)
#betaSamples1 <- betaBase + betaSamples1*rad/fac
#apply(betaSamples1,1,function(u){lines(s,X %*% u,col='green')})
#lines(s,accurate.y + barrier,type='l',col='red')
#lines(s,accurate.y - barrier,type='l',col='red')
#dev.new()

numCurves <- numSamples
betaSamples <- matrix(runif(numCurves*degreePlus,min=left,max=right),
	ncol=degreePlus)
 fac <- sqrt(apply(betaSamples,1,function(x){t(x) %*% M %*% x}))
betaBase <- matrix(rep(beta,numCurves),ncol=degreePlus,byrow=T)
betaSamples <- betaBase + betaSamples*rad/fac
#pdf(file='../../WriteUp/Graphics/Chapter3/chap_3_prob_2.pdf')
postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_2.eps", onefile=FALSE, horizontal=FALSE)
plot(s,accurate.y,type='l',ylim=c(miny-1,maxy+1))
apply(betaSamples,1,function(u){lines(s,X %*% u,col='green')})
lines(s,accurate.y + barrier,type='l',col='red')
lines(s,accurate.y - barrier,type='l',col='red')
lines(s,accurate.y,type='l')
#graphics.off()
dev.off()