Generate Mixed Non-Normal Quasi-Random Vectors

[digitalc6839]

Hello,

I need to be able to generate random numbers from a mixed non-normal set of vectors (e.g. 2 Normal, 1 Exponential, 2 Lognormal, etc.) from a Quasi-Random stream such as SOBOL. The documentation in vslnotes.pdf provides the code example (which I pasted in below) but this would bea 100x3 array of all Uniform. I need to be able to mix the vectors, so perhaps it is still 100 (rows) x 3 (columns)array but the first column is Normal, the second column is Uniform, and the third column is LogNormal.

I would appreciate any help with this anyone can provide.

include

include "mkl.h"

float r[100][3]; /* buffer for quasi-random numbers */

VSLStreamStatePtr stream;

/* Initializing */

vslNewStream( &stream, VSL_BRNG_SOBOL, 3 );

/* Generating */

vsRngUniform( VSL_METHOD_SUNIFORM_STD,

stream, 100*3, (float*)r, 2.0f, 3.0f );

/* Deleting the streams */

vslDeleteStream( &stream );

[Andrey_N_Intel]

Generation of distribution mixture using SOBOL QRNG can be organizedby means ofLeapfrogStream routine for individual components of quasi-random vector. Example below demonstrates this.

First component of the quasi-random vector is used to generate normally distrbuted random numbers.Second component of the quasi-random vector is for generation of uniform numbers. Finally, third component is for generation of lognormally distributed sequence.
Current implementation of VSL lognormal generator does not support Inverse Method for underlying Gaussian generator. So, direct call to RngLognormal function of the library is replaced with call to Gaussian generator which is followed by necessary transformations. MKL vector exponent function is used to transform normal numbers (low accuracy version of the function will speed up the application).

#include

#include "mkl.h"

#define SEED 3
#define N 100

main()
{
VSLStreamStatePtr stream, stream1, stream2, stream3;
int old_mode;

float r1, r2, r3;
int i;
float a, b, m1, m2, sigma1, sigma2, b1, beta;

/* Set low accuracy mode of calculations */
old_mode = vmlSetMode( VML_LA );
/* Initialize streams */
vslNewStream( &stream, VSL_BRNG_SOBOL, SEED );

vslCopyStream( &stream1, stream );
vslCopyStream( &stream2, stream );
vslCopyStream( &stream3, stream );

/* Leapfrog over individual components of quasi-random vector */
vslLeapfrogStream( stream1, 0, VSL_QRNG_LEAPFROG_COMPONENTS );
vslLeapfrogStream( stream2, 1, VSL_QRNG_LEAPFROG_COMPONENTS );
vslLeapfrogStream( stream3, 2, VSL_QRNG_LEAPFROG_COMPONENTS );

/* Call RNGs */
/* Uniform distribution parameters */
a = 0.0f; b = 1.0f;
/* Gaussian distribution parameters */
m1 = 10.0f; sigma1 = 3.0f;
/* Lognormal distribution parameters */
m2 = 0.0f; sigma2 = 1.0f; b1 = 3.0f; beta = 4.0f;

/* Generate normal numbers */
vsRngGaussian( VSL_METHOD_SGAUSSIAN_ICDF, stream1, N, r1, m1, sigma1 );
/* Generate uniform numbers */
vsRngUniform ( VSL_METHOD_SUNIFORM_STD, stream2, N, r2, a, b );
/* Generate lognormal numbers */
vsRngGaussian( VSL_METHOD_SGAUSSIAN_ICDF, stream3, N, r3, m2, sigma2 );
vsExp( N, r3, r3 );
for ( i = 0; i < N; i++ ) r3 = r3 * beta + b1;

/* Delete streams */
vslDeleteStream( &stream );
vslDeleteStream( &stream1 );
vslDeleteStream( &stream2 );
vslDeleteStream( &stream3 );

/* Restoring old mode */
vmlSetMode( old_mode );

return 0;
}

More details about MKL functions used in the ex ample are available in MKL Manual and VslNotes.Please, let us know how this scheme works for you.

[digitalc6839]

Thanks for the excellent answer. I appreciate it. I have read the documentation on the LeapfrogStream to include the vslnotes.pdf. I am a little confused though... In my understanding of theory, any use of Quasi-Random streams should be paired with the inverse cdf method of variant generation. I tested this with the Normal by pairing Niederr with BOXMULLER2 by running 10,000 vectors of N=100 and while calculating the mean of each 100. The results is that the mean is significantly biased from the true value even over 10,000 runs. When I replace BOXMULLER2 with ICDF then the bias completely goes away and the mean and Standard Deviation of the N=100 samples is much closer to theoretical.

