Today's Topics

• Some classical problems

• Inverse transform algorithm

• Rejection sampling

Problem

Given a probability distribution $\nu$, generate random variables $X_1,\dots ,X_n$ from $\nu$

$\nu$ can be specified in different forms

• Distribution function $F\left(x\right)$

• Density/Mass function $p\left(x\right)$

• Unnormalized density/mass function $q\left(x\right)\propto p\left(x\right)$

• Some general description of the distribution

Why We Need Sampling

Many quantities of interest in statistics can be written as expectations of random variables, expressed as high-dimensional integrals:

$$I=E_X\left[f\right(X\left)\right]=\int f\left(x\right)p\left(x\right)dx,X\sim p\left(x\right).$$

In many cases we are unable to evaluate $I$ exactly, but it is possible to find an unbiased estimator for $I$:

$$\hat{I}=\frac{1}{M}\sum _{i=1}^{M}f\left(X_i\right),X_i\sim p\left(x\right),i=1,\dots ,M.$$

The key challenge is to generate $X_1,\dots ,X_M\sim p\left(x\right)$.

Example - EM Algorithm

To maximize the marginal likelihood $$L(\theta ;X)=\int p(X|Z,\theta )p\left(Z|\theta \right)dZ,$$ EM algorithm consists of two steps:

• Expectation: Define $$Q\left(\theta |{\theta }^{\left(t\right)}\right)=E_{Z|X,{\theta }^{\left(t\right)}}\left[\log L\right(\theta ;X,Z\left)\right]$$

• Maximization: Solve $${\theta }^{(t+1)}=\underset{\theta }{\arg max}Q\left(\theta |{\theta }^{\left(t\right)}\right)$$

Example - EM Algorithm

When $Q\left(\theta |{\theta }^{\left(t\right)}\right)$ has no closed form, it needs to be approximated by Monte Carlo samples from $p(Z|X,{\theta }^{\left(t\right)})$.

This method is typically called Monte Carlo expectation-maximization (MCEM).

Example - Bayesian Statistics

At the core of Bayesian statistics is the sampling from posterior distribution $$p\left(\theta |x\right)\propto p(\theta ,x)=\pi \left(\theta \right)p\left(x|\theta \right).$$

In this case, $p\left(\theta |x\right)$ is given in the unnormalized form, since we only have access to $p(\theta ,x)=\pi \left(\theta \right)p\left(x|\theta \right)$.

Pseudo Random Number

• Strictly speaking, it is unlikely to get "real" random numbers in computers

• What existing algorithms generate are pseudo random numbers

• We omit these technical details

• Instead, we assume there is a "generator" that can produce independent random variables $U_1,U_2,\dots$ that follow the uniform distribution $Unif(0,1)$

• Our target is to generate other random numbers based on $U_1,U_2,\dots$

Classical Problems

• Random permutation

• One-pass random selection

• Generating normal random variables

Random Permutation

Problem: Generate a permutation of $1,\dots ,n$ such that all $n!$ possible orderings are equally likely.

• Widely used in sampling without replacement

• Can select only a subset

Algorithm

A simple and natural method:

• Generate a sequence of uniform random numbers $U_1,\dots ,U_n$

• Sort $U_1,\dots ,U_n$ and record their orders $I_1,\dots ,I_n$

• That is, $U_{I_1}\le \cdots \le U_{I_n}$

• Then $I_1,\dots ,I_n$ is a random permutation of $1,\dots ,n$

• The complexity of this algorithm is $O\left(n\log \right(n\left)\right)$

Algorithm

A better algorithm is as follows ^[1]:

• Here $Int\left(x\right)$ is the largest integer less than or equal to $x$

• This is known as the Fisher-Yates shuffling algorithm

• Complexity is $O\left(n\right)$

Proof Sketch

• For each $k$, $I$ follows a (discrete) uniform distribution on $\left[k\right]:=\{1,2,\dots ,k\}$

• When $k=n$, $P(P_n=i)=1/n$ for $i\in \left[n\right]$

• When $k=n-1$, $P(P_{n-1}=i|P_n=p_n)=1/(n-1)$ for $i\in \left[n\right]\setminus \left\{p_n\right\}$

• When $k=n-2$, $P(P_{n-2}=i|P_n=p_n,P_{n-1}=p_{n-1})=1/(n-2)$ for $i\in \left[n\right]\setminus \{p_n,p_{n-1}\}$

• By reduction, we have $P(P_1=p_1,\dots ,P_n=p_n)=1/n!$, where $\{p_1,\dots ,p_n\}$ is a permutation of $\left[n\right]$

One-Pass Random Selection

