Hello,
Inline.
Em 06-12-2012 03:52, Rebecca escreveu:
> 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)
You cannot expect this model, a first degree polynomial in X, to fit
your data.
> 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))
Get rid of cbind(), please. It's doing nothing and the construct
data.frame(cbind(...)) is potentially harmful.
> 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)
> }
Better could be
ConfidenceSet.Pointwise2 <- function(Boot.Result,alpha){
SetBounds <- matrix(ncol=2, nrow=n)
for(j in seq_len(ncol(Boot.Result)))
SetBounds[j, ] <- quantile(Result.Sort, probs = c(alpha/2, 1 -
alpha/2))
SetBounds
}
Simpler, no?
Hope this helps,
Rui Barradas
>
> 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.
[[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.
