Postscript version of this file

STAT 380 Lecture 10

Simulation

Monte Carlo computation of expected value:

To compute ${\rm E}(X)$: do experiment $n$ times to generate $n$ independent observations all with same distribution as $X$.

Get $X_1,\ldots,X_n$.

Use

$$\bar{X} = \frac{1}{n} \sum_1^n X_i$$

as an estimate of ${\rm E}(X)$

To estimate ${\rm E}(g(X))$: use same $X_i$ and compute

$$\frac{1}{n} \sum_1^n g(X_i) \, .$$

Random Number Generation

In practice: random quantities are not used.

Use: pseudo-random uniform random numbers.

They ``behave like'' iid Uniform(0,1) variables.

Many, many generators. One standard kind: linear congruential generators.

Start with $x_0$ an integer in range $0 \le x_0 < m$. Compute

$$x_{n+1} = a x_n+b \mod{m}$$

Here the mod means compute the remainder after division by $m$.

Integers $a$ and $b$ are chosen so that

• the sequence $x_0,\ldots,x_{m-1}$ runs through the set $\{0,\ldots,m-1\}$ (or as much of it is possible, sometimes).

• the sequence of pairs
  $$(x_0,x_1),(x_1,x_2),\ldots, (x_{m-1},x_m)$$
  fills up as much of the square $\{(i,j):1\le i,j \le m-1\}$.

• similarly for triples $(x_k,x_{k+1},x_{k+2})$.

Use

$$U_n = \frac{x_n}{m}$$

as Uniform(0,1) random numbers.

General Random Variables

We now pretend we can generate $U_1,U_2,\ldots$ which are iid Uniform(0,1).

Monte Carlo methods for generating a sample:

• Transformation:
  • Inverse Probability Integral Transform

  • Multivariate (e.g. Box Müller)

• Acceptance Rejection

• Markov Chain Monte Carlo

Inverse (Probability Integral) Transform

Fact: $F$ continuous CDF and $X,U$ satisfy

$$U=F(X)$$

then $U\sim$Uniform(0,1) if and only if $X\sim F$.

Proof: For simplicity: assume $F$ strictly increasing on inverval $(a,b)$ with $F(b)=1$ and $F(a)=0$.

If $U\sim$Uniform(0,1) then

$$P(X \le x) = P(F(X) \le F(x))$$

$$= P(U \le F(x))$$

$$= F(x)$$

Conversely: if $X\sim F$ and $0 < u < 1$ then there is a unique $x$ such that $F(x) = u$

$$P(U \le u) = P(F(X) \le F(x))$$

$$= P(X \le x)$$

$$= F(x)$$

$$= u$$

Application: generate $U$ and solve $F(X)=U$ for $X$ to get $X=F^{-1}(U)$ which has cdf $F$.

Example:For the exponential distribution

$$F(x) = 1-e^{-\lambda x}$$

Set $F(X)=U$ and solve to get

$$X=-\log(1-U)/\lambda$$

Observation: $U \sim 1-U$ so $X=-\log(U)/\lambda$ also has an Exponential($\lambda$) distribution.

Example:: for $F$ the standard normal cdf solving

$$F(X) = U$$

requires numerical approximation but such solutions are built in to many languages.

Special purpose transformations:

Example:: If $X$ and $Y$ are iid $N(0,1)$ define $R\ge 0$ and $\Theta\in[0,2\pi)$ by

$$X=R\cos(\Theta) \qquad Y=R\sin(\Theta)$$

NOTE: Book says

$$\Theta = \tan^{-1}(Y/X)$$

but this takes values in $(-\pi/2,\pi/2)$; notation

$$\tan^{-1}(x,y)$$

means angle $\theta$ in $[0,2\pi)$ so that

$$\frac{x}{\sqrt{x^2+y^2}} = \cos(\theta) \qquad \frac{y}{\sqrt{x^2+y^2}} = \sin(\theta) \, .$$

Then

$$\{(r,\theta) : r \le r^*, \theta \le \theta^*\}$$

is part of circle of radius $r^*$. (Start at $(r^*,0)$ and go clockwise to angle $\theta^*$.)

Compute joint cdf of $R,\Theta$ at $r^*,\theta^*$.

Define

$$\begin{multline*} C = \{(x,y) : x^2+y^2 \le (r^*)^2,\\ 0 \le \tan^{-1}(x,y) \le \theta^*\} \end{multline*}$$

Then

$$P(R\le r, \Theta \le\theta)$$

$$=\iint_C \frac{e^{-(x^2+y^2)/2}}{2\pi} dx dy$$

Do integral in polar co-ordinates to get

