presentation.Rnw

\documentclass[pdf]{beamer}

%\usepackage[italian]{babel}  % Date appears in italian and the syllables division changes to italian standard
\usepackage[utf8]{inputenc}   % It enables accented characters
\usepackage[T1]{fontenc}
\usepackage{hyperref}         % It enables hypertext links
\usepackage{url}			        % It enables url links
\usepackage{amsmath,amssymb}
\usepackage{multirow}         % It enables multirow in tabular

\usepackage{graphicx}         % It enables images
\usepackage{xcolor}
%\usepackage{subfigure}
\usepackage{caption}
\usepackage{subcaption}

%Disegni LaTeX
\usepackage{tikz}
\usetikzlibrary{shapes}
%\usetikzlibrary{shapes,snakes}
\usetikzlibrary{trees}
\usepackage{pgfplots,siunitx}
\pgfplotsset{compat=1.9}
\usepgfplotslibrary{units}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Nuovi comandi
\newcommand{\magnitude}[1]{\left\lVert #1 \right\rVert}
\DeclareMathOperator*{\argmin}{arg\,min}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usetheme{Pittsburgh}
%\mode<presentation>{\usetheme{Boadilla}}
%singapore
%\useoutertheme[right]{sidebar}
\usecolortheme{seahorse}
%\usecolortheme{dolphin}
\usefonttheme{professionalfonts}
\setbeamercovered{dynamic}
%\theoremstyle{definition}
%\newtheorem{definizione}{Definizione}
%\theoremstyle{plain}
%\newtheorem{teorema}{Teorema}

\definecolor{amber(sae/ece)}{rgb}{1.0, 0.49, 0.0}
\definecolor{ao(english)}{rgb}{0.0, 0.5, 0.0}
\definecolor{myblue}{RGB}{0,82,155}

\title{Local Polynomial Regression}
\subtitle{Statistical Machine Learning - individual project}
\author{Leonardo Stincone}

\date{18th July 2019}
\institute[units]{Università degli Studi di Trieste}

\logo{
\includegraphics[width=10mm, height=10mm]{img/logo-units_0.png}
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}


%% title frame
\begin{frame}
\titlepage
\end{frame}

<<setup, include = F, echo = F>>=
knitr::opts_chunk$set(echo = FALSE, fig.align = 'center')
@

<<libraries, include = F>>=
# Libraries ####
library(xtable)
library(mgcv)
library(splines)
library(MASS)
library(tidyverse)
library(gridExtra)
library(sm)
library(SemiPar)   # For lidar dataset

# Set plot themes
theme_set(theme_bw())

set.seed(42)
@

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \begin{frame}
% \frametitle{Table of contents}
% \tableofcontents
% \end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Problem statement: LIDAR dataset}

<<readData, cache = T>>=
data(lidar)
lidar <- lidar[sort.list(lidar$range), ]
@

\begin{columns}
\column{0.42\textwidth}
<<lidarPlot, fig.width = 4, fig.height = 4, cache = T>>=
ggplot(data = lidar,
       aes(x = range, y = logratio)) +
  geom_point() +
  labs(title = "LIDAR dataset")
@
\column{0.58\textwidth}
\small LIDAR = LIght Detection And Ranging
\begin{itemize}
\item it is a surveying method that measures distance to a target by illuminating the target with laser light and measuring the reflected light with a sensor
\item $x$: distance travelled before the light is reflected back to its source
\item $y$: logarithm of the ratio of received light from two laser sources
\end{itemize}
\end{columns}

% \bigskip
\vfill

\uncover<2->{
The objective is to estimate:
$$
f(x) = E \left[ Y \mid X=x \right]
$$
}

