# VB Quasi Random Example

← All NMath Code Examples

```Imports System
Imports System.Collections

Imports CenterSpace.NMath.Core

Namespace CenterSpace.NMath.Examples.VisualBasic

<summary>
A .NET example in Visual Basic showing how to generate a sequence of quasirandom points.
</summary>
Module QuasiRandomExample

Sub Main()

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).
Dim Dimensions As Integer = 1
Dim nqrg As New NiederreiterQuasiRandomGenerator(Dimensions)

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.
Dim NumPts As Integer = 100
Dim A As 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.
Dim H As 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.
Dim Seed As Integer = &H345
Dim Stream As New RandomNumberStream(Seed, RandomNumberStream.BasicRandGenType.MersenneTwister)
Dim UnifDist As New DoubleRandomUniformDistribution()
Dim pseudoRandom As New DoubleVector(NumPts, Stream, UnifDist)
H.Reset()
Console.WriteLine()
Console.WriteLine("{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.
Dimensions = 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 BitArrays containing the polynomial coefficients with the

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.
Dim primitivePolys(Dimensions - 1) As BitArray
primitivePolys(0) = New BitArray(New Boolean() {True, True}) x + 1
primitivePolys(1) = New BitArray(New Boolean() {True, True, True}) x^2 + x + 1
primitivePolys(2) = New BitArray(New Boolean() {True, False, True, True}) x^3 + x + 1
primitivePolys(3) = New BitArray(New Boolean() {True, True, False, True}) x^3 + x^2 + 1
primitivePolys(4) = New BitArray(New Boolean() {True, False, False, True, True}) x^4 + x + 1
primitivePolys(5) = New BitArray(New Boolean() {True, True, False, False, True}) x^4 + x^3 + 1

Dim Sobol As New SobolQuasiRandomGenerator(Dimensions, primitivePolys)
Dim NumPoints As Integer = 180000
Dim Points As New DoubleMatrix(Dimensions, NumPoints)
Sobol.Fill(Points)
Dim Sum As Double = 0

For I = 0 To NumPoints - 1
Sum += F(Points.Col(I))
Next

Sum /= NumPoints
Dim actualIntegralValue As Double = -0.022
Console.WriteLine("Actual integral value = " & actualIntegralValue)
Console.WriteLine()
Console.WriteLine("Monte-Carlo approximated integral value = " & Sum)

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

End Sub

<summary>
F(x) = cos(x(0))*2cos(2x(1))*...*6cos(6*x(5))
</summary>
<param name="x"></param>
<returns></returns>
Function F(ByVal X As DoubleVector) As Double

If (X.Length <> 6) Then
Throw New InvalidArgumentException("x must have length 6!")
End If

Dim Y As Double = 1

For i = 1 To X.Length
Y *= i * Math.Cos(i * X(i - 1))
Next
Return Y
End Function

End Module
End Namespace
```
← All NMath Code Examples
Top