---
title: 'Semiparametric regression: local methods, bootstrap'
output:
  html_document: default
  pdf_document: default
date: "2023-10-28"
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Data


```{r dati}
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))
```

```{r simulated data}
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


```{r smoothers definitions}
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)
## 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

```{r}
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}} \]

```{r}
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}} \]


```{r}
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"))
```
```{r}
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


```{r}
bootsmoother=function(xx,yy,...){
  d=data.frame(xx=xx,yy=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)
```

```{r}
B=100
bts1.sample=sapply(1:B,FUN=function(x) bootsmoother(xx=d$x,yy=d$y,h=0.1))
```

```{r}
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


```{r}
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)
```

```{r}
B=100
bts2.sample=sapply(1:B,FUN=function(x) bootsmoother.sp(xx=d$x,yy=d$y,h=0.1))
```


```{r}
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

```{r}
bts.sample=bts1.sample
conflev=0.95
```


Percentile ci

```{r}
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)
```

```{r}
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


```{r smootherforboot}
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])
```

```{r}
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)

```



```{r smootherforboot2}
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])
```

```{r}
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)

```
```{r message=FALSE, warning=FALSE}
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)

```