# Punktweise und simultane Konfidenzintervalle fuer
# polynomiale Regressionsfunktionen:
# Unterstellt wird, dass die Beobachtungen (x[i],y[i])
# dem Modell
#    y[i] = f_*(x[i]) + epsilon[i]
# mit einer ganzrationalen Funktion f_* gegebener Ordnung
# (degree) und unabhaengigen, nach N(0,sigma^2) verteilten
# Fehlern epsilon[i], 1 <= i <= n.
# Ausgegeben/gezeichnet werden punktweise und simultane
# (1-alpha)-Konfidenzintervalle fuer f_*(xx[i]) fuer alle
# Komponenten xx[i] eines Vektors xx.

PolRegrCB <- function(y,x,degree=1,
		xx=NULL,alpha=0.05,truefxx=NULL,trueth=NULL){
	if (is.null(xx)){
		xx <- sort(x)
	}
	n <- length(x)
	p <- degree+1
	df <- n - p
	Dx <- matrix(1,nrow=n,ncol=p)
	for (i in 2:p){
		Dx[,i] <- Dx[,(i-1)]*x
	}
	thhat <- qr.solve(Dx,y)
	yhat <- Dx %*% thhat
	shat <- sqrt(sum((y - yhat)^2)/df)
	Dxx <- matrix(rep(1,length(xx)*p),ncol=p)
	for (i in 2:p){
		Dxx[,i] <- Dxx[,(i-1)]*xx
	}
	fhatxx <- Dxx %*% thhat
	shatxx <- shat * sqrt( diag(Dxx %*% qr.solve(t(Dx)%*%Dx,t(Dxx))) )
	cv.pw  <- qt(1 - alpha/2, df=df)
	cv.sim <- sqrt(p * qf(1 - alpha, df1 = p, df2 = df))
	cb.low.pw  <- fhatxx - cv.pw*shatxx
	cb.upp.pw  <- fhatxx + cv.pw*shatxx
	cb.low.sim <- fhatxx - cv.sim*shatxx
	cb.upp.sim <- fhatxx + cv.sim*shatxx
	
	par(mai=c(0.4,0.4,0.05,0.05))
	plot(x,y,pch="o",xlab="",ylab="",xlim=c(min(xx),max(xx)))

	if (!is.null(truefxx)){
		points(xx,truefxx,type="l",lwd=2,lty=2)
	}
	if (!is.null(trueth)){
		truefxx <- Dxx %*% trueth
		points(xx,truefxx,type="l",lwd=2,lty=2)
	}

	points(xx,fhatxx,type="l",lwd=2,col="blue")
	points(xx,cb.low.pw,type="l",lwd=1,col="blue")
	points(xx,cb.upp.pw,type="l",lwd=1,col="blue")
	points(xx,cb.low.sim,type="l",lwd=2)
	points(xx,cb.upp.sim,type="l",lwd=2)
	return(list(fhat=fhatxx,
		lower.bounds.pw=cb.low.pw,upper.bounds.pw=cb.upp.pw,
		lower.bounds.sim=cb.low.sim,upper.bounds.sim=cb.upp.sim,
		sigmahat=shat))
}