Estatística Não-Paramétrica:
Projeto I

Arthur Silva1, Julyet Alves2, Luan Sousa3

20/09/2023

1 Introdução

O presente trabalho é parte dos processos avaliativos da disciplina Estatística Não-Paramétrica, ofertada pelo Departamento de Estatística e Matemática Aplicada (DEMA) da Universidade Federal do Ceará (UFC), no curso de Estatística, e ministrada pelo professor Manoel Santos4 (rpubs).

O intuito do trabalho é aplicar testes não-paramétricos a um conjunto de dados selecionado, a fim de inferir sobre de qual distribuição tais dados advêm.

2 Conhecendo os dados

A partir do seguinte vetor de dados (152.7, 172, 172.5, 173.3, 193, 204.7, 216.5, 234.9, 262.6, 422.6), de tamanho 10, buscamos encontrar qual distribuição de probabilidade melhor se ajusta a essa amostra, que é descrita em McCool(1974), tratando-se do tempo de vida em horas de 10 suportes de um certo tipo.

Podemos concluir, em um primeiro olhar, que tratam-se de dados contínuos, já que no caso discreto os números são inteiros não negativos. Dessa forma, já reduzimos bastante o espaço de distribuições possíveis.

Também podemos previamente, retirar alguns valores sumarizados dessa amostra, como média e mediana, e comparar tais valores com as propriedades das distribuições contínuas mais conhecidas. A tabela a seguir aponta alguns resultados.

m = density(dados)
moda = m$x[which.max(m$y)]
sumario = data.frame(rbind(summary(dados)),var(dados), moda, skewness(dados), kurtosis(dados))
colnames(sumario) = c('Mín.', '1º Q', 'Mediana', 'Média', '3º Q', 'Máx.', 'Var.', 'Moda', 'Assimetria', 'Curtose')
Mín. 1º Q Mediana Média 3º Q Máx. Var. Moda Assimetria Curtose
152.7 172.7 198.85 220.48 230.3 422.6 6147.44 182.95 1.59 1.52

Vemos que a média é maior do que a mediana, sendo a moda o menor valor dentre elas; e que a variância é relativamente alta. Também há a presença de um outlier, que é o valor máximo da amostra, já que este encontra-se longe da média.

A fim de confirmar se a amostra é simétrica ou assimétrica, podemos criar um histograma, mostrando a frequência dos valores e a curva da distribuição amostral, além de criar um boxplot para comprovar a presença do outlier:

Os gráficos nos apontam que estamos trabalhando com uma distribuição assimétrica à direita, limitando um pouco mais a quantidade de distribuições possíveis.

Finalmente, podemos começar a inferir com mais confiança de qual distribuição esses dados vêm, tendo em mente as seguintes observações:

  1. Distribuição contínua;
  2. \(\bar{x} > Md > Mo\);
  3. Assimetria à direita.

Podemos então, testar se nossos dados se ajustam adequadamente aos seguintes modelos, os quais podem atender às observações feitas anteriormente:

  • Exponencial;
  • Gama;
  • Log-Normal;
  • Weibull;
  • Birnbaum-Saunders.

3 Estimações e Testes

Inicialmente, para testar se nossos dados vêm ou não de cada distribuição acima utilizando métodos não paramétricos, precisamos ter estimativas dos parâmetros das distribuições que iremos testar, o que será encontrado com o pacote fitdistrplus, que nos dá estimativas de máxima verossimilhança por padrão.

Esse pacote também trás a função gofstat(), que testa a qualidade do ajuste das estimativas do parâmetro para a distribuição testada, além de performar testes não paramétricos (Kolmogorov-Smirnov, Cramer-von Mises e Anderson-Darling) para cada distribuição, usando os parâmetros estimados.

A Estimativa de Máxima Verossimilhança (EMV) é um método estatístico amplamente empregado para estimar os parâmetros de uma distribuição de probabilidade que melhor descreve um conjunto de dados. A ideia por trás da EMV é encontrar os valores dos parâmetros que tornam os dados observados mais prováveis sob a distribuição especificada.

3.1 Exponencial

Estimativa para \(\lambda\) caso \(\text{dados} \sim Exp(\lambda)\):

