Monte Carlo methods 

Simulation of random variables with a given distribution
An ideal random number generator returns a sequence of independent uniformly
distributed random variables 

       U1, U2, U3, ....... from Uniform (0,1)

A pseudo-random number generator returns a sequence of almost independent
and approximately uniformly distributed random variables. 

From this sequence,
how to obtain random variables with any given distribution function F(x)?

Example 1 (Bernoulli)
Let U be Uniform (0,1) variable. Define

	      X = 1 if U < p

		= 0 otherwise

Then X has a Bernoulli distribution with the parameter p.

Example 2 (Binomial)
Similarly, generate n independent Bernoulli variables X1, X2, ... Xn.
Then their sum X1 + X2 + ... + Xn has a Binomial (n,p) distribution.


Generating continuous random variables

Theorem:
Let X be a continuous random variable with a cumulative distribution function 

	      F(x) = P { X <= x }

Then: random variable Y = F(X) has a Uniform (0,1) distribution.

Method
This theorem is used for Monte Carlo simulations. In order to generate
a continuous random variable X with a given c.d.f. F(x),

     1. Generate a Uniform random variable U.
     2. Let X = F-1(U). Then X has the desired distribution F.

This is called the Inverse Transform method.

Example 3 (Exponential)
For an Exponential distribution with parameter a,
                       -ax
	   F(x) = 1 - e            

therefore, applying the Theorem above, one can generate

	   X = -(1/a) ln(1-U)

which has an Exponential distribution, if U is Uniform (0,1).

Certainly, 

	   X = -(1/a) ln(U) 
	   
also has Exponential distribution with parameter a, because
(1-U) is also a Uniform(0,1) variable.

Example 4 (Gamma)
It is difficult to invert a c.d.f. of a Gamma distribution, but one can
generate a Gamma distributed random variable by taking a sum of independent
Exponential variables.

Example 5 (Normal)
It is also difficult to invert a Normal c.d.f. There is an alternative way
of generating Normal random variables.

If U1 and U2 are independent Uniform(0,1) variables, then
               _________
	 X1 = V -2 ln U1 sin(2*Pi*U2)
and            _________
	 X2 = V -2 ln U1 cos(2*Pi*U2)

are independent Standard Normal variables. 

Non-standard Normal variables can be obtained by unstandardizing X1 or X2.

Generating discrete random variables

Now, let us extend the Inverse Transform method to discrete distributions.

Let X be a discrete random variable with a probability mass function

	      P(x) = P { X = x }

For simplicity, let us think that X is a nonnegative integer.
Again, start with a Uniform(0,1) variable U. Clearly,

  P {P(0)+P(1)+...+P(m-1)  <  U  <  P(0)+P(1)+...+P(m-1)+P(m)} = P(m)

because for any a and b, P( a < U < b ) = b-a.

Define X to be the smallest m for which P(0)+P(1)+...+P(m) > U.
Then X has the desired distribution P(x), because

    P{X=x} = P{P(0)+...+P(x-1) < U < P(0)+...+P(x)} = P(x).


Example 6 (Geometric)
Of course, one can generate a Geometric(p) random variable by taking
a sequence of Bernoulli(p) variables X1, X2, ... [see example 1]
and letting
		   X = min { n | Xn = 1 } - 1

However, it is easier to use the general algorithm. According to it,

		      X = [ ln(U)/ln(q) ] 
			    
has a Geometric distribution with parameter p.


Example 7 (Poisson)
For the Poisson distribution, there is an alternative method, which requires
a sequence of Standard Uniform variables U1, U2, ... .
                                   -a
     X = min { n | U1*U2*...*Un < e  }

This X will have a Poisson distribution with parameter a.

A Matlab program for generating various random variables
The lines starting with % are comments

% Generate a standard uniform random variable 

rand

% Generate a Bernoulli random variable with the parameter p=0.24, call it X

X = (rand < 0.24);

% ";" in the end suppresses the output.

% Generate a 200-by-1 vector of independent standard uniform random variables 
% Call it U.

U = rand(200,1);

% ";" in the end suppresses the output.

% Generate 200 independent Exponential random variables with the parameter 5

X = - 1/5 * log(U);

% Draw a histogram of X

hist(X)

% Compute the mean of X

mean(X)

% Pretend that X are interarrival times. Compute the vector of arrival times.

S = zeros(200,1);
for n=1:200
   S(n) = sum( X(1:n) );
end;

% Plot the arrival times with dots

plot(S,'.')

% Generate a Binomial random variable X with parameters n=10 and p=0.55

p = 0.55;
n = 10;
X = sum( rand(n,1) < p );

% Generate a Normal(0,1) random variable.

randn

% Here are Matlab programs used for various demonstrations in class:
% Poisson process of arrivals
% Bernoulli and Binomial processes
% Brownian motion
% Central Limit Theorem

% Here is a user-friendly Matlab tutorial. 
And here is a more detailed one, for beginners.