---
title: 'Semiparametric regression: local methods, comparison'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

In what follows we use the MC simulation strategy to compare two different smoothers with respect to their mean square errors as a function of the bandwith.

# Data

```{r simulated data}
n=200
#x=sort(runif(n,0,1))

#x=rbeta(n,1,6)
truef=function(x) sin(2*pi*x)
signaltonoise=2

simula=function(x=seq(0,1,length=n),signalnoise=signaltonoise){
  my=truef(x)
  y=rnorm(n,my,sd=sd(my)/signalnoise)
  return(data.frame(x=x,y=y,truef=my))
}

d=simula(x=sort(c(0,1,rbeta(n-2,3,3))))
xobs=d$x
plot(d$x,d$y)
lines(d$x,d$truef,col="darkred",lwd=2)
```



# 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
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))
}
binsmoother=function(x,xx,yy,h){
  seq(min(xx),max(xx),length=(max(xx)-min(xx))/h)
  X=model.matrix(yy~cut(xx,seq(min(xx),max(xx),length=(max(xx)-min(xx))/h),include.lowest = TRUE))
  Xnew=model.matrix(x~cut(x,seq(min(xx),max(xx),length=(max(xx)-min(xx))/h),include.lowest = TRUE))
  Xnew %*% solve(t(X)%*%X,t(X)%*%yy)
}
```


# Bias and variance via MC simulation

```{r dof computation}
smoother1=hneighmean
h.vec1=seq(20,(1+floor(n/2)),by=3)
#h.vec1=seq(0.1,0.3,length=50)

zeri=rep(0,length=length(d$y))
columnj=function(j,h) {
  zeri[j]=1
  smoother1(d$x,d$x,zeri,h=h)
}
computeDOF=function(h){
  L=sapply(1:n,columnj,h=h)
  sum(diag(L))
}
dof.vec1=sapply(h.vec1,computeDOF)


smoother2=kernelreg
h.vec2=seq(0.15,4,length=51)[-1]
columnj=function(j,h) {
  zeri[j]=1
  smoother2(d$x,d$x,zeri,h=h)
}
computeDOF=function(h){
  L=sapply(1:n,columnj,h=h)
  sum(diag(L))
}
dof.vec2=sapply(h.vec2,computeDOF)

#smoother=loessreg
#h.vec=seq(0.05,1,length=5)

dofforcomparison=3
h1=h.vec1[which.min(abs(dof.vec1-dofforcomparison))]
h2=h.vec2[which.min(abs(dof.vec2-dofforcomparison))]
```


```{r smoothers plot}
plot(d$x,d$y,pch=20,xlab="x",ylab="y")
curve(smoother1(x,d$x,d$y,h=h1),add=TRUE,lwd=2,col="red")
curve(smoother2(x,d$x,d$y,h=h2),add=TRUE,lwd=2,col="darkgreen")

```

Comparison of bias as a function of $x$

```{r}
M=1000
#ris=matrix...
#for (i in 1:M) {
#  d.i=simula()
#  rius[] = smoother(d.i$x,d.i$x,d.i$y,h=hfix)
#}
n=200 #sample size
xseq=xobs # values of x in which the bias and variance are calculated
truefseq=truef(xseq) # true values of f
ff.hat=sapply(rep(1,M),FUN=function(x) {d.i=simula(x=xobs)
  smoother1(xseq,d.i$x,d.i$y,h=h1)})
biasf.hat.1=(apply(ff.hat,1,mean)-truefseq)
varf.hat.1=apply(ff.hat,1,var)
ff.hat=sapply(rep(1,M),FUN=function(x) {d.i=simula(x=xobs)
  smoother2(xseq,d.i$x,d.i$y,h=h2)})
biasf.hat.2=(apply(ff.hat,1,mean)-truefseq)
varf.hat.2=apply(ff.hat,1,var)
#matplot(x,cbind(biasf.hat^2,varf.hat),col=c("blue","red"),pch=20)
matplot(xseq,cbind(biasf.hat.1,biasf.hat.2),col=c("darkred","darkgreen"),pch=20)
matplot(xseq,cbind(varf.hat.1,varf.hat.2),col=c("darkred","darkgreen"),pch=20)
matplot(xseq,cbind(biasf.hat.1^2+varf.hat.1,biasf.hat.2^2+varf.hat.1),col=c("darkred","darkgreen"),pch=20)
```


Bias and variance as a whole as a function of bandwidth

```{r}
M=1000

onesim=function(z){
  d.i=simula()
  sapply(h.vec1,FUN=function(x) smoother1(d.i$x,d.i$x,d.i$y,x))
}
#onesim(0.1)
allsim=sapply(1:M,onesim,simplify="array")

var2=function(x) mean(x^2)-mean(x)^2

my=d$truef

f.mean=apply(allsim,1:2,mean)
f.bias=apply(f.mean,2,FUN=function(x) x-my)
f.variance=apply(allsim,1:2,var2)

bias.overall=apply(abs(f.bias)^2,2,mean)
variance.overall=apply(f.variance,2,mean)
mse.overall=bias.overall+variance.overall
# # alternative calculation
# all.error=apply(allsim,2:3,FUN=function(x) x-my)
# f.mse.all=apply(all.error,1:2,FUN=function(x) mean(x^2))
# f.bias2=apply(all.error,1:2,mean)
# f.var2=apply(all.error,1:2,FUN=var2)

matplot(dof.vec1,cbind(bias.overall,variance.overall,mse.overall),col=c("blue","red","black"),lwd=c(1,1,2),type="l",lty=c(1,1,2))
abline(v=dof.vec1[which.min(mse.overall)])

```
```{r}
M=1000

onesim=function(z){
  d.i=simula()
  sapply(h.vec2,FUN=function(x) smoother2(d.i$x,d.i$x,d.i$y,x))
}
#onesim(0.1)
allsim=sapply(1:M,onesim,simplify="array")

var2=function(x) mean(x^2)-mean(x)^2

f.mean=apply(allsim,1:2,mean)
f.bias=apply(f.mean,2,FUN=function(x) x-my)
f.variance=apply(allsim,1:2,var2)

bias.overall2=apply(abs(f.bias)^2,2,mean)
variance.overall2=apply(f.variance,2,mean)
mse.overall2=bias.overall2+variance.overall2
# # alternative calculation
# all.error=apply(allsim,2:3,FUN=function(x) x-my)
# f.mse.all=apply(all.error,1:2,FUN=function(x) mean(x^2))
# f.bias2=apply(all.error,1:2,mean)
# f.var2=apply(all.error,1:2,FUN=var2)

matplot(dof.vec2,cbind(bias.overall2,variance.overall2,mse.overall2),col=c("blue","red","black"),lwd=c(1,1,2),type="l",lty=c(1,1,2))
abline(v=dof.vec2[which.min(mse.overall2)])
```


```{r}
matplot(cbind(dof.vec1,dof.vec2),cbind(mse.overall,mse.overall2),col=c("blue","red"),lwd=2,type="l",lty=1)

```