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 e1071
package 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:
- 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.
- 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).
- J'ai mal calculé le kurtosis de la population. (Aie)
- L'échantillonnage kurtosis est intrinsèquement biaisé et vous ne pouvez jamais le corriger pour des échantillons aussi petits.
r
unbiased-estimator
kurtosis
shabbychef
la source
la source
;
à la fin de vos déclarations. Vous avez bien fait de pré-allouer, mais pas besoin de remplirNA
,kvals <- numeric(length = n_trial)
aurait suffi. Avecrep
, vous n'avez besoin que des arguments 1 et 3 de votre appel (par exemplerep(NA, 10)
). Dans lafor
configuration de la boucle,1:n_trial
peut être dangereux si la programmation; est mieuxseq_along(kvals)
ouseq_len(n_trial)
dans ce cas. Enfin, si vous n'avez pas besoin de forcer l'impression, supprimez leprint()
toursummary()
- vous n'en aurez besoin que si vous ne travailliez pas de manière interactive avec R. HTH.print
. Les argumentsrep
étaient certainement erronés.Réponses:
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.
la source
e1071
donne 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 lemoments
package, 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.[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:Nous l'utilisons comme ceci:
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.la source
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)?source('foo.r')