\[ CV_{k} = \sum_{k = 1}^{K} \frac{n_k}{n} Err_k,\]
Onde \(Err_k = \sum_{i \in C_k} \frac{(y_i \neq \hat y_i)}{n_k}\), para quando a variável é discreta, representando a taxa de classificação incorreta.
Problemas que envolvem modelagem são aqueles nos quais temos interesse em estimar a relação entre uma ou mais variáveis resposta, denotadas por \(y\) e uma ou mais variáveis explicativas \(x\). Uma das abordagens possíveis na resolução deste tipo de problema é a suavização de uma curva por Kernel.
Em geral, um Kernel suavizador define um conjunto de pesos \({W_i(x)}_{i=1}^{n}\) para cada \(x\), e \[ \hat f(x) = \sum_{i = 1}^{n} W_i(x) y_i.\]
Para um Kernel suavizador, temos que a sequência de pesos \({W_i(x)}_{i=1}^{n}\) é representada ao descrever a forma da função de pesos \(W_i(x)\) por uma função de densidade, com um parâmetro de escala que ajusta o tamanho e forma dos pesos perto de \(x\). Assim, refere-se à essa função de forma como um kernel K, onde: \[ \int K(u)d(u) = 1.\]
O parâmetro de escala denominado por \(\lambda\) produz a sequência de pesos: \[ W_{hi}(x) = \frac{K \frac{(x - x_i)}{\lambda}} {\sum_{i = 1}^{n} K \frac{(x - x_i)}{\lambda}} \]
É crucial para a performance do estimador que a largura de banda, ou \(\lambda\), seja bem definida. Uma das abordagens mais utilizadas para encontrar este valor é a LOOCV (Leave-One-Out Cross Validation), de forma a minimizar: \[ EQ = \sum_{i = 1}^{n} (y_i - \hat f(x_i))^2 + \lambda \int \hat f''(x_i) dx\]
library(labestData)
library(lattice)
library(latticeExtra)
da <- data.frame(y = PaulaEx1.13.19$renda,
x = PaulaEx1.13.19$estud)
# Qual o valor ótimo para o lambda?
p1 <- xyplot(y ~ x, data = da, type = c("p", "g"),
col = "pink3", xlab = "x - renda",
ylab = "y - tempo de estudo",
pch = 16,
main = expression(lambda~"=1"))+
layer(panel.lines(ksmooth(da$x, da$y, "normal", 1),
lwd = 2,
col = "turquoise"))
p2 <- xyplot(y ~ x, data = da, type = c("p", "g"),
col = "pink3", xlab = "x - renda",
ylab = "y - tempo de estudo",
pch = 16,
main = expression(lambda~"=1.5"))+
layer(panel.lines(ksmooth(da$x, da$y, "normal", 1.5),
lwd = 2,
col = "turquoise"))
p3 <- xyplot(y ~ x, data = da, type = c("p", "g"),
col = "pink3", xlab = "x - renda",
ylab = "y - tempo de estudo",
pch = 16,
main = expression(lambda~"=2"))+
layer(panel.lines(ksmooth(da$x, da$y, "normal", 2),
lwd = 2,
col = "turquoise"))
print(p1, position = c(0, 0, 1/3, 1), more = TRUE)
print(p2, position = c(1/3, 0, 2/3, 1), more = TRUE)
print(p3, position = c(2/3, 0, 1, 1))
#---------------------
(n <- nrow(da))
[1] 50
k <- 5
da$i <- ceiling((1:n)/(n/k))
# Parte os dados em k conjuntos disjuntos.
das <- split(da, f = da$i)
cen <- expand.grid(fold = 1:k,
l = 3:250)
kfol <- mapply(f = cen$fold,
d = cen$l,
FUN = function(f, d) {
j <- da$i != f
b <- ksmooth(subset(da, j)$x,
subset(da, j)$y,
"normal", bandwidth = d)
# Erro de treinamento
cvt <- sum((subset(da, j)$y - b$y)^2,
na.rm = TRUE)/length(subset(da, j)$y)
#--------------------
# Erro de validação
p <- b$y[which.min(abs(b$x-das[[f]]$x))]
cvval <- crossprod((das[[f]]$y -
p)^2)/length(das[[f]]$y)
return(c(cvt = cvt, cvval = cvval))
})
kfol <- cbind(cen, as.data.frame(t(kfol)))
str(kfol)
'data.frame': 1240 obs. of 4 variables:
$ fold : int 1 2 3 4 5 1 2 3 4 5 ...
$ l : int 3 3 3 3 3 4 4 4 4 4 ...
$ cvt : num 1274371 1571846 1545997 1693674 1892323 ...
$ cvval: num 9.15e+12 3.67e+11 6.57e+11 7.58e+11 6.96e+10 ...
xyplot(cvt + cvval ~ l, groups = fold , data = kfol,
auto.key = list(corner = c(0.95, 0.95)),
type = c("p", "g", "o"),
scales = "free",
xlab = expression(lambda))
#-------------
which.min(kfol$cvval)
[1] 5
kfol[5,] # lambda = 3
fold l cvt cvval
5 5 3 1892323 69568605954
#-----------------------------------------------------------------------
# leave-one-out
# Cenários
cen <- expand.grid(fold = 1:n,
l = 3:15)
# Obtendo os erros para cada cenário
nfol <- mapply(f = cen$fold,
d = cen$l,
FUN = function(f, d) {
# browser()
b <- ksmooth(da[-f,]$x,
da[-f,]$y,
"normal", bandwidth = d)
# Erro de treinamento
cvt <- sum((da[-f]$y - b$y)^2,
na.rm = TRUE)/length(da[-f]$y)
#--------------------
# Erro de validação
p <- b$y[which.min(abs(b$x-da[f, ]$x))]
cvval <- (p - da$y[f])^2
return(c(cvt = cvt, cvval = cvval))
})
nfol <- cbind(cen, as.data.frame(t(nfol)))
which.min(nfol$cvval)
[1] 585
nfol[585,] # lambda = 14
fold l cvt cvval
585 35 14 1042973 5.752992
xyplot(cvt + cvval ~ l, groups = fold , data = nfol,
type = c("p", "g", "o"),
scales = "free",
xlab = expression(lambda))