O R para Ajuste de Distribuições


26 de novembro de 2007

A Estatística Descritiva possui ferramentas para que o analista quantitativo elabore modelos probabilísticos para as variáveis analisadas. O R possui um conjunto de funções para o ajuste de distribuições de probabilidade aos dados.

Abaixo, apresentamos uma lista de funções do R utilizadas para auxiliar o ajuste de distribuições de probabilidade.

hist(): histograma de frequência  
stem(): o ramo-e-folhas  
density(): ajusta densidade de forma não-paramétrica  
cut(): construção de faixas de valores para uma variável  
ecdf() : função de distribuição acumulada empírica  
fitdistr() : ajuste de distribuição univariada por máxima verossimilhança  
chisq.test(): teste qui-quadrado  
shapiro.test(): teste para a normalidade  
jarque.bera.test():teste para a normalidade  
ks.test(): teste de Kolmogorov-Smirnov  
lillie.test() : teste para normalidade  
ad.test(): teste para normalidade  
qqnorm(): o normal qqplot

Os gráficos exploratórios e as medidas-resumo podem funcionar como primeira abordagem para identificar uma distribuição de probabilidades aderente aos dados.

Histograma

O histograma é uma das formas mais simples de inferir sobre a função densidade de probabilidade de uma variável. Embora seja utilizado com freqüência, representa um ajuste pouco suave.

Considere o seguinte exemplo que gera dados da distribuição normal e explora o histograma como ferramenta de ajuste de distribuição.

  > #Iniciando a semente
  > set.seed(51007)
  > #Definindo o tamanho da amostra
  > n<-100
  > #Gerando dados da distribuição normal
  > x.norm<-rnorm(n)


  > #Construindo um histograma de frequencias
  > h<-hist(x.norm,prob=T)
  > #Visualizando a dispersão dos pontos nos intervalos.
  > rug(x.norm)
  > #Ajustando a curva normal aos dados
  > curve(dnorm(x),-3,3,add=T)

pict

Figura 1: Histograma, ajuste de curva normal teórica


As informações armazenadas no objeto h são:

  > names(h)
  [1] "breaks"      "counts"      "intensities" "density"     "mids"
  [6] "xname"       "equidist"
breaks: limites superiores dos intervalos de classe  
counts: frequência absoluta nos intervalos de classe  
intensities: densidade de frequência  
density: densidade de frequência  
mids: pontos médios dos intervalos  
xname: nome da variável  
equidist: variável lógica indicando que os intervalos possuem o mesmo comprimento

A partir das informações geradas pelo histograma, vamos remontar a tabela de freqüências em classes.

O comando apresentado na seqüência transforma a variável x.norm, originalmente quantitativa, em uma variável qualitativa cujos níveis representam as faixas de valores de uma tabela de freqüência em classes.

  > # Criando os intervalos de classe
  > classes<-cut(x.norm,breaks=h$breaks)
  > # Tabela de frequencias em classes
  > table(classes)
  classes
  (-3,-2.5] (-2.5,-2] (-2,-1.5] (-1.5,-1] (-1,-0.5]  (-0.5,0]   (0,0.5]   (0.5,1]
          1         2        10         9        16         7        24        15
    (1,1.5]   (1.5,2]   (2,2.5]
          9         5         2





ClassesFreq. Abs.



1 (-3,-2.5] 1
2 (-2.5,-2] 2
3 (-2,-1.5] 10
4 (-1.5,-1] 9
5 (-1,-0.5] 16
6 (-0.5,0] 7
7 (0,0.5] 24
8 (0.5,1] 15
9 (1,1.5] 9
10 (1.5,2] 5
11 (2,2.5] 2




Tabela 1: Tabela de Frequencia em Faixas de Valores

Na Tabela 1 é possível visualizar as freqüências absolutas nos intervalos gerados ao construir, de modo automático, o histograma. As informações nas tabelas podem ser mais detalhadas como acontece na Tabela 2.







ClassesPonto MédioFreq. Abs.Freq.Rel.





1 (-3,-2.5] -2.75 1 0.01
2 (-2.5,-2] -2.25 2 0.02
3 (-2,-1.5] -1.75 10 0.1
4 (-1.5,-1] -1.25 9 0.09
5 (-1,-0.5] -0.75 16 0.16
6 (-0.5,0] -0.25 7 0.07
7 (0,0.5] 0.25 24 0.24
8 (0.5,1] 0.75 15 0.15
9 (1,1.5] 1.25 9 0.09
10 (1.5,2] 1.75 5 0.05
11 (2,2.5] 2.25 2 0.02






Tabela 2: Tabela de Frequencia em Faixas de Valores

