# C# Quasi Random Example

← All NMath Code Examples

```using System;
using System.Collections;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.Core.Examples.CSharp
{
/// <summary>
/// A .NET example in C# showing how to generate a sequence of quasirandom points.
/// </summary>
class QuasiRandomExample
{
static void Main( string[] args )
{
// A quasirandom sequence is a sequence of n-tuples that fills n-space
// more uniformly than uncorrelated random points. NMath provides both
// Sobol and Niederreiter quasirandom number generators.

// Create a Niederreiter quasirandom number generator object with dimension 1
// to generate quasirandom numbers (not vectors).
int dim = 1;
var nqrg = new NiederreiterQuasiRandomGenerator( dim );

// You can fill an existing matrix or array. (The points are the columns of the matrix,
// so the number of rows in the given matrix must be equal to the Dimension of the
// quasirandom number generator.)

// Here we create and fill a 1 x numPts matrix with quasirandom numbers from a uniform
// (0,1) distribution. The resulting numbers will occupy the first row of A.
int numPts = 100;
var A = new DoubleMatrix( nqrg.Dimension, numPts );
nqrg.Fill( A );

// Create a histogram to verify that the quasirandom numbers are uniformly
// distributed in the interval (0,1). We create 10 bins expecting that each
// bin will contain approximately numPts/10 points.
var h = new Histogram( 10, 0, 1 );

Console.WriteLine();

Console.WriteLine( "{0} uniform (0,1) quasirandom points:", numPts );
Console.WriteLine( h.StemLeaf() );

// Compare the above results with those from a pseudo-random generator. One expects
// a not-so uniform distribution in space.
int seed = 0x345;
var stream = new RandomNumberStream( seed, RandomNumberStream.BasicRandGenType.MersenneTwister );
var unifDist = new DoubleRandomUniformDistribution();
var pseudoRandom = new DoubleVector( numPts, stream, unifDist );
h.Reset();
Console.WriteLine( "\n{0} uniform (0,1) pseudo-random points:", numPts );
Console.WriteLine( h.StemLeaf() );

// Quasi-random numbers are useful for evaluating the integral of a function in the
// unit n-dimensional cube - essentially calculating the average of the function
// at a set of randomly selected points.

// In the following example we use a Sobol generator to approximate the integral
// of a function F over the 6-dimensional unit cube. The function is defined by
// F(x) = 1*cos(1*x1) * 2*cos(2*x2) * 3*cos(3*x3) *...* 6*cos(6*x6)
// where x = (x1, x2,..., x6) is a point in 6-dimensional space.
dim = 6;

// In this example we provide our own primitive polynomials for initializing
// the generator. Primitive polynomials have coefficients in {0, 1} and are
// specified by BitArray's containing the polynomial coefficients with the
// leading coefficient at index 0.

// We can supply either dim polynomials, or dim - 1 polynomials. If we specify
// dim - 1 polynomials the primitive polynomial for the first dimension will
// be initialized with a default value.
var primitivePolys = new BitArray[dim];
primitivePolys[0] = new BitArray( new bool[] { true, true } ); // x + 1
primitivePolys[1] = new BitArray( new bool[] { true, true, true } ); // x^2 + x + 1
primitivePolys[2] = new BitArray( new bool[] { true, false, true, true } ); // x^3 + x + 1
primitivePolys[3] = new BitArray( new bool[] { true, true, false, true } ); // x^3 + x^2 + 1
primitivePolys[4] = new BitArray( new bool[] { true, false, false, true, true } ); // x^4 + x + 1
primitivePolys[5] = new BitArray( new bool[] { true, true, false, false, true } ); // x^4 + x^3 + 1

var sobol = new SobolQuasiRandomGenerator( dim, primitivePolys );
int numPoints = 180000;
var points = new DoubleMatrix( dim, numPoints );
sobol.Fill( points );
double sum = 0;
for ( int i = 0; i < numPoints; i++ )
{
sum += F( points.Col( i ) );
}
sum /= numPoints;
double actualIntegralValue = -0.022;
Console.WriteLine( "Actual integral value = " + actualIntegralValue );
Console.WriteLine( "\nMonte-Carlo approximated integral value = " + sum );

Console.WriteLine();
Console.WriteLine( "Press Enter Key" );
}

/// <summary>
/// F(x) = cos(x[0])*2cos(2x[1])*...*6cos(6*x[5])
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
static double F( DoubleVector x )
{
if ( x.Length != 6 )
throw new InvalidArgumentException( "x must have length 6!" );

double y = 1;
for ( int i = 1; i <= x.Length; i++ )
{
y *= i * Math.Cos( i * x[i - 1] );
}
return y;
}

} // class

} // namespace

```
← All NMath Code Examples
Top