\end{frame}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{What does local mean?}
\begin{columns}
\column[t]{0.52\textwidth}
\only<1->{
\small If we had enough points with $x=x_0$

<<local, fig.height = 4, fig.width = 4, cache = T>>=
par(mar = c(2,1,1,1))

n = 20
h = .1

# x = sort(runif(n, 0, 1))

x0 = 0.5
x = rep(x0, n)
ff = function(x){sqrt(3*x + .5)}
m = ff(x)
y = m + rnorm(n,0,0.1)

xlim = c(0,1)
# ylim = c(min(y) - .05 * (max(y) - min(y)),
         # max(y) + .05 * (max(y) - min(y)))
ylim = c(.5, 2.1)

plot(x, y, xlim = xlim, ylim = ylim, yaxt = "n", xaxt = "n",
     xlab = "", ylab = "", yaxs = "i",
     col = "blue", pch = 20)

# axis(1, at = c(x0-h, x0+h), labels = c("", ""))
axis(1, at = x0, labels = expression(x[0]))

# rect(x0 - h, ylim[1], x0 + h, ylim[2],
#      lwd=2, col = hsv(184/360, .57, .98, alpha = 0.2), border=NA)


curve(ff(x), add = TRUE, lwd = 1,
      xlim = c(-.1, 1.1))

# ind <- abs(x-x0) < h
# points(x[ind], y[ind], col = "blue", pch = 20)

segments(x0, -1, x0, mean(y), lwd = .5, lty = 2)
points(x0, mean(y), col = "red", pch = 20, cex = 2.5)

# lines(c(x0-h, x0+h), ylim[c(1,1)],
#       lwd = 2, col = "blue")
@
}
\column[t]{0.52\textwidth}
\only<2->{
\small We can consider points "close" to $x_0$

<<neighborhood, fig.height = 4, fig.width = 4, cache = T>>=
par(mar = c(2,1,1,1))

n = 100
h = .1

x = sort(runif(n, 0, 1))
x0 = 0.5
ff = function(x){sqrt(3*x + .5)}
m = ff(x)
y = m + rnorm(n,0,0.1)

xlim = c(0,1)
# ylim = c(min(y) - .05 * (max(y) - min(y)),
#          max(y) + .05 * (max(y) - min(y)))
ylim = c(.5, 2.1)

plot(x, y, xlim = xlim, ylim = ylim, yaxt = "n", xaxt = "n",
     xlab = "", ylab = "", yaxs = "i")

axis(1, at = c(x0-h, x0+h), labels = c("", ""))
axis(1, at = x0, labels = expression(x[0]))

rect(x0 - h, ylim[1], x0 + h, ylim[2],
     lwd=2, col = hsv(184/360, .57, .98, alpha = 0.2), border=NA)


curve(ff(x), add = TRUE, lwd = 1,
      xlim = c(-.1, 1.1))

ind <- abs(x-x0) < h
points(x[ind], y[ind], col = "blue", pch = 20)

segments(x0, -1, x0, mean(y[ind]), lwd = .5, lty = 2)
points(x0, mean(y[ind]), col = "red", pch = 20, cex = 2.5)

lines(c(x0-h, x0+h), ylim[c(1,1)],
      lwd = 2, col = "blue")
@
}
\end{columns}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{2 basic ideas}
\begin{columns}
\column[t]{0.52\textwidth}
\begin{block}{$k$ nearest neighbors}
$$ \hat{f}(x) = \frac{1}{k} \sum_{i=1}^n y_i I_{N_k(x)}(x_i) $$

\end{block}
\column[t]{0.52\textwidth}
\begin{block}{Neighborhood of radius $h$}
$$\hat{f}(x) = \frac{\sum_{i=1}^n y_i I_{[0,h]}(|x-x_i|)}{\sum_{i=1}^n I_{[0,h]}(|x-x_i|)}$$

\end{block}
\end{columns}


\begin{columns}
\column[t]{0.52\textwidth}
<<lidarKnn, fig.height = 4, fig.width = 6, cache = T>>=
knn=function(x,xx,yy,k){
  mean(yy[abs(xx-x)<=(sort(abs(xx-x))[k])])
}
knn2=function(x,xx,yy,k){  
  mapply(knn,x,MoreArgs=list(xx=xx,yy=yy,k=k))
}

data_legend <- tibble(x = rep(mean(lidar$range), 4),
                      y = rep(mean(lidar$logratio),4),
                      k = factor(str_c("k = ", c(6,6,30,30)),
                                 levels = str_c("k = ", c(6,30))))

ggplot(data = lidar,
       aes(x = range, y = logratio)) +
  geom_point() +
  stat_function(fun = function(x){knn2(x,lidar$range,lidar$logratio,k=6)},
                size = 1, col = "red") +
  stat_function(fun = function(x){knn2(x,lidar$range,lidar$logratio,k=30)},
                size = 1, col = "darkgreen") +
  geom_line(data = data_legend,
            aes(x = x, y = y, color = k), size = 1) +
  scale_color_manual(values = c("red", "darkgreen")) +
  labs(title = "LIDAR dataset",
       subtitle = "k nearest neighbors")
