VB.NET Independent Streams Example

Imports System

Imports CenterSpace.NMath.Core
Imports CenterSpace.NMath.Matrix

Namespace CenterSpace.NMath.Core.Examples.VisualBasic

  ' A .NET example in C# showing how to create streams of independent random numbers.
  Module IndependentStreamsExample

    Sub Main()

      ' NMath provides classes for generating several independent streams of random
      ' numbers using the leapfrop and skip-ahead (blocksplitting) methods.

      ' IMPORTANT NOTE:
      ' Skip-ahead and leapfrog streams are only supported for basic random number
      ' generator types which provide a more efficient algorithm than generation
      ' of that full sequence to pick out a required subsequence. The static variable
      ' SkipAheadStream.SupportedGeneratorTypes contains a list of the supported
      ' types for SkipAheadStream.

      Dim Seed As Integer = &H124
      Dim GenType As RandomNumberStream.BasicRandGenType = RandomNumberStream.BasicRandGenType.MultiplicativeCongruent31

      ' Create a SkipAheadRandomStreams object for generationg 10 independent random number streams, 
      ' each of length 100, using the Skip Ahead method.
      Dim NumberOfStreams As Integer = 5
      Dim StreamLength As Integer = 1000
      Dim IndStreams As New SkipAheadRandomStreams(Seed, GenType, NumberOfStreams, StreamLength)

      ' Create a streamLength x numberOfStreams matrix A and use the SkipAheadRandomStreams object 
      ' created above to fill it with random numbers. Each of the numberOfStreams columns 
      ' in the resulting matrix will contain streamLength random deviates distributed 
      ' uniformly in the interval (0, 1). Each column of random numbers will 
      ' be independent from the others.
      Dim Dist As New DoubleRandomUniformDistribution()
      Dim A As New DoubleMatrix(StreamLength, NumberOfStreams)
      IndStreams.Fill(Dist, A)

      ' Print out the correlation matrix for the sample to test
      ' independence. Correlations should be near zero for independent
      ' samples.
      Dim CorrelationMatrix As DoubleSymmetricMatrix = RankCorrelation(A, True)
      Console.WriteLine("Spearmans Rank Correlation for independent random samples:")
      Dim I, J As Integer
      For I = 0 To CorrelationMatrix.Rows - 1
        For J = 0 To I - 1
          Console.WriteLine("Correlation between column {0} and row {1} is {2}", I, J, CorrelationMatrix(I, J))
        Next
      Next

      ' You can also generate independent streams from different probability
      ' distributions. Here we generate 3 streams each with a different
      ' distribution using the leapfrog technique.
      NumberOfStreams = 3
      Dim LeapfrogIndStreams As New LeapfrogRandomStreams(Seed, GenType, NumberOfStreams, StreamLength)
      Dim Distributions(NumberOfStreams) As IRandomNumberDistribution(Of Double)
      Distributions(0) = New DoubleRandomBetaDistribution()
      Distributions(1) = New DoubleRandomExponentialDistribution()
      Distributions(2) = New DoubleRandomLogNormalDistribution()
      Dim C As DoubleMatrix = LeapfrogIndStreams.Next(Distributions)

      ' Print out the rank correlation matrix for the sample to test
      ' independence. Correlations should be near zero for independent
      ' samples.
      CorrelationMatrix = RankCorrelation(C, True)
      Console.WriteLine()
      Console.WriteLine("Spearmans Rank Correlation for samples from different distributions:")
      For I = 0 To CorrelationMatrix.Rows - 1
        For J = 0 To I - 1
          Console.WriteLine("Correlation between column {0} and row {1} is {2}", I, J, CorrelationMatrix(I, J))
        Next
      Next

      ' Independent streams from discrete distributions are also supported.
      StreamLength = 5
      GenType = RandomNumberStream.BasicRandGenType.WinchannHillCombined
      LeapfrogIndStreams = New LeapfrogRandomStreams(GenType, NumberOfStreams, StreamLength)
      Dim DiscreteDistributions(NumberOfStreams) As IRandomNumberDistribution(Of Integer)
      DiscreteDistributions(0) = New IntRandomBernoulliDistribution(0.6)
      DiscreteDistributions(1) = New IntRandomGeometricDistribution(0.2)
      DiscreteDistributions(2) = New IntRandomPoissonDistribution(0.8)
      Dim DiscreteRandomSamples()() As Integer = LeapfrogIndStreams.Next(Of Integer)(DiscreteDistributions)
      ' The ith independent stream is in the array discreteRandomSamples[i].
      For I = 0 To NumberOfStreams - 1
        Console.WriteLine()
        Console.WriteLine("Distribution: {0}", DiscreteDistributions(I).GetType())
        For J = 0 To DiscreteRandomSamples(I).Length - 1
          Console.Write("{0}{1}", DiscreteRandomSamples(I)(J), If(J = DiscreteRandomSamples(I).Length, "\n", ", "))
        Next
      Next

      Console.WriteLine()
      Console.WriteLine("Press Enter Key")
      Console.Read()
    End Sub

    ' <summary>
    ' Computes the Spearman rank correlation matrix for a set of random inputs. The random 
    ' inputs are taken to be the columns of the input matrix. The symmetric,
    ' positive definite matrix whose i,j entry is the Spearman correlation coefficient
    ' between the inputs in column i and column j is computed and returned.
    ' </summary>
    ' <param name="A">Matrix whose columns are the inputs whose Spearman rank correlation
    ' coefficient is computed.</param>
    ' <param name="useMidRankForTies">If true mid ranks are used for ranking ties.
    ' </param>
    ' <returns>The symmetric, positive definite matrix whose row i, column j 
    ' entry is  the correlation coefficient between the inputs in column i and column j
    ' </returns>
    Function RankCorrelation(ByVal A As DoubleMatrix, ByVal useMidRankForTies As Boolean) As DoubleSymmetricMatrix

      'Console.WriteLine("A.Rows " & A.Rows)
      'Console.WriteLine("A.Cols " & A.Cols)

      Dim N As Integer = A.Rows
      'Console.WriteLine(N)
      Dim SortedData(N - 1) As Double
      'Console.WriteLine(SortedData.Length)
      Dim Indices(N - 1) As Integer
      'Console.WriteLine(Indices.Length)

      Dim XRanks As New DoubleVector(N)
      Dim YRanks As New DoubleVector(N)

      Dim R As New DoubleSymmetricMatrix(A.Cols)

      For I = 0 To A.Cols - 1
        For J = 0 To A.Cols - 1
          If (I = J) Then
            R(I, J) = 1.0
          Else
            RankData(A.Col(I), SortedData, Indices, XRanks, useMidRankForTies)
            RankData(A.Col(J), SortedData, Indices, YRanks, useMidRankForTies)
            R(I, J) = RankCorrelationFromRanks(XRanks, YRanks)
          End If
        Next
      Next
      Return R
    End Function

    Sub RankData(ByVal Data As DoubleVector, ByVal SortedData() As Double, ByVal Indices() As Integer, ByVal Ranks As DoubleVector, ByVal useMidRankForTies As Boolean)

      'Console.WriteLine("data: " & Data.Length & " sorted " & SortedData.Length)

      If (Data.Length <> SortedData.Length) Then
        Throw New InvalidArgumentException("data and sortedData arrays do not have the same length in RankData")
      End If

      If (Data.Length <> Indices.Length) Then
        Throw New InvalidArgumentException("data and indicies arrays do not have the same length in RankData")
      End If

      For I = 0 To Indices.Length - 1
        SortedData(I) = Data(I)
        Indices(I) = I
      Next

      Array.Sort(Of Double, Integer)(SortedData, Indices)

      For I = 0 To SortedData.Length - 1
        Ranks(I) = Array.IndexOf(Of Integer)(Indices, I) + 1
      Next

      If (useMidRankForTies) Then
        ApplyMidranks(SortedData, Ranks)
      End If

    End Sub

    Sub ApplyMidranks(ByRef Data() As Double, ByRef Ranks As DoubleVector)

      For I = 1 To Data.Length - 1
        If (Data(I) = Data(I - 1)) Then
          Dim Start As Integer = I - 1
          Dim Count As Integer = 0
          Dim Ranksum As Double = 0

          While (Start + Count < Data.Length & Data(Start) = Data(Start + Count))
            Ranksum += Ranks(Start + Count)
            Count = Count + 1
          End While

          Dim MidRanks As Double = Ranksum / CType(Count, Double)
          For S = 0 To Count - 1
            Ranks(Start + S) = MidRanks
          Next
          I = Start + Count
        End If
      Next
    End Sub

    Function RankCorrelationFromRanks(ByRef xRanks As DoubleVector, ByRef yRanks As DoubleVector) As Double

      Dim N As Integer = xRanks.Length
      If (yRanks.Length <> N) Then
        Throw New InvalidArgumentException("All arrays must have same length in RankCorreleation")
      End If

      Dim Sum As Double = 0 ''

      For I = 0 To N - 1

        Dim D As Double = xRanks(I) - yRanks(I)
        Sum += (D * D)
      Next
      Return (1.0 - ((6 * Sum) / (N * (N * N - 1))))
    End Function

  End Module
End Namespace