Na Sessão vimos com usar distribuições de probabilidade no R. Estas distribuições tem expressões conhecidas e são indexadas por um ou mais parâmetros. Portanto, conhecer a distribuição e seu(s) parâmetro(s) é suficiente para caracterizar completamente o comportamento distribuição e extrair resultados de interesse.
Na prática em estatística em geral somente temos disponível uma amostra e não conhecemos o mecanismo (distribuição) que gerou os dados. Muitas vezes o que se faz é: (i) assumir que os dados são provenientes de certa distribuição, (ii) estimar o(s) parâmetro(s) a partir dos dados. Depois disto procura-se verificar se o ajuste foi ``bom o suficiente'', caso contrário tenta-se usar uma outra distribuição e recomeça-se o processo.
A necessidade de estudar fenômenos cada vez mais complexos levou ao desenvolvimento de métodos estatísticos que às vezes requerem um flexibilidade maior do que a fornecida pelas distribuições de probabilidade de forma conhecida. Em particular, métodos estatísticos baseados em simulação podem gerar amostras de quantidades de interesse que não seguem uma distribuição de probabilidade de forma conhecida. Isto ocorre com frequência em métodos de inferência Bayesiana e métodos computacionalmente intensivos como bootstrap, testes Monte Carlo, dentre outros.
Nesta sessão vamos ver como podemos, a partir de um conjunto de dados explorar os possíveis formatos da distribuição geradora sem impor nenhuma forma paramétrica para função de densidade.
A estimação de densidades é implementada no R pela função density
.
O resultado desta função é bem simples e claro: ela produz uma função de densidade obtida a partir dos dados sem forma paramétrica
conhecida. Veja este primeiro exemplo que utiliza o conjunto de dados precip
que já vem com o R
e contém valores médios de precipitação em 70 cidades americanas.
Nos comandos a seguir vamos carregar o conjunto de dados, fazer um histograma de frequências relativas e depois adicionar a este histograma a linha de densidade estimada, conforma mostra a Figura .
> data(precip) > hist(precip, prob=T) > precip.d <- density(precip) > lines(precip.d)
Portanto podemos ver que a função density
``suaviza'' o histograma, capturando e
concentrando-se nos principais aspectos dos dados disponíveis.
Vamos ver na Figura
uma outra forma de visualizar os dados e sua densidade estimada, agora sem fazer o histograma.
> plot(precip.d) > rug(precip)
Embora os resultados mostrados acima seja simples e fáceis de entender, há muita coisa por trás deles!
Não vamos aqui estudar com detalhes esta função e os fundamentos teóricos nos quais se baseiam esta implementação computacional pois isto estaria muito além dos objetivos e escopo deste curso.
Vamos nos ater às informações principais que nos permitam compreender o básico necessário sobre o uso da função.
Para maiores detalhes veja as referências na documentação da função, que pode ser vista digitando help(density)
Basicamente, a função density
produz o resultado visto anteriormente
criando uma sequência de valores no eixo-X e estimando a densidade em cada ponto
usando os dados ao redor deste ponto.
Podem ser dados pesos aos dados vizinhos de acordo com sua proximidade ao ponto a ser estimado.
Vamos examinar os argumentos da função.
> args(density) function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"), window = kernel, width, give.Rkern = FALSE, n = 512, from, to, cut = 3, na.rm = FALSE)Os dois argumentos chave são portanto
bw
e kernel
que controlam a distância na qual se procuram vizinhos e o peso a ser dado a cada vizinho, respectivamente.
Para ilustrar isto vamos experimentar a função com diferentes valores para o argumento bw
.
Os resultados estão na Figura .
Podemos notar que o grau de suavização aumenta a medida de aumentamos os valores deste argumento
e as densidades estimadas podem ser bastante diferentes!
plot(density(precip, bw=1), main='') rug(precip) lines(density(precip, bw=5), lty=2) lines(density(precip, bw=10), lty=3) legend(5, 0.045, c('bw=1', 'bw=5', 'bw=10'), lty=1:3)
O outro argumento importante é tipo de função de pesos, ao que chamamos de núcleo (kernel). O R implementa vários núcleos diferentes cujos formatos são mostrados na Figura .
(kernels <- eval(formals(density)$kernel)) plot (density(0, bw = 1), xlab = "", main="kernels com bw = 1") for(i in 2:length(kernels)) lines(density(0, bw = 1, kern = kernels[i]), col = i) legend(1.5,.4, legend = kernels, col = seq(kernels), lty = 1, cex = .8, y.int = 1)
Utilizando diferentes núcleos no conjunto de dados precip
obtemos os resultados mostrados
na Figura .
Note que as densidades estimadas utilizando os diferentes núcleos são bastante similares!
plot(density(precip), main='') rug(precip) lines(density(precip, ker='epa'), lty=2) lines(density(precip, ker='rec'), col=2) lines(density(precip, ker='tri'), lty=2, col=2) lines(density(precip, ker='biw'), col=3) lines(density(precip, ker='cos'), lty=3, col=3) legend(0, 0.035, legend = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine"), lty = rep(1:2,3), col = rep(1:3, each=2))
Portanto, inspecionando os resultados anteriores podemos concluir que
a largura de banda (bandwidth - bw
)
é o que mais influencia a estimação de densidade, isto é, é o argumento mais importante.
O tipo de núcleo (kernel) é de importância secundária.
Bem, a esta altura voce deve estar se perguntando:
mas como saber qual a largura de banda adequada?
A princípio podemos tentar diferentes valores no argumento bw
e inspecionar os resultados.
O problema é que esta escolha é subjetiva.
Felizmente para nós vários autores se debruçaram sobre este problema e descobriram
métodos automáticos de seleção que que comportam bem na maioria das situações práticas.
Estes métodos podem ser especificados no mesmo argumento bw
, passando agora para este argumento
caracteres que identificam o valor, ao invés de um valor numérico.
No comando usado no início desta sessão onde não especificamos o argumento bw
foi utilizado o valor ``default'' que é o método nrd0
que implementa a regra prática de Silverman.
Se quisermos mudar isto para o método de Sheather & Jones podemos fazer como nos comandos abaixo
que produzem o resultado mostrado na Figura .
precip.dSJ <- density(precip, bw='sj') plot(precip.dSJ) rug(precip)
Os detalhes sobre os diferentes métodos implementados estão na documentação da função bw.nrd
.
Na Figura ilustramos resultados obtidos com os diferentes métodos.
data(precip) plot(density(precip, n=1000)) rug(precip) lines(density(precip, bw="nrd"), col = 2) lines(density(precip, bw="ucv"), col = 3) lines(density(precip, bw="bcv"), col = 4) lines(density(precip, bw="SJ-ste"), col = 5) lines(density(precip, bw="SJ-dpi"), col = 6) legend(55, 0.035, legend = c("nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi"), col = 1:6, lty = 1)
Faithful
e obter estimação de densidade para as variáveis 'tempo de erupção' e 'duração da erupção.
airquality
e densidades estimadas para as 4 variáveis medidas neste conjunto de dados.