fit_exp = fitdist(dados,'exp'); fit_exp
## Fitting of the distribution ' exp ' by maximum likelihood 
## Parameters:
##         estimate Std. Error
## rate 0.004535559 0.00136012

Testes para a distribuição, levando em conta a hipótese \(H_{0}: \text{dados} \sim Exp(\hat{\lambda})\):

gof = gofstat(fit_exp); gof
## Goodness-of-fit statistics
##                              1-mle-exp
## Kolmogorov-Smirnov statistic 0.4997162
## Cramer-von Mises statistic   0.5508352
## Anderson-Darling statistic   2.5899552
## 
## Goodness-of-fit criteria
##                                1-mle-exp
## Akaike's Information Criterion  129.9161
## Bayesian Information Criterion  130.2187
# testes manuais, para ver o valor-p
ks.test(dados,'pexp',fit_exp$estimate)
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dados
## D = 0.49972, p-value = 0.007828
## alternative hypothesis: two-sided
cvm.test(dados,'pexp',fit_exp$estimate)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## omega2 = 0.55084, p-value = 0.02755
ad.test(dados,'pexp',fit_exp$estimate)
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## An = 2.59, p-value = 0.04563

Assim, analisando com os níveis de confiança mais usuais no teste de Kolmogorov-Smirnov, como 5%, 2% e 1%, rejeitamos a hipótese de que \(\text{dados} \sim Exp(\hat{\lambda})\).

Podemos ainda, comparar graficamente a distribuição empírica dos dados com a distribuição teórica, e também o histograma com a curva teórica, como mostra a figura a seguir, a fim de ilustrar nossa decisão:

3.2 Gama

Estimativa para \(\alpha\; \text{e}\; \lambda\) caso \(\text{dados} \sim Gama(\alpha;\lambda)\):

fit_gama = fitdist(dados,'gamma'); fit_gama
## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## shape 11.56190158 5.07021762
## rate   0.05244191 0.02349471

Testes para a distribuição, levando em conta a hipótese \(H_{0}: \text{dados} \sim Gama\left(\hat{\alpha};\hat{\lambda}\right)\):

gof_g = gofstat(fit_gama); gof_g
## Goodness-of-fit statistics
##                              1-mle-gamma
## Kolmogorov-Smirnov statistic  0.18533689
## Cramer-von Mises statistic    0.09871694
## Anderson-Darling statistic    0.68133410
## 
## Goodness-of-fit criteria
##                                1-mle-gamma
## Akaike's Information Criterion    115.2275
## Bayesian Information Criterion    115.8327
# testes manuais, para ver o valor-p
ks.test(dados,'pgamma', shape=fit_gama$estimate[1], rate=fit_gama$estimate[2])
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dados
## D = 0.18534, p-value = 0.8223
## alternative hypothesis: two-sided
cvm.test(dados,'pgamma', fit_gama$estimate[1], fit_gama$estimate[2])
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## omega2 = 0.098717, p-value = 0.6004
ad.test(dados,'pgamma', fit_gama$estimate[1], fit_gama$estimate[2])
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Gamma distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## An = 0.68133, p-value = 0.5705

Com isso, há fortes evidências para suspeitar que nossos dados seguem uma distribuição Gama com os parâmetros estimados. Iremos agora, plotar os gráficos de comparação entre distribuição empírica vs acumulada teórica; e o histograma dos dados vs curva teórica:

Logo, comprovamos que a distribuição Gama possui um ajuste relativamente bom para modelar tais dados, visto que as curvas teóricas ficaram bem próximas às amostrais, apesar do histograma revelar certa discrepância.

3.3 Log-Normal

Partindo agora para o modelo Log-Normal, queremos testar a hipótese de que nossos dados seguem essa lei. Porém, vamos primeiro estimar os parâmetros:

fit_lnorm = fitdist(dados,'lnorm'); fit_lnorm
## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate Std. Error
## meanlog 5.3519435 0.08814779
## sdlog   0.2787478 0.06232629

Com nossas estimativas em mãos, podemos então realizar nossos testes, bem como verificar a qualidade do ajuste, como mostram as saídas seguintes:

gof_ln = gofstat(fit_lnorm); gof_ln
## Goodness-of-fit statistics
##                              1-mle-lnorm
## Kolmogorov-Smirnov statistic  0.16334586
## Cramer-von Mises statistic    0.07630523
## Anderson-Darling statistic    0.54921679
## 
## Goodness-of-fit criteria
##                                1-mle-lnorm
## Akaike's Information Criterion    113.8687
## Bayesian Information Criterion    114.4739
# testes manuais, para ver o valor-p
ks.test(dados,'plnorm', fit_lnorm$estimate[1], fit_lnorm$estimate[2])
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dados
## D = 0.16335, p-value = 0.9146
## alternative hypothesis: two-sided
cvm.test(dados,'plnorm', fit_lnorm$estimate[1], fit_lnorm$estimate[2])
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## omega2 = 0.076305, p-value = 0.7244
ad.test(dados,'plnorm', fit_lnorm$estimate[1], fit_lnorm$estimate[2])
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: log-normal distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## An = 0.54922, p-value = 0.6929

Depreendemos dos testes acima, que há evidências muito fortes para não rejeitar a hipótese nula, pois temos, em todos os testes, os valores-p altos, inclusive, mais altos que os dos testes anteriores para as outras distribuições.

Para checar, iremos também construir os gráficos comparativos já vistos anteriormente:

É fácil ver que tal distribuição tem melhor ajuste que as anteriores, tanto por meio dos métodos numéricos quanto dos gráficos, pois nota-se uma aproximação maior entre a distribuição empírica e a distribuição acumulada, bem como entre a frequência relativa dos dados junto à densidade da função teórica.

3.4 Weibull

Iremos repetir os mesmos passos anteriores para chechar se nossos dados seguem distribuição Weibull:

fit_wei = fitdist(dados,'weibull'); fit_wei
## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters:
##         estimate Std. Error
## shape   2.935658  0.6335469
## scale 246.395648 28.3155506

Então, realizamos os testes com as estimativas:

gof_w = gofstat(fit_wei); gof_w
## Goodness-of-fit statistics
##                              1-mle-weibull
## Kolmogorov-Smirnov statistic     0.2193159
## Cramer-von Mises statistic       0.1546648
## Anderson-Darling statistic       0.9481572
## 
## Goodness-of-fit criteria
##                                1-mle-weibull
## Akaike's Information Criterion      118.6026
## Bayesian Information Criterion      119.2078
# testes manuais, para ver o valor-p
ks.test(dados, 'pweibull', fit_wei$estimate[1], fit_wei$estimate[2])
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dados
## D = 0.21932, p-value = 0.6462
## alternative hypothesis: two-sided
cvm.test(dados, 'pweibull', fit_wei$estimate[1], fit_wei$estimate[2])
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## omega2 = 0.15466, p-value = 0.3805
ad.test(dados, 'pweibull', fit_wei$estimate[1], fit_wei$estimate[2])
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  Parameters assumed to be fixed
## 
## data:  dados
## An = 0.94816, p-value = 0.3831

É possível ver que os testes nos apontam que a distribuição Weibull apenas se ajusta melhor quando comparada à distribuição Exponencial, pois tem valores-p menores do que os dos outros modelos.

Mostramos esse resultado graficamente, como mostra a seguir:

3.5 Birnbaum-Saunders

A distribuição Birnbaum-Saunders não está inclusa nos pacotes nativos de distribuições do R, e então, para poder estimá-la, usaremos a biblioteca extraDistr, que possui as funções de que precisamos.

Começamos estimando, manualmente, os parâmetros da distribuição, já que a função fitdist() necessita de valores iniciais quando usada para distribuições não listadas no pacote (ver Delignette-Muller, Dutang, e Siberchicot (2023)).

Utilizamos como base, os valores iniciais apropriados para este modelo, vistos em sala de aula:

xb = mean(dados)
hb = 1/mean(1/dados)

alfa_s = sqrt(2*sqrt(xb/hb)-1)
beta_s = sqrt(xb*hb)

xb; hb; alfa_s; beta_s
## [1] 220.48
## [1] 203.8853
## [1] 1.039134
## [1] 212.0204

Agora, podemos estimar normalmente:

