In this posting, we will be
- getting familiar with NAG .NET library, specifically its class G05 for generators for creating pseudo-random numbers and low-discrepancy sequences.
- creating easy-to-use C# IRandom interface and its two implementations, which are effectively wrapping specific NAG methods for generating random numbers.
- interfacing our C# program with Excel by using Excel-DNA.
COOLED BY RANDOMNESS
Generating random numbers is a serious business. A fundamental difference in random number generation is the one between pseudorandom numbers and low-discrepancy sequence. Statistical properties of pseudo-random numbers are very close to "true" random numbers. However, low-discrepancy sequences are designed to give a more even distribution in space, because each number in sequence is designed to be maximally avoiding of the others. For this reason, they are preferred choice for Quasi Monte Carlo methods. This method may use less sample paths to converge to a true value faster than method using pseudo-random numbers.
NAG
Numerical Algorithms Group has been developing impressive collection of numerical libraries with already proven industrial strength. The company itself has over forty years of experience and staffed with experts from the fields of mathematics and statistics. The documentation of class methods is, without a doubt, the best I have seen. Also, there are a lot of useful "copy-paste-run" examples, just to get you started. Also, NAG customer support is one of the best I have experienced. Support people were studying and solving my problem, instead of just trying to shake me off.
As a teaser, NagLibrary namespace for .NET is presented in the picture below.
Unfortunately, as with all the best things in our life, also NAG is not free. However, the company is offering a trial period, before you make any purchasing decision for the product.
NAG random number generators
At this posting, we are concentrating on NAG .NET library class G05 for generating random numbers. Class G05 contains numerous different methods for this purpose, but we are using only two specific methods in this posting. We are going to create easy C# wrapper classes for the following methods.
- g05sk which generates a vector of pseudorandom numbers taken from a normal distribution. We will use Mersenne Twister algorithm for creating pseudorandom numbers.
- g05ym which generates a uniformly distributed low-discrepancy sequence. We will use Sobol sequence for creating low-discrepancy sequence.
STEP ONE: C# program
Create a new C# class project (RandomGenerator). First, we are going to create interface (IRandom) for random class implementations as shown in the picture below. This interface has only one public method - getRandom. For this method, we define the size of random matrix (two dimensional array) we would like to receive (number of random numbers and number of dimensions). This method is then returning required matrix, filled with random numbers.
using System; // namespace RandomGenerator { public interface IRandom { double[,] getRandom(int rows, int cols); } }
Next, we are going to create two different interface implementations for our IRandom interface. The first implementation is using Mersenne Twister algorithm for creating a set of pseudo-random numbers. Class implementation is shown in the picture below. Insert a new class (NAG_mersenneTwister) into your project and copyPaste the code into it. Remember also to add reference to NagLibrary32.dll file which you will find in your NAG folder.
using System; using NagLibrary; // namespace RandomGenerator { // wrapper for NAG mersenne twister pseudorandom numbers generator public class NAG_mersenneTwister : IRandom { private const int generator = 3; // hard-coded mersenne twister generator type private const int subGenerator = 1; private const double mean = 0.0; private const double var = 1.0; // public double[,] getRandom(int nRows, int nCols) { double[,] rnd = new double[nCols, nRows]; int errorState = 0; G05.G05State g05State = new G05.G05State(generator, subGenerator, out errorState); if (errorState != 0) throw new Exception("g05State error"); // double[] temp = new double[nRows]; for (int i = 0; i < nCols; i++) { G05.g05sk(nRows, mean, var, g05State, temp, out errorState); if (errorState != 0) throw new Exception("g05sk error"); // for (int j = 0; j < nRows; j++) { rnd[i, j] = temp[j]; } } return rnd; } } }
The second implementation is using Sobol sequence for creating a set of low-discrepancy numbers. Class implementation is shown in the picture below. Insert a new class (NAG_sobolSequence) into your project and copyPaste the code into it.
using System; using NagLibrary; // namespace RandomGenerator { // wrapper for NAG sobol sequence generator public class NAG_sobolSequence : IRandom { private const int generatorType = 1; // hard-coded sobol sequence generator type private const int skipFirstItems = 2500; private const int returnOrder = 1; // // estimated coefficients for rational approximations // used when transforming uniform variates to normal private const double a1 = -39.6968302866538; private const double a2 = 220.946098424521; private const double a3 = -275.928510446969; private const double a4 = 138.357751867269; private const double a5 = -30.6647980661472; private const double a6 = 2.50662827745924; // private const double b1 = -54.4760987982241; private const double b2 = 161.585836858041; private const double b3 = -155.698979859887; private const double b4 = 66.8013118877197; private const double b5 = -13.2806815528857; // private const double c1 = -0.00778489400243029; private const double c2 = -0.322396458041136; private const double c3 = -2.40075827716184; private const double c4 = -2.54973253934373; private const double c5 = 4.37466414146497; private const double c6 = 2.93816398269878; // private const double d1 = 0.00778469570904146; private const double d2 = 0.32246712907004; private const double d3 = 2.445134137143; private const double d4 = 3.75440866190742; // // break points for transformation function const double p_low = 0.02425; const double p_high = 1 - p_low; // public double[,] getRandom(int nRows, int nCols) { int errorStatus = 0; int n = 32 * nCols + 7; int[] initializationInformation = new int[n]; double[,] rnd = new double[nCols, nRows]; // // initialize sobol quasi-random numbers generator G05.g05yl(generatorType, nCols, initializationInformation, skipFirstItems, out errorStatus); if (errorStatus != 0) throw new Exception("g05yl error"); // // generate uniformly-distributed quasi-random numbers G05.g05ym(nRows, returnOrder, rnd, initializationInformation, out errorStatus); if (errorStatus != 0) throw new Exception("g05ym error"); // // transformation into normally-distributed quasi-random numbers for (int i = 0; i < nCols; i++) { for (int j = 0; j < nRows; j++) { rnd[i, j] = normsinv(rnd[i, j]); } } // return result matrix return rnd; } // private double normsinv(double u) { // transform uniformly-distributed number into normal double q, r; double n = 0.0; // // throw an error if a given uniform number is out of bounds if ((u <= 0.0) || (u >= 1.0)) throw new Exception("given uniform number is out of bounds"); // if (u < p_low) { // rational approximation for lower region q = Math.Sqrt(-2.0 * Math.Log(u)); n = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); goto exitPoint; } // if (u <= p_high) { // rational approximation for mid region q = u - 0.5; r = q * q; n = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1); goto exitPoint; } // if (u < 1.0) { // rational approximation for upper region q = Math.Sqrt(-2.0 * Math.Log(1.0 - u)); n = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } exitPoint: return n; } } }
The numbers received from NAG Sobol generator are mapped to uniform space. For normalizing these numbers, I have used an algorithm for computing the inverse normal cumulative distribution function, developed by Peter Acklam. At this point, we have created the parts of the program, in which are creating the actual random numbers. Next, we are going to interface our C# program with Excel by using Excel-DNA.
STEP TWO: Interfacing C# to Excel with Excel-DNA
Interfacing is done by using Excel-DNA. Insert a new class (RandomGenerator) into your project and copyPaste the code below into it. The complete example for this particular scheme has already been fully covered and explained in here.using System; using ExcelDna.Integration; using System.Windows.Forms; // namespace RandomGenerator { public static class RandomGenerator { public static void execute() { try { // create Excel application object dynamic Excel = ExcelDnaUtil.Application; // // use mersenne twister algorithm NAG_mersenneTwister mt = new NAG_mersenneTwister(); double[,] MT_random = mt.getRandom(1000, 2); Excel.Range["_mersenneTwister"] = Excel.WorksheetFunction.Transpose(MT_random); // // use sobol sequence NAG_sobolSequence sobol = new NAG_sobolSequence(); double[,] MT_sobol = sobol.getRandom(1000, 2); Excel.Range["_sobolSequence"] = Excel.WorksheetFunction.Transpose(MT_sobol); } catch (Exception e) { MessageBox.Show(e.Message.ToString()); } } } }
At this point, we have been creating NAG random number generator wrapper classes as implementations of IRandom interface. Moreover, we have created interfacing program, which connects C# program with Excel. In this project, we are practically using Excel only as a platform for data output for C#.
STEP THREE: Excel and VBA
Open a new Excel workbook and set the following named ranges into worksheet (Sheet1).
- Range "_mersenneTwister" (A1:B1000)
- Range "_sobolSequence" (D1:E1000)
Option Explicit ' Private Sub CommandButton1_Click() Application.Run ("execute") End Sub '
Now, while this workbook is still open, doubleClick RandomGenerator.xll file in your \\RandomGenerator\bin\Release folder. After this, xll can be used by Excel and our C# function (execute) is available to be called from VBA program (Application.Run). VBA event handler program will call and start C# program execute, which will create and send the both sets of random numbers into named Excel ranges.
STEP FOUR: Scatter plots
Next, I have plotted the both sets of random numbers into a two-dimensional space. Scatter plots are shown in the picture below. On the left side, we have random numbers created with NAG pseudo-random numbers generator (Mersenne Twister). On the right side, we have random numbers created with NAG Low-discrepancy numbers generator (Sobol sequence). The upper part is showing normalized data and the lower part is showing data mapped back to uniform plane by using Excel NORMSDIST worksheet function.
The picture clearly shows the difference between the two types of generators. The both generators are filling two-dimensional space randomly, but low-discrepancy sequence (on the right side) fills the space more uniformly, avoiding that undesired clustering effect what is appearing on the left side (pseudo-random generator).
Thanks for reading again.
-Mike
ReplyDeleteYou really have an awesome blog. You doing great and I really love it. Thanks for posting. God bless.
Zea
www.imarksweb.org