¿Cómo puedo ajustar un conjunto de datos a una distribución de Pareto en R?

Question

Tenga, digamos, los siguientes datos:

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

Quiero una forma sencilla de ajustar esto (y varios otros conjuntos de datos) a una distribución de Pareto. Lo ideal sería que los valores teóricos de salida, menos idealmente los parámetros.

Answer 1

Bueno, si tienes una muestra $X_1, ..., X_n$ a partir de una distribución de pareto con parámetros $m>0$ y $\alpha>0$ (donde $m$ es el parámetro de límite inferior y $\alpha$ es el parámetro de forma) la log-verosimilitud de esa muestra es:

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

esto es un aumento monotónico en $m$ por lo que el maximizador es el mayor valor que es consistente con los datos observados. Dado que el parámetro $m$ define el límite inferior del soporte de la distribución de Pareto, el óptimo es

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

que no depende de $\alpha$ . A continuación, utilizando trucos de cálculo ordinario, la MLE para $\alpha$ debe satisfacer

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

un poco de álgebra simple nos dice que la MLE de $\alpha$ es

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

En muchos sentidos importantes (por ejemplo, la eficiencia asintótica óptima en el sentido de que logra el límite inferior de Cramer-Rao), esta es la mejor manera de ajustar los datos a una distribución de Pareto. El código R siguiente calcula la MLE para un conjunto de datos dado, 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

Editar: Basándonos en el comentario de @cardinal y en el mío que aparece a continuación, también podemos observar que $\hat{\alpha}$ es el recíproco de la media muestral del $\log(X_i /\hat{m})$ que resultan tener una distribución exponencial. Por lo tanto, si tenemos acceso a un software que pueda ajustar una distribución exponencial (lo que es más probable, ya que parece surgir en muchos problemas estadísticos), entonces el ajuste de una distribución de Pareto se puede lograr transformando el conjunto de datos de esta manera y ajustándolo a una distribución exponencial en la escala transformada.

Answer 2

Puede utilizar el fitdist proporcionada en fitdistrplus paquete:

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