A Simulation Study of Skewness and Kurtosis using Splus

1998-09-28


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.


Back to Statistics 743