Data

d=read.table("lidar.dat",header=TRUE)
names(d)=c("x","y")
d=d[sort.list(d$x),]
d$x=(d$x-min(d$x))/(max(d$x)-min(d$x))
n=200
x=sort(runif(n,0,1))
x=seq(0,1,length=n)
truef=function(x) sin(2*pi*x)
my=truef(x)
signaltonoise=2

simula=function(signalnoise=signaltonoise){
  y=rnorm(n,my,sd=sd(my)/signalnoise)
  return(data.frame(x=x,y=y,truef=my))
}
d=simula()
plot(d$x,d$y)

The smoothers

runmean0=function(x,xx,yy,h)
  mean(yy[abs(xx-x)<h])
runmean=function(x,xx,yy,h)  
  mapply(runmean0,x,MoreArgs=list(xx=xx,yy=yy,h=h))
runmean(c(0.5,0.6),d$x,d$y,h=0.1)
## [1]  0.04881761 -0.49564193
## Mean of nearest neighbours
hneighmean0=function(x,xx,yy,h){
  mean(yy[abs(xx-x)<sort(abs(xx-x))[h]])
}
hneighmean=function(x,xx,yy,h)  
  mapply(hneighmean0,x,MoreArgs=list(xx=xx,yy=yy,h=h))
## Kernel regression
kernelreg=function(x,xx,yy,h,kernel="box") ksmooth(xx,yy,kernel=kernel,bandwidth=h,x.points=x)$y
## Loess
loessreg = function(x,xx,yy,h,degree=2,family="gaussian") {
  fit=loess(yy~xx,span=h,surface="direct")
  predict(fit,newdata=data.frame(xx=x))
}

smoother=kernelreg

Obtain the smoothing matrix

zeri=rep(0,length=length(d$y))
columnj=function(j,h,smoother="smoother") {
  zeri[j]=1
  smoother=get(smoother)
  smoother(x,d$x,zeri,h=h)
}
L=sapply(1:n,columnj,h=0.2,smoother="smoother")

Confidence intervals and bands

We obtain pointwise confidence intervals for \(f(x_i)\) based on \[ \hat{f}-f \thicksim_{as} \mathcal N(0,\sigma^2LL^T) \] as \[ \hat{f}(x_i) \pm z_{1-\alpha/2}\sigma \sqrt{[LL^T]_{ii}} \]

conflev=0.95
f.hat=L%*%d$y
sigma2.hat=sum((d$y-f.hat)^2)/(length(d$y)-sum(diag(L)))
f.low=f.hat+qnorm((1-conflev)/2)*sqrt(sigma2.hat*diag(L%*%t(L)))
f.up=f.hat-qnorm((1-conflev)/2)*sqrt(sigma2.hat*diag(L%*%t(L)))
plot(d$x,d$y)
matlines(d$x,cbind(f.hat,f.low,f.up),lty=c(1,2,2),lwd=2,col="red")

A confidence band, that is \(f_{low}\), \(f_{up}\) such that \[ P(f_{low}(x)\leq f(x) \leq f_{up}(x)\mbox{ for all }x\in\{x_1,\ldots,x_r\}) \geq 1-\alpha \] the Bonferroni correction as \[ \hat{f}(x_i) \pm z_{1-\alpha/(2n)}\sigma \sqrt{[LL^T]_{ii}} \]

conflev1=1-(1-conflev)/n
f.hat=L%*%d$y
sigma2.hat=sum((d$y-f.hat)^2)/(length(d$y)-sum(diag(L)))
f.low.sim=f.hat+qnorm((1-conflev1)/2)*sqrt(sigma2.hat*diag(L%*%t(L)))
f.up.sim=f.hat-qnorm((1-conflev1)/2)*sqrt(sigma2.hat*diag(L%*%t(L)))
plot(d$x,d$y)
matlines(d$x,cbind(f.hat,f.low,f.up,f.low.sim,f.up.sim),lty=c(1,2,2,3,3),lwd=2,col=c("black","red","red","blue","blue"))

