Simple Skew-Normal PRNG in JavaScript

The basic JavaScript Math.Random returns values with an approximately uniform distribution. While this is often sufficient, many situations require normal distributions. Some situations may be even better served by an approximately normal, but asymmetric, distribution. These are known as skew-normal distributions.

Skew-Normal Distributions

A normal distribution has the same mean, median, and mode at its center, while a skew-normal distribution has the mean, median, and mode spread across the distribution. A shape parameter (α) describes the skewness of the distribution. The plot below (from Skbkekas) shows a couple of example distributions with various shape parameters.

Skew normal densities

There are some real-life situations that behave in a skew-normal fashion. For instance, income distribution is skew-normal, with the vast majority of households having a fairly standard income and a long tail of extreme values. Another example is a game that may want to assign values or award prizes in a skew-normal manner.

The Code

A random value from a given distribution is known as a variate. Given the uniform distribution of JavaScript’s Math.Random, getting a skew-normal variate requires a couple of steps.

Box-Muller Transform

We first want to get normal variates from our uniform variates. Although there are various methods for this, the Box-Muller transform is considered efficient and easy to implement. The Box-Muller transform produces two independent normal variates which we’ll need to get the skew-normal variate, so it works out well.


const randomNormals = (rng) => {
    let u1 = 0, u2 = 0;
    //Convert [0,1) to (0,1)
    while (u1 === 0) u1 = rng();
    while (u2 === 0) u2 = rng();
    const R = Math.sqrt(-2.0 * Math.log(u1));
    const Θ = 2.0 * Math.PI * u2;
    return [R * Math.cos(Θ), R * Math.sin(Θ)];
};

In essence, the uniform variates are used to establish a magnitude (R) and direction (Θ) in polar coordinates. Those are translated into a Cartesian system, and the result describes two independent, normally distributed random variates.

The argument rng represents any uniform random number generator. For a basic case, you can simply pass in Math.random.

Skew-normal transform

Next, we convert the two independent normal variates generated above into a skew-normal variate based off work by Azzalini and Dalla Valle (1996).


const randomSkewNormal = (rng, ξ, ω, α = 0) => {
    const [u0, v] = randomNormals(rng);
    if (α === 0) {
        return ξ + ω * u0;
    }
    const 𝛿 = α / Math.sqrt(1 + α * α);
    const u1 = 𝛿 * u0 + Math.sqrt(1 - 𝛿 * 𝛿) * v;
    const z = u0 >= 0 ? u1 : -u1;
    return ξ + ω * z;
};

The ξ, ω, and α refer to the location (mean), scale (standard deviation), and shape (skewness). You can see that when skewness is 0, the return value is simply the normal variate transformed in a standard way to the given mean with the given standard deviation.

When skewness is not 0, we first determine the correlation coefficient (𝛿) from the shape parameter (α). This particular relationship is described by Azzalini (1985) and corresponds to line 6 above.

Then we can use the correlation coefficient to produce a correlated pair of variates via the formula on line 7. This method of correlating independent random variables is well documented in various statistics textbooks.

Azzalini (1996) gives us an efficient way to determine the skew-normal variate for random number generation that avoids rejection sampling, namely line 8. We can then transform this value to the given location and scale, giving us a skew-normal variate!

References

  • Azzalini, A. (1985). A class of distributions which inclues the normal ones. Scand. J. Statist.
  • Azzalini, A. & Dalla Valle, A. (1996). The multivariate skew-normal distribution. Biometrika.