# R commands to accompany the binned and kernel regression
# section

wmap = read.table("E:/Vrontos/mathimata/A-mathimata-2017-2018/MSc-Statistics-Non-papametric-Statistics[EARINO]/Non-Parameteric-Regression/DataAndRcommands/wmap.dat",header=T)
hvec = c(15, 50, 75)

#----------------------------------------------------------
# A function to find the binned estimate

binned = function(x,y,h,xlab="x",ylab="y",estonly=FALSE,
                  linecol=2)
{
  bin = pmax(ceiling((x-min(x))/h),1)
  fhat = tapply(y,bin,mean)
  xlo = (0:(length(fhat)-1))*h+min(x)
  xhi = (1:length(fhat))*h+min(x)
  if(!estonly)
  {
    plot(x,y,xlab=xlab,ylab=ylab,pch=".",cex=3,
         cex.axis=1.3,cex.lab=1.3)
  }
  lines(rbind(xlo,xhi),rbind(fhat,fhat),lwd=3,col=linecol)
}
# Plots
binned(wmap$ell[1:700],wmap$Cl[1:700],15,"l","Cl")
binned(wmap$ell[1:700],wmap$Cl[1:700],50,"l","Cl")
binned(wmap$ell[1:700],wmap$Cl[1:700],75,"l","Cl")

for(h in hvec)
{
  postscript(file=paste("binnedwmap",h,".eps",sep=""),
             width=10,height=8)
  binned(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl")
  dev.off()
}


#----------------------------------------------------------
# Now, the local averaging example

localave = function(x,y,h,xlab="x",ylab="y",estonly=FALSE,
                    linecol=2)
{
  xg = seq(min(x),max(x),length=100)
  yhat = rep(0,length(xg))
  for(i in 1:length(xg))
  {
    yhat[i] = mean(y[x >= (xg[i]-(h/2)) & x <= (xg[i]+(h/2))])
  }
  if(!estonly)
  {
    plot(x,y,xlab=xlab,ylab=ylab,pch=".",cex=3,
         cex.axis=1.3,cex.lab=1.3)
  }
  lines(xg,yhat,lwd=3,col=linecol)
}
# Plots
binned(wmap$ell[1:700],wmap$Cl[1:700],15,"l","Cl")
localave(wmap$ell[1:700],wmap$Cl[1:700],15,"l","Cl",
         linecol=4,estonly=TRUE)
legend(0,-2000,c("Binned","Local Averaging"),col=c(2,4),
       lwd=3,cex=1.3)

for(h in hvec)
{
  postscript(file=paste("localavewmap",h,".eps",sep=""),
             width=10,height=8)
  binned(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl")
  localave(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl",
           linecol=4,estonly=TRUE)
  legend(0,-2000,c("Binned","Local Averaging"),col=c(2,4),
         lwd=3,cex=1.3)
  dev.off()
}


#----------------------------------------------------------
# Now, the kernel smoother

kernsmooth = function(x,y,h,kernel,xlab="x",ylab="y",
                      estonly=FALSE,linecol=2)
{
  xg = seq(min(x),max(x),length=100)
  yhat = rep(0,length(xg))
  for(i in 1:length(xg))
  {
    yhat[i] = sum(kernel((x-xg[i])/h)*y)/
      sum(kernel((x-xg[i])/h))
  }
  if(!estonly)
  {
    plot(x,y,xlab=xlab,ylab=ylab,pch=".",cex=3,
         cex.axis=1.3,cex.lab=1.3)
  }
  lines(xg,yhat,lwd=3,col=linecol)
}

tricube = function(x)
{
  pmax((1-abs(x)^3)^3,0)
}

for(h in hvec)
{
  postscript(file=paste("gausskernwmap",h,".eps",sep=""),
             width=10,height=8)
  localave(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl",linecol=2)
  kernsmooth(wmap$ell[1:700],wmap$Cl[1:700],h,dnorm,"l","Cl",
             estonly=TRUE,linecol=4)
  kernsmooth(wmap$ell[1:700],wmap$Cl[1:700],h,tricube,"l","Cl",
             estonly=TRUE,linecol=5)
  
  
  #   binned(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl",estonly=TRUE)
  legend(0,-2000,c("boxcar kernel","Gaussian kernel","tricube kernel"),col=c(2,4,5),
         lwd=3,cex=1.3)
  dev.off()
}


# To stress boundary bias
h = 50

postscript(file=paste("boundbiaswmap",h,".eps",sep=""),
           width=10,height=8)
plot(wmap$ell[1:100],wmap$Cl[1:100],xlab="l",ylab="Cl",
     pch=".", cex=5, cex.axis=1.3, cex.lab=1.3)
localave(wmap$ell[1:700],wmap$Cl[1:700],h,"l","Cl",
         linecol=2,estonly=TRUE)
kernsmooth(wmap$ell[1:700],wmap$Cl[1:700],h,dnorm,"l","Cl",
           estonly=TRUE,linecol=4)
kernsmooth(wmap$ell[1:700],wmap$Cl[1:700],h,tricube,"l","Cl",
           estonly=TRUE,linecol=5)
legend(0,3000,c("boxcar kernel","Gaussian kernel","tricube kernel"),col=c(2,4,5),
       lwd=3,cex=1.3)
dev.off()