C++ Coding - Random Numbers (Update c++11 VS2012)

The example question for these solutions can be found on my website (click here).

2.1 Random Numbers

To use the new random number generator we need to include the random library, the cmath library for any calculations and also iostream to show results onscreen. We first create a new project with an empty program with the correct libraries, and then declare a variable of type mt19937. This declares a new random number generator, which generates pseudo random sequence of integers defined by the Mersenne Twister algorithm. A computer can only generate a random sequence of integers, but of course we can then take that sequence of integers and convert it to any required distribution. Some of the conversions are simple but others are more complex, luckily we now have inbuilt c++ conversion to all standard distributions (more on this later). This means you will always have to create a generator to pass as an argument to the probability distribution you want to generate.

Once you have declared the generator, the next number in the sequence is called with the operator() function, or as can be seen in the program below.

#include <iostream>
#include <cmath>
#include <random>
using namespace std;

int main()
{
  // declare a random number generator
  // here we choose the merton twister engine (32bit)
  mt19937 rng;
  cout << rng() << endl;
  // output is 3499211612
}

download

We have called this function once here to return the value 3499211612. You will notice that the number is the same every time the program is run. This is because a random number generator always follows a predefined sequence, and the default start point in the sequence is always the same. We can call the next few numbers by writing a loop or by simple copy/pasting the cout line a few times to get

  cout << rng() << endl;
  cout << rng() << endl;
  cout << rng() << endl;
  // output is 
  //   3499211612
  //   581869302
  //   3890346734

download

This repetition in the sequence is in fact extremely useful when bug checking your code as you can eliminate the randomness as the reason your results change. When doing numerical computations it is hardly ever required to start the sequence of with a random start point, since you are only ever interested in what happens on average over a large number of simulations. If you wish to change the starting point of sequence to get a different result or to rerun calculations on the same set of random numbers, then we need to set a seed. You set the seed with the function seed(int), or as follows

  rng.seed(123);
  cout << rng() << endl;

download

Note that the 123 here does not refer the 123rd number in the sequence, but rather this refers to the number 123 in the sequence, and the next result will be generated with that number as input. This means that inputting 124 as the seed will put you in a completely unrelated position in the sequence. Seeds will normally only need to be set once (and only once) in your programs, unless you want to compare results in some way.

There maybe times where you want to set a really random sequence (to generate demonstrations, different results etc) then the best way is to use the type random_device. Declare a variable of that type and you can set any number of random seeds with the operator() function, the next code demonstrates this

  random_device randomSeed; // this uses hardware parameters (clocks etc) to generate a random seed
  rng.seed(randomSeed()); // set the seed from the random device
  cout << rng() << endl;
  // my output this time is 1776944864
  // but it's different for each run!!!

download

From here on in we will be using the default seed which should be consistent across platforms.

For the uniform distribution we use the syntax as shown in this program

#include <iostream>
#include <cmath>
#include <random>
using namespace std;

int main()
{
  // declare a random number generator
  // here we choose the merton twister engine (32bit)
  mt19937 rng;
  // a uniform distribution
  uniform_real_distribution<double> U(0.,1.);
  // get a number from the uniform distribution
  cout << U(rng) << endl;
  return 0;
}

download

To calculate the probability using simulations, note that if u is a random draw from the uniform distribution on (0, 1) then probability and expected value (average) are linked as follows

So to calculate this with simulations we need to count how many times u lands in the interval, and then divide by the total simulations. We simply need then to store N the total number of simulations as an integer, and store the running sum of the payoff which in this case is 1 if 0.25 < u < 0.5 and 0 otherwise. Your code might look a bit like this

#include <iostream>
#include <cmath>
#include <random>
using namespace std;

int main()
{
  // declare a random number generator
  // here we choose the merton twister engine (32bit)
  mt19937 rng;
  // a uniform distribution
  uniform_real_distribution<double> U(0.,1.);
  
  // total number of simulations
  int N=1000;
  // keep track of payoff
  double sum=0;
  for(int i=0;i<N;i++)
  {
    double u=U(rng); // generate the random number 
    if(0.25 < u && u < 0.5) // split the interval condition
    {
      sum += 1;// add in payoff
    }
  }
  
  cout << " P(0.25 < u < 0.5) = " << sum/N << endl;
  return 0;
}

download

