Aprimora sobre validação cruzada.

parent 2a227488
---
title: "Validação Cruzada"
author: Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira
date: 2017-08-24
date: 2018-08-24
#bibliography: ../config/Refs.bib
#csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
......@@ -16,6 +16,20 @@ opts_chunk$set(
# Definindo o modelo
Para demonstrar o uso funcionamento da validação cruzada, será usado um
conjunto de dados simulado apenas pela conveniência de serem conhecidos
parâmetros.
As observações serão provenientes do modelo
$$
y = \beta_{0} + \beta_{1} x + \beta_{3} \sin(x/\beta_{4}) + \epsilon,
$$
em que $\beta = [\beta_0, \beta_1, \beta_2, \beta_3]$ são os parâmetros
de regressão e $\epsilon$ um erro normal de média 0 e variância
$\sigma^2$. Os reais valores usados para esses parâmetros na simulação
estão expostos dentro do código.
```{r, messages = FALSE}
#-----------------------------------------------------------------------
# Pacotes.
......@@ -48,10 +62,22 @@ da$y <- eta(da$x) + epsilon(da$x)
# Exibindo os dados.
plot(y ~ x, data = da)
curve(eta, add = TRUE, col = "orange", lwd = 3)
legend("bottomright",
bty = "n",
legend = c("Sinal verdadeiro", "Valores observados"),
lty = c(1, 0),
lwd = c(3, 1),
col = c("orange", "black"),
pch = c(NA, 1))
```
# Validação cruzada *holdout*
Na validação cruzada *holdout* é deixado de fora, para avaliação do
desempenho de predição, uma porporção dos dados. Com a porção
complementar, um modelo é ajustado aos dados, ocasião na qual são
estimadas as quantidade livres do modelo, ou seja, os parâmetros.
```{r}
#-----------------------------------------------------------------------
# Partindo os dados em treino e validação.
......@@ -69,11 +95,13 @@ str(das)
# Criando uma lista de fórmulas para repetir os dados nos paineis.
dmax <- 15
form <- replicate(dmax, y ~ x)
names(form) <- 1:dmax
names(form) <- sprintf("pol. grau: %d", 1:dmax)
# Diagramas de dispersão dos dados um ajuste variando o grau.
xyplot.list(form,
data = das[["tra"]],
cex = 0.6,
layout = c(3, 5),
as.table = TRUE) +
layer(panel.smoother(form = y ~ poly(x, degree = panel.number()),
method = lm)) +
......@@ -81,6 +109,7 @@ xyplot.list(form,
layer(panel.points(x = x,
y = y,
pch = 19,
cex = 0.6,
col = "red"),
data = das[["val"]])
......@@ -89,17 +118,16 @@ deg <- 1:dmax
cv <- sapply(deg,
FUN = function(d) {
# Ajuste com os dados de treino.
m0 <- lm(y ~ poly(x, degree = d),
data = das[["tra"]])
# Predição para os dados de validação.
yhat <- predict(m0,
newdata = das[["val"]])
# Erro de treino.
# crossprod(fitted(m0) - das[["tra"]]$y)/ntra
Etra <- deviance(m0)/ntra
# Erro de validação.
Eval <- crossprod(yhat - das[["val"]]$y)/nval
c(ErrTra = Etra, ErrVal = Eval)
m0 <- lm(y ~ poly(x, degree = d), data = das[["tra"]])
# Predição com os dados de treino e de teste.
yhat_inter <- predict(m0, newdata = das[["tra"]]) # Dados de treino.
yhat_exter <- predict(m0, newdata = das[["val"]]) # Dados de teste.
# Erros de predição
Etra <- crossprod(das[["tra"]]$y - yhat_inter)/ntra
Eval <- crossprod(das[["val"]]$y - yhat_exter)/ntra
# Retorno.
c(ErrTra = Etra,
ErrVal = Eval)
})
cv <- cbind(deg, as.data.frame(t(cv)))
......@@ -113,12 +141,21 @@ xyplot(ErrTra + ErrVal ~ deg,
# Validação cruzada com $k$-*folds*
Na validação cruzada do tipo $k$-folds, os dados são dividos nos
registros em 5 conjuntos de tamanho aproximadamente igual (igual quando
$n$ for múltiplo de $k$).
Ao todo são $k$ situações no qual se faz exatamente o que foi feito no
procedimento *holdout*. Todavia, com essa divisão, cada registro ficará
como observação no conjunto de teste pelo menos uma vez. Isso faz com
que a avaliação do modelo tenha maior cobertura em termos de condições.
```{r}
#-----------------------------------------------------------------------
# Método k-fold.
# Gerando um conjunto maior de observações do mesmo modelo.
da <- data.frame(x = runif(200, 0, 100))
da <- data.frame(x = runif(207, 0, 100))
da$y <- eta(da$x) + epsilon(da$x)
nrow(da)
......@@ -127,7 +164,10 @@ n <- nrow(da)
k <- 10
da$i <- ceiling((1:n)/(n/k))
# Parte os dados em k conjuntos disjuntos.
# Número de registros por fold.
table(da$i)
# Parte os dados em k conjuntos disjuntos (gera uma lista).
das <- split(da, f = da$i)
str(das)
......@@ -138,6 +178,9 @@ names(form) <- sprintf("fold %d", 1:k)
xyplot.list(form,
data = da,
type = "n",
xlab = "Grau do polinômio",
ylab = "Erro",
cex = 0.6,
as.table = TRUE) +
# Curva do modelo verdadeiro.
layer(panel.curve(eta, col = "orange", lwd = 3)) +
......@@ -148,18 +191,18 @@ xyplot.list(form,
layer(panel.points(x = x[da$i == packet.number()],
y = y[da$i == packet.number()],
pch = 19,
cex = 0.6,
col = "red")) +
# Ajuste do modelo aos dados de treinamento.
layer(panel.smoother(x = x[da$i != packet.number()],
y = y[da$i != packet.number()],
form = y ~ poly(x, degree = 5),
method = lm)) +
# # Ajuste do modelo aos dados de treinamento.
# layer(panel.smoother(x = x[da$i != packet.number()],
# y = y[da$i != packet.number()],
# form = y ~ poly(x, degree = 5),
# method = lm)) +
# Linhas de referência.
layer(panel.grid(), under = TRUE)
# Avaliando cenários varrendo os fold e variando o grau do polinômio.
cen <- expand.grid(fold = 1:k,
deg = 1:15)
cen <- expand.grid(fold = 1:k, deg = 1:15)
# Avaliando cada cenário.
kfol <- mapply(f = cen$fold,
......@@ -167,27 +210,32 @@ kfol <- mapply(f = cen$fold,
FUN = function(f, d) {
# Ajuste do modelo aos dados de treinamento.
j <- da$i != f
m0 <- lm(y ~ poly(x, degree = d),
data = subset(da, j))
# Erro de treinamento.
ErrTra <- deviance(m0)/sum(j)
# Erro de validação.
yhat <- predict(m0, newdata = das[[f]])
ErrVal <- crossprod(yhat - das[[f]]$y)/length(yhat)
return(c(ErrTra = ErrTra, ErrVal = ErrVal))
m0 <- lm(y ~ poly(x, degree = d), data = subset(da, j))
# Predição com os dados de treino e de teste.
yhat_inter <- predict(m0, newdata = subset(da, j)) # Dados de treino.
yhat_exter <- predict(m0, newdata = das[[f]]) # Dados de teste.
# Erros de predição
Etra <- crossprod(subset(da, j)$y - yhat_inter)/length(yhat_inter)
Eval <- crossprod(das[[f]]$y - yhat_exter)/length(yhat_exter)
# Retorno.
c(ErrTra = Etra, ErrVal = Eval)
})
kfol <- cbind(cen, as.data.frame(t(kfol)))
# Curva de erro em cada fold.
xyplot(ErrTra + ErrVal ~ deg | fold,
xyplot(ErrTra + ErrVal ~ deg | factor(fold),
data = kfol,
xlab = "Grau do polinômio",
ylab = "Erro",
as.table = TRUE,
auto.key = list(corner = c(0.95, 0.95)),
type = "o")
# Os folds juntos para um mesmo tipo de erro.
xyplot(ErrTra + ErrVal ~ deg,
xlab = "Grau do polinômio",
ylab = "Erro",
groups = fold,
data = kfol,
type = "o") +
......@@ -205,8 +253,7 @@ xyplot(ErrTra + ErrVal ~ deg,
# k-fold ou leave-one-out (n-fold)?
# Cenários.
cen <- expand.grid(fold = 1:n,
deg = 1:15)
cen <- expand.grid(fold = 1:n, deg = 1:15)
# Obtendo os erros para cada cenário.
nfol <- mapply(f = cen$fold,
......@@ -239,6 +286,8 @@ xyplot(ErrTra + ErrVal ~ deg,
# Leave-one-out vs k-fold.
xyplot(ErrVal ~ deg,
data = nfol,
xlab = "Grau do polinômio",
ylab = "Erro",
type = c("p", "a")) +
as.layer(xyplot(ErrVal ~ (deg + 0.15),
data = kfol,
......@@ -348,8 +397,15 @@ resul
xyplot(eqm + Vy + Vh + Bh ~ deg,
data = resul,
auto.key = list(corner = c(0.95, 0.95))) +
glayer(panel.smoother(..., se = FALSE))
type = c("p", "o"),
xlab = "Grau do polinômio",
ylab = "Componentes do erro quadrático médio",
auto.key = list(corner = c(0.95, 0.95),
text = c(expression(EQM(hat(y))),
expression(Var(y)),
expression(Var(hat(y))),
expression(Bias(hat(y))))))# +
# glayer(panel.smoother(..., se = FALSE))
```
```{r, eval = FALSE, include = FALSE}
......@@ -395,7 +451,7 @@ da$y <- eta(da$x) + epsilon(da$x)
# Cria lista com fórmulas.
dmax <- 15
form <- replicate(dmax, y ~ x)
names(form) <- 1:dmax
names(form) <- sprintf("Pol. grau: %d", 1:dmax)
xyplot.list(form,
data = da,
......@@ -425,11 +481,11 @@ xyplot.list(form,
return(yfit)
})
# Ajuste médio ponto a ponto.
yfitm <- rowMeans(resul)
panel.lines(x = xseq,
y = yfitm,
col = "green2",
lwd = 3)
# yfitm <- rowMeans(resul)
# panel.lines(x = xseq,
# y = yfitm,
# col = "green2",
# lwd = 3)
# O verdadeiro modelo.
panel.curve(eta,
col = "orange",
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment