The Report Requires the Doing of Certain Tasks and the Answer to Certain Questions. What

The report requires the doing of certain tasks and the answer to certain questions. What is required is normally in one of three forms:

a)Output from running Matlab code (normally Figures (plots) but sometimes numerical data) which you should paste into your document.

b)Provision of Matlab code. When including code based on the provided code please remove all comments in the original code (for brevity and clarity) but include any comments you need to make on the modifications you have made to the code.

c)Written answers

Remember that the report should be your own work, as an individual, and suspected plagiarism/collusion will be investigated.

Introduction to Monte Carlo Simulation

Monte Carlo simulation is a computational approach thatusesrepetitive random sampling to obtain numerical results. It is often used for engineering, physics and mathematics problems but also beyond those worlds and into, for example, financial modelling. It is particularly helpful when using other mathematical methods is either difficult or actually impossible. In this software laboratory we use Monte Carlo simulation to generate draws from input probability distributions in order to draw from or obtain a more complicated output probability distribution/probability density function (pdf). It also has uses in both numerical integration and in optimization (e.g. random search). The Monte Carlo simulation approach was developed in the mid-late 1940s by Stanislaw Ulam, whilst working on the US nuclear weapons programme.The name was coined by Nicholas Metropolis, one of his colleagues, as a secret project code name,as he saw similarities with games of chance in casinos like the Monte Carlo casino in Monaco[1].

Random Numbers and Computing

Truly random numbers obviously cannot be produced by computer programs. They must be provided by something external like radioactive decay. Such sequences are available but are clumsy to use and often not adequate in terms of speed and number. Pseudorandom numbers are a sequence of numbers generated by an algorithm in a way that the numbers look statistically independent and uniformly distributed. This is the main method used in random number generators. Sometimes these are further tailored to meet specific needs, leading to quasi-random numbers.

MATLAB’s Probability and Statistics Capabilities

MATLAB has a phenomenal number and variation of probability and statistics tools (as well as (pseudo)random number generators etc. (for some specific distributions e.g. normal (Gaussian), Poisson) it has an entire Statistics Toolbox). We are deliberately ignoring most of that in order to focus on understanding the Monte Carlo approach properly. In particular we will look at a general approach, originating from a uniform random variable on the interval (0,1), to generating random numbers that follow various distributions (which we can define). This is despite MATLAB having built in generators for some random variables. We will also not particularly exploit MATLAB’s matrix capability on this occasion but will use a for loop and generate our random trials sequentially.

MATLAB Code Provided

You are given code “H62NUM_Monte_Carlo_Script.m” and function “cdfcalc.m”.

• Even though this remains available to you on Moodle you are well advised to keep an unaltered copy of this code where you keep your other MATLAB code, so that you can easily return to it (and make further copies) in order to approach later Tasks.

• There are parts of the main script file provided where you are encouraged to (or need to) edit the code and there are places where you should not make changes. These are indicated by EDITABLE CODE and DO NOT EDIT comments.
• As a general rule what you do in the main loop for Monte Carlo simulation requires some care – you can easily accidentally slow your code down and end up sitting around waiting for an excessive time for a result.
• Remember semi colons (;) to suppress output to screen, unless you are running a very small number of test trials otherwise the continual outputting to screen will really slow down execution
• Do not put anything in the loop that does not require repeated evaluation (so doing a one time calculation for a value that stays constant should be done before the loop)
• When you are not given explicit instructions otherwise it is sensible to try a small number of trials first and build up to a large number – this is useful in debugging and also helps you gauge (very approximately) how long the code will take to run for a larger number of trials
• You don’t have to store the precise results of each trial in memory (except if you are running a very small number of trials and are just checking the operation of the code).
• Remember that Ctrl-C can be used to terminate execution of MATLAB code.

What does the provided code do?

This code is set up to allow the generation of continuous random variables, as long as the appropriate pdf is used in the code. NOTE: It can easily be adapted to deal with discrete random variables, but our focus will be with continuous random variables.

Furthermore it allows transformation of those random variables and the combination of multiple random variables (by various operations including addition and multiplication). So we can generate samples of complicated random variables for which we would struggle mathematically to generate formulae for their pdfs. In fact we can generate lots of these random samples and generate histograms that lead us to numerical approximations of the pdf. In other words we can use this code as an approach to generate pdfs (and, if we want, other information) for complicated random variables without having to do difficult (or sometimes impossible) mathematics.

The samples of our chosen input random variables are obtained by generating uniform random variables on the interval (0,1) which map onto the cdf of our chosen random variable (recall that the cdf is a function that starts at 0 for – and rises to 1 by + (if not before)) and from which we read off an equivalent random variable from our chosen distribution. We do this “reading off” either by a formula (when the formula for the target cdf can be inverted to make the target random variable the subject) or by using interpolation when we have the cdf and the corresponding arguments in the form of 1D arrays (“vectors”). We use linear interpolation for simplicity. More about how a uniform random variable can be used to generate our target random variable is given in the Appendix.

You should look at the code and try to figure out what it is doing. In essence there are the following aspects:

• Set up the number of trials and the number of bins
• We set up the pdfs of any input random variables (the formulae and any parameter values).
• If it is impossible, or difficult, to generate a simple formula that can transform a uniform random variable into the desired random variable, then we instead integrate the desired pdf (or pdfs) to generate a set of points in two 1D arrays that represent the cdf and its arguments, from which we can effectively replace the transformation formula with an interpolation. The function cdfcalc (in file cdfcalc.m)provides this (by doing a very basic trapezoidal numerical integration for simplicity). We will need to specify the range over which the integration is done.
• We set up the range of possible output values for our target random variable and define the bins into which they will be placed, including initializing counters for each bin (count(1) for bin 1 etc.).

…. / ….
1 / 2 / 3 / …. / nobins

lowhigh

binwidth

Notice that

binwidth=(high-low)/nobins(1)

and

bin number = ((bin centre value – bin 1 centre value)/binwidth) + 1 (2)

and a direct consequence of this is that rounding to the nearest integer allows us to calculate

binto increment = round(((random value- bin 1 centre value)/binwidth) + 1)(3)

• The main trial loop is repeated many times. Inside it we generate one or more uniformly distributed random samples[2](using MATLAB’s rand function) and use these (either using a specific formula or by interpolation of the cdf vector) to generate samples from one or more target pdfs. We then operate on and combine these to generate a sample of the output random variable we are attempting to characterize and understand.
• We collect these samples in the bins. Each bin has a counter associated with it and when a sample lands within a bin the counter is incremented. With a very small loss of accuracy due to having a finite bin width you are storing all the information you need by incrementing the relevant bin’s counter. There is no real need to store the results of specific trials in memory though MATLAB’s hist function, which we are not using, plots histograms directly from the full data set held in memory for a given bin number. In some cases (with very large trial numbers – which are desirable for greatest accuracy) this can really slow or sometimes crash your code. Instead, here, the precise values for each trial are overwritten by the values from the next one.
• After the loop has finished we can plot the output “histogram” (effectively a plot of the counter values for the different bin numbers). Strictly this would be a bar graph (which MATLAB’s hist would generate) but we will plot it as a line graph to better see its relationship to the pdf which we generate next.
• We then scale the bin counters (divide by the total number of trials to get the probability per bin, and then divide by the binwidth to get the pdf given that probability in a bin = pdf x binwidth) and return the horizontal axis to the sample space (convert the bin numbers to the appropriate points on the output range i.e. use (2)) in order to plot an approximation of the pdf. This pdf (in the form of the “vectors” (arrays) for its value and for its argument) can be used to obtain probability estimates and values e.g. mean and variance.

We may also do other things during the loop if we are looking for specific output – e.g. we may keep a running total of all the samples, in order to find the mean directly at the end, and/or we may track how often the sample is above or below a threshold value, in order to directly approximate the probability of this happening. However, if we don’t do these things, we can still acquire good estimates of them (and other interesting statistics) from the output “histogram” and corresponding(approximate) pdf. Note that this (approximate) pdf will always be noisy but will look smoother when more trials are used[3].

Continuous Distribution Reminders

For a continuous random variableX we have the expectation of f(X) to be

wherep(x) is the probability density function (pdf) of X.

The cumulative distribution function (cdf) is

The mean of X is given by the expectation of X, i.e.

and the variance of X is given by

where the final step requires a little algebra to see.

Probability density function for popular continuous random variables

Uniform:foraxb (0 otherwise)E(X)=(a+b)/2, V(X)=(b-a)2 /12

(Negative) Exponential: for x0 (0 otherwise) E(X)=-1, V(X)=-2

Gaussian (normal):-x E(X)=µ, V(X)=2

We will also find the cdf for the exponential random variable useful for calculating the necessary transformation from a uniform random variable. For the exponential distribution the cdf is . It is possible to write a formula for the Gaussian (normal) distribution cdf in terms of the error function (erf()) but we will use the interpolation method with Gaussians.

Specific Tasks and Questions

When doing the tasks you may find it helpful to have Microsoft Word (or similar) open in order to paste in your output plots, relevant code, and first attempts at answering questions. For any plots (or other output data) please state your simulation parameters (parameters for the input pdfs, what the target random variable is as a function of the input random variables, and the numbers of trials and bins)

Task 1

Use the code to generate an approximation to a pdf fora uniform random variable, taking values between 0 and 1.

(Use x1=u1; in the code)

Try the following four cases: i) 10 bins, 1000 trials ii) 10 bins, 10000 trials iii) 100 bins, 1000 trials iv) 100 bins, 10000 trials. This may seem futile AS YOU ALREADY KNOW THE FORMULA but this is a necessary preparation – we need to trust that we really do generate the desired pdf before we can consider using it to generate more complicated pdfs. How the plots compare to a plot of the formula (generate separately, e.g. from Command Line)? What is the impact of changing both the number of bins and the number of trials?

Task 2

Use the code to generate an approximation to a pdf for an exponential random variable with a mean (-1) of 1. Choosingthe range requires a little care as this random variable has no maximum value. You will need to impose one by having a range for your bins between 0 and an upper limit of 10.

There is a simple formula that transforms a uniform random variable into an exponential random variable.means you need the lines like x1=-log(1-u1)/lambda; in the code (where log in MATLAB is a natural logarithm).

Try the following four cases, again: i) 10 bins, 1000 trials ii) 10 bins, 10000 trials iii) 100 bins, 1000 trials iv) 100 bins, 10000 trials. Again, this may seem futile but, again, this is a necessary preparation – we need to trust that we really do generate the desired pdf before we can consider using it to generate more complicated pdfs. Again, how do the plots compare to a plot of the formula (generate this separately, e.g. from the Command Line) and what is the impact of changing both the number of bins and the number of trials?

Task 3

Use the code to generate an approximation to a pdf for a standard (mean (µ) of 0, variance (2) of 1) normal (Gaussian) random variable. Again, choosingthe range requires a little care as this random variable has no maximum or minimum value but you will need to impose them by having a range for your bins between some lower and some upper limit. Generally it is sensible to make this symmetric about the mean. Use -6 and +6 respectively here.

There is no simple formula that transforms a uniform random variable into a Gaussian random variable (without using additional functions like the error function) so here is a good opportunity to use the cdf interpolation approach instead. You will need to use the cdfcalc function in the preparation before the loop and you will need to use the interp1 function in obtaining x1 from u1 in the loop. So, for example, x1=interp1(cdf1vec,x1vec,u1);

Try just the following case of 100 bins, 1000000 trials.Again, this may have seemed futile but, again, this is a necessary preparation. Again, how does the pdf plot compare to a plot of the formula (generate this separately, e.g. from the Command Line). You will have a little wait for this one – but there is a point to having many trials…

Task 4a

When we add random variables the mean of the resulting random variable is the sum of the means (i.e. E(X1+X2)=E(X1)+E(X2)) and the resulting variance is the sum of the variances (i.e. V(X1+X2)=V(X1)+V(X2)).

You are going to add the random variablesobtainable from Tasks 1 to 3 (in some cases with small changes) but now you will have to change the range of the output random variable to take into account the addition outcomes.

Provide the following, using 100 bins and 100,000 trials, and stating the ranges used on your input random variables.

i)the (approximate) pdf for adding two uniform random variables (on (0,1))

ii)the (approximate) pdf for adding two standard Gaussian random variable (mean (µ)of 0, variance (2) of 1) [you can use the same preliminary cdf integration for both of these as identical (though independent) random variables]

iii)the (approximate) pdf for adding a uniform random variable (on (0,1)) and an exponential random variable with mean (-1) of 1

iv)the (approximate) pdf for adding a uniform random variable (on (0,1)) and a Gaussian random variable (mean (µ) of 0, variance (2) of 0.0625 (=1/16))

v)the (approximate) pdf for adding a uniform random variable (on (0,1)) and a Gaussian random variable (mean (µ) of 0, variance(2) of 16)

Calculate(after the lab!) the theoretical pdfs mathematically in these relatively simple cases. For bothi) and iii) perform the necessary convolution (mathematically) of the input pdfs to show this.

Task4b

Explain, with mathematics and the provision of a code snippet (i.e. only the relevant code), how you could use the pdf information in pdfvecand bincentrevecto calculate a mean. Perform this calculation (for case iii in Task 4a) and compare with the number expected from your understanding of the theory.(Note: you couldwrite a separate script that does this and run it after running the main script to generate the data- in which case do not clear the workspace (memory) before (or in) running your new script!)

Also enable (uncomment) the running total in the loop (and relevant lines before and after the loop) that allows the mean to be found directly and compare the value found that way.

Task 5

i)Provide the (approximate) pdfs for adding 10 uniform random variables (each on (0,1)) [you have already added 2 in Task 4a]

ii)Provide the (approximate) pdfs for adding a)2, b) 10 exponential random variables each with mean of 1

In each case where you add 10 calculate the mean and variance you would expect the sum random variable to take. Separately plot (e.g. from the Command Line) the Gaussian pdf of equivalent variance and mean to the summed random variable, and discuss how it compares.

Say (without necessarily doing it – optional!) what would happen if you added many Gaussian random variables.