Park Miller distributions

I have posted previously about the Park Miller method for generating uniform random numbers. This post is about visualising the distribution and any modifications of the .random() function. The histographer is used to visualise the distribution and estimate the descriptive statistics. The default ParkMiller.random() uniform distribution is shown first. Thereafter a Marsaglia transform is applied to generate the standard normal distribution. Lastly, the standard normal is modified to give a guassian distribution.

The transforms to the .random() method are implemented in the inherited abstract random number generator class. This means that they can be applied to any implementation of the .random() method. I am interested to eventually see how they respond to other random number generators in the future.

Code

package net.avdw.gui
{
	import flash.display.Sprite;
	import flash.events.Event;
	import flash.utils.getTimer;
	import net.avdw.align.alignHorizontalCenterVerticalMiddleTo;
	import net.avdw.demo.ADemo;
	import net.avdw.math.Range;
	import net.avdw.random.generator.ParkMillerRng;
	
	public class ParkMillerRandomStreamingHistograph extends ADemo
	{
		private var histogram:StreamingHistograph;
		
		public function ParkMillerRandomStreamingHistograph()
		{
			if (stage)
				setup();
			else
				addEventListener(Event.ADDED_TO_STAGE, setup);
		}
		
		private function setup(e:Event = null):void
		{
			histogram = new StreamingHistograph("ParkMiller.random()", Range.closedOpen(0, 1));
			removeEventListener(Event.ADDED_TO_STAGE, setup);
			
			addChild(histogram);
			alignHorizontalCenterVerticalMiddleTo(stage, histogram);
			
			addEventListener(Event.ENTER_FRAME, update);
		}
		
		private const rng:ParkMillerRng = new ParkMillerRng();
		private function update(e:Event):void
		{
			var ticksPerSec:uint = 1000 / 50;
			var time:uint = getTimer();
			while (getTimer() - time < ticksPerSec)
				histogram.add(rng.random());
			
			histogram.refresh();
		}
	
	}

}

Uniform

Park Miller RNG applied with a Marsaglia Transform

Code

package net.avdw.random.generator
{
	
	/**
	 * ...
	 * @author Andrew van der Westhuizen
	 */
	public class ParkMillerRng extends ARng
	{
		/* The original seed used by this number generator */
		protected var seed:uint;
		protected var walkingNumber:uint;
		
		/**
		 * Setups the random number generator given a seed.
		 * If no seed is provided then a random seed is selected.
		 * @param	seed
		 */
		public function ParkMillerRng(seed:uint = 0)
		{
			if (seed == 0)
				seed = uint.MAX_VALUE * Math.random();
			
			walkingNumber = this.seed = seed;
		}
		
		override public function random():Number
		{
			return (walkingNumber = (walkingNumber * 16807) % 2147483647) / 0x7FFFFFFF + 0.000000000233;
		}
		
		override public function reset():ARng
		{
			walkingNumber = seed;
			return this;
		}
	}
}

Standard Normal

Park Miller RNG applied with a Marsaglia Transform

Code

/**
 * Marsaglia transform
 * http://en.wikipedia.org/wiki/Marsaglia_polar_method
 * http://blog.controul.com/2009/04/standard-normal-distribution-in-as3/
 * @return A number from the standard normal distribution
 */
private var useCache:Boolean = true;
private var cache:Number = 0;
final public function stdNormal():Number
{
	useCache = !useCache;
	if (useCache)
		return cache;
		
	var x:Number, y:Number, w:Number;
	do 
	{
		x = float( -1, 1);
		y = float( -1, 1);
		w = x * x + y * y;
	}
	while (w >= 1 || w == 0);
	
	w = Math.sqrt( -2 * Math.log(w) / w);
	cache = x * w;
	return y * w;
}

Gaussian

Park Miller RNG applied with a Marsaglia Transform

Code

/**
 * Returns a number that fits the gaussian parameters
 * @param	mean
 * @param	stddev
 * @return 	A number fitting the gaussian parameters
 */
final public function gaussian(mean:Number, stddev:Number):Number {
	return mean + stdNormal() * stddev;
}