VB Quasi Random Example

← All NMath Code Examples

 

Imports System
Imports System.Collections

Imports CenterSpace.NMath.Core

Namespace CenterSpace.NMath.Core.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)
      H.AddData(A.Row(0))

      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()
      H.AddData(pseudoRandom)
      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 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.
      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")
      Console.Read()

    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