# Is sample kurtosis hopelessly biased?

I am looking at the sample kurtosis of a fairly skewed random variable, and the results seem inconsistent. To simply illustrate the problem, I looked at the sample kurtosis of a log-normal RV. In R (which I am slowly learning):

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


The summary I get is

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


According to Wikipedia, the kurtosis for this log-normal RV should be around 114. Clearly the sample kurtosis is biased.

Doing some research I found that sample kurtosis is biased for small sample sizes. I used the 'G2' estimator as provided by the e1071 package in CRAN, and got very similar results for this sample size.

The question: which of the following characterize what is going on:

1. The standard error of the sample kurtosis is simply very large for this RV (even though the hand-wavey common estimate of the standard error is of order $1/\sqrt{n}$). Alternatively, I used too few samples (2048) in this study.
2. These implementations of sample kurtosis suffer from numerical problems which might be corrected by e.g. Terriberry's method (in much the same way that Welford's method gives better results than the naive method for sample variance).
3. I computed the population kurtosis incorrectly. (ouch)
4. Sample kurtosis is inherently biased and you can never fix it for such small sample sizes.
-
BTW, since I am a beginner in R, I would appreciate any comments regarding my code, however brief, both on points of style and substance. In particular, I was hoping there might be a more elegant way to phrase the for loop. –  shabbychef Dec 8 '10 at 3:05
on the R style, you don't need ; on the ends of your statements. You did right to pre-allocate, but no need to fill with NA, kvals <- numeric(length = n_trial) would have sufficed. With rep, you only need arguments 1 and 3 from your call (e.g. rep(NA, 10)). In the for loop set-up, 1:n_trial can be dangerous if programming; better is seq_along(kvals) or seq_len(n_trial) in this case. Finally, if you don't need to force printing, drop the print() round summary() - you'd only need it if you weren't working interactively with R. HTH. –  Gavin Simpson Dec 8 '10 at 8:28
thanks! definitely what I was looking for. I was running this from a file, so I needed the print. The args to rep were certainly erroneous. –  shabbychef Dec 8 '10 at 17:46

There's a bias correction. It's not huge. I believe the sampling variance of the kurtosis is proportional to the eighth (!) central moment, which can be enormous for a lognormal distribution. You would need millions of trials (or far more) in a simulation to detect bias unless the CV is tiny. (Plot a histogram of kvals to see how extraordinarily skewed they are.)

The correct kurtosis is indeed about 113.9364.

As far as R style goes, it can be convenient to encapsulate the simulation in a function so you can easily modify the sample size or number of trials.

-
The G2 estimator from e1071 gives the 'standard' bias correction; see cran.r-project.org/web/packages/e1071/e1071.pdf . Using this estimator instead of g2, implemented by the moments package, had little effect, as I noted in the Q. The dependence on the eighth moment would indeed imply the sample size is too small here. –  shabbychef Dec 8 '10 at 17:55

[Just on the R Style - @whuber has answered the Kurtsosis Q]

This was a bit too complicated to stick into a comment. For such simple loops like the one you use, we can combine @whuber's suggestion of encapsulating the simulation in a function with the replicate() function. replicate() takes care of allocation and running the loop for you. An example is given below:

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


We use it like this:

> 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


Note that I use the rlnorm() function to generate the log-normal random variable. It is equivalent to exp(rnorm()) in your loop but uses the correct tool, and we allow our function to pass on user-specified parameters of the log-normal distribution.

> set.seed(123)
> exp(rnorm(1))
[1] 0.5709374
> set.seed(123)
> rlnorm(1)
[1] 0.5709374

-
+1 Nice example. –  whuber Dec 8 '10 at 15:35
+1 for set.seed, which would help in examples like this. Is there a substantial reason for encapsulating in a function (e.g. the R interpreter will pre-compile functions, thus there is some speedup), or is stylistic (e.g. encapsulating functions, like broccoli, is good for you), or somewhere in between (e.g. there are, say, lots of operators in R which act on functions, so one should get used to functional programming)? –  shabbychef Dec 8 '10 at 17:49
@shabbychef: I think the main one is effort. You could run your code with the loop etc repeatedly and it would be fine, but you have to keep running all the lines of your code. By encapsulating it, you run 1 line of R code for each simulation you want to run. There is no speed up as R doesn't compile anything ahead of time IIRC. –  Gavin Simpson Dec 8 '10 at 18:07
thanks for the clarification. Since this is all in a small file, it is one line anyway: source('foo.r') ;) –  shabbychef Dec 8 '10 at 18:10