Woher wissen Sie, ob meine Daten zur Pareto-Verteilung passen?

Ich habe eine Probe, die ein Vektor mit 220 Zahlen ist. Hier ist ein Link zu einem Histogramm meiner Daten. . Und ich möchte überprüfen, ob meine Daten zu einer Pareto-Verteilung passen, aber ich möchte keine QQ-Diagramme mit dieser Verteilung sehen, aber ich brauche eine genaue Antwort mit p-Wert in R, wie zum Beispiel den Anderson-Darling-Test für Normalität ( ad.test) . Wie könnte ich das machen? Bitte seien Sie so genau wie möglich.

(PS) Zunächst denke ich, dass Glen_b in seinen obigen Kommentaren zur Nützlichkeit eines solchen Tests Recht hat: Reale Daten sind sicherlich nicht genau Pareto-verteilt, und für die meisten praktischen Anwendungen wäre die Frage: "Wie gut ist die Pareto-Näherung?" - und das QQ-Diagramm ist ein guter Weg, um die Qualität einer solchen Annäherung zu zeigen.

Wie auch immer, Sie können Ihren Test mit der Kolmogorov-Smirnov-Statistik durchführen, nachdem Sie die Parameter mit maximaler Wahrscheinlichkeit geschätzt haben. Diese Parameterschätzung verhindert, dass der Wert von verwendet wird , sodass Sie einen parametrischen Bootstrap durchführen können, um ihn zu schätzen. Wie Glen_b im Kommentar mitteilt, kann dies mit dem Lilliefors-Test verbunden werden .$p$ks.test

Hier sind einige Zeilen R-Code.

Definieren Sie zunächst die Grundfunktionen für Pareto-Verteilungen.

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

Die folgende Funktion berechnet die MLE der Parameter (Begründungen in Wikipedia ).

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

Diese Funktionen berechnen die KS-Statistik und verwenden den parametrischen Bootstrap, um den Wert zu schätzen .$p$

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

Nun zum Beispiel eine Stichprobe aus einer Pareto-Distribution:

> # 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

... und aus a :$\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

Beachten Sie, dass ich nicht behaupte, dass dieser Test unvoreingenommen ist: Wenn die Stichprobe klein ist, kann eine gewisse Verzerrung bestehen. Der parametrische Bootstrap berücksichtigt die Unsicherheit bei der Parameterschätzung nicht gut (überlegen Sie, was passieren würde, wenn Sie mit dieser Strategie naiv testen, ob der Mittelwert einer normalen Variablen mit unbekannter Varianz Null ist).

PS Wikipedia sagt ein paar Worte dazu. Hier sind zwei weitere Fragen, für die eine ähnliche Strategie vorgeschlagen wurde: Anpassungstest für eine Mischung , Anpassungstest für eine Gammaverteilung .