Como ajusto um conjunto de dados a uma distribuição de Pareto no R?

22

Tenha, digamos, os seguintes dados:

8232302  684531  116857   89724   82267   75988   63871   
  23718    1696     436     439     248     235

Deseja uma maneira simples de ajustar isso (e vários outros conjuntos de dados) a uma distribuição de Pareto. Idealmente, produziria os valores teóricos correspondentes, menos idealmente os parâmetros.

Felix
fonte
O que se pretende "combinando valores teóricos"? As expectativas das estatísticas de pedidos, dadas as estimativas de parâmetros? Ou alguma outra coisa?
Glen_b -Reinstala Monica

Respostas:

33

Bem, se você tem uma amostra X1,...,Xn de uma distribuição de pareto com os parâmetros m>0 e α>0 (onde m é o parâmetro do limite inferior e α é o parâmetro de forma), a probabilidade logarítmica dessa amostra é:

nlog(α)+nαlog(m)(α+1)i=1nlog(Xi)

este é um aumento monotônico em , portanto, o maximizador é o maior valor que é consistente com os dados observados. Como o parâmetro m define o limite inferior do suporte para a distribuição de Pareto, o ideal émm

m^=miniXi

o que não depende de . Em seguida, usando truques comuns de cálculo, o MLE para deve satisfazerααα

nα+nlog(m^)i=1nlog(Xi)=0

alguma álgebra simples nos diz que o MLE de éα

α^=ni=1nlog(Xi/m^)

Em muitos sentidos importantes (por exemplo, eficiência assintótica ideal, pois atinge o limite inferior de Cramer-Rao), esta é a melhor maneira de ajustar dados a uma distribuição de Pareto. O código R abaixo calcula o MLE para um determinado conjunto de dados X,.

pareto.MLE <- function(X)
{
   n <- length(X)
   m <- min(X)
   a <- n/sum(log(X)-log(m))
   return( c(m,a) ) 
}

# example. 
library(VGAM)
set.seed(1)
z = rpareto(1000, 1, 5) 
pareto.MLE(z)
[1] 1.000014 5.065213

Edit: Com base nos comentários de @cardinal e I abaixo, também podemos observar que é o inverso da média da amostra do 's, que acontece com tem uma distribuição exponencial. Portanto, se tivermos acesso a software que possa se ajustar a uma distribuição exponencial (o que é mais provável, pois parece surgir em muitos problemas estatísticos), o ajuste de uma distribuição de Pareto poderá ser realizado transformando os dados dessa maneira e ajustando-os a uma distribuição exponencial na escala transformada. log(Xi/ m )α^log(Xi/m^)

Macro
fonte
3
(+1) Podemos escrever as coisas um pouco mais sugestivamente, observando que é distribuído exponencialmente com a taxa . disso e pela invariância dos MLEs em transformação, concluímos imediatamente que , onde substituímos por na última expressão. Isso também sugere como podemos usar o software padrão para ajustar um Pareto, mesmo que nenhuma opção explícita esteja disponível. ct ct = 1 / ˉ Y m mYi=log(Xi/m)αα^=1/Y¯mm^
cardeal
@ cardinal - Então, é o inverso da média da amostra dos 's, que possuem uma distribuição exponencial. Como isso nos ajuda? log(Xi/ m )α^log(Xi/m^)
Macro
2
Olá, Macro. O que eu estava tentando enfatizar era que o problema de estimar os parâmetros de um Pareto pode ser (essencialmente) reduzido ao da estimativa da taxa de um exponencial: Através da transformação acima, podemos converter nossos dados e problemas em um (talvez) mais familiar e extraia imediatamente a resposta (assumindo que nós, ou nosso software, já sabemos o que fazer com uma amostra de exponenciais).
cardeal
Como posso medir o erro desse tipo de ajuste?
Emanuele
@emanuele, a variação aproximada de um MLE é o inverso da matriz de informações de Fisher, que exigirá que você calcule pelo menos uma derivada da probabilidade de log. Ou você pode usar um tipo de reamostragem de autoinicialização para estimar o erro padrão.
Macro
8

Você pode usar a fitdistfunção fornecida no fitdistrpluspacote:

library(MASS)
library(fitdistrplus)
library(actuar)

# suppose data is in dataPar list
fp <- fitdist(dataPar, "pareto", start=list(shape = 1, scale = 500))
#the mle parameters will be stored in fp$estimate
akashrajkn
fonte
Deveria ser isso library(fitdistrplus)?
Sean
1
Sim @Sean, edição de resposta em conformidade
Kevin L chaves
Observe que a chamada para library(actuar)é necessária para que isso funcione.
jsta
O que fp $ estima ["shape"] representa neste caso? Talvez seja o alfa estimado? Ou beta?
Albert Hendriks