R Code to Simulate the Sampling Distribution of the Coefficient of Excess for Exponential Data

> excess
function (x) 
{
    mean((x - mean(x))^4)/(mean((x - mean(x))^2)^2) - 3
}
 
> rexpexcess
function (n, B, ...) 
{
    ex <- apply(matrix(rexp(n * B), nrow = B), 1, excess)
    hist(ex, xlab = paste("Coefficient of excess with n =", n), 
        ...)
    invisible(ex)
}
 
> mean(rexpexcess(50,1000))
[1] 3.168413