Problem: Suppose that there is a sequence of nonnegative numbers $a_1,\dots ,a_n$, and one can access them sequentially. The task is to randomly select an index $i^{\ast }$ such that $P(i^{\ast }=i)=a_i/{\sum }_{j=1}^n{a}_j$, with only one pass of data access.

• Useful when accessing data is expensive

• For example, reading data from hard disks

• Can be extended to multiple independent selections

Algorithm ^[2]

We need to show that $P(i^{\ast }=i)=a_i/S_n$ for $i=1,\dots ,n$, where $S_k={\sum }_{j=1}^k{a}_j$.

Proof Sketch

• Now we attempt to prove a stronger conclusion: after $k$ iterations, we have $P(i^{\ast }=i)=a_i/S_k$ for $i=1,\dots ,k$

• Clearly this holds for $k=1$

• Suppose that the result holds for $k=l$, and then in the $(l+1)$-th iteration:

  • For $i=1,\dots ,l$, $i^{\ast }=i$ means that in the new iteration $i^{\ast }$ is not updated
  • For $i=l+1$, $i^{\ast }=l+1$ means $i^{\ast }$ is updated
  • The probability of updating is $a_{l+1}/S_{l+1}$, which is exactly $P(i^{\ast }=l+1)$
  • Finally, for $i=1,\dots ,l$, $P(i^{\ast }=i)=\left(a_i/S_l\right)\cdot (1-a_{l+1}/S_{l+1})=a_i/S_{l+1}$

Generating Normal Random Variables

Problem: Given a uniform random number generator, simulate random variables $X_1,\dots ,X_n\stackrel{iid}{\sim }N(0,1)$.

• Foundation of many simulation algorithms

• Extensively studied

• Many different algorithms

Box-Muller Transform

1. Generate two independent uniform random numbers $U_1$ and $U_2$

2. Let $Z_1=\sqrt{-2\log \left(U_1\right)}\cos \left(2\pi U_2\right)$ and $Z_2=\sqrt{-2\log \left(U_1\right)}\sin \left(2\pi U_2\right)$

3. Then $Z_1$ and $Z_2$ are two independent $N(0,1)$ random variables

• Requires evaluating functions $\log \left(x\right)$, $\sqrt{x}$, $\sin \left(x\right)$, and $\cos \left(x\right)$

• May be inefficient when lots of random numbers are required

Inverse Transform Algorithm

Using the general inverse transform algorithm (introduced later):

1. Generate uniform random number $U$

2. Set $Z={\Phi }^{-1}\left(U\right)$, where $\Phi \left(x\right)$ is the c.d.f. of $N(0,1)$

• However, computing ${\Phi }^{-1}(\cdot )$ may be difficult and inefficient

Ziggurat Method

Essentially a carefully designed rejection method (introduced later). From [3]:

• Practically one of the most efficient methods to generate normal random numbers

Other Methods

Many other methods exist, see [3]

Implementation

We do a simple comparison between R's built-in rnorm() and various normal random number generators based on the Ziggurat method.

library(dplyr)
library(bench)
library(dqrng)
library(RcppZiggurat)
bench::mark(
    rnorm(1e7), dqrnorm(1e7), zrnorm(1e7),
    min_iterations = 10, max_iterations = 10, check = FALSE
) %>% select(expression, min, median)

## # A tibble: 3 × 3
##   expression          min   median
##   <bch:expr>     <bch:tm> <bch:tm>
## 1 rnorm(1e+07)    417.7ms  420.4ms
## 2 dqrnorm(1e+07)   94.9ms   99.9ms
## 3 zrnorm(1e+07)    73.1ms   75.6ms

Univariate random number generation

• Inverse transform algorithm

• Rejection sampling

Discrete Distribution

Suppose we have a discrete distribution with the following probability mass function:

$$p\left(x_i\right)\equiv P(X=x_i)=p_i,i=0,1,\dots ,\sum _j{p}_j=1$$

The following algorithm ^[1] can be used to generate $X\sim p\left(x\right)$:

Discrete Distribution

This procedure can be compactly expressed as $X=F^{-1}\left(U\right)$, where $F(\cdot )$ is the c.d.f. of $X$, and $F^{-1}(\cdot )$ is the (generalized) inverse c.d.f.:

$$F^{-1}\left(p\right)=inf\{x\in R:p\le F(x\left)\right\}.$$

This is why this method is called the inverse transform algorithm for discrete distributions.

Continuous Distribution

For continuous distributions, the algorithm is more straightforward. To generate $X\sim p\left(x\right)$:

1. Generate uniform random number $U$

