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

Bueno, si tiene una muestra de una distribución pareto con los parámetros y (donde es el parámetro de límite inferior y es el parámetro de forma), la probabilidad logarítmica de eso muestra es: $X_1, ..., X_n$ $m>0$ $\alpha>0$ $m$ $\alpha$

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 es un aumento monotónico en , por lo que el maximizador es el valor más grande que es consistente con los datos observados. Como el parámetro define el límite inferior del soporte para la distribución de Pareto, lo óptimo es $m$ $m$

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

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

que no depende de . Luego, usando trucos de cálculo ordinarios, el MLE para debe satisfacer $\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$

un poco de álgebra simple nos dice que el MLE de es $\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})}$

En muchos sentidos importantes (por ejemplo, eficiencia asintótica óptima en el sentido de que alcanza 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 de R a continuación calcula el 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: según el comentario de @cardinal e I a continuación, también podemos observar que es el recíproco de la media muestral de los , que suceden a tener una distribución exponencial. Por lo tanto, si tenemos acceso a un software que puede ajustarse a una distribución exponencial (lo que es más probable, ya que parece surgir en muchos problemas estadísticos), entonces se puede ajustar una distribución de Pareto transformando el conjunto de datos de esta manera y ajustándolo a una distribución exponencial en la escala transformada. $\hat{\alpha}$ $\log(X_i /\hat{m})$

Macro
fuente

(+1) Podemos escribir cosas un poco más sugestivas al notar que se distribuye exponencialmente con rate . De esto y de la invariancia de los MLEs bajo transformación, concluimos de inmediato que , donde reemplazamos por en la última expresión. Esto también sugiere cómo podríamos usar un software estándar para ajustar un Pareto, incluso si no hay una opción explícita disponible.

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$

cardenal

@cardinal - Entonces, es el recíproco de la media muestral de los , que tienen una distribución exponencial. ¿Cómo nos ayuda esto?

\hat{α}

$\hat{\alpha}$

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

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

Macro

Hola macro El punto que intentaba señalar era que el problema de estimar los parámetros de un Pareto puede reducirse (esencialmente) al de la estimación de la tasa de un exponencial: a través de la transformación anterior, podemos convertir nuestros datos y problemas en un (quizás) más familiar e inmediatamente extraiga la respuesta (suponiendo que nosotros, o nuestro software, ya sepamos qué hacer con una muestra de exponenciales).

cardenal

¿Cómo puedo medir el error de este tipo de ajuste?

emanuele

@emanuele, la varianza aproximada de un MLE es la inversa de la matriz de información del pescador, que requerirá que calcules al menos una derivada de la probabilidad logarítmica. O bien, podría usar una especie de remuestreo de arranque para estimar el error estándar.

Macro

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

Respuestas: