¿Cómo saber si mis datos se ajustan a la distribución de Pareto?

Question

Tengo una muestra que es un vector con 220 números. Aquí hay un enlace a un histograma de mis datos. . Y deseo comprobar si mis datos se ajustan a una distribución de Pareto, pero no quiero ver gráficos QQ con esa distribución, sino que necesito una respuesta exacta con valor p en R, como la prueba de Anderson-Darling para la normalidad ( ad.test ). ¿Cómo podría hacerlo? Por favor, sea tan específico como pueda.

Answer 1

(PS) En primer lugar, creo que Glen_b tiene razón en sus comentarios anteriores sobre la utilidad de una prueba de este tipo: los datos reales seguramente no son exactamente distribuida de Pareto, y para la mayoría de las aplicaciones prácticas la pregunta sería "¿qué tan buena es la aproximación de Pareto?" - y el gráfico QQ es una buena manera de mostrar la calidad de dicha aproximación.

De todos modos, puede hacer su prueba con el estadístico de Kolmogorov-Smirnov, después de estimar los parámetros por máxima verosimilitud. Esta estimación de los parámetros impide utilizar el $p$ -valor de ks.test , por lo que se puede hacer bootstrap paramétrico para estimarlo. Como dice Glen_b en el comentario, esto se puede conectar a Prueba de Lilliefors .

Aquí hay algunas líneas de código R.

En primer lugar, defina las funciones básicas para tratar las distribuciones 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)

La siguiente función calcula el MLE de los parámetros (justificaciones en Wikipedia ).

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

Y esta función calcula el estadístico KS, y utiliza el bootstrap paramétrico para estimar la $p$ -valor.

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))
}

Ahora, por ejemplo, una muestra procedente de una distribución 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

...y de un $\chi^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

Tenga en cuenta que no pretendo que esta prueba sea insesgada: cuando la muestra es pequeña, puede existir cierto sesgo. El bootstrap paramétrico no tiene bien en cuenta la incertidumbre en la estimación de los parámetros (piense en lo que ocurriría al utilizar esta estrategia para probar ingenuamente si la media de alguna variable normal con varianza desconocida es cero).

PS Wikipedia dice unas palabras sobre esto. Aquí hay otras dos preguntas para las que se sugirió una estrategia similar: Prueba de bondad de ajuste para una mezcla , prueba de bondad de ajuste para una distribución gamma .