2. Set $X=F^{-1}\left(U\right)$, where $F\left(x\right)$ is the c.d.f. of $X$

Proof Sketch

To prove that $X\sim p\left(x\right)$, we only need to show $P(X\le x)=F\left(x\right)$:

$$\begin{array}{cccc}P(X\le x)& =P\left(F^{-1}\right(U)\le x)& & \text{(by algorithm)}\\ & =P(U\le F(x\left)\right)& & \text{(monotonicity of}F(\cdot )\text{)}\\ & =F\left(x\right)& & \text{(c.d.f of uniform distribution)}\end{array}$$

Challenge

• In general, the inverse transform algorithm only applies to univariate distributions (but we have a generalized method called measure transport)

• Also, evaluating $F^{-1}(\cdot )$ may be difficult

Example

• Consider the exponential distribution with mean 1, whose distribution function is given by $F\left(x\right)=1-e^{-x}$

• Using the inverse transform algorithm we have $X=-\log (1-U)$

• Since $1-U$ also follows $Unif(0,1)$, we can simply do $X=-\log \left(U\right)$

• To generate exponential random variable with mean $\theta$, we have $X=-\theta \log \left(U\right)$

Rejection Sampling

Rejection sampling is a general technique for exact sampling. It applies to both discrete and continuous distributions, and is not limited to univariate distributions.

Suppose we want to generate $X\sim f\left(x\right)$, and we have an existing method to sample from another distribution $X\sim g\left(x\right)$. $f\left(x\right)$ and $g\left(x\right)$ can be interpreted as probability mass functions or density functions.

Rejection Sampling

Let $c>0$ be a constant such that $f\left(x\right)/g\left(x\right)\le c$ for all $x\in X$, where $X$ stands for the support of $g\left(x\right)$.

Then the rejection sampling method proceeds as follows:

1. Generate $Y\sim g\left(y\right)$

2. Generate a uniform random number $U$

3. If $U\le f\left(Y\right)/\left[cg\right(Y\left)\right]$, set $X=Y$. Otherwise return to Step 1.

Diagram

From [1]:

Theorem ^[1]

1. The generated random variable $X$ follows $f\left(x\right)$.

2. The number of iterations of the algorithm is a geometric random variable with mean $c$.

Question: What is the best possible value of $c$ ?

Proof Sketch

• $P(\text{accepted}|Y=y)=f\left(y\right)/\left[cg\right(y\left)\right]$

• $P\left(\text{accepted}\right)=\int P(\text{accepted}|Y=y)g\left(y\right)dy=1/c$

• Let $X=(-\infty ,x_1]\times \cdots \times (-\infty ,x_p]$, and then

$$\begin{array}{cc}P(X\in X)& =P(Y\in X,\text{accepted})+(1-1/c)P(X\in X)\\ & ={\int }_{X}P(\text{accepted}|Y=y)g\left(y\right)dy+(1-1/c)P(X\in X)\\ & =\left(1/c\right){\int }_{X}f\left(y\right)dy+(1-1/c)P(X\in X)\end{array}$$

• This shows that $P(X\in X)={\int }_{X}f\left(y\right)dy$

Geometric Interpretation

Consider the set

$$S_M\left(h\right)=\left\{\right(x,z):0\le z\le Mh(x),x\in R^p,z\in R\},$$

where $h$ is a density function, and $M>0$ is a constant.

We show that if $(X,Z)\sim Unif\left(S_M\right(f\left)\right)$, then $X\sim f\left(x\right)$.

Geometric Interpretation

Geometric Interpretation

Let $S=S_M\left(f\right)$ and $X=(-\infty ,x_1]\times \cdots \times (-\infty ,x_p]$, and then by definition,

$$\begin{array}{cc}P(X\in X)& =\frac{vol(S\cap X\times [0,+\infty \left)\right)}{vol\left(S\right)}=\frac{{\int }_X{\int }_0^{Mf\left(x\right)}dzdx}{{\int }_{R^p}{\int }_0^{Mf\left(x\right)}dzdx}\\ & ={\int }_{X}f\left(x\right)dx,\end{array}$$

which proves that $X$ has density $f\left(x\right)$.

Conversely, if $X\sim h\left(x\right)$ and $Z|\{X=x\}\sim Unif(0,Mh(x\left)\right)$, then $(X,Z)\sim Unif\left(S_M\right(h\left)\right)$.

Challenge

• Finding a good proposal distribution $g\left(x\right)$ can be hard

• Finding the constant $c$ is even harder

• Rejection sampling may "waste" a lot of random numbers

• If $f\left(x\right)$ is high-dimensional, the problem would be much more challenging

References