Simulation

Recall the strong theorem of large numbers: let be a sequence of independent random variables having a common distribution and let . Then with probability 1,

.

It follows that one way to estimate is to simulate from and form the arithmetic mean of the x’s as the estimate.

In general can estimate

by

where are drawn independently from .

Example:

Suppose a random variable X has a beta distribution with parameters a, b:

if

where is called a Beta function.

We want to evaluate . This is an incomplete beta function which has been tabulated.

Alternatively define

Then

So estimate by simulating from and letting . 

This general procedure is called Monte Carlo Integration.

How do you draw random variables from ?

---

• Underlying many of the methods we will discuss for simulation will be the requirement to generate random number e.g. a random number between 0 and 1.
• Random number generation and assessing random number generators is a major industry in itself. We will just assume that we have some method for doing it.

e.g. then is approximately .

Note the examples in the book

– generating a random permutation

– estimating the number of distinct entries in a large set

Simulating continuous random variables

The Inverse Transformation Method (method 1)

Proposition:

Consider a random variable with CDF . Define a new random variable . Then ~ .

Proof:



Corollary:

Let . Define a random variable where F is a CDF. Then F is the CDF of X.

Proof:

. 

Thus, if we know how to invert , we can simulate X by generating a random u ~ U(0, 1) and letting .

Example:

Simulating a Weibull. For a Weibull random variable X have:

Then is that value of X such that:

X generated in this fashion will have Weibull distribution.

Note:

This method only works when is invertible which will not often be the case.

The rejection method (method 2)

Suppose there is a density which is “close” to the density that we wish to simulate from but it is much easier to simulate from than (e.g. might be gamma and Weibull). Then provided c such that

for all x, we can use to get simulations from .

Note:

must have support at least as big as .

Here is how it works:

Step 1	Simulate Y having density and simulate a random number U.
Step 2	If set . Otherwise return to step 1.

Claim: that the value X has density function .

Proof:

The denominator does not involve x.

But .

Note:

Since , the number of iterations until an acceptance will be geometric with mean c. Thus it is important to choose so that c is small.

Note:

However, a difficulty with the rejection method is that in many applications, c is hard to compute. Often end up choosing c very conservatively large thereby resulting in very high computational costs.

Sampling Importance Resampling (method 3) (Rubin, 1987)

Again, assume there is a density which is close to the density f that we want to simulate from. Then to generate a sample of size n from f, proceed as follows:

1. draw from
2. sample a value x from the set where the probability of sampling each is proportional to
3. sample a second value x using the same procedure, but excluding the already sampled value from the set
4. repeatedly sample without replacement n – 2 more times.

“Proof”:

Each of the n x’s is drawn with possibility:

as



Note:

SIR is subject to the same problems as rejection sampling but has the advantage that it converges provided without explicitly requiring .

Example: (beta distribution)

Suppose want to simulate from , . We chose as our candidate (also called the “envelope” distribution) , .

Then is bounded iff

Use either rejection sampling or SIR.

If then use

Can simulate from by the inverse transformation method:

So if

Methods for simulating Normal Random Variables

1. Sum of 12 uniforms (approximate algorithm)

,

and X is approximately normal by Central Limit Theorem. This is a fairly crude approximation but may be adequate for some applications.

2. Rejection Sampling

First note that if (i.e. Z ~ N(0, 1))

Then letting

,

( )

Try candidate distribution (i.e. exp(1))

Then

( since if which is impossible )

So, the procedure is:

1. generate a realization from (Inverse transform)
2. generate
3. if then accept x. Else go to 1.
4. generate . If .

Note:

This is quite efficient requiring an average of 1.32 iterations per acceptance.

3. Box-Muller

Let X, Y ~ N(0, 1), X, Y independent

Consider a transformation to polar coordinates:

To get the joint distribution of and need Jacobian of the transformation

Since

for

So and are independent. Furthermore .

So proceed by generating and and then setting

resulting in two independent standard normal random variables.

Simulating a beta random variable

Consider the n-th smallest of n + m – 1 random 0-1 numbers.

has a beta distribution with params n of m.

A problem with this approach is that finding the n-th smallest of n + m – 1 random numbers is expensive if n, m are large.