Ist die Probenkurtose hoffnungslos voreingenommen?

8

Ich betrachte die Stichproben-Kurtosis einer ziemlich verzerrten Zufallsvariablen, und die Ergebnisse scheinen inkonsistent zu sein. Um das Problem einfach zu veranschaulichen, habe ich mir die Beispielkurtose eines logarithmisch normalen Wohnmobils angesehen. In R (was ich langsam lerne):

library(moments); 

samp_size = 2048;
n_trial = 4096;

kvals <- rep(NA,1,n_trial); #preallocate
for (iii in 1:n_trial) {
    kvals[iii] <- kurtosis(exp(rnorm(samp_size)));
}
print(summary(kvals));

Die Zusammenfassung, die ich bekomme, ist

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11.87   28.66   39.32   59.17   61.70 1302.00 

Laut Wikipedia sollte die Kurtosis für dieses logarithmisch normale Wohnmobil bei 114 liegen. Die Kurtosis der Probe ist eindeutig voreingenommen.

Bei einigen Nachforschungen stellte ich fest, dass die Probenkurtose bei kleinen Probengrößen voreingenommen ist. Ich habe den 'G2'-Schätzer verwendet, wie er im e1071Paket in CRAN bereitgestellt wird , und habe für diese Stichprobengröße sehr ähnliche Ergebnisse erhalten.

Die Frage : Welche der folgenden Aussagen charakterisieren, was los ist:

  1. 1/n
  2. Diese Implementierungen von Probe Kurtosis leiden unter numerischen Problemen , die durch korrigiert werden könnten , zB Terriberry Methode (in die gleichen Weise , dass Welford Methode bessere Ergebnisse als die naive Methode zur Stichprobenvarianz gibt).
  3. Ich habe die Populationskurtose falsch berechnet. (Autsch)
  4. Die Probenkurtose ist von Natur aus voreingenommen und kann bei so kleinen Probengrößen niemals behoben werden.
shabbychef
quelle
Übrigens, da ich ein Anfänger in R bin, würde ich mich über Kommentare zu meinem Code freuen, auch wenn diese kurz sind, sowohl in Bezug auf Stil als auch in Bezug auf Inhalt. Insbesondere hatte ich gehofft, dass es eine elegantere Möglichkeit geben könnte, die for-Schleife zu formulieren.
Shabbychef
1
Auf dem R-Stil brauchen Sie nicht ;an den Enden Ihrer Aussagen. Sie haben richtig vorab zugewiesen, aber keine Notwendigkeit zu füllen NA, kvals <- numeric(length = n_trial)hätte genügt. Mit repbenötigen Sie nur die Argumente 1 und 3 aus Ihrem Aufruf (zB rep(NA, 10)). In der forSchleifeneinstellung 1:n_trialkann beim Programmieren gefährlich sein; besser ist seq_along(kvals)oder seq_len(n_trial)in diesem Fall. Wenn Sie das Drucken nicht erzwingen müssen, lassen Sie die print()Runde fallen summary()- Sie brauchen sie nur, wenn Sie nicht interaktiv mit R. HTH arbeiten.
Gavin Simpson
Vielen Dank! definitiv das, wonach ich gesucht habe. Ich habe dies aus einer Datei ausgeführt, also brauchte ich die print. Die Argumente repwaren sicherlich falsch.
Shabbychef

Antworten:

6

Es gibt eine Bias-Korrektur . Es ist nicht riesig. Ich glaube, dass die Stichprobenvarianz der Kurtosis proportional zum achten (!) Zentralmoment ist, was für eine logarithmische Normalverteilung enorm sein kann. Sie würden Millionen von Versuchen (oder weit mehr) in einer Simulation benötigen, um Verzerrungen zu erkennen, es sei denn, der Lebenslauf ist winzig. (Zeichnen Sie ein Histogramm von kvals, um zu sehen, wie außergewöhnlich schief sie sind.)

Die richtige Kurtosis liegt in der Tat bei 113.9364.

In Bezug auf den R-Stil kann es praktisch sein, die Simulation in eine Funktion zu kapseln, damit Sie die Stichprobengröße oder die Anzahl der Versuche leicht ändern können.

whuber
quelle
2
Der G2-Schätzer von e1071gibt die 'Standard'-Vorspannungskorrektur an; Siehe cran.r-project.org/web/packages/e1071/e1071.pdf . Die Verwendung dieses Schätzers anstelle von g2, der vom momentsPaket implementiert wurde , hatte wenig Wirkung, wie ich im Q festgestellt habe. Die Abhängigkeit vom achten Moment würde tatsächlich bedeuten, dass die Stichprobengröße hier zu klein ist.
Shabbychef
5

[Nur auf dem R-Stil - @whuber hat die Kurtsosis Q beantwortet]

Dies war etwas zu kompliziert, um in einen Kommentar einzusteigen. Für solche einfachen Schleifen wie die von Ihnen verwendete können wir den Vorschlag von @ whuber kombinieren, die Simulation in einer Funktion mit der replicate()Funktion zu kapseln . replicate()kümmert sich um die Zuordnung und führt die Schleife für Sie aus. Ein Beispiel ist unten angegeben:

require(moments)
foo <- function(size, trials, meanlog = 0, sdlog = 1) {
    replicate(trials,
              kurtosis(rlnorm(size, meanlog = meanlog, 
                              sdlog = sdlog)))
}

Wir benutzen es so:

> set.seed(1)
> out <- foo(2048, 10000)
> summary(out)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10.93   28.77   39.99   62.53   62.58 1557.00

Beachten Sie, dass ich die rlnorm()Funktion verwende, um die logarithmische normale Zufallsvariable zu generieren. Es entspricht exp(rnorm())in Ihrer Schleife, verwendet jedoch das richtige Werkzeug, und wir erlauben unserer Funktion, benutzerdefinierte Parameter der Protokollnormalverteilung weiterzugeben.

> set.seed(123)
> exp(rnorm(1))
[1] 0.5709374
> set.seed(123)
> rlnorm(1)
[1] 0.5709374
Gavin Simpson
quelle
+1 für set.seed, was in solchen Beispielen helfen würde. Gibt es einen wesentlichen Grund für die Kapselung in einer Funktion ( z. B. kompiliert der R-Interpreter Funktionen vor, daher gibt es eine gewisse Beschleunigung) oder ist er stilistisch ( z. B. ist die Kapselung von Funktionen wie Brokkoli gut für Sie) oder irgendwo dazwischen ( zB gibt es zum Beispiel viele Operatoren in R, die auf Funktionen einwirken, also sollte man sich an die funktionale Programmierung gewöhnen)?
Shabbychef
@ Shabbychef: Ich denke, die wichtigste ist die Anstrengung. Sie könnten Ihren Code wiederholt mit der Schleife usw. ausführen, und es wäre in Ordnung, aber Sie müssen weiterhin alle Zeilen Ihres Codes ausführen. Durch die Kapselung führen Sie für jede Simulation, die Sie ausführen möchten, eine Zeile R-Code aus. Es gibt keine Beschleunigung, da R vor IIRC nichts kompiliert.
Gavin Simpson
1
Danke für die Klarstellung. Da dies alles in einer kleinen Datei ist, ist es sowieso eine Zeile : source('foo.r');)
shabbychef