ESTATÍSTICA NÃO PARAMÉTRICA
PRÁTICA 9: ESTIMAÇÃO DE DENSIDADES E MÉTODOS DE SUAVIZAÇÃO
\[\\[0.05in]\]1 PACOTES NECESSÁRIOS
Para esta nona aula prática de Estatística Não Paramétrica, precisaremos dos seguintes pacotes:
require(ggplot2) # Gráficos
require(dplyr) # Manipulação de dados
require(MASS) # Para a base de dados mcycle2 MOTIVAÇÕES PRÁTICAS
Inicialmente, vamos ver quais são as motivações práticas desta aula e uma rápida análise descritiva das variáveis pertinentes para as análises de hoje.
2.1 MODELOS DE CARROS
A base de dados mtcars é um
conjunto de dados clássico incluído no R,
frequentemente utilizado para exemplos de análise estatística e
aprendizado de máquina. Ele contém informações sobre 32
modelos de carros extraídas da edição de 1974 da revista
Motor Trend. Os dados foram projetados para explorar as
relações entre o desempenho de veículos,
consumo de combustível e outras características
relacionadas. As variáveis desta base de dados
são:
mpg: Milhas por galão de combustível.cyl: Número de cilindros.disp: Deslocamento do motor (em polegadas cúbicas).hp: Potência bruta (em cavalos de força).drat: Razão do diferencial traseiro.wt: Peso do veículo (em milhares de libras).qsec: Tempo de 1/4 de milha (em segundos).vs: Tipo de motor (0 = motor em V, 1 = motor em linha).am: Tipo de transmissão (0 = automática, 1 = manual).gear: Número de marchas.carb: Número de carburadores.
O comando a seguir traz a base de dados:
mtcars## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
## Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
## Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
## Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
## Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
## Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
## Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
## Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
## Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
## Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
## Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
## AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
## Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2Para a aula de hoje, precisaremos das variáveis
mpg e ph.
2.1.1 KERNEL PARA DENSIDADES
Para ilustrar a estimação de densidades a partir do kernel,
vamos utilizar os dados da variável mpg.
Inicialmente, vamos visualizar o gráfico histograma destes dados
considerando a regra de Sturges para o cálculo do
número de intervalos ideal k:
\[ k = \lceil 1 + \log_2(n) \rceil \]
em que \(n\) é o tamanho da amostra.
O histograma é então feito com os seguintes comandos:
# Regra de Sturges:
n = length(mtcars$mpg)
k = ceiling(1 + log2(n))
# Histograma.
ggplot(mtcars, aes(x = mpg)) +
geom_histogram(aes(y = ..density..), bins = k, fill = "lightblue", color = "black") +
labs(x = "Milhas por Galão (mpg)", y = "Densidade") +
theme_minimal()
2.1.2 MÉTODOS DE SUAVIZAÇÃO
Para os métodos de suavização, estaremos
interessados em suavizar a relação entre as variáveis
hp (variável x) e mpg (variável y).
Inicialmente, vamos visualizar a associação destas duas
variáveis em um gráfico de dispersão:
# Gráfico de dispersão
ggplot(mtcars, aes(x = hp, y = mpg)) +
geom_point(color = "blue", size = 3) + # Pontos de dispersão
labs(
x = "Potência do Motor (hp)",
y = "Milhas por Galão (mpg)"
) +
theme_minimal()
Utilizaremos diferentes métodos não paramétricos para suavizar esta relação.
2.2 IMPACTOS EM CAPACETE
A base de dados mcycle, disponível no
pacote MASS do R, consiste em medições
coletadas durante um experimento que investigou os efeitos de impactos
em capacetes usados por motociclistas. Esses dados são
amplamente utilizados para demonstrar métodos
de suavização não paramétrica devido à sua natureza não linear
e à relação complexa entre as variáveis. As variáveis
coletadas em 133 mensurações foram:
times: Tempo decorrido após o impacto em Milisegundosaccel: Aceleração registrada no capacete (medido em g)
O comando a seguir traz a base de dados
mcycle## times accel
## 1 2.4 0.0
## 2 2.6 -1.3
## 3 3.2 -2.7
## 4 3.6 0.0
## 5 4.0 -2.7
## 6 6.2 -2.7
## 7 6.6 -2.7
## 8 6.8 -1.3
## 9 7.8 -2.7
## 10 8.2 -2.7
## 11 8.8 -1.3
## 12 8.8 -2.7
## 13 9.6 -2.7
## 14 10.0 -2.7
## 15 10.2 -5.4
## 16 10.6 -2.7
## 17 11.0 -5.4
## 18 11.4 0.0
## 19 13.2 -2.7
## 20 13.6 -2.7
## 21 13.8 0.0
## 22 14.6 -13.3
## 23 14.6 -5.4
## 24 14.6 -5.4
## 25 14.6 -9.3
## 26 14.6 -16.0
## 27 14.6 -22.8
## 28 14.8 -2.7
## 29 15.4 -22.8
## 30 15.4 -32.1
## 31 15.4 -53.5
## 32 15.4 -54.9
## 33 15.6 -40.2
## 34 15.6 -21.5
## 35 15.8 -21.5
## 36 15.8 -50.8
## 37 16.0 -42.9
## 38 16.0 -26.8
## 39 16.2 -21.5
## 40 16.2 -50.8
## 41 16.2 -61.7
## 42 16.4 -5.4
## 43 16.4 -80.4
## 44 16.6 -59.0
## 45 16.8 -71.0
## 46 16.8 -91.1
## 47 16.8 -77.7
## 48 17.6 -37.5
## 49 17.6 -85.6
## 50 17.6 -123.1
## 51 17.6 -101.9
## 52 17.8 -99.1
## 53 17.8 -104.4
## 54 18.6 -112.5
## 55 18.6 -50.8
## 56 19.2 -123.1
## 57 19.4 -85.6
## 58 19.4 -72.3
## 59 19.6 -127.2
## 60 20.2 -123.1
## 61 20.4 -117.9
## 62 21.2 -134.0
## 63 21.4 -101.9
## 64 21.8 -108.4
## 65 22.0 -123.1
## 66 23.2 -123.1
## 67 23.4 -128.5
## 68 24.0 -112.5
## 69 24.2 -95.1
## 70 24.2 -81.8
## 71 24.6 -53.5
## 72 25.0 -64.4
## 73 25.0 -57.6
## 74 25.4 -72.3
## 75 25.4 -44.3
## 76 25.6 -26.8
## 77 26.0 -5.4
## 78 26.2 -107.1
## 79 26.2 -21.5
## 80 26.4 -65.6
## 81 27.0 -16.0
## 82 27.2 -45.6
## 83 27.2 -24.2
## 84 27.2 9.5
## 85 27.6 4.0
## 86 28.2 12.0
## 87 28.4 -21.5
## 88 28.4 37.5
## 89 28.6 46.9
## 90 29.4 -17.4
## 91 30.2 36.2
## 92 31.0 75.0
## 93 31.2 8.1
## 94 32.0 54.9
## 95 32.0 48.2
## 96 32.8 46.9
## 97 33.4 16.0
## 98 33.8 45.6
## 99 34.4 1.3
## 100 34.8 75.0
## 101 35.2 -16.0
## 102 35.2 -54.9
## 103 35.4 69.6
## 104 35.6 34.8
## 105 35.6 32.1
## 106 36.2 -37.5
## 107 36.2 22.8
## 108 38.0 46.9
## 109 38.0 10.7
## 110 39.2 5.4
## 111 39.4 -1.3
## 112 40.0 -21.5
## 113 40.4 -13.3
## 114 41.6 30.8
## 115 41.6 -10.7
## 116 42.4 29.4
## 117 42.8 0.0
## 118 42.8 -10.7
## 119 43.0 14.7
## 120 44.0 -1.3
## 121 44.4 0.0
## 122 45.0 10.7
## 123 46.6 10.7
## 124 47.8 -26.8
## 125 47.8 -14.7
## 126 48.8 -13.3
## 127 50.6 0.0
## 128 52.0 10.7
## 129 53.2 -14.7
## 130 55.0 -2.7
## 131 55.0 10.7
## 132 55.4 -2.7
## 133 57.6 10.7Esta base de dados será utilizada para os
exercícios onde vocês trabalharão com as
variáveis times e accel
3 KERNEL PARA DENSIDADES
Para estimar as densidades a partir do Kernel,
utilizaremos a função density do pacote
base do R. Essa função tem os
seguintes argumentos:
x= dados da variável xbw= largura de banda (bandwidth)kernel= “gaussian” (normal), “epanechnikov” e “rectangular” (uniforme)
Importante destacar que o argumento bw, além de
aceitar números, também aceita os caracteres:
nrd0= É a regra de Silverman. Default doRe ideal para distribuições unimodais ou aproximadamente normais.SJ-ste= Métodos plug-in baseados no minimizador do erro quadrático integrado (MISE). Ideal para distribuições multimodais e assimétricas.ucv= Unbiased Cross-Validation. Ideal para distribuições muito complexas ou quando é necessária uma alta precisão.bcv= Biased Cross-Validation. Resulta em uma suavização ligeiramente maior, útil se você deseja evitar superajustes aos dados.
Para descobrir o valor de fato de \(h\) com cada um destes argumentos,
basta usar as funções bw.nrd0(x),
bw.SJ(x), bw.ucv(x) e bw.bcv(x),
respectivamente.
Para inspecionar, inicialmente, o efeito de \(h\) na estimação da densidade, vamos considerar o Kernel gaussiano e diferentes valores de \(h\) = 2.5, 3 e 3.5:
# Kernel Gaussiano com diferentes valores de h
dens_gaussian1 <- density(mtcars$mpg, kernel = "gaussian",bw = 2.5)
dens_gaussian2 <- density(mtcars$mpg, kernel = "gaussian",bw = 3)
dens_gaussian3 <- density(mtcars$mpg, kernel = "gaussian",bw = 3.5)Agora vamos incorporar a densidade estimada no histograma:
ggplot(mtcars, aes(x = mpg)) +
geom_histogram(aes(y = ..density..),bins = k, fill = "lightblue", color = "black", alpha = 0.7) +
geom_line(data = data.frame(x = dens_gaussian1$x, y = dens_gaussian1$y), aes(x = x, y = y, color = "h=2.5"), size = 1) +
geom_line(data = data.frame(x = dens_gaussian2$x, y = dens_gaussian2$y), aes(x = x, y = y, color = "h=3"), size = 1) +
geom_line(data = data.frame(x = dens_gaussian3$x, y = dens_gaussian3$y), aes(x = x, y = y, color = "h=3.5"), size = 1) +
labs(title = "Histograma e Densidades com Diferentes Kernels", x = "MPG", y = "Densidade") +
scale_color_manual(name = "Kernel", values = c("h=2.5" = "red", "h=3" = "green", "h=3.5" = "blue")) +
theme_minimal() +
theme(legend.position = "right") +
ggtitle("Histograma de 'MPG' com Densidade e Diferentes Kernels") +
theme(plot.title = element_text(hjust = 0.5))
Observa-se que quanto menor o valor de h, menor a suavização. Por outro lado, quanto maior o valor de h, maior a suavização. Como os dados são assimétricos, vamos utilizar o argumento “SJ-ste” e alternar entre as funções kernel (uniforme, gaussiano e epanechnikov):
# Tamanho h segundo o argumento "SJ-ste":
bw.SJ(mtcars$mpg)## [1] 2.801755# Kernels uniforme, gaussiano e epanechnikov:
dens_uniform <- density(mtcars$mpg, kernel = "rectangular",bw = "SJ-ste")
dens_gaussian <- density(mtcars$mpg, kernel = "gaussian",bw = "SJ-ste")
dens_epanechnikov <- density(mtcars$mpg, kernel = "epanechnikov",bw = "SJ-ste")Agora vamos incorporar a densidade estimada no histograma:
ggplot(mtcars, aes(x = mpg)) +
geom_histogram(aes(y = ..density..),bins=k, fill = "lightblue", color = "black", alpha = 0.7) +
geom_line(data = data.frame(x = dens_uniform$x, y = dens_uniform$y), aes(x = x, y = y, color = "Uniforme"), size = 1) +
geom_line(data = data.frame(x = dens_gaussian$x, y = dens_gaussian$y), aes(x = x, y = y, color = "Gaussiano"), size = 1) +
geom_line(data = data.frame(x = dens_epanechnikov$x, y = dens_epanechnikov$y), aes(x = x, y = y, color = "Epanechnikov"), size = 1) +
labs(title = "Histograma e Densidades com Diferentes Kernels", x = "MPG", y = "Densidade") +
scale_color_manual(name = "Kernel", values = c("Uniforme" = "red", "Gaussiano" = "green", "Epanechnikov" = "blue")) +
theme_minimal() +
theme(legend.position = "right") +
ggtitle("Histograma de 'MPG' com Densidade e Diferentes Kernels") +
theme(plot.title = element_text(hjust = 0.5))
Ambas as funções Kernel Epanechnikov e Gaussiano se adaptaram muito bem aos dados e são muito similares.
3.1 EXERCÍCIO: IMPACTOS EM CAPACETE
Para a variável times:
1- Faça um histograma com a regra de Sturges.
2- Investigue o efeito do kernel Gaussiano com h=6, 8 e 10.
3- Faça a sua estimativa final de densidade utilizando kernel nestes dados.
4 SUAVIZADOR KERNEL
Para fazer o suavizador Kernel no R,
utilizaremos a função ksmooth do pacote base do
R. Essa função possui os seguintes
argumentos:
x: os valores da variavel x.y: os valores da variável y.kernel: kernel utilizado (box= uniforme enormal= gaussiano).bandwidth: largura de banda.x.points: valores x para suavização.
Por simplicidade, vamos explorar somente o
kernel gaussiano. Para avaliar a associação
entre as variáveis hp e mpg, vamos
investigar este kernel para diferentes larguras de banda (h=10,
h=20 e h=50):
# As variáveis x e y
x <- mtcars$hp
y <- mtcars$mpg
# Valores de x de interesse para suavização
x.points <- seq(min(x),max(x),1)
# Suavizador kernel para h = 10
suaviza10 <- ksmooth(x, y, kernel = "normal", bandwidth = 10, x.points=x.points)
data_h10 <- data.frame(hp = suaviza10$x, mpg = suaviza10$y, h = "h = 10")
# Suavizador kernel para h = 20
suaviza20 <- ksmooth(x, y, kernel = "normal", bandwidth = 20, x.points=x.points)
data_h20 <- data.frame(hp = suaviza20$x, mpg = suaviza20$y, h = "h = 20")
# Suavizador kernel para h = 50
suaviza50 <- ksmooth(x, y, kernel = "normal", bandwidth = 50, x.points=x.points)
data_h50 <- data.frame(hp = suaviza50$x, mpg = suaviza50$y, h = "h = 50")
suaviza <- rbind(data_h10, data_h20, data_h50)Agora vamos fazer o gráfico que possui os dados originais e as respectivas suavizações:
# Gráfico de dispersões com as suavizações gaussianas
ggplot() +
geom_point(data = mtcars, aes(x = hp, y = mpg), color = "blue", size = 2, alpha = 0.7) +
geom_line(data = suaviza, aes(x = hp, y = mpg, color = h), size = 1) +
scale_color_manual(values = c("red", "green", "orange")) +
labs(title = "Kernel Smoothing Normal com Diferentes Larguras de Banda",
x = "Horsepower (hp)", y = "Miles per Gallon (mpg)",
color = "Bandwidth") +
theme_minimal()
Conforme pode ser visto neste gráfico, quanto maior a largura de banda h, maior a suavização. Entretanto, é de interesse encontrar o valor ótimo de \(h\). Utilizaremos os métodos da validação cruzada e empírico para tal finalidade.
4.1 CRITÉRIO DE SILVERMAN
Vimos em sala o critério de Silverman para a escolha do valor
de \(h\). No R, ele pode
ser implementado com o seguinte comando:
x <- mtcars$hp
silverman_h <- 0.9 * min(sd(x), IQR(x) / 1.34) * length(x)^(-1/5)
silverman_h## [1] 28.04104Assim, obtemos um valor de h = 28.04
4.2 VALIDAÇÃO CRUZADA
Para cada valor fixado de h, removemos uma observação da base de dados, por vez. Após a remoção, o kernel é aplicado nas n-1 observações restantes e estima-se o valor da resposta para a observação faltante. Dessa forma, podemos calcular medidas de acurácia como o erro quadrático e, ao final, obter uma média destes erros para todas as observações removidas. O valor de h que fornecer o menor erro quadrático médio será escolhido.
Para os possíveis valores de \(h\), utilizaremos o intervalo:
\[ h \in [\frac{amplitude}{n};amplitude] \]
Tem-se os seguintes comandos para que seja efetuada a validação cruzada:
# Dados de exemplo
x <- mtcars$hp
y <- mtcars$mpg
n=length(mtcars$mpg)
# Intervalo de valores para h
inf = (max(x)-min(x))/n
sup = max(x)-min(x)
h_values <- seq(ceiling(inf),ceiling(sup), by = 1)
# Função para cálculo de erro usando Validação Cruzada
cv_error <- function(h) {
n <- length(x)
errors <- numeric(n)
for (i in 1:n) {
# Deixar um ponto fora
x_novo <- x[-i]
y_novo <- y[-i]
x_test <- x[i]
# Ajuste do kernel smoothing
y_pred <- ksmooth(x_novo, y_novo, kernel = "normal", bandwidth = h, x.points = x_test)$y
# Calcular o erro
errors[i] <- (y[i] - y_pred)^2
}
# Retornar o erro médio e removendo indeterminações
mean(errors,na.rm=T)
}
# Aplicar CV para diferentes valores de h
errors <- sapply(h_values, cv_error)
# Encontrar o h que minimiza o erro
best_h <- h_values[which.min(errors)]
best_h## [1] 13Assim, obtemos um valor h de 13. Para finalizar, vamos comparar ambos os valores de h:
# Suavizador gaussiano com h=Silverman
suaviza_sil <- ksmooth(x, y, kernel = "normal", bandwidth = silverman_h, x.points=x.points)
data_sil <- data.frame(hp = suaviza_sil$x, mpg = suaviza_sil$y, h = "h = Silverman")
# Suavizador gaussiano com h=VC
suaviza_vc <- ksmooth(x, y, kernel = "normal", bandwidth = best_h, x.points=x.points)
data_vc <- data.frame(hp = suaviza_vc$x, mpg = suaviza_vc$y, h = "h = Validação Cruzada")
suaviza <- rbind(data_sil, data_vc)
# Gráfico de dispersões com as suavizações gaussianas
ggplot() +
geom_point(data = mtcars, aes(x = hp, y = mpg), color = "blue", size = 2, alpha = 0.7) +
geom_line(data = suaviza, aes(x = hp, y = mpg, color = h), size = 1) +
scale_color_manual(values = c("red", "green")) +
labs(title = "Kernel Smoothing Normal com Diferentes Larguras de Banda",
x = "Horsepower (hp)", y = "Miles per Gallon (mpg)",
color = "Bandwidth") +
theme_minimal()
A suavização oriunda pelo Silverman é mais suave e se demonstra mais interessante.
4.2.1 EXERCÍCIOS: IMPACTOS EM CAPACETE
Para as variáveis accel(y) e times (x) da
base de dados mcycle:
1- Investigue a suavização kernel com diferentes valores de h para suavizar a relação entre essas variáveis.
2- Utilize os critérios de Silverman e Validação Cruzada para guiar a sua escolha sobre a melhor largura da janela h.