% Monte Carlo Simulation
% Monte Carlo simulations are used in a variety of areas, including risk modeling for securities, particle transport, etc. The basic idea is to use distribution functions for random input variables to generate a single set of inputs. These can then be used to provide a single random output variable. If this process is repeated many times, then one can generate a large set of output samples, which can then be used to produce statistical results for this output variable.
% Random number generators almost always produce numbes with auniform distribution between 0 and 1. That is, repeated calls produce a set of random numbers that are evenly distributed between 0 and 1. If Ru is a number produced by such an algorithm, then a uniform distribution between a and b can be produced with the following algorithm:
%  N=a+R_u \;\left(b-a\right)
%  Many random distributions are modeled as "normal" distributions, which represent the ubiquitous "bell" curve. In this case, we specify both a mean and standard deviation for the distribution. Many random number generators will produce a normal distribution with a mean of 0 and a standard deviation of 1. Using this type of distribution, we can generate a normal distribution with a different standard deviation and mean from:
% N=\sigma R_n +M
% where R_n is a normally distributed random number, M is the desired mean and \sigma is the desired standard deviation.
% In Matlab, random numbers are generated using the rand or randn functions. The rand function produces uniformly distributed random numbers between 0 and 1 and the randn function produces normally distributed random numbes with a mean of 0 and a standard deviation of 1.
% As an example, the commands below calculate the average of 100 random numbers with a uniform distribution between –1 and 1. The correct answer, of course, is zero. The more points we sample in our simulation, the more accurate the result is likely to be. You can run the simulation a few times to see that the result will change each time we sample new values. You can also get a sense how small the answer is for 1000 samples. Try running more samples to see the reduction in error.

list=-1+2*rand(1000,1);
plot(list,'o')
mean(list)
hist(list)

% The first plot shows the randomly samples values and the second plot shows a histogram.Given that our random variable is uniformly distributed, this histogram should be flat.
% As a simple demonstration of the use of Monte Carlo methods for evaluating uncertainty, let's assume that we are interested in some parameter S, which is a function of an uncertain input variable z, according to:
% S=\frac{1}{\sqrt{1+z^2 }}
% If we assume that z varies randomly between 2 and 5, then S will vary randomly between 0.19 and 0.45 (computing these by evaluating S at the extreme values of z), but the distribution will not be uniform. To simulate the probability distribution for S, given the assumed distribution for z, we can do the following:
% Generate a set of random values for z, varying between 2 and 5.
% For each random z value, generate a corresponding value for S using the above formula.
% Generate a probability distribution for S using a histogram.
%
% A simulation of this random output variable, using 10,000 samples, would be as follows:

n=10000;
z=3*rand(n,1)+2;
s=1./sqrt(1+z.^2);
[ncount,x]=hist(s,n/200);
plot(ncount,x,'o')

% Notice here that I generate a list of n random numbers and store them in a vector z. Then I generate a vector of outputs in a single command. This could also be done in a for loop, but that approach is much less efficient than the approach shown above. Also, I have chosen to use 10,000 samples for this simulation. Using more samples will provide more accurate statistics, but at the cost of additional compute time.
% Once you have a set of random outputs, you can analyze them in a variety of ways, as shown below. Once again, Matlab will be substantially faster if you avoid loops.

%average value of S
mean(s)
%variance of S
var(s)
%standard deviation of S
std(s)
%probability that S is greater than 0.3
sum(s>0.3)/n

% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu