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

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

n \log (α) + n α \log (m) - (α + 1) \sum_{i = 1}^{n} \log (X_{i})

$n \log(\alpha) + n \alpha \log(m) - (\alpha+1) \sum_{i=1}^{n} \log(X_i)$

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

\hat{m} = min_{i} X_{i}

$\hat{m} = \min_{i} X_i$

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

\frac{n}{α} + n \log (\hat{m}) - \sum_{i = 1}^{n} \log (X_{i}) = 0

$\frac{n}{\alpha} + n \log( \hat{m} ) - \sum_{i=1}^{n} \log(X_i) = 0$

alguma álgebra simples nos diz que o MLE de é $\alpha$

\hat{α} = \frac{n}{\sum_{i = 1}^{n} \log (X_{i} / \hat{m})}

$\hat{\alpha} = \frac{n}{\sum_{i=1}^{n} \log(X_i/\hat{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. $\hat{\alpha}$ $\log(X_i /\hat{m})$

Macro
fonte

(+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.

Y_{i} = \log (X_{i} / m)

$Y_i = \log(X_i/m)$

α

$\alpha$

\hat{α} = 1 / \bar{Y}

$\hat\alpha = 1/\bar Y$

m

$m$

\hat{m}

$\hat m$

cardeal

@ cardinal - Então, é o inverso da média da amostra dos 's, que possuem uma distribuição exponencial. Como isso nos ajuda?

\hat{α}

$\hat{\alpha}$

\log (X_{i} / \hat{m})

$\log(X_i/\hat{m})$

Macro

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

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

Respostas: