Thursday, November 24, 2011

Normally Distributed Random Numbers

So you have a random generator that produces uniformly distributed random numbers from 0 to 1 and you want to create a normal distributed random number generator.

        /// <summary>
        /// Returns a normally distributed random number
        /// </summary>
        /// <param name="mean">The mean value of the distribution</param>
        /// <param name="stdDev">The standard deviation of the distribution</param>
        /// <returns>A normally distributed random number</returns>
        public Double NextDouble(double mean, double stdDev)
        {
            Random rand = new Random();
            double u1 = rand.NextDouble(), u2 = rand.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }

The only problem is that if you want to create lots of random numbers all the logarithms and square roots just bog down the machine especially if you want real time performance. So what do you do ?

Simple, trade memory for computation time. I created a simple class that is used like the Random class of the .NET framework but it internally caches normally distributed numbers upon initialization.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Emgu_Tracker_Advanced_V2.Source
{
    public class RandomNormalCached : Random
    {
        #region Variables
        protected int noCached = 101;
        /// <summary>
        /// Array of precomputed number with 0 mean and 1 variance
        /// </summary>
        protected double[] precomputedNumbers;
        /// <summary>
        /// Internal state counter
        /// </summary>
        protected int counter;
        #endregion

        #region Constructor
        public RandomNormalCached()
            : base()
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        /// <summary>
        /// Constructor with custom number of cached random numbers
        /// </summary>
        /// <param name="_noCached">Number of cached random numbers</param>
        /// <param name="_seed">Random seed initializer</param>
        public RandomNormalCached(int _noCached, int _seed)
            : base(_seed)
        {
            counter = 0;
            noCached = _noCached;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        public RandomNormalCached(int _seed)
            : base(_seed)
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }
        #endregion

        #region Functions
        public Double NextDouble(double mean, double stdDev)
        {
            double result = mean + stdDev * precomputedNumbers[counter];
            counter = (counter + 1) % noCached;
            return result;
        }

        /// <summary>
        /// Returns a nornal distributed random number
        /// </summary>
        /// <param name="mean">The mean/expected value</param>
        /// <param name="stdDev">The standard deviation</param>
        /// <returns>Normal distributed normal number</returns>
        protected Double ComputeNormalDistribution(double mean, double stdDev)
        {
            double u1 = base.NextDouble(), u2 = base.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }
        #endregion
    }
}