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

10

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.

stjudent
quelle
1
Das Ergebnis eines statistischen Tests sagt Ihnen nicht, dass Ihre Daten eine Pareto-Verteilung haben . Tatsächlich können Sie ziemlich sicher sein, dass echte Daten keine Pareto-Verteilung haben. Ein Test zeigt Ihnen lediglich, ob Sie über genügend Daten verfügen, um die Abweichung von Pareto zu ermitteln. Das heißt, wenn es alles ablehnt, was es sagt, ist "Ja, die Stichprobengröße war groß genug, um Ihnen zu sagen, was Sie bereits wussten". Warum sollten Sie eine solche Übung durchführen, die Ihre eigentliche Frage nicht beantworten kann?
Glen_b -State Monica
Ist Ihre Frage wirklich nichts anderes als "Welche Codezeilen schreibe ich, damit Programm R Prozedur X ausführt"? Dann ist es hier kein Thema. Es könnte sich um eine Programmierfrage handeln. Wenn Ihre Frage einen statistischen Aspekt hat (z. B. "Ist dies sinnvoll?"),
Sollten Sie
1
Nun zum Anderson-Darling-Test (oder zu dem Kolmogorov-Smirnov, den @Zen oben vorgeschlagen hat). Dies sind Tests für vollständig spezifizierte Verteilungen. Das heißt, damit die Tests die gewünschten Eigenschaften haben, müssen Sie alle Parameter a priori ( NICHT schätzen ) angeben . Daher können Sie für diese Übung keine von beiden verwenden, da Sie keine festgelegten Parameter haben. (Vermutlich tun Sie dies auf Vorschlag einer anderen Person. Es ist sehr schwierig, jemandem über einen Vermittler Missverständnisse zu erklären.)
Glen_b - Monica
Wofür machst du diese Tests? zB welche Aktionen ändern sich je nachdem, ob Sie ablehnen oder nicht ablehnen?
Glen_b -State Monica
Sie sollten sich immer ein QQ-Diagramm ansehen, unabhängig von Ihrem Motiv. Und Sie sollten keinen "exakten" P-Wert fetischisieren. Ein anderer Test würde Ihnen einen anderen "exakten" P-Wert geben.
Nick Cox

Antworten:

12

(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 .pks.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 :χ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 .

Elvis
quelle
3
Wenn Sie die Verteilung der Teststatistik für die Schätzung von Parametern auf diese Weise anpassen, handelt es sich nicht um einen KS-Test (obwohl er auf einer KS-Statistik basiert), sondern um eine bestimmte Art von Lilliefors-Test . Dies ist nicht mehr nichtparametrisch, aber man kann eine durch Simulation für jede gegebene Verteilung konstruieren. Lilliefors tat dies speziell für das Normale und Exponentielle ... in den 1960er Jahren.
Glen_b -State Monica
Danke für diesen Kommentar @Glen_b Das wusste ich nicht.
Elvis
Kein Problem; Es ändert nichts am Inhalt Ihrer Arbeit (was in Ordnung ist), nur an dem, was es heißen sollte.
Glen_b -Rate State Monica
@Glen_b Ich habe einige wesentliche Änderungen in meiner Antwort vorgenommen, nochmals vielen Dank!
Elvis