L'échantillonnage kurtosis est-il désespérément biaisé?

8

Je regarde l'échantillon kurtosis d'une variable aléatoire assez asymétrique, et les résultats semblent incohérents. Pour illustrer simplement le problème, j'ai regardé l'échantillon kurtosis d'un VR log-normal. En R (que j'apprends lentement):

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

Le résumé que je reçois est

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

Selon Wikipédia , le kurtosis pour ce VR log-normal devrait être d'environ 114. Il est clair que l'échantillon de kurtosis est biaisé.

En faisant des recherches, j'ai découvert que l'échantillonnage kurtosis est biaisé pour les petits échantillons. J'ai utilisé l'estimateur «G2» fourni par le e1071package dans CRAN, et j'ai obtenu des résultats très similaires pour cette taille d'échantillon.

La question : laquelle des caractéristiques suivantes caractérise ce qui se passe:

  1. L'erreur standard de l'échantillon kurtosis est simplement très grande pour ce VR (même si l'estimation commune ondulée à la main de l'erreur standard est de l'ordre ). Alternativement, j'ai utilisé trop peu d'échantillons (2048) dans cette étude.1/n
  2. Ces implémentations de l'échantillonnage kurtosis souffrent de problèmes numériques qui pourraient être corrigés par exemple par la méthode de Terriberry (de la même manière que la méthode de Welford donne de meilleurs résultats que la méthode naïve pour la variance de l'échantillon).
  3. J'ai mal calculé le kurtosis de la population. (Aie)
  4. L'échantillonnage kurtosis est intrinsèquement biaisé et vous ne pouvez jamais le corriger pour des échantillons aussi petits.
shabbychef
la source
BTW, puisque je suis un débutant en R, j'apprécierais tout commentaire concernant mon code, aussi bref soit-il, tant sur les points de style que sur le fond. En particulier, j'espérais qu'il pourrait y avoir une manière plus élégante d'exprimer la boucle for.
shabbychef
1
sur le style R, vous n'avez pas besoin ;à la fin de vos déclarations. Vous avez bien fait de pré-allouer, mais pas besoin de remplir NA, kvals <- numeric(length = n_trial)aurait suffi. Avec rep, vous n'avez besoin que des arguments 1 et 3 de votre appel (par exemple rep(NA, 10)). Dans la forconfiguration de la boucle, 1:n_trialpeut être dangereux si la programmation; est mieux seq_along(kvals)ou seq_len(n_trial)dans ce cas. Enfin, si vous n'avez pas besoin de forcer l'impression, supprimez le print()tour summary()- vous n'en aurez besoin que si vous ne travailliez pas de manière interactive avec R. HTH.
Gavin Simpson
Merci! certainement ce que je cherchais. J'exécutais cela à partir d'un fichier, donc j'avais besoin de print. Les arguments repétaient certainement erronés.
shabbychef

Réponses:

6

Il y a une correction de biais . Ce n'est pas énorme. Je pense que la variance d'échantillonnage de la kurtosis est proportionnelle au huitième (!) Moment central, ce qui peut être énorme pour une distribution log-normale. Vous auriez besoin de millions d'essais (ou bien plus) dans une simulation pour détecter les biais, sauf si le CV est minuscule. (Tracez un histogramme des kvals pour voir à quel point ils sont asymétriques.)

Le kurtosis correct est en effet d'environ 113,9364.

En ce qui concerne le style R, il peut être pratique d'encapsuler la simulation dans une fonction afin que vous puissiez facilement modifier la taille de l'échantillon ou le nombre d'essais.

whuber
la source
2
L'estimateur G2 de e1071donne la correction de biais «standard»; voir cran.r-project.org/web/packages/e1071/e1071.pdf . L'utilisation de cet estimateur au lieu de g2, implémenté par le momentspackage, a eu peu d'effet, comme je l'ai noté dans le Q. La dépendance au huitième moment impliquerait en effet que la taille de l'échantillon est trop petite ici.
shabbychef
5

[Juste sur le style R - @whuber a répondu au Kurtsosis Q]

C'était un peu trop compliqué pour s'en tenir à un commentaire. Pour des boucles aussi simples que celle que vous utilisez, nous pouvons combiner la suggestion de @ whuber d'encapsuler la simulation dans une fonction avec la replicate()fonction. replicate()s'occupe de l'allocation et de l'exécution de la boucle pour vous. Un exemple est donné ci-dessous:

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

Nous l'utilisons comme ceci:

> 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

Notez que j'utilise la rlnorm()fonction pour générer la variable aléatoire log-normale. Il est équivalent à exp(rnorm())dans votre boucle mais utilise le bon outil, et nous permettons à notre fonction de transmettre les paramètres spécifiés par l'utilisateur de la distribution log-normale.

> set.seed(123)
> exp(rnorm(1))
[1] 0.5709374
> set.seed(123)
> rlnorm(1)
[1] 0.5709374
Gavin Simpson
la source
+1 pour set.seed, ce qui aiderait dans des exemples comme celui-ci. Y a-t-il une raison substantielle pour encapsuler dans une fonction ( par exemple, l'interpréteur R précompilera les fonctions, donc il y a un certain accélération), ou est stylistique ( par exemple , les fonctions d'encapsulation, comme le brocoli, sont bonnes pour vous), ou quelque part entre les deux ( par exemple, il y a, par exemple, beaucoup d'opérateurs dans R qui agissent sur les fonctions, donc on devrait s'habituer à la programmation fonctionnelle)?
shabbychef
@shabbychef: Je pense que le principal est l'effort. Vous pouvez exécuter votre code avec la boucle, etc. à plusieurs reprises et ce serait bien, mais vous devez continuer à exécuter toutes les lignes de votre code. En l'encapsulant, vous exécutez 1 ligne de code R pour chaque simulation que vous souhaitez exécuter. Il n'y a aucune accélération car R ne compile rien à l'avance IIRC.
Gavin Simpson
1
Merci pour la clarification. Comme tout est dans un petit fichier, il s'agit quand même d'une seule ligne source('foo.r')
:;