Commit 06eebccd authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Adiciona tutorial de métodos de gradiente descendente.

parent 49382ee8
......@@ -25,6 +25,8 @@ navbar:
menu:
- text: "Validação cruzada"
href: tutorials/01-cross-validation.html
- text: "Gradiente descendente"
href: tutorials/02-gradient-methods.html
- text: "Scripts"
icon: fa-file-text
href: scripts/
......
......@@ -36,7 +36,7 @@ eta <- function(x, beta = c(10, 0.05, 2)) {
beta[1] + beta[2] * x + beta[3] * sin(x/10)
}
# Retorna um rúido normal padrão.
# Retorna um rdo normal padrão.
epsilon <- function(x) {
rnorm(n = length(x), mean = 0, sd = 1)
}
......
---
title: "Métodos de Gradiente Descendente"
author: Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira
date: 2017-09-05
#bibliography: ../config/Refs.bib
#csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
```{r, include = FALSE}
source("../config/setup.R")
opts_chunk$set(
cache = FALSE,
message = FALSE,
warning = FALSE)
```
# Gradiente Descendente
Para materializar o algoritmo de gradiente descendente (GD), vamos
considerar a estimação de parâmetros em modelos de regressão não linear.
Primeiramente, modelos não lineares são estimados considerando
algoritmos baseados em Newton-Raphson que fazem uso da primeira e
segunda derivadas do modelo em relação aos parâmetros. O GD só
utiliza-se da primeira derivada.
## Simulação do modelo
```{r}
# Carregando pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)
```
```{r}
#-----------------------------------------------------------------------
# Simulando dados do modelo Michaelis Menten.
# Retorna os valores determinados pelo modelo para a média de Y|x.
eta <- function(x, theta) {
theta[1] * x/(theta[2] + x)
}
# Valores para os parâmetros.
pars <- c(tha = 10,
thv = 2)
# Gráfico da função média.
curve(eta(x, pars),
from = 0,
to = 20)
# Retorna um rúido normal padrão.
epsilon <- function(x, ...) {
rnorm(n = length(x), mean = 0, ...)
}
# Simulando do modelo.
set.seed(321)
da <- data.frame(x = sort(runif(n = 50, min = 3, max = 20)))
da$y <- eta(x = da$x, theta = pars) + epsilon(x = da$x, sd = 0.1)
plot(y ~ x, data = da)
curve(eta(x, theta = pars),
add = TRUE,
col = "orange")
```
## Inspeção da função objetivo
Vamos avaliar a função objetivo (função perda ou custo) para um grid de
valores dos parâmetros.
```{r}
# Função perda (custo ou objetivo).
loss <- function(theta, y, x) {
sum((y - eta(x, theta = theta))^2)
}
loss(theta = pars, y = da$y, x = da$x)
# Grid de valores dos parâmetros.
grid <- expand.grid(tha = seq(0, 20, length.out = 41),
thv = seq(-3, 8, length.out = 41))
head(grid)
# Avalia a função objetivo em cada ponto do grid.
grid$rss <- mapply(tha = grid$tha,
thv = grid$thv,
FUN = function(tha, thv) {
loss(theta = c(tha, thv),
y = da$y,
x = da$x)
})
# Define escala de cores para visualização.
colr <- colorRampPalette(rev(brewer.pal(11, "Spectral")),
space = "rgb")
# Exibe a superfície da função objetivo (com log para mais detalhes).
levelplot(log(rss) ~ tha + thv,
data = grid,
col.regions = colr,
contour = TRUE) +
layer(panel.abline(v = pars["tha"],
h = pars["thv"],
lty = 2))
wireframe(log(rss) ~ tha + thv,
data = grid,
drape = TRUE,
col.regions = colr(101),
panel.3d.wireframe = wzRfun::panel.3d.contour,
type = c("bottom"))
```
## O algoritmo
A função objetivo (custo) para o problema em mãos é
$$
L(\theta) = \sum_{i = 1}^{n} (y_{i} - f(x_{i}, \theta))^2,
$$
em que $f$ é o modelo considerado para a média de $Y|x$, ou seja,
$$
f(x_{i}, \theta) = \frac{\theta_a x_i }{\theta_v + x_i},
\quad \theta = (\theta_a, \theta_v)'.
$$
A derivada primeira da função $L$ com relação a cada um dos seus
argumentos fornece o vetor gradiente. Você pode usar o Wolfran para
ganhar tempo no computo dessas derivadas. Por exemplo, passe a expressão
para o motor do Wolfran
["derivative of (y - a * x/(v + x))^2 wrt a"](https://www.wolframalpha.com/input/?i=derivative+of+(y+-+a*x%2F(v%2Bx))%5E2+wrt+a).
Você também pode usar o R para fazer essas derivadas de forma analítica
ou numérica. Outra opção é usar o
[wxMaxima](http://andrejv.github.io/wxmaxima/).
```{r}
#-----------------------------------------------------------------------
# Gradiente.
# Derivadas analíticas com o R.
loss_expr <- expression((y - tha * x/(thv + x))^2)
D(expr = loss_expr, name = "tha")
D(expr = loss_expr, name = "thv")
# Função gradiente.
grad <- function(theta, y, x) {
tha <- theta[1]
thv <- theta[2]
part <-
cbind(partial_tha =
-(2 * (x/(thv + x) * (y - tha * x/(thv + x)))),
partial_thv =
2 * (tha * x/(thv + x)^2 * (y - tha * x/(thv + x))))
# Derivada da soma = soma das derivadas.
colSums(part)
}
# Avaliando a função objetivo.
loss(theta = c(10, 2), y = da$y, x = da$x)
# Avaliando a função gradiente.
grad(theta = c(10, 2),
y = da$y,
x = da$x)
#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente.
# Tolerância: parar quando a mudança for irrelevante ou exceder número
# máximo de iterações.
tol <- 1E-5
niter <- 10000
# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.001
# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)
# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(tha = 18, thv = -1)
repeat {
delta <- alpha * grad(theta = th[k, ],
y = da$y,
x = da$x)
th[k + 1, ] <- th[k, ] - delta
# print(th[k, ])
if (k > niter | all(abs(delta) < tol)) break
k <- k + 1
}
th <- th[1:k, ]
#-----------------------------------------------------------------------
# Traço do algoritmo.
levelplot(rss ~ tha + thv,
data = grid,
contour = TRUE) +
layer(panel.abline(v = pars["tha"],
h = pars["thv"],
lty = 2)) +
layer(panel.points(x = th[, 1],
y = th[, 2],
type = "o",
cex = 0.1))
# Número de iterações.
k
# Estimativas.
th[k, ]
# A descida.
lth <- apply(th,
MARGIN = 1,
FUN = loss,
y = da$y,
x = da$x)
plot(log(lth) ~ seq_along(lth), type = "o", cex = 0.3)
#-----------------------------------------------------------------------
# Faça testes.
# Repita as instruções acima com os valores abaixo e veja o algorítmo
# divergir.
th[k, ] <- c(tha = 10, thv = 3)
alpha <- 0.1
#-----------------------------------------------------------------------
# Ajustando o modelo usando um Gauss-Newton.
n0 <- nls(y ~ tha * x/(thv + x),
data = da,
start = c(tha = 18, thv = -1),
trace = TRUE)
summary(n0)
```
# Gradiente Descendente Estocástico
Em termos simples, o gradiente descendente estocástico difere apenas no
fato de não usar todos os dados mas sim usar uma fração dos dados. Esse
subconjunto dos dados é obtido com reamostragem dos dados originais. A
reamostragem é o componente estocástico do algorítmo. Com isso, o traço
da função objetivo deixa de ser suave, o que impossibilita o uso de
critérios de parada baseados em diferenças de valores consecutivos.
Então o mais comum é executar o algoritmo até exceder o número máximo de
iterações. O número máximo de iterações deve ser escolhido de forma a
garantir suficiente proximidade com o ponto de ótimo da função.
```{r}
#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente Estocástico.
# Para quando exceder número máximo de iterações.
niter <- 30000
# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.01
# Número de elementos amostrados do conjunto de treinamento.
index <- 1:nrow(da)
nele <- 3
# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)
# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(tha = 18, thv = -1)
while (k < niter) {
db <- da[sample(index, size = nele), ]
delta <- alpha * grad(theta = th[k, ],
y = db$y,
x = db$x)
th[k + 1, ] <- th[k, ] - delta
# print(th[k, ])
k <- k + 1
}
#-----------------------------------------------------------------------
# O traço do algorítmo.
th <- th[1:k, ]
levelplot(rss ~ tha + thv,
data = grid,
contour = TRUE) +
layer(panel.abline(v = pars["tha"],
h = pars["thv"],
lty = 2)) +
layer(panel.points(x = th[, 1],
y = th[, 2],
type = "l"))
# Estimativas.
th[k, ]
# A descida.
lth <- apply(th,
MARGIN = 1,
FUN = loss,
y = da$y,
x = da$x)
plot(log(lth), type = "o", cex = 0.3)
# Vendo apenas as últimas iterações.
plot((tail(lth, n = 500)), type = "o", cex = 0.3)
```
# Gradiente Descendente Boosting
```{r}
#-----------------------------------------------------------------------
# Gerando dados.
curve(10 + 0.05 * x + 2 * sin(x/10), 0, 100)
# Retorna os valores determinados pelo modelo.
eta <- function(x, beta = c(10, 0.05, 2)) {
beta[1] + beta[2] * x + beta[3] * sin(x/10)
}
# Simulando do modelo.
set.seed(567)
da <- data.frame(x = sort(runif(200, 0, 100)))
da$y <- eta(da$x) + epsilon(da$x, sd = 1)
# Exibindo os dados.
plot(y ~ x, data = da)
curve(eta, add = TRUE, col = "orange", lwd = 3)
```
## Exemplo com função "patamar"
A função patamar é uma função por partes definidas sobre o domínio de
$x$ dividido em $p$ subconjuntos justapostos em que para cada parte a
função é de valor constante, ou seja
$$
f(x) = \begin{cases}
y_1, & x_0 < x \leq x_1 \\
y_2, & x_1 < x \leq x_2 \\
\vdots & \vdots \\
y_p, & x_{p-1} < x \leq x_p.
\end{cases}
$$
Vamos considerar a versão com $p = 2$. Dessa forma, a função retornará o
valor $y_1$ para valores $x < \beta$ e $y_2$ para valores de $x \geq
\beta$, em que $\beta$ é o ponto de corte do domínio.
```{r}
#-----------------------------------------------------------------------
# A função h(x) será uma função "patamar".
# Modelo inicial com predito sendo a média das observações.
m0 <- lm(y ~ 1, data = da)
plot(y ~ x, data = da)
lines(da$x, fitted(m0), type = "o", col = 2)
# Função patamar. Retorna ajustado e resíduo.
h <- function(beta, y, x) {
dico <- ifelse(x < beta, 0, 1)
fit <- ave(y, dico, FUN = mean)
list(dico = dico, fit = fit, res = y - fit)
}
# Função custo.
loss <- function(beta, y, x) {
fit <- h(beta, y, x)$fit
sum((y - fit)^2)
}
# Verifica o comportamento da função.
# bseq <- seq(0, 100, 1)
# lseq <- sapply(bseq, FUN = loss, y = da$y, x = da$x)
# plot(lseq ~ bseq, type = "o")
# Usando a optim().
# optim(par = c(50),
# fn = loss,
# y = da$y,
# x = da$x)$par
#-----------------------------------------------------------------------
# O algorítmo.
# Cria uma lista para armazenar os artefatos produzidos.
n <- nrow(da)
k <- 12
resul <- replicate(k,
data.frame(fit = numeric(n),
res = numeric(n)),
simplify = FALSE)
# Inicia a lista com resultados do primeiro ajuste.
resul[[1]]$fit <- fitted(m0)
resul[[1]]$res <- residuals(m0)
bval <- c(NA, numeric(k - 1))
# Executa o gradiente boosting.
for (i in 2:k) {
# Encontra o beta minimizando a função perda.
beta <- optimize(f = loss,
y = resul[[i - 1]]$res,
x = da$x,
interval = c(0, 100))$minimum
# Exibe o progresso.
cat(i, "\t", beta, "\n")
bval[i] <- beta
# Extrai valores ajustados e resíduos.
hfit <- h(beta,
y = resul[[i - 1]]$res,
x = da$x)
# Atualiza os valores ajustados e resíduos.
resul[[i]]$fit <- resul[[i - 1]]$fit + hfit$fit
resul[[i]]$res <- da$y - resul[[i]]$fit
}
#-----------------------------------------------------------------------
# Exibe os resultados.
# Lista replicando a mesma fórmula.
form <- replicate(k, y ~ x)
names(form) <- 1:k
xyplot.list(form,
data = da,
type = "n",
as.table = TRUE) +
layer({
i <- which.packet()
fit <- resul[[i]]$fit
panel.points(x = x,
y = y,
col = ifelse(y > fit, "blue", "red"))
panel.lines(da$x,
fit,
col = "green2",
lwd = 3)
panel.abline(v = bval[i], lty = 2)
})
# Converte lista em tabela.
names(resul) <- as.character(1:k)
result <- plyr::ldply(resul, .id = "k")
result$k <- factor(result$k, levels = 1:k)
result$x <- da$x
# O comportamento dos resíduos com o acréscimo de funções.
xyplot(res ~ x | k,
data = result,
type = c("n", "smooth"),
span = 0.3,
as.table = TRUE) +
layer(panel.points(x = x, y = y,
col = ifelse(y > 0, "blue", "red"))) +
layer(panel.abline(h = 0, lty = 2))
```
## Exemplo com funções segmentadas lineares
```{r}
#-----------------------------------------------------------------------
# Ajustar modelos segmentados.
library(segmented)
# Função que ajusta o modelo segmentado com um ponto de quebra. Já está
# com a função de perda quadrática embutida na segmented().
h <- function(y, x) {
m <- lm(y ~ x)
s <- segmented(m, seg.Z = ~x)
fit <- fitted(s)
brk <- summary.segmented(s)$psi[1, 2]
list(brk = brk, fit = fit, res = y - fit)
}
# Verifica o que a função retorna.
str(h(y = da$y, x = da$x))
# Cria uma lista para armazenar os artefatos produzidos.
n <- nrow(da)
k <- 6
resul <- replicate(k,
data.frame(fit = numeric(n),
res = numeric(n)),
simplify = FALSE)
# Inicia a lista com resultados do primeiro ajuste.
resul[[1]]$fit <- fitted(m0)
resul[[1]]$res <- residuals(m0)
bval <- c(NA, numeric(k - 1))
# Executa o gradiente boosting.
for (i in 2:k) {
# Extrai valores ajustados e resíduos.
hfit <- h(y = resul[[i - 1]]$res,
x = da$x)
# Extrai o ponto de quebra.
bval[i] <- hfit$brk
# Atualiza os valores ajustados e resíduos.
resul[[i]]$fit <- resul[[i - 1]]$fit + hfit$fit
resul[[i]]$res <- da$y - resul[[i]]$fit
}
#-----------------------------------------------------------------------
# Exibe os resultados.
# Lista replicando a mesma fórmula.
form <- replicate(k, y ~ x)
names(form) <- 1:k
xyplot.list(form,
data = da,
type = "n",
as.table = TRUE) +
layer({
i <- which.packet()
fit <- resul[[i]]$fit
panel.points(x = x,
y = y,
col = ifelse(y > fit, "blue", "red"))
panel.lines(da$x,
fit,
col = "green2",
lwd = 3)
panel.abline(v = bval[i], lty = 2)
})
# Converte lista em tabela.
names(resul) <- as.character(1:k)
result <- plyr::ldply(resul, .id = "k")
result$k <- factor(result$k, levels = 1:k)
result$x <- da$x
# O comportamento dos resíduos com o acréscimo de funções.
xyplot(res ~ x | k,
data = result,
type = c("n", "smooth"),
span = 0.3,
as.table = TRUE) +
layer(panel.points(x = x, y = y,
col = ifelse(y > 0, "blue", "red"))) +
layer(panel.abline(h = 0, lty = 2))
```
# Cuidados com o tunning
```{r}
#-----------------------------------------------------------------------
# Debugando com problema em uma dimensão.
m0 <- lm(dist ~ speed, data = cars)
coef(m0)
loss <- function(beta, y, xx) {
mean((y - (coef(m0)[1] + beta * xx))^2)
}
loss(y = cars$dist, xx = cars$speed, beta = coef(m0)[2])
deviance(m0)
D(expression((y - (beta0 + beta * x))^2),
name = "beta")
mygrad <- function(y, xx, beta) {
beta0 <- coef(m0)[1]
mean(-(2 * (xx * (y - (beta0 + beta * xx)))))
}
library(numDeriv)
numDeriv::grad(func = loss,
x = coef(m0)[2] + 1,
y = cars$dist,
xx = cars$speed)
mygrad(y = cars$dist,
xx = cars$speed,
beta = coef(m0)[2] + 1)
# bseq <- seq(0, 8, length.out = 50)
bseq <- seq(-10, 20, length.out = 50)
lseq <- sapply(bseq, FUN = loss, y = cars$dist, x = cars$speed)
plot(lseq ~ bseq, type = "l")
# alpha de 0.001 OK
# alpha de 0.003 OK ZIG-ZAG
# alpha de 0.005 OK
alpha <- 0.0038
n <- 100
bt <- numeric(n)
k <- 1
bt[k] <- -5
while (k < n) {
gr1 <- numDeriv::grad(func = loss,
x = bt[k],
y = cars$dist,
xx = cars$speed)
gr2 <- mygrad(y = cars$dist,
xx = cars$speed,
beta = bt[k])
l <- loss(y = cars$dist, x = cars$speed, beta = bt[k])
cat(gr1, "\t", gr2, "\t", bt[k], "\t", l, "\n")
bt[k + 1] <- bt[k] -
alpha * gr2
k <- k + 1
}
plot(bseq, lseq, type = "l")
lt <- sapply(bt, FUN = loss, y = cars$dist, x = cars$speed)
lines(lt ~ bt, type = "o", col = 2)
abline(v = coef(m0)[2])
head(bt)
#-----------------------------------------------------------------------
```
# Links úteis
* <http://www.ccs.neu.edu/home/vip/teach/MLcourse/4_boosting/slides/gradient_boosting.pdf>;
* <https://cran.r-project.org/web/packages/bst/>;
* <https://cran.r-project.org/web/packages/mboost/>;
* <https://cran.r-project.org/web/packages/xgboost/>;