@
\column[t]{0.52\textwidth}
<<lidarRadiusH, fig.height = 4, fig.width = 6, cache = T>>=
runmean=function(x,xx,yy,h){
  mean(yy[abs(xx-x)<=h])
}
runmean2=function(x,xx,yy,h){  
  mapply(runmean,x,MoreArgs=list(xx=xx,yy=yy,h=h))
}

data_legend <- tibble(x = rep(mean(lidar$range), 4),
                      y = rep(mean(lidar$logratio),4),
                      h = factor(str_c("h = ", c(4,4,20,20)),
                                 levels = str_c("h = ", c(4,20))))

ggplot(data = lidar,
       aes(x = range, y = logratio)) +
  geom_point() +
  stat_function(fun = function(x){runmean2(x,lidar$range,lidar$logratio,h=4)},
                size = 1, col = "red") +
  stat_function(fun = function(x){runmean2(x,lidar$range,lidar$logratio,h=20)},
                size = 1, col = "darkgreen") +
  geom_line(data = data_legend,
            aes(x = x, y = y, color = h), size = 1) +
  scale_color_manual(values = c("red", "darkgreen")) +
  labs(title = "LIDAR dataset",
       subtitle = "Neighborhood of radius h")
@
\end{columns}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Nadaraya-Watson kernel regression}

$$
\hat{f}(x) = \sum_{i=1}^n \ell_i(x) y_i
$$

with:
$$
\ell_i(x) = \frac{K\left(\frac{x-x_i}{h}\right)}{\sum_{j=1}^n K\left(\frac{x-x_j}{h}\right)}
$$

where $K(\cdot)$ is a kernel function that satisfies:
\begin{itemize}
\item $ K(x)\geq 0$
\item $ \int K(x)dx=1$
\item $ \int xK(x)dx=0$
\item $ \int x^2K(x)dx>0$
\end{itemize}


% $$
% I_{[0,h]}(|x-x_i|) = I_{[-1,1]}\left(\frac{x-x_i}{h}\right)
% $$

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Some proposed kernels}

<<plotKernel, cache = T>>=
ylim = c(0, 1.2)

plot_kernel = function(ff){
  plot(0, 0, type = "n",
       xlim = c(-2,2), ylim = ylim)
  segments(-3, 0, 3, 0, lwd = 1, lty = 2)
  curve(ff(x), add = T, lwd = 3, col = "red")
}
@

\resizebox{\textwidth}{!}{
\begin{tabular}{cccc}
Uniform &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/uniform-1}
}%
& 
Triangle &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/triangle-1}
}%
\\
$\frac{1}{2}I_{[-1,1]}(u)$ & &  $(1-|u|)I_{[-1,1]}(u)$ & \\ & & & \\ & & & \\
Triweight &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/triweight-1}
}%
& 
Quartic &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/quartic-1}
}%
\\
 $\frac{35}{32}(1-u^2)^3I_{[-1,1]}(u)$ & &   $\frac{15}{16}(1-u^2)^2I_{[-1,1]}(u)$ & \\ & & & \\ & & & \\
Cosine &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/cosine-1}
}%
& 
Epanechnikov &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/epanechnikov-1}
}%
\\
$\frac{\pi}{4}\cos\left(\frac{\pi}{2}u\right)I_{[-1,1]}(u)$ & &    $\frac{3}{4}(1-u^2)I_{[-1,1]}(u)$  & \\ & & & \\ & & & \\
Gaussian &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/gaussian-1}
}%
& 
Tricube &
\multirow{4}{*}{
\includegraphics[width=0.18\textwidth]{figure/tricube-1}
}%
\\
$\frac{1}{\sqrt{2\pi}} e^{-u^2/2}$ & & $\frac{70}{81}(1-|u|^3)^3I_{[-1,1]}(u)$ & \\ & & & \\ & & & \\

\end{tabular}
}


\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Design bias, boundary bias and concavity bias}

<<nw_bias, fig.height = 4, fig.width = 8>>=
par(mar = c(2, .5, 1, .5))
laymat = matrix(c(1:8), byrow = FALSE, nrow = 2)
laymat = laymat[c(1, 1, 2), ]
layout(laymat)

