The function skew(x) computes the skewness coefficient g1 for an array of observations x.
> skew function(x) { mean((x - mean(x))^3)/((mean((x - mean(x))^2))^1.5) }
The function excess(x) computes the coefficient of excess g2 for an array of observations x.
> excess function(x) { (mean((x - mean(x))^4)/((mean((x - mean(x))^2))^2)) - 3 }
The function normex(n, nboot, theta) generates nboot independent random samples, each with n observations, from a standard normal distribution, and places them in a matrix. A user-specified function theta() is applied to each row of the matrix, returning a vector of theta() values for the nboot samples.
> normex function(n, nboot, theta) { data <- matrix(rnorm(n * nboot), nrow = nboot) apply(data, 1, theta) }
The function gammex(n, nboot, theta, al) generates nboot independent random samples, each with n observations, from a gamma distribution with shape parameter al, and places them in a matrix. A user-specified function theta() is applied to each row of the matrix, returning a vector of theta() values for the nboot samples.
> gammex function(n, nboot, theta, al) { data <- matrix(rgamma(n * nboot, al), nrow = nboot) apply(data, 1, theta) }
These functions can be used to study the sampling distributions of g1 and g2 for different sample sizes n and, in the case of the gamma distribution, for different values of the shape parameter al. A histogram is a useful way to display the results. Pick nboot large enough to give a smooth histogram but not so large that the computation takes too long.
> hist(normex(5,1000,skew))
The simulated values of g1 are distributed more or less symmetrically about the true value, which is 0.
> hist(normex(5,1000,excess))
Most of the simulated values of g2 are much less than the true value, which is 0.
> hist(gammex(5,1000,skew,1.5))
All of the simulated values of g1 are less than the true value, which is 2 / sqrt(1.5) = 1.63.
> hist(gammex(5,1000,excess,1.5))
All of the simulated values of g2 are much less than the true value, which is 6 / 1.5 = 4.