## Source: Chapter11.r (with modifications) in ## http://www.stat.colostate.edu/computationalstatistics/CompStatR.zip ############################################################################ ### EXAMPLE 11.9 - 11.11 CONFIDENCE REGIONS ############################################################################ # smoothdat.x, smoothdat.y = data for smoothing # S = smoothing matrix # yhat = original smooth estimate # boot.yhat = bootstrap smooths ############################################################################ ##INITIAL VALUES cbsmooth=read.table("confidenceband.dat",header=T) smoothdat.y = cbsmooth$y smoothdat.x = cbsmooth$x plot(smoothdat.x,smoothdat.y,xlab="Predictor",ylab="Response") s = function(x){x^2;} plot(s, min(smoothdat.x), max(smoothdat.x), add=T) ##CONSTRUCT RUNNING-LINE SMOOTHING MATRIX ##Span of 15 chosen by cross-validation of scatter.smooth() s = 7 #span = 2*s+1 S = matrix(0,50,50) for (i in 1:50) { lo=max(1,i-s) ; hi=min(50,i+s) design.mat=cbind(1,smoothdat.x[lo:hi]) h=design.mat%*%solve(t(design.mat)%*%design.mat)%*%t(design.mat) where=ifelse(i<26-s-1,dim(h)[2]-s,s+1) S[i,lo:hi]=h[where,] } yhat = S%*%cbind(smoothdat.y) resids = smoothdat.y-yhat #Note: qr(S%*%t(S))$rank is 47 so watch out for inverting ##BOOTSTRAP set.seed(323) boot.yhat=matrix(0,nrow=50,ncol=1000) for (i in 1:1000) { newfit=yhat+sample(resids,50,replace=T) boot.yhat[,i]=S%*%cbind(newfit) } pointwise.lower=apply(boot.yhat,1,quantile,.025) pointwise.upper=apply(boot.yhat,1,quantile,.975) ## plot pointwise confidence band plot(smoothdat.x,smoothdat.y) lines(smoothdat.x,yhat,lwd=2) polygon(c(smoothdat.x,rev(smoothdat.x)),c(pointwise.lower,rev(pointwise.upper)), density=-1,col=grey(.7),border=NA) lines(c(min(smoothdat.x),max(smoothdat.x)),-c(.00,.00),lty=2) points(smoothdat.x,smoothdat.y,pch=16) ##INFLATION METHOD (Example 11.10) #Used bisection to search for correct omega. Result below. omega=1.66 Lhat=yhat-pointwise.lower Uhat=pointwise.upper-yhat adj.pointwise.lower=yhat-omega*Lhat adj.pointwise.upper=yhat+omega*Uhat is.completely.within=function(yhat,lo,up) { all(yhat>=lo & yhat<=up) } number.within=apply(boot.yhat,2,is.completely.within, lo=adj.pointwise.lower,up=adj.pointwise.upper) sum(number.within)/1000 # plot graph plot(smoothdat.x,smoothdat.y) lines(smoothdat.x,yhat,lwd=2) polygon(c(smoothdat.x,rev(smoothdat.x)),c(pointwise.lower,rev(pointwise.upper)), density=-1,col=grey(.7),border=NA) lines(c(min(smoothdat.x),max(smoothdat.x)),-c(.00,.00),lty=2) points(smoothdat.x,smoothdat.y,pch=16) lines(smoothdat.x,adj.pointwise.upper,lty=2) lines(smoothdat.x,adj.pointwise.lower,lty=2) ##SUPERIMPOSITION METHOD (Example 11.11) ##Using same span as above s.hat.star=matrix(0,nrow=50,ncol=1000) v.star=rep(0,1000) eigen.decomp=eigen(S%*%t(S)) SST.inv=eigen.decomp$vectors%*%diag(1/eigen.decomp$values)%*% t(eigen.decomp$vectors) denom=50-2*sum(diag(S))+sum(diag(S%*%t(S))) set.seed(382) for (i in 1:1000) { ystar=yhat+sample(resids,50,replace=T) s.hat.star[,i]=S%*%cbind(ystar) sigma2.hat=sum((s.hat.star[,i]-ystar)^2)/denom v.star[i]=t(cbind(s.hat.star[,i]-ystar))%*%SST.inv%*% rbind(s.hat.star[,i]-ystar)/sigma2.hat } sorting=(1:1000)[order(v.star)] s.hat.star=s.hat.star[,sorting] v.star=sort(v.star) ###OUTPUT GRAPHS plot(smoothdat.x,smoothdat.y, xlab="X", ylab="Y", main="Superimposition") lines(smoothdat.x,yhat,lwd=2) polygon(c(smoothdat.x,rev(smoothdat.x)),c(pointwise.lower,rev(pointwise.upper)), density=-1,col=grey(.7),border=NA) lines(c(min(smoothdat.x),max(smoothdat.x)),-c(.00,.00),lty=2, lwd=2, col="red") points(smoothdat.x,smoothdat.y,pch=16) for (i in seq(25,975,by=50)) { lines(smoothdat.x,s.hat.star[,i]) } # a simplified graph plot(smoothdat.x,smoothdat.y, xlab="X", ylab="Y", main="Superimposition") lines(c(min(smoothdat.x),max(smoothdat.x)),-c(.00,.00),lty=2, lwd=2, col="red") for (i in seq(25,975,by=50)) { lines(smoothdat.x,s.hat.star[,i]) }