Estou observando a curtose de amostra de uma variável aleatória bastante distorcida e os resultados parecem inconsistentes. Para ilustrar o problema, observei a amostra de curtose de um RV log-normal. Em R (que estou aprendendo lentamente):
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));
O resumo que recebo é
Min. 1st Qu. Median Mean 3rd Qu. Max.
11.87 28.66 39.32 59.17 61.70 1302.00
Segundo a Wikipedia , a curtose para esse RV log-normal deve ser em torno de 114. Claramente a curtose da amostra é tendenciosa.
Ao fazer algumas pesquisas, descobri que a curtose da amostra é tendenciosa para amostras pequenas. Usei o estimador 'G2', conforme fornecido pelo e1071
pacote no CRAN, e obtive resultados muito semelhantes para esse tamanho de amostra.
A questão : qual das seguintes opções caracteriza o que está acontecendo:
- O erro padrão da curtose da amostra é simplesmente muito grande para este RV (mesmo que a estimativa comum do erro padrão da onda manual seja da ordem ). Como alternativa, usei poucas amostras (2048) neste estudo.
- Essas implementações de curtose amostra sofrem de problemas numéricos que podem ser corrigidos pelo ex método de Terriberry (da mesma maneira que o método de Welford dá melhores resultados do que o método ingênuo de variância da amostra).
- Eu calculei a curtose populacional incorretamente. (ai)
- A curtose da amostra é inerentemente tendenciosa e você nunca pode corrigi-la para tamanhos de amostra tão pequenos.
r
unbiased-estimator
kurtosis
shabbychef
fonte
fonte
;
no final de suas declarações. Você fez o direito de pré-alocar, mas não precisa preencherNA
,kvals <- numeric(length = n_trial)
seria suficiente. Comrep
, você só precisa dos argumentos 1 e 3 da sua chamada (por exemplorep(NA, 10)
). Nafor
configuração do loop,1:n_trial
pode ser perigoso se estiver programando; melhor éseq_along(kvals)
ouseq_len(n_trial)
neste caso. Por fim, se você não precisar forçar a impressão, solte aprint()
rodadasummary()
- você só precisaria se não estivesse trabalhando de maneira interativa com a R. HTH.print
. Os argumentosrep
foram certamente errôneos.Respostas:
Há uma correção de viés . Não é enorme. Acredito que a variância amostral da curtose é proporcional ao oitavo (!) Momento central, que pode ser enorme para uma distribuição lognormal. Você precisaria de milhões de tentativas (ou muito mais) em uma simulação para detectar viés, a menos que o CV seja pequeno. (Faça um histograma de kvals para ver como eles são extraordinariamente distorcidos.)
A curtose correta é realmente sobre 113.9364.
No que diz respeito ao estilo R, pode ser conveniente encapsular a simulação em uma função para que você possa modificar facilmente o tamanho da amostra ou o número de tentativas.
fonte
e1071
fornece a correção de viés 'padrão'; consulte cran.r-project.org/web/packages/e1071/e1071.pdf . O uso desse estimador em vez de g2, implementado pelomoments
pacote, teve pouco efeito, como observei no Q. A dependência no oitavo momento implicaria de fato que o tamanho da amostra é muito pequeno aqui.[Apenas no estilo R - @whuber respondeu à Kurtsosis Q]
Isso foi um pouco complicado demais para ser comentado. Para loops tão simples como o que você usa, podemos combinar a sugestão do @ whuber de encapsular a simulação em uma função com a
replicate()
funçãoreplicate()
cuida da alocação e executa o loop para você. Um exemplo é dado abaixo:Usamos assim:
Observe que eu uso a
rlnorm()
função para gerar a variável aleatória log-normal. É equivalente aoexp(rnorm())
seu loop, mas usa a ferramenta correta, e permitimos que nossa função repasse parâmetros especificados pelo usuário da distribuição log-normal.fonte
set.seed
, o que ajudaria em exemplos como este. Existe uma razão substancial para encapsular uma função ( por exemplo, o intérprete R pré-compila funções, portanto, existe alguma aceleração) ou é estilístico ( por exemplo , funções encapsuladoras, como brócolis, são boas para você) ou em algum lugar intermediário ( por exemplo, existem, digamos, muitos operadores em R que atuam em funções; portanto, deve-se acostumar à programação funcional)?source('foo.r')
;)