ungraf2 = function(x, y, ff, bw=1,
                   miny = min(y)-0.05*(max(y)-min(y)),
                   maxy = max(y)+0.05*(max(y)-min(y))){
  plot(x, y, ylim = c(miny, maxy), yaxs = "i", yaxt = "n",
       xlab = "", ylab = "", col = gray(0.7))
  rug(x, line = 0.5)
  
  curve(ff(x), add = TRUE, lwd = 2, n = 100, col = "black") # True regression
  a = ksmooth(x, y, bandwidth = bw)
  lines(a$x, a$y, lwd = 2, col = "darkgreen") # Nadaraya–Watson
  # rect(min(x1),-0.1,max(x1),1,lwd=2,col=hsv(184/360,.57,.98,alpha=0.2),border=NA)
  # lines(range(x1),-0.1*c(1,1),lwd=2,col="blue")
  # sig=0.25*bw/qnorm(0.75)
  # curve(miny+0.1*(maxy-miny)*dnorm(x,0.5,sig)/dnorm(0.5,0.5,sig),add=TRUE,lwd=2,col="lightblue")
  # hist(x,freq=FALSE,main="",yaxt="n",col=gray(0.7),border="white")
}

# set.seed(42)
n = 200

# Computing linear

x1_unif = seq(0, 1, length = n) #sort(runif(n,0,1))
ff1 = function(x){0.5 + x}
m = ff1(x1_unif)
y1_unif = m + rnorm(n, 0, 0.025)

x = qnorm(seq(0, 1, length = (n + 2))[2:(n + 1)]) #sort(rnorm(n,0,1))
min = min(x)
max = max(x)
x1_norm = (x - min(x)) / (max(x) - min(x))
# ff1 = function(x){0.5 + x}
m = ff1(x1_norm)
y1_norm = m + rnorm(n, 0, 0.025)



# Plotting linear

miny1 = min(y1_unif, y1_norm) - 0.05*(max(y1_unif, y1_norm) - min(y1_unif, y1_norm))
maxy1 = max(y1_unif, y1_norm) + 0.05*(max(y1_unif, y1_norm) - min(y1_unif, y1_norm))

ungraf2(x = x1_unif, y = y1_unif, ff = ff1, bw = 0.5,
        miny = miny1, maxy = maxy1)
curve(dunif(x),
      yaxt = "n")

ungraf2(x = x1_norm, y = y1_norm, ff = ff1, bw = 0.5,
        miny = miny1, maxy = maxy1)
curve(dnorm(x, mean = .5, sd = 1/(max - min)),
      yaxt = "n")



# Computing quadratic

x2_unif = seq(0, 1, length = n) #sort(runif(n,0,1))
ff2 = function(x){(x - 0.5)^2}
m = ff2(x2_unif)
y2_unif = m + rnorm(n,0,0.025)

x2_norm = qnorm(seq(0, 1, length = (n + 2))[2:(n + 1)]) #sort(rnorm(n,0,1))
min = min(x)
max = max(x)
x2_norm = (x - min(x)) / (max(x) - min(x))
# ff2 = function(x){(x - 0.5)^2}
m = ff2(x2_norm)
y2_norm = m + rnorm(n, 0, 0.025)


# Plotting quadratic

miny2 = min(y2_unif, y2_norm) - 0.05*(max(y2_unif, y2_norm) - min(y2_unif, y2_norm))
maxy2 = max(y2_unif, y2_norm) + 0.05*(max(y2_unif, y2_norm) - min(y2_unif, y2_norm))

ungraf2(x = x2_unif, y = y2_unif, ff = ff2, bw = 0.4,
        miny = miny2, maxy = maxy2)
curve(dunif(x),
      yaxt = "n")

ungraf2(x = x2_norm, y = y2_norm, ff = ff2, bw = 0.4,
        miny = miny2, maxy = maxy2)
curve(dnorm(x, mean = .5, sd = 1/(max - min)),
      yaxt = "n")

@



\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Design bias: what's happening?}

<<designBias, out.width='\\textwidth', fig.width = 6, fig.height = 5, cache = T>>=
n = 50
h = .1

# x = sort(runif(n, 0, 1))
mean = .4
sd = .04
x = sort(rnorm(n, mean, sd))
x0 = 0.5
ff = function(x){sqrt(3*x + .5)}
m = ff(x)
y = m + rnorm(n,0,0.025)

xlim = c(0,1)
# ylim = c(min(y) - .05 * (max(y) - min(y)),
#          max(y) + .05 * (max(y) - min(y)))
ylim = c(.6, 2)