fit_bs = fitdist(dados, 'fatigue', start = list(alpha = alfa_s, beta = beta_s)); fit_bs
## Fitting of the distribution ' fatigue ' by maximum likelihood 
## Parameters:
##          estimate  Std. Error
## alpha   0.2824708  0.06315268
## beta  212.0496397 18.75181333

E mais uma vez, iremos realizar os testes de qualidade de ajuste:

gof_bs = gofstat(fit_bs); gof_bs
## Goodness-of-fit statistics
##                              1-mle-fatigue
## Kolmogorov-Smirnov statistic    0.17069155
## Cramer-von Mises statistic      0.08289995
## Anderson-Darling statistic      0.57742322
## 
## Goodness-of-fit criteria
##                                1-mle-fatigue
## Akaike's Information Criterion      113.9435
## Bayesian Information Criterion      114.5487
# testes manuais
ks.test(dados, 'pfatigue', fit_bs$estimate[1], fit_bs$estimate[2])
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dados
## D = 0.17069, p-value = 0.8872
## alternative hypothesis: two-sided
cvm.test(dados, 'pfatigue', fit_bs$estimate[1], fit_bs$estimate[2])
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: distribution 'pfatigue'
##  Parameters assumed to be fixed
## 
## data:  dados
## omega2 = 0.0829, p-value = 0.6858
ad.test(dados, 'pfatigue', fit_bs$estimate[1], fit_bs$estimate[2])
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: distribution 'pfatigue'
##  Parameters assumed to be fixed
## 
## data:  dados
## An = 0.57742, p-value = 0.6652

Para elucidar nossas estatísticas, como é de praxe, usaremos métodos gráficos:

4 Comparação geral

Para melhor evidenciar e auxiliar na escolha correta da nossa distribuição, os gráficos e tabela a seguir (contendo os valores-p) trazem um condensando do que foi apanhado nas análises anteriores, a fim de facilitar a visualização e validar nossas decisões.

Exponencial Gama Log_Normal Weibull Birnbaum_Saunders
Kolmogorov-Smirnov 0.0078 0.8223 0.9146 0.6462 0.8872
Cramer-von Mises 0.0276 0.6004 0.7244 0.3805 0.6858
Anderson-Darling 0.0456 0.5705 0.6929 0.3831 0.6652

5 Conclusão

Partindo de uma amostra pequena — de apenas 10 elementos — e realizando os testes não-paramétricos, bem como utilizando-se de métodos gráficos, concluímos que a distribuição que mais se ajusta aos dados é o modelo Log-Normal, que performou melhor nos testes, bem como apresentou curvas teóricas mais próximas à distribuição empírica e ao histograma apresentados, além de ter estimadores com baixo erro padrão, apesar de que a distribuição proposta em (Lemonte 2006), para o mesmo conjunto, tenha sido o modelo Birnbaum-Saunders.

Referências

Delignette-Muller, Marie-Laure, Christophe Dutang, e Aurelie Siberchicot. 2023. "fitdistrplus: Help to Fit of a Parametric Distribution to Non-Censored or Censored Data". https://cran.r-project.org/web/packages/fitdistrplus/fitdistrplus.pdf.
Faraway, Julian, George Marsaglia, John Marsaglia, e Adrian Baddeley. 2022. "goftest: Classical Goodness-of-Fit Tests for Univariate Distributions". https://cran.r-project.org/web/packages/goftest/goftest.pdf.
Larson, Ron, Betsy Farber, e Cyro; traducão técnica Patarra. 2004. Estatı́stica aplicada. Prentice Hall.
Lemonte, Artur José. 2006. "Inferência sobre os parâmetros da distribuição Birnbaum-Saunders bi-paramétrica". https://repositorio.ufpe.br/bitstream/123456789/6339/1/arquivo7199_1.pdf.
Moore, David S. 2005. "A estatıstica básica e sua prática". LTC 6rd ed., Rio de Janeiro, Brazil.
Morettin, Pedro A, e Wilton O Bussab. 2017. Estatı́stica básica. Saraiva Educação SA.
R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org.
Xie, Yihui. 2016. bookdown: Authoring Books and Technical Documents with R Markdown. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/bookdown.

  1. ↩︎

  2. ↩︎

  3. ↩︎

  4. ↩︎