> 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