plot(x, y, xlim = xlim, ylim = ylim, yaxt = "n", xaxt = "n",
     xlab = "", ylab = "", yaxs = "i")

axis(1, at = c(x0-h, x0+h), labels = c("", ""))
axis(1, at = x0, labels = expression(x[0]))

rect(x0 - h, ylim[1], x0 + h, ylim[2],
     lwd=2, col = hsv(184/360, .57, .98, alpha = 0.2), border=NA)

curve(ff(x), add = TRUE, lwd = 1,
      xlim = c(-.1, 1.1))

ind <- abs(x-x0) < h
points(x[ind], y[ind], col = "blue", pch = 20)

segments(x0, -1, x0, mean(y[ind]), lwd = .5, lty = 2)
points(x0, mean(y[ind]), col = "red", pch = 20, cex = 3)


curve(ylim[1] + 0.2 * (ylim[2] - ylim[1]) * dnorm(x, mean, sd) / dnorm(0,0,sd),
      add = TRUE, lwd = 2, col = "lightblue")

lines(c(x0-h, x0+h), ylim[c(1,1)],
      lwd = 2, col = "blue")
@

\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: idea}

Locally the Nadaraya-Watson estimator is a Weighted Least Square Estimator:
$$
\hat{f}_{NW}(x_0) = \underset{a}{\mbox{argmin}} \sum_{i=1}^n K\left(\frac{x_i-x_0}{h}\right)(y_i - a)^2
$$

\hfill

\uncover<2->{
Idea: instead of approximating $f(x_0)$ with a constant value $a$, we could approximate it with a polynomial $p_{x_0}(u, \boldsymbol{a})$.

\hfill

Taylor polynomial approximation:
$$
p_{x_0}(u; \boldsymbol{a}) = a_0 + a_1(u-x_0) + \frac{a_2}{2!}(u-x_0)^2 +\ldots + \frac{a_d}{d!}(u-x_0)^d
$$
}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: estimation}

We can estimate the coefficients of $p_{x_0}(u;\boldsymbol{a})$ as:
$$
\hat{\boldsymbol{a}}(x_0) = \underset{\boldsymbol{a}}{\mbox{argmin}}
\sum_{i=1}^n K\left(\frac{x_i-x_0}{h}\right)\left(y_i - p_{x_0}(x_i;{\boldsymbol{a}})\right)^2
$$

Thus, we can define the estimator for $f(x)$ in $x_0$ just computing $p_{x_0}(u;\hat{\boldsymbol{a}})$ in $x_0$:

$$
\hat{f}(x_0) = p_{x_0}(x_0; \hat{\boldsymbol{a}})
$$

Then we can repeat the process for each value of $x$ in a grid and obtain $\hat{f}(x)$.

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: example}
<<fitCustomLoess_1, cache = T>>=
n <- 100
x <- sort(runif(n,0,1))
truefun <- function(x){sin(x*2*pi) * (x<0.5) + 0}
y <- truefun(x) + rnorm(n,0,0.1)


grid = seq(from = 0, to = 1, by = .001)
# grid = seq(from = 0, to = 1, by = .1)


# Linear loess
my_loess <- function(x, y, grid, span){
  stima <- rep(NA, length(grid))
  
  a0 <- rep(NA, length(grid))
  a1 <- rep(NA, length(grid))
  
  # Computing
  for (i in 1:length(grid)){
    fit <- lm(y~x, weights = dnorm(x, grid[i], span))
    a0[i] <- coef(fit)[1]
    a1[i] <- coef(fit)[2]
    stima[i] = predict(fit, data.frame(x = grid[i]))
  }
  return(list(stima = stima,
              a0 = a0,
              a1 = a1))
}


# Fit models

span1 <- .1
span2 <- .04

interesting_points <- c(.1, .25, .4, .8)


fit_myloess_1 <- my_loess(x, y, grid, span = span1)
fit_myloess_2 <- my_loess(x, y, grid, span = span2)



# Fit 1
stima <- fit_myloess_1$stima
a0 <- fit_myloess_1$a0
a1 <- fit_myloess_1$a1


# Create dataframes
data_fit1_grid <- tibble(grid,
                    a0, a1,
                    stima)

data_fit1_grid2 <- data_fit1_grid %>% 
  select(x = grid, y = stima)


data <- tibble(x, y)

data_fit1_grid_small <- data_fit1_grid %>% 
  filter(grid %in% interesting_points)

data_fit1_cross_small <- crossing(data_fit1_grid_small, data)