So, my fundamental understanding would be that any Quasi-Random stream should always be paired with inverse CDF. When running any non-Quiasi-Random stream (such as MT2203) it is then faster to use something like BOXMULLER2 instead of inverse CDF. Is this correct?

My confusion stems around your statement "Current implementation of VSL lognormal generator does not support Inverse Method for underlying Gaussian generator". The help system list four different methods for vdRngLognormal and they all appear to be inverse CDFs. They are listed as below. How can I tell from looking at the help system which of the non-uniform methods are available and which are not?

VSL_METHOD_SLOGNORMAL_ICDF
VSL_METHOD_DLOGNORMAL_ICDF
VSL_METHOD_SLOGNORMAL_ICDF_ACCURATE
VSL_METHOD_DLOGNORMAL_ICDF_ACCURATE

By the way, the on page 31 of vslnotes.pdf the word sentence "VSL implements the leapfrog method trough the function vslLeapfrogStream:" has the word "trough" which should probably be "through".

[Andrey_N_Intel]

Thanks for your answer.
We know that method name for Lognormal distribution generator is not accurately set; we work to fix it. Also, thank you for communicating the misprint in the documentation.

Inverse CDF method for Gaussian generator is faster than Box-Muller methodin 32 bit and EM64T versions of the library.

[digitalc6839]

Thank you.

I ran a quick test comparing ICDF and BOXMULLER2, and found BOXMULLER2 to be about twice as fast as ICDF. Of course, they are both blazingly fast but I would like to know if I am doing something wrong as in your post you indicate that ICDF is faster. I generated vectors of both of size 10,000 and iterated on this 10,000 times (for a total of 100 Million numbers of each). The results on my computer are ...

ICDF = 4.25 seconds

BOXMULLER2 = 1.859 seconds

Of course BOXMULLER2 is faster than BOXMULLER (without the two aka single variant instead of pairs) but BOXMULLER still comes in at 3.328 seconds. The code I used is below.

Am I doing something wrong?

double r1,r2;
 
 VSLStreamStatePtr stream,stream1,stream2;
int errcode,i;
 DWORD startTime,endTime,deltaTime;

 
 //setup the stream
 // VSL_BRNG_NIEDERR SL_BRNG_SOBOL // VSL_BRNG_MT19937 // VSL_BRNG_NIEDERR 
errcode = vslNewStream( &stream, VSL_BRNG_NIEDERR, 2 );
 CheckVslError( errcode );

//two streas from Quasi-Random
//number of streams must match the last argument in vslNewStream
vslCopyStream( &stream1, stream );
 vslCopyStream( &stream2, stream );

/* Leapfrog over individual components of quasi-random vector */
 vslLeapfrogStream( stream1, 0, VSL_QRNG_LEAPFROG_COMPONENTS );
 vslLeapfrogStream( stream2, 1, VSL_QRNG_LEAPFROG_COMPONENTS );

//Timing for ICDF Loop through 10,000 iterations of vectors of size N
startTime = GetTickCount();
for (i=0;i<10000;i++)
{

//Normal (Gaussian) Double Precision 
errcode = vdRngGaussian(VSL_METHOD_DGAUSSIAN_ICDF,stream1, N, r1, 10.0f, 1.0f );
CheckVslError( errcode );

}
endTime = GetTickCount();
deltaTime = endTime - startTime;
 printf ("Completetion Time ICDF:%ld
",deltaTime);

 //Timing for BoxMuller2 Loop through 10,000 iterations of vectors of size N
startTime = GetTickCount();
for (i=0;i<10000;i++)
{

//Normal (Gaussian) Double Precision 
errcode = vdRngGaussian(VSL_METHOD_DGAUSSIAN_BOXMULLER2,stream2, N, r2, 10.0f, 1.0f );
CheckVslError( errcode );

}
endTime = GetTickCount();
deltaTime = endTime - startTime;
 printf ("Completetion Time BoxMuller2:%ld
",deltaTime);

 vslDeleteStream( &stream );
 vslDeleteStream( &stream1 );
 vslDeleteStream( &stream2 );

[Andrey_N_Intel]

What version of the library do you use? Speed-up of ICDF method vs. Box-Muller method is expected in MKL 10. Thanks, Andrey