# Beispiel 13.10 zu Eichbaendern:

n <- 200
m <- c(4,Inf)
sigma <- 0.1
alpha <- 0.05

X <- (1:n)/n
Y <- X + sigma*rnorm(n)
yy <- seq(-0.2,1.2,length.out=201)

xxl <- xxu <- matrix(rep(0,length(m)*length(yy)),ncol=length(m))
res <- lm(Y ~ X)
ah <- res$coefficients[1]
bh <- res$coefficients[2]
sh2 <- sum(res$residuals^2)/(n-2)

mX <- mean(X)
Q <- sum((X - mX)^2)
dh <- sh * sqrt(qf(1-alpha,1,n-2)) / sqrt(Q) / abs(bh)
kh1 <- (sh2/bh^2)*(1/m + 1/n)*qf(1-alpha,1,n-2)
kh2 <- (sh2/bh^2)*(1/Q)      *qf(1-alpha,1,n-2)

xxh <- (yy - ah)/bh
for (i in 1:length(m))
{
	xxl[,i] <- xxh + (kh2*(xxh-mX) - sqrt((1-kh2)*kh1[i] + kh2*(xxh-mX)^2))/(1-kh2)
	xxu[,i] <- xxh + (kh2*(xxh-mX) + sqrt((1-kh2)*kh1[i] + kh2*(xxh-mX)^2))/(1-kh2)
}

par(cex=1.2,mai=c(0.5,0.5,0.1,0.1))
plot(X,Y, xlim=c(-0.1,1.1),ylim=c(-0.2,1.2),pch=16,xlab="",ylab="")
abline(a=ah,b=bh,lwd=1,lty=2)
for (i in 1:length(m))
{
	lines(xxl[,i],yy,lty=1,lwd=2)
	lines(xxu[,i],yy,lty=1,lwd=2)
}
j7 <- sum(yy <= 0.7)
abline(h=yy[j7],lty=1)
for (i in 1:length(m))
{
	abline(v=xxl[j7,i],lty=1)
	abline(v=xxu[j7,i],lty=1)
}