I'm trying to find a bootstrap based confidence band for a linear model.
I have created a data set with X and Y
X=runif(n,-1.25,1.25)
e=rnorm(n,0,1)
Y=exp(3*X)+5*sin((30*X)/(2*pi))+2*e
fit=lm(Y~X)
summary(fit)
I define a bootstrap function named PairedBootstrap which is not listed here.
Than I try many ways to find the confidence band. One way is to predict Y using
the model I get above for 1000 times.
Data.Orig=data.frame(cbind(X,Y))
B=1000
Boot.Result=matrix(nrow=B,ncol=n)
for ( b in 1:B){
Data.Orig.Boot=PairedBootstrap(Data.Orig)
fit.Boot=predict(fit,newdata=Data.Orig.Boot,type="response")
Boot.Result[b,]=fit.Boot
}
And try to find 95% confidence interval for 1000 copies of y corresponding to
each x.
ConfidenceSet.Pointwise=function(Boot.Result,alpha){
n=ncol(Boot.Result)
B=nrow(Boot.Result)
SetBounds=matrix(ncol=2,nrow=n)
for(j in 1:n){
Result.Sort=sort(Boot.Result[,j])
SetBounds[j,1]=Result.Sort[floor(0.5*B*alpha)]
SetBounds[j,2]=Result.Sort[ceiling(B*(1-0.5*alpha))]
}
return(SetBounds)
}
And then try to line them up. But the result is not what I want.
alpha=0.05
Boot.Pointwise=ConfidenceSet.Pointwise(Boot.Result,alpha)
lines(Data.Orig[,1],Boot.Pointwise[,1], col=2, lwd=3)
Could someone tell where I'm wrong, or is there another better way to do it?
Thank you so much!
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