sim.dist=replicate(1000,max(abs(rnorm(n,0,sqrt(sigma2.hat*diag(L%*%t(L))))/sqrt(sigma2.hat*diag(L%*%t(L))))))
mquant=quantile(sim.dist,conflev)
f.low.sim2=f.hat-mquant*sqrt(sigma2.hat*diag(L%*%t(L)))
f.up.sim2=f.hat+mquant*sqrt(sigma2.hat*diag(L%*%t(L)))
plot(d$x,d$y)
matlines(d$x,cbind(f.hat,f.low,f.up,f.low.sim,f.up.sim,f.low.sim2,f.up.sim2),lty=c(1,2,2,3,3,2,2),lwd=2,col=c("black","red","red","blue","blue","green","green"))

Bootstrap strategy one: non parametric

bootsmoother=function(xx,yy,...){
  d.i=d[sample(1:nrow(d),replace=TRUE),]
  smoother(d$x,d.i$x,d.i$y,...)
}
res=bootsmoother(d$x,d$y,h=0.1)
B=100
bts1.sample=sapply(1:B,FUN=function(x) bootsmoother(xx=d$x,yy=d$y,h=0.1))
plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
matlines(d$x,t(apply(bts1.sample,1,quantile,prob=c(0.025,0.975))),col="red",lty=1,lwd=2)

Bootstrap strategy two: semi parametric

y.hat=smoother(d$x,d$x,d$y,h=0.1)
resid=d$y-y.hat
bootsmoother.sp=function(xx,yy,...){
  d$ynew=y.hat+sample(resid,replace=FALSE)
  smoother(d$x,d$x,d$ynew,...)
}
res=bootsmoother.sp(d$x,d$y,h=0.1)
B=100
bts2.sample=sapply(1:B,FUN=function(x) bootsmoother.sp(xx=d$x,yy=d$y,h=0.1))
plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
matlines(d$x,t(apply(bts2.sample,1,quantile,prob=c(0.025,0.975))),col="red",lty=1,lwd=2)

Bootstrap ci

bts.sample=bts1.sample
conflev=0.95

Percentile ci

plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
matlines(d$x,t(apply(bts.sample,1,quantile,prob=c((1-conflev)/2,1-(1-conflev)/2))),col="red",lty=1,lwd=2)

plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
boot.summary=apply(bts.sample,1,FUN=function(x) c(mean(x),sd(x)))
matlines(d$x,cbind(boot.summary[1,]+qnorm((1-conflev)/2)*boot.summary[2,],
                   boot.summary[1,]-qnorm((1-conflev)/2)*boot.summary[2,]),col="red",lty=1,lwd=2)

Using boot

library(boot)
smootherforboot=function(d,inds,...){
  d.b=d[inds,]
  smoother(d$x,d.b$x,d.b$y,...)
}
d.boot=boot(d,smootherforboot,R=1000,h=0.1)
out.ci=sapply(1:nrow(d),FUN=function(x) boot.ci(d.boot,index=x,type="norm")$normal[-1])
plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
matlines(d$x,t(out.ci),col="red",lty=1,lwd=2)

smootherforboot2=function(d,inds,...){
  d$ynew=y.hat+resid[inds]
  smoother(d$x,d$x,d$ynew,...)
}
d.boot=boot(d,smootherforboot2,R=1000,h=0.1)
out.ci=sapply(1:nrow(d),FUN=function(x) boot.ci(d.boot,index=x,type="norm")$normal[-1])
plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
matlines(d$x,t(out.ci),col="red",lty=1,lwd=2)

out.ci=sapply(1:nrow(d),FUN=function(x) boot.ci(d.boot,index=x,type="bca")$bca[-(1:3)])
plot(d$x,d$y)
lines(d$x,smoother(d$x,d$x,d$y,h=0.1))
boot.summary=apply(bts.sample,1,FUN=function(x) c(mean(x),sd(x)))
matlines(d$x,t(out.ci),col="red",lty=1,lwd=2)