Quasi-Random Sequences in IronPython QuickStart Sample

Illustrates how to generate quasi-random sequences like Fauré and Sobol sequences using classes in the Numerics.NET.Statistics.Random namespace in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This sample demonstrates how to use quasi-random sequences to compute multi-dimensional integrals efficiently and accurately.

The sample shows how to use both Halton and Sobol sequences to evaluate a specific 5-dimensional integral. It demonstrates the key differences between these sequence types and provides practical examples of their usage. The program:

• Creates a Halton sequence and uses it to estimate a 5-dimensional integral
• Shows how to iterate through the sequence points and evaluate the integrand
• Demonstrates progress monitoring by displaying intermediate results
• Illustrates the use of the more sophisticated Sobol sequences
• Shows how to skip points in a Sobol sequence efficiently
• Compares the results to the known exact value of the integral

The sample serves as a practical introduction to quasi-random number generation and shows how these sequences can be more effective than pseudo-random numbers for certain numerical integration tasks. It includes examples of proper sequence initialization and demonstrates recommended usage patterns for both sequence types.

The code

import numerics

from Extreme.Statistics.Random import *
from Extreme.Mathematics import *
from Extreme.Mathematics.LinearAlgebra import *

# Illustrates the use of quasi-random sequences by computing
# a multi-dimensional integral.

# This QuickStart Sample demonstrates the use of
# quasi-random sequences by computing
# a multi-dimensional integral.

# We will use one million points.
capacity = 10000
# The number of dimensions:
dimension = 5

# We will evaluate the function
#
#    Product(i = 1 -> # dimensions) |4 x[i] - 2|
#
# over the hypercube 0 <= x[i] <= 1. The value of this integral
# is exactly 1.

# This vector will hold the current vector in the sequence:
point = Vector.Create(dimension)
# Create the sequence:
sequence = HaltonSequence(point, capacity)

print "# iter.  Estimate"
# Compute the integral by summing over all points:
sum = 0.0
for i in range(1, capacity):
	if i % 1000 == 0:
		print "{0:6}  {1:8.4f}".format(i, sum / i)

	# Evaluate the integrand:
	functionValue = 1.0
	for j in range(dimension):
		functionValue *= abs(4.0*point[j]-2.0)
	sum += functionValue

	sequence.MoveNext()

# Print the final result.
print "Final estimate: {0:8.4f}".format(sum / capacity)
print "Exact value: 1.0000"