# Fit 2
stima <- fit_myloess_2$stima
a0 <- fit_myloess_2$a0
a1 <- fit_myloess_2$a1

# Create dataframes
data_fit2_grid <- tibble(grid,
                         a0, a1,
                         stima)

data_fit2_grid2 <- data_fit2_grid %>% 
  select(x = grid, y = stima)


# data <- tibble(x, y)

data_fit2_grid_small <- data_fit2_grid %>% 
  filter(grid %in% interesting_points)

data_fit2_cross_small <- crossing(data_fit2_grid_small, data)



# Union of datasets
data_grid_small <- rbind(data_fit1_grid_small %>% 
                           mutate(span = span1),
                         data_fit2_grid_small %>% 
                           mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))

data_cross_small <- rbind(data_fit1_cross_small %>%
                            mutate(span = span1),
                          data_fit2_cross_small %>% 
                            mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))

data_grid2 <- rbind(data_fit1_grid2 %>% 
                      mutate(span = span1),
                    data_fit2_grid2 %>% 
                      mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))
@

<<fitCustomLoess_2, cache = T>>=
# QUadratic loess
my_loess_quad <- function(x, y, grid, span){
  stima <- rep(NA, length(grid))
  
  a0 <- rep(NA, length(grid))
  a1 <- rep(NA, length(grid))
  a2 <- rep(NA, length(grid))
  
  # Computing
  for (i in 1:length(grid)){
    fit <- lm(y ~ x + I(x^2), weights = dnorm(x, grid[i], span))
    a0[i] <- coef(fit)[1]
    a1[i] <- coef(fit)[2]
    a2[i] <- coef(fit)[3]
    stima[i] = predict(fit, data.frame(x = grid[i]))
  }
  return(list(stima = stima,
              a0 = a0,
              a1 = a1,
              a2 = a2))
}


# Fit models

fit_myloess_quad_1 <- my_loess_quad(x, y, grid, span = span1)
fit_myloess_quad_2 <- my_loess_quad(x, y, grid, span = span2)



# Fit 1
stima <- fit_myloess_quad_1$stima
a0 <- fit_myloess_quad_1$a0
a1 <- fit_myloess_quad_1$a1
a2 <- fit_myloess_quad_1$a2


# Create dataframes
data_fit_quad1_grid <- tibble(grid,
                         a0, a1, a2,
                         stima)

data_fit_quad1_grid2 <- data_fit_quad1_grid %>% 
  select(x = grid, y = stima)


data <- tibble(x, y)

data_fit_quad1_grid_small <- data_fit_quad1_grid %>% 
  filter(grid %in% interesting_points)

data_fit_quad1_cross_small <- crossing(data_fit_quad1_grid_small, data)



# Fit 2
stima <- fit_myloess_quad_2$stima
a0 <- fit_myloess_quad_2$a0
a1 <- fit_myloess_quad_2$a1
a2 <- fit_myloess_quad_2$a2

# Create dataframes
data_fit_quad2_grid <- tibble(grid,
                         a0, a1, a2,
                         stima)

data_fit_quad2_grid2 <- data_fit_quad2_grid %>% 
  select(x = grid, y = stima)


# data <- tibble(x, y)

data_fit_quad2_grid_small <- data_fit_quad2_grid %>% 
  filter(grid %in% interesting_points)

data_fit_quad2_cross_small <- crossing(data_fit_quad2_grid_small, data)



# Union of datasets
data_quad_grid_small <- rbind(data_fit_quad1_grid_small %>% 
                           mutate(span = span1),
                         data_fit_quad2_grid_small %>% 
                           mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))

data_quad_cross_small <- rbind(data_fit_quad1_cross_small %>%
                            mutate(span = span1),
                          data_fit_quad2_cross_small %>% 
                            mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))

data_quad_grid2 <- rbind(data_fit_quad1_grid2 %>% 
                      mutate(span = span1),
                    data_fit_quad2_grid2 %>% 
                      mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))

data_parabola_1 <- crossing(data_fit_quad1_grid_small, tibble(x = seq(from = -.1, to = 1.1, by = 0.001))) %>% 
  mutate(y = a0 + a1*x + a2*x^2)

data_parabola_2 <- crossing(data_fit_quad2_grid_small, tibble(x = seq(from = -.1, to = 1.1, by = 0.001))) %>% 
  mutate(y = a0 + a1*x + a2*x^2)

