Search code examples
javascriptnormal-distribution

How to get a random value from truncated normal distribution with mean and std in JavaScript?


What is a JS alternative to the same Python implementation?

import matplotlib.pyplot as plt
from scipy.stats import truncnorm
import numpy as np
mean = 1
std = 2
clip_a = -4
clip_b = 3

a, b = (clip_a - mean) / std, (clip_b - mean) / std
x_range = np.linspace(-3 * std, 3 * std, 1000)
plt.plot(x_range, truncnorm.pdf(x_range, a, b, loc = mean, scale = std));

enter image description here

I'd like to get a random value according to the distribution (in JS the same code with size=1):

dist = truncnorm.rvs(a, b, loc = mean, scale = std, size=1000000)
plt.hist(dist);

enter image description here


Solution

  • Here is a JS function that implements a truncated, skew-normal pseudo random number generator (PRNG). It is based on this blog post by Tom Liao and has been extended to consider lower and upper bounds (truncation).

    Essentially, the function is called recursively until a variate within the desired bounds is found.
    You can pass your own random number generator using the rng property, though Math.random will be used by default. Also, as you didn't ask for a skew-normal distribution, you can just ignore the skew property as it defaults to 0. This will leave you with a truncated normal PRNG, just as you asked for.

    function randomTruncSkewNormal({
      rng = Math.random,
      range = [-Infinity, Infinity],
      mean,
      stdDev,
      skew = 0
    }) {
      // Box-Muller transform
      function 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(Θ)];
      }
    
      // Skew-normal transform
      // If a variate is either below or above the desired range,
      // we recursively call the randomSkewNormal function until
      // a value within the desired range is drawn
      function randomSkewNormal(rng, mean, stdDev, skew = 0) {
        const [u0, v] = randomNormals(rng);
        if (skew === 0) {
          const value = mean + stdDev * u0;
          if (value < range[0] || value > range[1])
            return randomSkewNormal(rng, mean, stdDev, skew);
          return value;
        }
        const sig = skew / Math.sqrt(1 + skew * skew);
        const u1 = sig * u0 + Math.sqrt(1 - sig * sig) * v;
        const z = u0 >= 0 ? u1 : -u1;
        const value = mean + stdDev * z;
        if (value < range[0] || value > range[1])
          return randomSkewNormal(rng, mean, stdDev, skew);
        return value;
      }
    
      return randomSkewNormal(rng, mean, stdDev, skew);
    }
    

    Calling this function in the following manner

    const data = [];
    for (let i = 0; i < 50000; i++) {
      data.push({
        x: i,
        y: randomTruncSkewNormal({ 
          range: [-4,3], 
          mean: 1, 
          stdDev: 2
        })
      });
    }
    

    and plotting the data using your charting library of choice should give your the desired output. vega plot of truncated normal distribution

    I also made a small Observable notebook interactively demonstrating the function which you might want to check out as well.