14 Explorando distribuições de probabilidade empíricas

Na Sessão 12 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.

14.1 Estimação de densidades

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 27.

  > data(precip)
  > hist(precip, prob = T, main = "")
  > precip.d <- density(precip)
  > lines(precip.d)


PIC

Figura 27: Histograma para os dados precip e a densidade estimada usando a função density.


Portanto podemos ver que density() “suaviza” o histograma, capturando e concentrando-se nos principais aspectos dos dados disponíveis. Vamos ver na Figura 28 uma outra forma de visualizar os dados e sua densidade estimada, agora sem fazer o histograma.

  > plot(precip.d)
  > rug(precip)


PIC

Figura 28: Dados precip e a densidade estimada usando a função density.


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, 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, ...)
  NULL

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 29. 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)


PIC

Figura 29: Densidade estimada usando a função density com diferentes valores para o argumento bw.


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 30.

  > (kernels <- eval(formals(density.default)$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, 0.4, legend = kernels, col = seq(kernels), lty = 1,
  +     cex = 0.8, y.int = 1)


  [1] "gaussian"     "epanechnikov" "rectangular"  "triangular"   "biweight"
  [6] "cosine"       "optcosine"

PIC

Figura 30: Diferentes núcleos implementados pela função density.


Utilizando diferentes núcleos no conjunto de dados precip obtemos os resultados mostrados na Figura 31. 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))


PIC

Figura 31: Densidade estimada usando a função density com diferentes valores para o argumento kernel.


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 32.

  > precip.dSJ <- density(precip, bw = "sj")
  > plot(precip.dSJ)
  > rug(precip)


PIC
Figura 32: Densidade estimada para os dados precip usando a função density com critério de Sheather & Jones para seleção da largura de banda.


Os detalhes sobre os diferentes métodos implementados estão na documentação de bw.nrd(). Na Figura 33 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)


PIC
Figura 33: Diferentes métodos para largura de banda implementados pela função density.


14.2 Exercícios

1.
Carregar o conjunto de dados faithful e obter estimação de densidade para as variáveis tempo de erupção e duração da erupção.
2.
Carregar o conjunto airquality e densidades estimadas para as 4 variáveis medidas neste conjunto de dados.
3.
Rodar e estudar os exemplos da sessão examples da documentação da função density.