# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Exemplo:
#     Intervalo de predicao no modelo classico de regressao linear.
# Modelo: 
#     y = beta0 + beta1 x + epsilon
#==============================================================================
set.seed(12345)
n          <-  10      # Tamanho de cada amostra
alpha      <-  0.5     # Intercepto
beta       <-  0.5     # Coeficiente angular
sigma      <-  0.1     # Desvio padrao do erro
epsilon    <-  rnorm(n,0,sigma)
x          <-  rnorm(n,0,1)
y          <- alpha + beta*x + epsilon
xseq       <- seq(-6,6,length=500)
yseq       <- alpha + beta*xseq
yfit       <- lm(y ~ x)
x0         <- cbind(rep(1,500),xseq)
ypred      <- c(yfit$coef%*%t(x0))
s2pred     <- t(yfit$res)%*%yfit$res/yfit$df.residual
var.ypred  <- rep(NA,500)
for(i in 1:500){
    var.ypred[i] <- s2pred + t(x0[i,])%*%vcov(yfit)%*%x0[i,]
}
ep.ypred   <- sqrt( var.ypred )
talpha     <- qnorm(0.975)
liminf     <- ypred-talpha*ep.ypred
limsup     <- ypred+talpha*ep.ypred

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(x,y,pch=15,xlab=expression(x),ylab=expression(y),
    xlim=c(-7,7),ylim=c(-3.5,3.5))
lines(xseq,yseq,col="red",lwd=2)
lines(xseq,ypred,col="blue",lwd=2)
lines(xseq,liminf,col="black",lwd=2,lty=2)
lines(xseq,limsup,col="black",lwd=2,lty=2)
text(-5,3,"n = 10")
legend(-7,3,legend=c("Verdadeira","Ajustada","Int. Pred."),
    col=c("red","blue","black"),lty=c(1,1,2),bty="n")