General Sampling Methods

Computing expectation empirically requires methods to generate iid random variables having a given distribution $F$ . The starting assumption will be that there is a technique to produce iid uniform(0,1) random variables.

Inverse Transform Method

An operator $F$ on real numbers is a distribution function if the following conditions are satisfied:

$F$ is non-negative.
$\lim_{x \rightarrow -\infty} F(x) = 0$
$\lim_{x \rightarrow \infty} F(x) = 1$
$F$ is non-decreasing
$F$ is right-continuous.

When $U$ is a uniform random variable, note that $P[U \le F(x)] = F(x)$ . If $F$ is strictly increasing, then $F^{-1}$ exists and is also strictly increasing. Furthermore, $P[F^{-1} U \le x] = P[U \le F(x)] = F(x)$ . Hence $X = F^{-1} U$ is a random variable with distribution $F$ .

More generally, it is not always the case that the distribution function is invertible. For example, a function with compact support is not invertible. One can define a pseudo-inverse

$I(u) = \inf\{x : u \le F(x)\}$

Since $P[U <= F(x)] = P[I(U) \le x]$ then $I(IU) \sim F$ :

Suppose $u <= F(x)$ . Then $I(u) \le x$ since $I(u)$ is a lower bound. Suppose $I(u) \le x$ . Then for any $\epsilon > 0$ there is a $p < x + \epsilon$ such that $u \le F(p)$ since $I(u)$ is the greatest lower bound. But since $F$ is non-decreasing, $u \le F(p) \le F(x + \epsilon)$ . Since this is true for all $\epsilon > 0$ , it follows that $u \le F(x+)$ and since $F$ is right-continuous, $F(x+) = F(x)$ . Hence $u \le F(x)$ .

Note that $I$ is non-decreasing:

Let $u < v$ . Then when $v \le F(x)$ it follows that $u \le F(x)$ and so $I(u) \le I(v)$ .

and $I = F^{-1}$ when $F$ is strictly increasing:

Since $F$ is strictly increasing, $F^{-1}$ exists and is strictly increasing. The smallest value of $x$ for which the bound $u \le F(x)$ is satisfied is $x = F^{-1}(u)$ .

And so the pseudo-inverse corresponds with the inverse when the inverse exists.

Example : simulating an exponential random variable

Here the density is $f(x) = \lambda e^-{\lambda x} \mathbb{I}(x \ge 0)$ and the distribution is $F(x) = e^{-\lambda x} \mathbb{I}(x \ge 0)$ . The inverse is $x = -\lambda^{-1} \ln(u)$

Example : simulating a random variable with density $f(x) = \sum_{i \ge 0} a_i \delta(x - a_i)$ .

In principle, find the partial sum $\sum_{i <= n} a_i \le u < \sum_{i <= n+1} a_i$ , and set $I(u) = x_n$ . Unless the partial sums have a closed form representation that can be inverted the procedure has an unbounded processing time.

Box-Mueller Method

Consider the joint density of independent standard normal random variables:

$f(x,y) = f(x)f(y) = (2\pi)^{-1} e^{-\frac{1}{2} (x^2 + y^2)}$

Switch from a rectangular to a polar coordinate system:

$x = r \cos \theta, y = r \sin \theta$

$f_{r,\theta}(r,\theta) = f(r\cos\theta,r \sin \theta) \frac{\partial (x,y)}{\partial (r,\theta)} = (2\pi)^{-1} e^{-\frac{1}{2} r^2} r$

Switch from $r$ to $z = r^2$ coordinates:

$f_{z,\theta}(z,\theta) = (2\pi)^{-1} \cdot (1/2) e^{-1/2 d}$

and recognize that this is the product of densities of independent uni(0, $2\pi$ ) and exp(1/2) random variables. Hence for

$X = Z^{1/2} \cos \Theta$

$Y = Z^{1/2} \sin \Theta$

$X,Y\sim$ normal(0,1) when $Z \sim$ exp(1/2) and $\Theta \sim$ uni(0, $2\pi$ )

Since the starting basis of simulation are uni(0,1) random variables, the inverse transform method gives

$X = (-2 \log U_1)^{1/2} \cos 2\pi U_2$

$Y = (-2 \log U_1)^{1/2} \sin 2\pi U_2$

where $U_1,U_2 \sim$ uni(0,1)

This technique is known as Box-Mueller approach. As mentioned in Sheldon, the efficiency can be improved by avoiding explicit computation of trigonometric functions. (to be completed)

Acceptance-Rejection Method

This technique was invented by Von Neumann and samples from a target distribution $f$ by generating samples from another distribution $g$ and then accepting or rejecting the sample if it passes some test.