data_parabola <- rbind(data_parabola_1 %>% 
                         mutate(span = span1),
                       data_parabola_2 %>% 
                         mutate(span = span2)) %>% 
  mutate(span_str = str_c("span = ", span))
@


<<LoessExplanation_1, fig.height = 4, fig.width = 8, cache = T>>=
data_fit1_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  # stat_function(fun = truefun, size = 1) +
  geom_point(data = data_fit1_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span1) / dnorm(0, mean = 0, sd = span1))) +
  # geom_line(data = data_fit1_grid2,
  #           aes(x = x, y = y),
  #           col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  geom_abline(aes(intercept = a0, slope = a1),
              color = "red") +
  facet_grid(~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "")
@



\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: example}
<<LoessExplanation_2, fig.height = 4, fig.width = 8, cache = T>>=
data_fit1_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  # stat_function(fun = truefun, size = 1) +
  geom_point(data = data_fit1_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span1) / dnorm(0, mean = 0, sd = span1))) +
  geom_line(data = data_fit1_grid2,
            aes(x = x, y = y),
            col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  # geom_abline(aes(intercept = a0, slope = a1),
  # color = "red") +
  facet_grid(~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "")
@
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: example}
<<LoessExplanation_3, fig.height = 4, fig.width = 8, cache = T>>=
data_fit1_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  geom_point(data = data_fit1_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span1) / dnorm(0, mean = 0, sd = span1))) +
  stat_function(fun = truefun, size = 1) +
  geom_line(data = data_fit1_grid2,
            aes(x = x, y = y),
            col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  # geom_abline(aes(intercept = a0, slope = a1),
  #             color = "red") +
  facet_grid(~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "")
@

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: example}
<<LoessExplanation_4, fig.height = 4, fig.width = 8, cache = T>>=
data_fit1_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  geom_point(data = data_fit1_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span1) / dnorm(0, mean = 0, sd = span1))) +
  stat_function(fun = truefun, size = 1) +
  geom_line(data = data_fit1_grid2,
            aes(x = x, y = y),
            col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  geom_abline(aes(intercept = a0, slope = a1),
              color = "red") +
  facet_grid(~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "")
@

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Design bias, boundary bias and concavity bias}
<<loess_bias, fig.height = 4, fig.width = 8, cache = T>>=
par(mar = c(2, .5, 1, .5))
laymat = matrix(c(1:8), byrow = FALSE, nrow = 2)
laymat = laymat[c(1, 1, 2), ]
layout(laymat)

ungraf3 = function(x, y, ff, bw=1,
                   miny = min(y)-0.05*(max(y)-min(y)),
                   maxy = max(y)+0.05*(max(y)-min(y))){
  plot(x, y, ylim = c(miny, maxy), yaxs = "i", yaxt = "n",
       xlab = "", ylab = "", col = gray(0.7))
  rug(x, line = 0.5)
  
  curve(ff(x), add = TRUE, lwd = 2, n = 100, col = "black") # True regression
  
  # The quartile are at +-0.25*bw
  a = ksmooth(x, y, bandwidth = bw)
  lines(a$x, a$y, lwd = 2, col = "darkgreen") # Nadaraya–Watson
  
  sig = 0.25 * bw / qnorm(0.75) # Standard deviation of the normal kernel
  a = sm.regression(x, y, h = sig, display = 'none')
  lines(a$eval.points, a$estimate, lwd = 2, col = "blue") # LOESS
  
  # rect(min(x1),-0.1,max(x1),1,lwd=2,col=hsv(184/360,.57,.98,alpha=0.2),border=NA)
  # lines(range(x1),-0.1*c(1,1),lwd=2,col="blue")
  # sig=0.25*bw/qnorm(0.75)
  # curve(miny+0.1*(maxy-miny)*dnorm(x,0.5,sig)/dnorm(0.5,0.5,sig),add=TRUE,lwd=2,col="lightblue")
  # hist(x,freq=FALSE,main="",yaxt="n",col=gray(0.7),border="white")
}


# Plotting linear

miny1 = min(y1_unif, y1_norm) - 0.05*(max(y1_unif, y1_norm) - min(y1_unif, y1_norm))
maxy1 = max(y1_unif, y1_norm) + 0.05*(max(y1_unif, y1_norm) - min(y1_unif, y1_norm))

ungraf3(x = x1_unif, y = y1_unif, ff = ff1, bw = 0.5,
        miny = miny1, maxy = maxy1)
