aula10
lsj
June 22, 2021
regressão:
1.1. otimização de um modelo por gradiente descendente
1.2. modelos e generalização
1.3. regularização: LASSO e ridge regression
1.4. regressão robusta
1.5. regressão com kernel
0. inicialização: bibliotecas para a aula de hoje
# para monitorar o tempo de processamento:
t00 = Sys.time()
# reprodutibilidade
set.seed(123)
# bibliotecas
library(glmnet)## Loading required package: Matrix## Loaded glmnet 3.0-2#require(doMC)
#library(carData)
library(olsrr)## This version of bslib is designed to work with shiny version 1.6.0 or higher.##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:olsrr':
##
## cementlibrary(quantreg)## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolverequire(robustbase)## Loading required package: robustbaselibrary(ggplot2)
library(kernlab)##
## Attaching package: 'kernlab'## The following object is masked from 'package:ggplot2':
##
## alpha1. Regressão
1.1 otimização de um modelo por gradiente descendente
Vamos ilustrar o algoritmo de otimização por gradiente descendente encontrando os parâmetros de um modelo linear simples: \(y=w_0+w_1 x + \epsilon\):
# https://www.r-bloggers.com/linear-regression-by-gradient-descent/
# simulação de um modelo linear: y=1+2x, com x entre -1 e 1 e ruÃdo gaussiano com sd = 0.2
n=100
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.2)
# vamos primeiro ajustar um modelo linear:
modelo_linear <- lm( y ~ x )
print(modelo_linear)##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 0.9892 1.9910# visualização:
plot(x,y, col=rgb(0.2,0.4,0.6,0.4), main='regressão linear por gradiente descendente')
abline(modelo_linear, col='blue')
# função custo: erro quadrático médio/2
custo <- function(X, y, w) {
sum( (X %*% w - y)^2 ) / (2*length(y))
}
# taxa de aprendizado e número de iterações
alpha <- 0.01
num_iters <- 1000
# guarda a história:
custo_historia <- double(num_iters)
w_historia <- list(num_iters)
# inicializa os parâmetros (w1,w2)=(0,0)
w <- matrix(c(0,0), nrow=2)
# adiciona a coluna de '1's para w1
X <- cbind(1, matrix(x))
# gradiente descendente:
for (i in 1:num_iters) {
error <- (X %*% w - y)
delta <- t(X) %*% error / length(y)
w <- w - alpha * delta
custo_historia[i] <- custo(X, y, w)
w_historia[[i]] <- w
}
# resultado
print(w)## [,1]
## [1,] 0.9888444
## [2,] 1.9114226# visualização do resultado:
par(mfrow=c(1,1))
plot(x,y, col=rgb(0.2,0.4,0.6,0.4), main='regressão linear por gradiente descendente')
for (i in c(1,3,6,10,14,seq(20,num_iters,by=10))) {
abline(coef=w_historia[[i]], col=rgb(0.8,0,0,0.3))
}
abline(coef=w, col='blue')
# visualização da convergência:
plot(custo_historia, type='line', col='blue', lwd=3, main='convergência da função custo', ylab='custo', xlab='Iterações')## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to first
## character
1.2 modelos e generalização
Para examinar como a complexidade afeta a generalização, vamos gerar um 'conjunto de treinamento' e um 'conjunto de teste' com o mesmo polinômio de segundo grau e vamos comparar o ajuste de polinômios de graus variados a esses dois conjuntos:
# reprodutibilidade
set.seed(123)
#modelo: y = 0.6-0.5*x+0.8*x^2
x = seq(0,1,0.01)
y = 0.6-0.5*x+0.8*x^2
#conjunto de treinamento: dados com erros de 0.02
ntreino = 10
xtreino = runif(ntreino,0,1)
ytreino = 0.6-0.5*xtreino+0.8*xtreino^2+rnorm(ntreino,0,sd = 0.02)
# conjunto de teste
ntest = 10
xtest = runif(ntest,0,1)
ytest = 0.6-0.5*xtest+0.8*xtest^2+rnorm(ntest,0,sd = 0.02)
plot(x,y,col='blue',type='l',ylim=c(0,1), main='vermelho: treino; verde: teste',lwd=2)
points(xtreino,ytreino,col='red',pch=20)
points(xtest,ytest,col='green',pch=20)
# vamos examinar modelos polinomiais de várias complexidades
# m: grau do polinômio
# custo: soma do quadrado dos resÃduos
custo = rep(0,ntreino-1)
custo_test = rep(0,ntreino-1)
grau = rep(0,ntreino-1)
for(i in 1:(ntreino-1)){
grau[i]=i
# ajuste de polinômio com lm()
yf=lm(ytreino ~ poly(xtreino, i, raw=TRUE))
model_treino = rep(0,ntreino)
for(j in 1:(i+1)){model_treino = model_treino + yf$coefficients[j]*xtreino^(j-1)}
custo[i] = sum((ytreino-model_treino)^2)
model_test = rep(0,ntreino)
for(j in 1:(i+1)){model_test = model_test + yf$coefficients[j]*xtest^(j-1)}
custo_test[i] = sum((ytest-model_test)^2)
}
custo## [1] 4.669985e-02 2.136856e-03 2.106459e-03 1.564326e-03 6.450436e-04
## [6] 3.614415e-04 3.128014e-04 1.378176e-04 2.080280e-20custo_test## [1] 0.038768122 0.005321278 0.005350283 0.005131396 0.013097903
## [6] 0.027322413 0.187435975 12.365744952 91.947179410plot(grau,custo,col='red',type='l',main='custo dos conjuntos de treinamento e teste')
points(grau,custo,col='red',pch=19)
lines(grau,custo_test,col='green')
points(grau,custo_test,col='green',pch=19)
legend("topright",legend=c("treino","teste"),pch=19,col=c("red","green"))
# polinômio de grau 9:
yf=lm(ytreino ~ poly(xtreino, 9, raw=TRUE))
xt = seq(0,1,0.001)
yt = rep(0,length(xt))
for(j in 1:10){yt = yt + yf$coefficients[j]*xt^(j-1)}
plot(x,y,col='blue',type='l',ylim=c(0,1), main='preto: polinômio de grau 9',lwd=2)
points(xtreino,ytreino,col='red',pch=20)
points(xtest,ytest,col='green',pch=20)
lines(xt,yt,lwd=3)
1.3 regularização: LASSO e ridge regression
Nesses métodos um termo de regularização torna os parâmetros mais estáveis
Embora os modelos ajustem os dados pior que mÃnimos quadrados comum (OLS, ordinary least squares), eles têm melhor capacidade de generalização
Vamos usar o pacote glmnet: ajusta um modelo linear com regularização
esta função tem 2 parâmetros principais: lambda é um parâmetro de regularização (controla a magnitude dos parâmetros da regressão)
alpha controla o método- alpha = 1 corresponde a LASSO, alpha = 0 a ridge regression e alpha entre 0 e 1 interpola entre esses dois métodos
# https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
# reprodutibilidade
set.seed(123)
# dados: vamos gerar dados com um modelo linear e introduzir alguns outliers
# modelo linear: y=1+2x, com x entre -1 e 1 e ruÃdo gaussiano com sd = 0.1
n=15
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.1)
# outliers: idem com sd = 1
x[n+1]=-0.6
y[n+1]=-2
x[n+2]=-0.4
y[n+2]=-1.5
# vamos começar plotando os dados e o modelo:
ymodelo = 1+2*x
plot(x,y,pch=20,main='vermelho: modelo, verde: lm com os outliers')
lines(x,ymodelo,type='l',col='red',lwd=2)
# modelo linear normal (OLS):
lm.mod = lm(y~x)
summary(lm.mod)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3846 -0.1214 0.1485 0.3178 0.5787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7952 0.1421 5.594 5.12e-05 ***
## x 2.3510 0.2496 9.419 1.09e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5824 on 15 degrees of freedom
## Multiple R-squared: 0.8554, Adjusted R-squared: 0.8457
## F-statistic: 88.71 on 1 and 15 DF, p-value: 1.092e-07abline(lm.mod,col='green',lwd=2)
# ridge regression:
# vamos considerar inicialmente vários valores para o parâmetro de regularização lambda
lambda <- 10^seq(10, -2, length = 100)
#ridge regression: alpha = 0
# como x só tem 1 dimensão e o glmnet precisa que x tenha ao menos duas, usamos o truque abaixo:
xx=cbind(x,0)
ridge.mod <- glmnet(xx, y, alpha = 0)
# coeficientes do ajuste
predict(ridge.mod, s = 0, type = 'coefficients')[1:2,]## (Intercept) x
## 0.8079507 2.1519292# plot
plot(x,y,pch=20)
lines(x,ymodelo,type='l',col='red',lwd=2)
abline(lm.mod,col='green',lwd=2)
a=predict(ridge.mod, s = 0, type = 'coefficients')[1]
b=predict(ridge.mod, s = 0, type = 'coefficients')[2]
lines(x,a+b*x,col='blue',lwd=2)
# valor ótimo de lambda para ridge por validação cruzada
lambdas <- 10^seq(3, -2, by = -.1)
cv_fit <- cv.glmnet(xx, y, alpha = 0, lambda = lambdas, nfolds = 10)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold# valor ótimo:
opt_lambda <- cv_fit$lambda.min
opt_lambda## [1] 0.01584893# modelos ajustados na CV
fit <- cv_fit$glmnet.fit
# coeficientes do ajuste
predict(fit, s = 0, type = 'coefficients')[1:2,]## (Intercept) x
## 0.7962139 2.3347219# agora LASSO alpha = 1
lasso.mod <- glmnet(xx, y, alpha = 1)
# coeficientes do ajuste
predict(lasso.mod, s = 0, type = 'coefficients')[1:2,]## (Intercept) x
## 0.7961648 2.3354854# plot
plot(x,y,pch=20)
lines(x,ymodelo,type='l',col='red')
abline(lm.mod,col='green')
lines(x,a+b*x,col='blue')
a2=predict(lasso.mod, s = 0, type = 'coefficients')[1]
b2=predict(lasso.mod, s = 0, type = 'coefficients')[2]
lines(x,a2+b2*x,col='black')
legend("topright",legend=c("modelo","lm",'ridge','lasso'),pch=19,col=c("red","green",'blue','black'))
1.4 regressão robusta
O R tem várias bibliotecas com métodos robustos, como a MASS, quantreg e robustbase
Os métodos de regressão robusta usam funções de custo menos sensÃvel a outliers
exemplo: função de Huber \[ L(a) = L_\delta(a) = \left\{ \begin{array}{lr} \frac{1}{2}a^2 & {\rm se}~ |a| \le \delta\\ \delta(|a|- \frac{1}{2}\delta) & {\rm se}~ |a| > \delta \end{array} \right. \] onde \(a = y_i -\mathbf x_i.\mathbf w\) é o resÃduo de uma medida
esta função é quadrática para valores pequenos de \(a\) e linear para valores grandes
Vamos dar uma examinada nos vários métodos usando os dados (x,y) do exemplo acima
# https://rpubs.com/dvallslanaquera/robust_regression
# reprodutibilidade
set.seed(123)
# dados:
n=100
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.2)
# outliers:
nout = 20
for(i in 1:nout){
x[n+i] = runif(1, -0.8, 0)
y[n+i] = 2.5 + rnorm(1,mean=2,sd=0.5)
}
# modelo linear
fitLS <- lm(y~x)
summary(fitLS)$r.squared## [1] 0.1628212p = predict(fitLS, newx = x)
sd_LS = sd(p-y)
mad_LS = mad(p-y)
# Outliers usando a figura de Cook
# a distância de Cook mede o efeito de se remover um certo ponto da análise e, assim, é útil para identificar outliers
ols_plot_cooksd_bar(fitLS)
# Huber loss
fitH <- rlm(y~x, k2 = 1.345)
p = predict(fitH, newx = x)
sd_H = sd(p-y)
mad_H = mad(p-y)
# LQS - faz um ajuste aos pontos "bons", sem os outliers
fitLQS <- lqs(y~x, method = "lms")
p = predict(fitLQS, newx = x)
sd_LQS = sd(p-y)
mad_LQS = mad(p-y)
# LTS: least trimmed squares - minimiza a soma dos quadrados dos resÃduos para uma subamostra dos dados
fitLTS <- lqs(y~x, method = "lts")
p = predict(fitLTS, newx = x)
sd_LTS = sd(p-y)
mad_LTS = mad(p-y)
# LAD: least absolute devition
fitLAD <- rq(y~x, tau = 0.5)
p = predict(fitLAD, newx = x)
sd_LAD = sd(p-y)
mad_LAD = mad(p-y)
# S-estimator
fitS <- lqs(y~x, method = "S")
p = predict(fitS, newx = x)
sd_S = sd(p-y)
mad_S = mad(p-y)
# MM-estimator: combinação dos estimadores M e S
fitMM <- rlm(y~x, method = "MM")
p = predict(fitMM, newx = x)
sd_MM = sd(p-y)
mad_MM = mad(p-y)
# desvio padrão e mÃnimo desvio absoluto dos vários estimadores
c(sd_LS,sd_H,sd_LQS,sd_LTS,sd_LAD,sd_S,sd_MM)## [1] 1.524960 1.567875 1.589155 1.582387 1.568842 1.579226 1.579755c(mad_LS,mad_H,mad_LQS,mad_LTS,mad_LAD,mad_S,mad_MM)## [1] 0.6686189 0.2624243 0.2457071 0.2476601 0.2615887 0.2406093 0.2419468# visualização
plot(x,y,
xlab = "x", ylab = "y", type = "p",
pch = 20, cex = .8)
abline(fitLS, col = 1)
abline(fitH, col = 2)
abline(fitLAD, col = 3)
abline(fitLTS, col = 4)
abline(fitLQS, col = 5)
abline(fitS, col = 6)
abline(fitMM, col = 7)
legend('topright', c("LS", "Huber", "LAD","LTS","LQS",
"S-estimador","MM-estimador" ),
lty = rep(1, 7), bty = "n",
col = c(1, 2, 3, 4, 5, 6, 7))
1.5 regressão com kernel
este tipo de regressão é muito usado em interpolação
# Nonparametric Regression and Cross-Validation, Yen-Chi Chen
# reprodutibilidade
set.seed(123)
# simulação dos dados
X = runif(200, min=0, max=4*pi)
Y = sin(X) + rnorm(200, sd=0.3)
plot(X,Y, pch=20)
# vamos usar a função ksmooth para esta regressão
# dois parâmetros principais: o tipo de kernel e a largura de banda (h)
# kernel gaussiano
Kreg = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = 0.9)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="purple")
# kernel box
Kreg = ksmooth(x=X,y=Y,kernel = "box",bandwidth = 0.9)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="orange")
# efeito de variar h
Kreg1 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =0.1)
Kreg2 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =0.9)
Kreg3 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =3.0)
plot(X,Y,pch=20)
lines(Kreg1, lwd=4, col="orange")
lines(Kreg2, lwd=4, col="purple")
lines(Kreg3, lwd=4, col="limegreen")
legend("topright", c("h=0.1","h=0.9","h=3.0"), lwd=6,col=c("orange","purple","limegreen"))
# escolha de h por validação cruzada
# leave-one-out-cross-validation (LOOCV)
n = length(X)
# n: sample size
h_seq = seq(from=0.2,to=2.0, by=0.1)
# smoothing bandwidths we are using
CV_err_h = rep(NA,length(h_seq))
for(j in 1:length(h_seq)){
h_using = h_seq[j]
CV_err = rep(NA, n)
for(i in 1:n){
X_val = X[i]
Y_val = Y[i]
# validation set
X_tr = X[-i]
Y_tr = Y[-i]
# training set
Y_val_predict = ksmooth(x=X_tr,y=Y_tr,kernel = "normal",bandwidth=h_using,x.points = X_val)
CV_err[i] = (Y_val - Y_val_predict$y)^2
# we measure the error in terms of difference square
}
CV_err_h[j] = mean(CV_err)
}
# erro medio da validação cruzada para os vários valores de h
CV_err_h## [1] 0.10183202 0.09692202 0.09471400 0.09323028 0.09230145 0.09174197
## [7] 0.09144213 0.09141794 0.09171727 0.09240507 0.09351454 0.09506037
## [13] 0.09704216 0.09944046 0.10225141 0.10545900 0.10905926 0.11305207
## [19] 0.11744657plot(x=h_seq, y=CV_err_h, type="b", lwd=3, col="blue",xlab="bandwidth h", ylab="LOOCV prediction error")
which(CV_err_h == min(CV_err_h))## [1] 8h_seq[which(CV_err_h == min(CV_err_h))]## [1] 0.9#exemplo com 5-CV:
n = length(X)
N_cv = 100
k = 5
cv_lab = sample(n,n,replace=F) %% k
## randomly split all the indices into k numbers
h_seq = seq(from=0.2,to=2.0, by=0.1)
CV_err_h = rep(0,length(h_seq))
for(i_tmp in 1:N_cv){
CV_err_h_tmp = rep(0, length(h_seq))
cv_lab = sample(n,n,replace=F) %% k
for(i in 1:length(h_seq)){
h0 = h_seq[i]
CV_err =0
for(i_cv in 1:k){
w_val = which(cv_lab==(i_cv-1))
X_tr = X[-w_val]
Y_tr = Y[-w_val]
X_val = X[w_val]
Y_val = Y[w_val]
kernel_reg = ksmooth(x = X_tr,y=Y_tr,kernel = "normal",bandwidth=h0,
x.points=X_val)
# WARNING! The ksmooth() function will order the x.points from
# the smallest to the largest!
CV_err = CV_err+mean((Y_val[order(X_val)]-kernel_reg$y)^2,na.rm=T)
# na.rm = T: remove the case of 'NA'
}
CV_err_h_tmp[i] = CV_err/k
}
CV_err_h = CV_err_h+CV_err_h_tmp
}
CV_err_h = CV_err_h/N_cv
plot(h_seq,CV_err_h, type="b", lwd=4, col="blue", xlab="Smoothing Bandwidth",
ylab="5-CV Error")
h_opt = h_seq[which(CV_err_h==min(CV_err_h))]
h_opt## [1] 0.9# tempo de processamento deste script
Sys.time()-t00## Time difference of 6.263521 secs