Vamos calcular algumas estatísticas com base na tabela de freqüência produzida.

  > #Frequencias Absoluntas nas classes
  > ni<-h$counts
  > ni
   [1]  1  2 10  9 16  7 24 15  9  5  2
  > # Frequencias Relativas nas classes
  > fi<-ni/n
  > fi
   [1] 0.01 0.02 0.10 0.09 0.16 0.07 0.24 0.15 0.09 0.05 0.02
  > # Média com base na tabela de frequência
  > xbarra<-h$mids%*%fi #produto de vetores
  > xbarra
        [,1]
  [1,] -0.08
  > # Desvios dos pontos médios em relação à média
  > desv<-h$mids-xbarra
  > desv
   [1] -2.67 -2.17 -1.67 -1.17 -0.67 -0.17  0.33  0.83  1.33  1.83  2.33
  > # Desvios quadráticos dos pontos médios em relação à média
  > desv2<-desv^2
  > desv2
   [1] 7.1289 4.7089 2.7889 1.3689 0.4489 0.0289 0.1089 0.6889 1.7689 3.3489
  [11] 5.4289
  > # Variância
  > varianc<-desv2%*%fi
  > varianc
         [,1]
  [1,] 1.2061

Conforme a Tabela 3 é possível constatar a diferença entre estatísticas calculadas em dados brutos e dados em classes.





dados em classesdados brutos



média -0.08 -0.10
variância 1.21 1.16




Tabela 3: Estatísticas Descritivas da Tabela em Intervalos e Dados Brutos

Ramo-e-Folhas

O Ramo-e-Folhas é um gráfico que possui as mesmas características do histograma, com a grande vantagem de reproduzir os valores numéricos e, consequentemente, não há perda de informação.

  > # grafico de ramo-e-folhas
  > stem(x.norm)

    The decimal point is at the |
  
    -2 | 5
    -2 | 410
    -1 | 9987766655
    -1 | 3221111000
    -0 | 9999888877766655
    -0 | 44220
     0 | 000001112222333344444
     0 | 55555556677778889
     1 | 0111112233
     1 | 66689
     2 | 02

Repare que, na figura acima, o formato da distribuição é semelhante ao exibido pelo histograma e o número de intervalos é exatamente o mesmo.

Densidade Não-Paramétrica

Uma forma mais sofisticada de estimar a densidade de probabilidade dos dados é feita através do comando density.


  > hist(x.norm,prob=T)
  > rug(x.norm)
  > curve(dnorm(x),-3,3,add=T,col='red')
  > lines(density(x.norm),col='blue')

pict

Figura 2: Função de Densidade Não-Paramétrica


O resultado pode ser visto como uma variante do histograma, porém suavizado.

Função de Distribuição Acumulada Empírica

Através do comando ecdf é possível encontrar, e plotar, Fe, a Função de Distribuição Acumulada Empírica ,ou seja, as probabilidades acumuladas com base nos dados.

Veja o exemplo da amostra abaixo

  > x<-c(3,5,7,9,11,13,17)
  > x

  [1]  3  5  7  9 11 13 17

O comportamento de Fe para este conjunto de dados é obtido da seguinte forma:

  > # Tamanho da amostra
  > n<-length(x)
  > # Função de Distribuição Empírica
  > Fe<-seq(1:n)/n
  > Fe

  [1] 0.1428571 0.2857143 0.4285714 0.5714286 0.7142857 0.8571429 1.0000000

O gráfico desta função é obtido por :


  > # Gráfico da Função de Distribuição Acumulada Empírica
  > plot(ecdf(x),ylab="Probabilidade Acumulada")

pict

Figura 3: Função de Distribuição Acumulada Empírica


  > #Função de Distribuição Acumulada Empírica Suavizada
  > Fe.s<-((seq(1:n)-0.5)/n)
  > Fe.s<-c(0,Fe.s,1)
  > Fe.s
  [1] 0.00000000 0.07142857 0.21428571 0.35714286 0.50000000 0.64285714 0.78571429
  [8] 0.92857143 1.00000000


  > plot(ecdf(x),ylab="Probabilidade Acumulada")
  > lines(c(min(x)-0.5,x,max(x)+0.5),Fe.s,col='red')

pict

Figura 4: Suavização da Função de Distribuição Acumulada Empírica


Com a Função de Distribuição Acumulada Suavizada, podemos estimar qualquer quantil que for de interesse.

Exercícios

