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>usingnamespacestd; intmain() {// declare a random number generator// here we choose the merton twister engine (32bit)mt19937 rng; cout <<rng() << endl;// output is 3499211612}

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

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;

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 seedrng.seed(randomSeed());// set the seed from the random devicecout <<rng() << endl;// my output this time is 1776944864// but it's different for each run!!!

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>usingnamespacestd; intmain() {// declare a random number generator// here we choose the merton twister engine (32bit)mt19937 rng;// a uniform distributionuniform_real_distribution<double>U(0.,1.);// get a number from the uniform distributioncout <<U(rng) << endl;return0; }

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>usingnamespacestd; intmain() {// declare a random number generator// here we choose the merton twister engine (32bit)mt19937 rng;// a uniform distributionuniform_real_distribution<double>U(0.,1.);// total number of simulationsint N=1000;// keep track of payoffdouble sum=0;for(int i=0;i<N;i++) { double u=U(rng);// generate the random numberif(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;return0; }

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>usingnamespacestd; intmain() {// declare a random number generator// here we choose the merton twister engine (32bit)mt19937 rng;// a normal distribution with mean 0 and variance 1normal_distribution<double>Phi(0.,1.);// get a number from the normal distributioncout <<Phi(rng) << endl;return0; }

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

// input parametersint totalRuns=1e6;// total number of simulations// interval specificationconstint n=21;// number of intervalsdouble a=-4.2,b=4.2;// start/end pointsdouble h=(b-a)/n;// fixed interval widthdouble x[n];// center of intervals// local storageint counter[n];// number of times landing in interval

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 distributionint jStar =floor( (phi-a)/h );// use the floor function to get j* for the nearest point to x_j* to phicout << phi << " " << jStar << " " << a+jStar*h << " " << a+(jStar+1)*h << endl;// output is// 0.13453 10 -0.2 0.2

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>usingnamespacestd; intmain() {// declare a random number generator// here we choose the merton twister engine (32bit)mt19937 rng;// a uniform distributionnormal_distribution<double>Phi(0,1.);// get a number from the normal distribution// create a histogram and output to screen// number of intervalsconstint n=21; int totalRuns=1e6;// start/end pointsdouble a=-4.2,b=4.2;// fixed interval widthdouble h=(b-a)/n;// center of intervalsdouble x[n];// number of times landing in intervalint counter[n];for(int j=0;j<n;j++) { x[j] = a+h/2.+j*h;// setup the position of the central pointscounter[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 distributionint jStar =floor( (phi-a)/h );// use the floor function to get j* for the nearest point to x_j* to phiif(jStar>=n || jStar<0)continue;// if you are outside the grid continuecounter[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 50return0; }

## No comments:

## Post a Comment