Como saber se meus dados se encaixam na distribuição Pareto?

10

Eu tenho uma amostra que é um vetor com 220 números. Aqui está um link para um histograma dos meus dados. . E desejo verificar se meus dados se encaixam em uma distribuição de Pareto, mas não quero ver gráficos de QQ com essa distribuição, mas preciso de uma resposta exata com valor-p em R, como o teste de Anderson-Darling para normalidade ( ad.test) . Como eu pude fazer isso? Por favor, seja o mais específico possível.

stjudent
fonte
1
O resultado de um teste estatístico não informa que seus dados têm uma distribuição Pareto . De fato, você pode ter certeza de que, se forem dados reais, eles não terão uma distribuição Pareto. Tudo o que um teste mostra é se você possui dados suficientes para captar o desvio de valor de Pareto que você possui. Ou seja, se rejeitar tudo o que diz é "sim, o tamanho da amostra foi grande o suficiente para lhe dizer o que você já sabia". Por que você faria esse exercício, que não pode responder à pergunta real que você tem?
Glen_b -Reinstate Monica
Sua pergunta é realmente nada além de 'que linhas de código eu escrevo para fazer o programa R executar o procedimento X'? Então está fora de tópico aqui. Pode ser qualificado como uma questão de programação. Se há um aspecto estatístico para sua pergunta (como '? Isso faz algum sentido fazer'), então você deve esclarecer e enfatizar os aspectos dela
Glen_b -Reinstate Monica
1
Agora, para o teste de Anderson-Darling (ou, nesse caso, o Kolmogorov-Smirnov que o @Zen sugeriu acima). Esses são testes para distribuições completamente especificadas . Ou seja, para que os testes tenham as propriedades desejadas, você deve especificar a priori ( NÃO estimar ) todos os parâmetros. Portanto, você não pode usar nenhum deles para este exercício porque não possui parâmetros pré-especificados. (Presumivelmente, você está fazendo isso por sugestão de outra pessoa. É muito difícil explicar conceitos errôneos a alguém por meio de um intermediário.)
Glen_b -Reinstala Monica 30/11/13
Para que você está fazendo esse teste? por exemplo, que ações serão alteradas dependendo de você rejeitar ou deixar de rejeitar?
Glen_b -Reinstate Monica
Você deve sempre olhar para um gráfico de QQ, independentemente do seu motivo. E você não deve fetichizar um valor P "exato". Um teste diferente forneceria um valor P "exato" diferente.
Nick Cox

Respostas:

12

(PS) Antes de tudo, acho que Glen_b está certo nos comentários acima sobre a utilidade de um teste: dados reais certamente não são exatamente distribuídos por Pareto, e para aplicações mais práticas a pergunta seria "quão boa é a aproximação de Pareto?" - e o gráfico QQ é uma boa maneira de mostrar a qualidade dessa aproximação.

De qualquer forma, você pode fazer seu teste com a estatística Kolmogorov-Smirnov, depois de estimar os parâmetros por máxima verossimilhança. Essa estimativa de parâmetro impede o uso do valor- de , portanto, você pode executar o bootstrap paramétrico para estimar. Como Glen_b diz no comentário, isso pode ser conectado ao teste de Lilliefors .pks.test

Aqui estão algumas linhas de código R.

Primeiro defina as funções básicas para lidar com as distribuições de Pareto.

# distribution, cdf, quantile and random functions for Pareto distributions
dpareto <- function(x, xm, alpha) ifelse(x > xm , alpha*xm**alpha/(x**(alpha+1)), 0)
ppareto <- function(q, xm, alpha) ifelse(q > xm , 1 - (xm/q)**alpha, 0 )
qpareto <- function(p, xm, alpha) ifelse(p < 0 | p > 1, NaN, xm*(1-p)**(-1/alpha))
rpareto <- function(n, xm, alpha) qpareto(runif(n), xm, alpha)

A função a seguir calcula o MLE dos parâmetros (justificativas na Wikipedia ).

pareto.mle <- function(x)
{
  xm <- min(x)
  alpha <- length(x)/(sum(log(x))-length(x)*log(xm))
  return( list(xm = xm, alpha = alpha))
}

p

pareto.test <- function(x, B = 1e3)
{
  a <- pareto.mle(x)

  # KS statistic
  D <- ks.test(x, function(q) ppareto(q, a$xm, a$alpha))$statistic

  # estimating p value with parametric bootstrap
  B <- 1e5
  n <- length(x)
  emp.D <- numeric(B)
  for(b in 1:B)
  {
    xx <- rpareto(n, a$xm, a$alpha);
    aa <- pareto.mle(xx)
    emp.D[b] <- ks.test(xx, function(q) ppareto(q, aa$xm, aa$alpha))$statistic
  }

  return(list(xm = a$xm, alpha = a$alpha, D = D, p = sum(emp.D > D)/B))
}

Agora, por exemplo, uma amostra proveniente de uma distribuição de Pareto:

> # generating 100 values from Pareto distribution
> x <- rpareto(100, 0.5, 2)
> pareto.test(x)
$xm
[1] 0.5007593

$alpha
[1] 2.080203

$D
         D 
0.06020594 

$p
[1] 0.69787

χ2(2)

> # generating 100 values from chi square distribution
> x <- rchisq(100, df=2)
> pareto.test(x)
$xm
[1] 0.01015107

$alpha
[1] 0.2116619

$D
        D 
0.4002694 

$p
[1] 0

Note que não afirmo que este teste seja imparcial: quando a amostra é pequena, pode existir algum viés. O bootstrap paramétrico não leva bem em conta a incerteza na estimativa do parâmetro (pense no que aconteceria ao usar essa estratégia para testar ingenuamente se a média de alguma variável normal com variação desconhecida é zero).

PS Wikipedia diz algumas palavras sobre isso. Aqui estão duas outras perguntas para as quais uma estratégia semelhante foi sugerida: Teste de qualidade do ajuste para uma mistura , teste de qualidade do ajuste para uma distribuição gama .

Elvis
fonte
3
Quando você ajusta a distribuição da estatística de teste para a estimativa de parâmetros dessa maneira, não é um teste KS (embora baseado em uma estatística KS) - é um tipo específico de teste de Lilliefors . Isso não é mais não paramétrico, mas é possível construir um via simulação para qualquer distribuição. Lilliefors fez isso especificamente para o normal e exponencial ... na década de 1960.
Glen_b -Reinstate Monica
Obrigado por este comentário @Glen_b Eu não sabia disso.
Elvis
Sem problemas; não muda nada sobre o conteúdo do que você está fazendo (o que é bom), apenas o que deveria ser chamado.
Glen_b -Reinstala Monica
@Glen_b Fiz algumas alterações substanciais na minha resposta, obrigado novamente!
Elvis