$$\int_0^{r^*}\int_0^{\theta^*} \frac{e^{-r^2/2}}{2\pi} r dr d\theta$$

The $\theta$ integral just gives

$$\theta^*/(2\pi)$$

The $r$ integral then gives

$$1-e^{-(r^*)^2/2}$$

So

$$P(R\le r^*,\Theta \le\theta^*) = \left(1-e^{-(r^*)^2/2}\right) \times \frac{\theta^*}{2\pi}$$

Product of two cdfs so: $R$ and $\Theta$ are independent.

Moreover $R^2$ has cdf

$$P(R^2 \le x) = 1-e^{-x/2}$$

which is Exponential with rate $1/2$.

$\Theta$ is Uniform$(0,2\pi)$.

So: generate $U_1, U_2$ iid Uniform(0,1).

Define

$$\Theta = 2\pi U_2$$

and

$$R=\sqrt{-2\log(U_1)}$$

Put

$$X=R\cos(\Theta) \qquad Y=R\sin(\Theta)$$

You have generated two independent $N(0,1)$ variables.

Acceptance Rejection

If $F(X)=U$ difficult to solve (eg $F$ hard to compute) sometimes use acceptance rejection.

Goal: generate $X\sim F$ where $F^\prime=f$.

Tactic: find density $g$ and constant $c$ such that

$$f(x) \le c g(x)$$

and s.t. can generate $Y$ from density $g$.

Algorithm:

1. Generate $Y$ which has density $g$.

2. Generate $U\sim$Uniform(0,1) independent of $Y$.

3. If $U \le f(Y)/\{cg(Y)\}$ let $X=Y$.

4. If not reject $Y$ and go back to step 1; repeat, generating new $Y$ and $U$ (independently) till you accept a $Y$.

Facts:

• Since you repeat the same experiment each trial: number of trials, $N$, success is Geometric (and sure to be finite).

• The probability of success on a individual trial is

  $$P[U \le f(Y)/\{cg(Y)\}]$$

  $$= \int P[U \le f(y)/\{cg(y)\}\vert Y=y]g(y) dy$$

  $$= \int \frac{f(y)}{cg(y)} g(y) dy$$

  $$= \frac{1}{c} \int f(y) dy$$

  $$= 1/c$$

  
  
  so
  $${\rm E}(N) = 1/c$$

Most important fact: $X\sim f$:

Proof: this is like the old craps example:

We compute

$$P( X \le x) = P[ Y \le x\vert U \le f(Y)/\{cg(Y)\}]$$

(condition on the first iteration where the condition is met).

Condition on $Y$ to get

$$P( X \le x)$$

$$= \frac{\int_{-\infty}^x P[U \le f(Y)/\{cg(Y)\}\vert Y=y]g(y)dy}{ P[U \le f(Y)/\{cg(Y)\}]}$$

$$= \frac{\int_{-\infty}^x \frac{ f(y)}{cg(y)} g(y)dy}{ 1/c}$$

$$= \int_{-\infty}^x f(y) dy$$

$$= F(x)$$

as desired.

Example:Half normal distribution: $X$ has density

$$f(x) = \frac{2e^{-x^2/2}}{\sqrt{2\pi}}1(x>0)$$

(Density of absolute value of $N(0,1)$.)

Use $g(x) = ae^{-ax}1(x>0$. To find $c$ maximize

$$\frac{f(x)}{g(x)} = \frac{e^{-x^2/2+ax}}{a\sqrt{2\pi}}$$

Rewrite as

$$\frac{e^{a^2/2} e^{-(x-a)^2/2}}{a\sqrt{2\pi}}$$

Maximum is at $x=a$ giving

$$c=\frac{e^{a^2/2}}{a\sqrt{2\pi}}$$

Choose $a$ to minimize $c$ (or $\log(c)$); get $a=1$ so

$$c=\frac{e}{\sqrt{2\pi}}$$

Algorithm is then: generate $Y=-\log(U_1)$ and then compute

$$e^{-(Y-1)^2}$$

Generate $U_2$ and if

$$U_2 \le e^{-(Y-1)^2}$$

accept $Y$ otherwise try again.

To generate $N(0,1)$: use a third $U$ to pick a sign at random: negative if $u < 1/2$ otherwise positive.

Discrete Distributions

The inverse transformation method works for discrete distributions.

$X$ has possible values $\{x_1,x_2,\ldots\}$ with probabilities $\{p_1,p_2,\ldots\}$ compute cumulative probabilities

$$P_k = p_1+\cdots+p_k$$

with $P_0=0$.

Generate $U\sim$Uniform(0,1).

Find $k$ such that

$$P_{k-1} < U \le P_k$$

Put $X=x_k$.