Once we are satisfied with the uniform random generator we can move onto the normal random generator. Again this is simply a matter now of declaring the type normal_distribution<double> to enable us to convert random integers into draws from a normal distribution. For the code we get

#include <iostream>
#include <cmath>
#include <random>
using namespace std;

int main()
{
  // declare a random number generator
  // here we choose the merton twister engine (32bit)
  mt19937 rng;
  // a normal distribution with mean 0 and variance 1
  normal_distribution<double> Phi(0.,1.);
  // get a number from the normal distribution
  cout << Phi(rng) << endl;
  return 0;
}

download

In the next part we are asked to build up a histogram plot of the normal distribution over the range of intervals

with i = −10,−9,…, 9, 10 and h = 0.4. Here we have 21 intervals, so we are going to have to store the central position of each of the intervals as well as keep count of the frequency at which a random draw lands in the interval. For simplicity we shall use arrays for storage, with a double for the interval position and integer for counting. The variable declarations should be, noting here that n must be constant to use with arrays

  // input parameters
  int totalRuns=1e6; // total number of simulations
  // interval specification
  const int n=21;  // number of intervals
  double a=-4.2,b=4.2;// start/end points
  double h=(b-a)/n;// fixed interval width
  double x[n];// center of intervals
  
  // local storage
  int counter[n];// number of times landing in interval

download

and you should then initialise the values. Because array indices work starting from 0 it will be better if we index the interval by j=0,1,...,20. Then the centre of the interval is given by the formula

(1)

It is also good practice to reset the value of the counter to zero.

Next we want to generate a random number and then find which interval it sits in, and to do this we need the floor function from the math library. We can rearrange (1) to find that a number phi sits in the j^∗th interval given by

We can test this with some code

  double phi=Phi(rng); // get a number from the normal distribution
  int jStar = floor( (phi-a)/h ); // use the floor function to get j* for the nearest point to x_j* to phi
  cout  << phi << " " << jStar << " " << a+jStar*h << " " << a+(jStar+1)*h << endl;
// output is
//   0.13453 10 -0.2 0.2

download

The interval selected is the 10th interval which corresponds to [-0.2,0.2], and clearly phi=0.13453 is in this interval. When you are coding this up you will need to check that you are within the bounds.

The final solution can be plotted (with gnuplot for example as I have used)

and the code to generate this plot is

#include <iostream>
#include <cmath>
#include <random>
using namespace std;

int main()
{
  // declare a random number generator
  // here we choose the merton twister engine (32bit)
  mt19937 rng;
  // a uniform distribution
  normal_distribution<double> Phi(0,1.);
  // get a number from the normal distribution
  
  // create a histogram and output to screen
  // number of intervals
  const int n=21;
  int totalRuns=1e6;
  // start/end points
  double a=-4.2,b=4.2;
  // fixed interval width
  double h=(b-a)/n;
  // center of intervals
  double x[n];
  // number of times landing in interval
  int counter[n];

  for(int j=0;j<n;j++)
  {
    x[j] = a+h/2.+j*h;  // setup the position of the central points
    counter[j]=0; // reset the counter
  }
  
  for(int run=0;run<totalRuns;run++) // run over all simulations
  {
    double phi=Phi(rng); // get a number from the normal distribution
    int jStar = floor( (phi-a)/h ); // use the floor function to get j* for the nearest point to x_j* to phi
    if(jStar>=n || jStar<0)continue; // if you are outside the grid continue
    counter[jStar]++; // add one to the correct counter
  }
  
  for(int j=0;j<n;j++) // output results
  {
    // output i, midpoint, frequency ( phi \in [ x_j-0.5*h , x_j+0.5*h ] )
    cout << j << " " << x[j] << " " << counter[j] << endl;
  }
// output looks like
// 0 -4 66
// 1 -3.6 291
// 2 -3.2 1022
// 3 -2.8 3340
// 4 -2.4 9043
// 5 -2 22022
// 6 -1.6 44645
// 7 -1.2 77928
// 8 -0.8 115622
// 9 -0.4 146577
// 10 0 158419
// 11 0.4 147059
// 12 0.8 115190
// 13 1.2 77761
// 14 1.6 44979
// 15 2 22082
// 16 2.4 9372
// 17 2.8 3252
// 18 3.2 1015
// 19 3.6 245
// 20 4 50
  return 0;
}

download