Non-Uniform Random Numbers QuickStart Sample (IronPython)

import numerics

from System import Array

from Extreme.Statistics.Distributions import *
from Extreme.Statistics.Random import *

# Illustrates generating non-uniform random numbers
# using the classes in the Extreme.Statistics.Random
# namespace.

# Random number generators and the generation
# of uniform pseudo-random numbers are illustrated
# in the UniformRandomNumbers QuickStart Sample.

# In this sample, we will generate numbers from
# an exponential distribution, and compare summary
# results to what would be expected from 
# the corresponding Poisson distribution.

meanTimeBetweenEvents = 0.42

# We will use the exponential distribution to generate the time 
# between events. The number of events per unit time follows
# a Poisson distribution.

# The parameter of the exponential distribution is the time between events.
exponential = ExponentialDistribution(meanTimeBetweenEvents) 
# The parameter of the Poisson distribution is the mean number of events
# per unit time, which is the reciprocal of the time between events:
poisson = PoissonDistribution(1 / meanTimeBetweenEvents)

# We use a MersenneTwister to generate the random numbers:
random = MersenneTwister()

# The totals array will track the number of events per time unit.
totals = Array.CreateInstance(int, 15)

currentTime = 0
endOfCurrentTimeUnit = 1
eventsInUnit = 0

totalTime = 0
count = 0

while currentTime < 100000:
    timeBetween = exponential.Sample(random)
    totalTime += timeBetween 
    count = count + 1
    # Alternatively, we could have written
    #   timeBetween = random.NextDouble(exponential)
    # which would give an identical result.
    currentTime += timeBetween
    while currentTime > endOfCurrentTimeUnit:
        if eventsInUnit >= totals.Length:
            eventsInUnit = totals.Length-1
        totals[eventsInUnit] = totals[eventsInUnit] + 1
        eventsInUnit = 0
        endOfCurrentTimeUnit = endOfCurrentTimeUnit + 1
    eventsInUnit = eventsInUnit + 1

print "{0}", totalTime / count
# Now print the totals
print "# Events    Actual  Expected"
for i in range(0, totals.Length):
	expected = 100000 * poisson.Probability(i)
	print "{0:8}  {1:8}  {2:8.1f}".format(i, totals[i], expected)