curve(dunif(x),
      yaxt = "n")

ungraf3(x = x1_norm, y = y1_norm, ff = ff1, bw = 0.5,
        miny = miny1, maxy = maxy1)
curve(dnorm(x, mean = .5, sd = 1/(max - min)),
      yaxt = "n")


# Plotting quadratic

miny2 = min(y2_unif, y2_norm) - 0.05*(max(y2_unif, y2_norm) - min(y2_unif, y2_norm))
maxy2 = max(y2_unif, y2_norm) + 0.05*(max(y2_unif, y2_norm) - min(y2_unif, y2_norm))

ungraf3(x = x2_unif, y = y2_unif, ff = ff2, bw = 0.4,
        miny = miny2, maxy = maxy2)
curve(dunif(x),
      yaxt = "n")

ungraf3(x = x2_norm, y = y2_norm, ff = ff2, bw = 0.4,
        miny = miny2, maxy = maxy2)
curve(dnorm(x, mean = .5, sd = 1/(max - min)),
      yaxt = "n")
@
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local Polynomial Regression on LIDAR dataset}
<<lidarLoess, fig.height = 6, fig.width = 9, cache = T>>=
data_legend <- tibble(x = rep(mean(lidar$range), 6),
                      y = rep(mean(lidar$logratio),6),
                      degree = as.character(c(0,0,1,1,2,2)))

ggplot(data = lidar,
       aes(x = range, y = logratio)) +
  geom_point() +
  geom_smooth(method = "loess", se = F, span = .7,
              method.args = list(degree = 2), color = "red", size = 1) +
  geom_smooth(method = "loess", se = F, span = .3,
              method.args = list(degree = 1), color = "blue", size = 1) +
  geom_smooth(method = "loess", se = F, span = .3,
              method.args = list(degree = 0), color = "darkgreen", size = 1) +
  geom_line(data = data_legend,
            aes(x = x, y = y, color = degree), size = 1) +
  scale_color_manual(values = c("darkgreen", "blue", "red")) +
  labs(title = "LIDAR dataset",
       subtitle = "Local Polynomial Regression")
@
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Local polynomial regression: linear vs quadratic}

<<LoessComparison, fig.height = 8, fig.width = 12, cache = T>>=
p_loess_linear <- data_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  geom_point(data = data_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span) / dnorm(0, mean = 0, sd = span))) +
  stat_function(fun = truefun, size = 1) +
  geom_line(data = data_grid2,
            aes(x = x, y = y),
            col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  geom_abline(aes(intercept = a0, slope = a1),
              color = "red") +
  facet_grid(span_str~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "", title = "Linear loess")


p_loess_quadratic <- data_quad_grid_small %>% 
  ggplot(aes(x = grid, y = stima)) +
  geom_point(data = data_quad_cross_small,
             aes(x = x, y = y,
                 color = dnorm(x, mean = grid, sd = span) / dnorm(0, mean = 0, sd = span))) +
  stat_function(fun = truefun, size = 1) +
  geom_line(data = data_quad_grid2,
            aes(x = x, y = y),
            col = "brown", size = 1) +
  geom_point(color = "red", size = 2) +
  geom_line(data = data_parabola,
            aes(x = x, y = y),
            color = "red") +
  facet_grid(span_str~grid) +
  theme(legend.position = "none") +
  scale_color_gradient(low = "gray90", high = "black") +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.3, 1.2)) +
  labs(x = "", y = "", title = "Quadratic loess")


grid.arrange(p_loess_linear,
             p_loess_quadratic,
             ncol = 1)
@

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Bibliography}

\begin{columns}
\column[t]{0.5\textwidth}
\small
\textbf{Fan, Gijbels} \\
\textit{Local Polynomial Modelling and its Applications} \\
Springer (1996)
\column[t]{0.5\textwidth}
\small
\textbf{Ruppert, Wand, Carroll} \\
\textit{Semiparametric Regression} \\
Cambridge University Press (2003)
\end{columns}

\begin{columns}
\column{0.5\textwidth}
\begin{center}
\includegraphics[width=0.6\textwidth]{img/cover_localPolynomialModellingAndItsApplication.png}
\end{center}
\column{0.5\textwidth}
\begin{center}
\includegraphics[width=0.6\textwidth]{img/cover_semiparametriRegression.png}
\end{center}
\end{columns}

\end{frame}


\end{document}