Os dados no arquivo a ser carregado contém o desempenho dos alunos no Processo Seletivo Estendido no ano de 2007.

  > # Leitura do conjunto de dados
  > pse<-read.csv('http://www.leg.ufpr.br/~joel/dados/pse2007.csv')
  > # Exibição das 6 primeiras linhas
  > head(pse)

    X CURSO mat bio qui geo fis por lit hist LEM
  1 1 EST N   0   0   5   3   2  10   3    1   5
  2 2 EST N   3   3   4   7   4   8   2    4   1
  3 3 EST N   5   2   2   3   3   7   4    2   2
  4 4 EST N   2   2   3   4   6   7   1    4   2
  5 5 EST N   3   4   4   9   3   9   3    7   8
  6 6 EST N   2   3   1   3   5   8   1    5   1

  > # Exibição das 6 últimas linhas
  > tail(pse)

        X     CURSO mat bio qui geo fis por lit hist LEM
  459 459 Mat Ind T   2   3   3   2   4   3   1    3   3
  460 460 Mat Ind T   1   1   2   2   4   5   1    3   2
  461 461 Mat Ind T   3   1   3   3   3   8   4    3   5
  462 462 Mat Ind T   4   3   3   3   5   7   0    4   2
  463 463 Mat Ind T   3   4   3   6   6   8   2    4   4
  464 464 Mat Ind T   4   3   3   5   7   8   1    2   3

  1. Calcule as principais medidas resumo para estes dados

      > # Resumindo as informacoes de todas as variaveis.
      > summary(pse)

             X                     CURSO          mat             bio
       Min.   :  1.0   EST N          :161   Min.   :0.000   Min.   :0.000
       1st Qu.:116.8   MAT (Bac Lic) T:104   1st Qu.:1.000   1st Qu.:1.000
       Median :232.5   Mat Ind T      : 92   Median :2.000   Median :2.000
       Mean   :232.5   MAT (Lic) N    :107   Mean   :2.634   Mean   :2.422
       3rd Qu.:348.2                         3rd Qu.:4.000   3rd Qu.:3.000
       Max.   :464.0                         Max.   :8.000   Max.   :7.000
            qui             geo             fis             por            lit
       Min.   :1.000   Min.   :0.000   Min.   :0.000   Min.   : 2.0   Min.   :0.000
       1st Qu.:2.000   1st Qu.:3.000   1st Qu.:3.000   1st Qu.: 5.0   1st Qu.:1.000
       Median :3.000   Median :4.000   Median :3.000   Median : 7.0   Median :2.000
       Mean   :3.069   Mean   :3.836   Mean   :3.554   Mean   : 6.7   Mean   :1.922
       3rd Qu.:4.000   3rd Qu.:5.000   3rd Qu.:5.000   3rd Qu.: 8.0   3rd Qu.:3.000
       Max.   :8.000   Max.   :9.000   Max.   :8.000   Max.   :12.0   Max.   :6.000
            hist            LEM
       Min.   :0.000   Min.   :0.000
       1st Qu.:3.000   1st Qu.:2.000
       Median :4.000   Median :3.000
       Mean   :3.950   Mean   :3.237
       3rd Qu.:5.000   3rd Qu.:5.000
       Max.   :8.000   Max.   :8.000
  2. Ilustre graficamente a diferença entre os cursos

    Neste item, vamos primeiramente criar uma nova variavel que soma as notas de todas as provas realizadas pelo candidato.

      > # Criando uma nova variavel e anexando ao conjunto pse
      > pse$notafinal<-with(pse,mat+bio+qui+geo+fis+por+lit+hist+LEM)
      > # Resumindo as informacoes desta variavel
      > summary(pse$notafinal)

         Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
        11.00   26.00   30.00   31.33   35.00   62.00

    Utilizaremos o box-plot como ferramenta para comparação.


      > # Construindo o boxplot com as notas finais
      > boxplot(pse$notafinal~pse$CURSO)

    pict


  3. Examine a correlação entre os desempenhos nas diferentes disciplinas
      > # Selecionando todas as linhas e as colunas de 3 ate 11
      > discip<-pse[,3:11]
      > # Calculando a matriz de correlacao entre as notas
      > correlacoes<-cor(discip)
      > correlacoes

                  mat       bio         qui       geo        fis        por
      mat  1.00000000 0.2579604  0.13570161 0.2789203 0.34611361 0.16180803
      bio  0.25796043 1.0000000  0.15974264 0.2527182 0.23427131 0.15794839
      qui  0.13570161 0.1597426  1.00000000 0.1237300 0.14197504 0.09533673
      geo  0.27892025 0.2527182  0.12372997 1.0000000 0.21084688 0.40538627
      fis  0.34611361 0.2342713  0.14197504 0.2108469 1.00000000 0.13704096
      por  0.16180803 0.1579484  0.09533673 0.4053863 0.13704096 1.00000000
      lit  0.07527067 0.1160417  0.15972211 0.1640895 0.07056438 0.16296187
      hist 0.13068988 0.2051098  0.14710358 0.3606197 0.13940618 0.29867822
      LEM  0.15263761 0.1433847 -0.03162456 0.2150677 0.11234875 0.24979889
                  lit      hist         LEM
      mat  0.07527067 0.1306899  0.15263761
      bio  0.11604169 0.2051098  0.14338466
      qui  0.15972211 0.1471036 -0.03162456
      geo  0.16408949 0.3606197  0.21506770
      fis  0.07056438 0.1394062  0.11234875
      por  0.16296187 0.2986782  0.24979889
      lit  1.00000000 0.2250333  0.17445784
      hist 0.22503329 1.0000000  0.22367091
      LEM  0.17445784 0.2236709  1.00000000