############################################################################ # # # REGRESSÃO NÃO PARAMÉTRICA - EXEMPLO FARAWAY # # # # Cybelli Barbosa e Thiago Brandão # # # ############################################################################ require (faraway) ## dados do exemplo A (simulado) # função: f(x)=sin³(2*pi*x³) data(exa) plot (y ~ x, exa, main="Exemplo A", pch=".") lines (m ~ x, exa) ## dados do exemplo B (simulado): # função: f(x)=0 data(exb) plot (y ~ x, exb, main="Exemplo B", pch=".") lines (m ~ x, exb) ## Dados reais Old Faithful: # O tempo de espera entre as erupções e a duração da erupção do gêiser Old Faithful # em Yellowstone National Park, Wyoming, EUA. # eruptions: duração da erupção em minutos. # waiting: tempo de espera até a proxima erupção, em minutos. data(faithful) plot (waiting ~ eruptions, faithful, main="Old Faithful", pch=".") #### MÉTODO DO NUCLEO: ### Estimador Kernel: ## lambda = 0.1: plot(waiting ~ eruptions, faithful, main="bandwidth=0.1", pch=".") lines(ksmooth(faithful$eruptions, faithful$waiting, "normal", 0.1)) ## lambda = 0.5: plot(waiting ~ eruptions, faithful, main="bandwidth=0.5", pch=".") lines(ksmooth(faithful$eruptions, faithful$waiting, "normal", 0.5)) ## lambda = 2: plot(waiting ~ eruptions, faithful, main="bandwidth=2", pch=".") lines(ksmooth(faithful$eruptions, faithful$waiting, "normal", 2)) # A bandwidth=0.5 é a melhor escolha entre as três opções. #### Validação cruzada: ## Old Faithful library (sm) hm <- hcv(faithful$eruptions, faithful$waiting, display="lines") # Critério de validação cruzada como uma função da suavidade, o mínimo ocorre em 0.424. sm.regression(faithful$eruptions, faithful$waiting, h=hm, xlab="eruptions", ylab="waiting") # Estimador kernel com o valor estimado de suavização. ## Ex A hm <- hcv(exa$x, exa$y, display="lines") sm.regression(exa$x, exa$y, h=hm, xlab="x", ylab="y") ## Ex B hm <- hcv(exb$x, exb$y, display="lines") sm.regression(exb$x, exb$y, h=0.005, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.0141, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.018, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.022, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.026, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.282, xlab="x", ylab="y") sm.regression(exb$x, exb$y, h=0.5, xlab="x", ylab="y") ## alterando o parametro h no gráfico: library (manipulate) manipulate({ sm.regression(exb$x, exb$y, h=h.range, xlab="x", ylab="y") }, h.range=slider(0.005,0.5) ) #### MÉTODO DAS SPLINES: ### Suavização ## Old Faithful plot (waiting ~ eruptions, faithful, pch=".") lines(smooth.spline(faithful$eruptions, faithful$waiting)) ## Exemplo A plot (y ~ x, exa, pch=".") lines(exa$x, exa$m) lines(smooth.spline(exa$x, exa$y), lty=2) ## Exemplo B plot (y ~ x, exb, pch=".") lines(exb$x, exb$m) lines(smooth.spline(exb$x, exb$y), lty=2) ### Regressão ## (aplicado ao exemplo A) rhs <- function (x,c) ifelse (x > c, x - c, 0) curve (rhs(x,0.5), 0, 1) knots <- 0:9/10; knots # nós ou janelas igualmente espaçados dm <- outer (exa$x, knots, rhs) matplot (exa$x, dm, type="l", col=1) ## ajuste de regressão g <- lm(exa$y ~ dm) plot (y ~ x, exa, pch=".", xlab="x", ylab="y") lines (exa$x, predict(g)) ## adensamento dos nós newknots <- c(0, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95) dmn <- outer(exa$x, newknots, rhs) gn <- lm(exa$y ~ dmn) plot (y ~ x, exa, pch=".", xlab="x", ylab="y") lines (exa$x, predict(gn)) ## Fino ajuste: library (splines) matplot(bs(seq(0,1,length=1000), df=12), type="l", ylab=" ", col=1) sm1 <- lm(y ~ bs(x,12), exa) plot (y ~ x, exa, pch=".") lines(m ~ x, exa) lines(predict(sm1) ~ x, exa, lty=2) # poderia ser melhorado incluindo nós na região de maior curvatura e menos na região mais plana #### MÉTODO DOS POLINOMIOS ## Old Faithful plot (waiting ~ eruptions, faithful, pch=".") f <- loess(waiting ~ eruptions, faithful) i <- order (faithful ~ eruptions) lines(f$x[i], f$fitted[i]) ## Exemplo A plot(y ~ x, exa, pch=".") lines(exa$x, exa$m, lty=1) f <- loess(y ~ x, exa) lines(f$x, f$fitted, lty=2) f <- loess(y ~ x, exa, span=0.22) lines(f$x, f$fitted, lty=5) ## Exemplo B plot (y ~ x, exb, pch=".") f <- loess(y ~ x, exb) lines(f$x, f$fitted, lty=2) f <- loess(y ~ x, exb, span=1) lines(f$x, f$fitted, lty=5) lines(exb$x, exb$m) #---------------------------------------------------------------------