Today's Topics

• Importance sampling

• Measure transport sampler

Importance sampling

Strictly speaking, importance sampling (IS) is not a method to obtain sample points X1,…,Xn that follow the target distribution p(x).

Instead, it is a technique to estimate expectations related to p(x):

μ=EX[f(X)]=∫f(x)p(x)dx,X∼p(x).

A Direct Solution

Of course, one direct method to approximate μ is to generate X1,…,XM∼p(x), and then an unbiased estimator for μ is given by

ˆμ=1MM∑i=1f(Xi),Xi∼p(x), i=1,…,M.

Suppose that we use rejection sampling to get Xi based on a proposal distribution q(x). There are two issues here:

• Rejection sampling discards sample points
• f(x) may be close to zero outside a region A for which P(X∈A) is small

Example

We want to estimate p=P(X>π) and μ=E(X|X>π), where X∼N(0,1).

Naive solution: sample X1,…,XMiid∼N(0,1), and get ˆp=1MM∑i=1I{Xi>π},ˆμ=∑Mi=1Xi⋅I{Xi>π}∑Mi=1I{Xi>π}.

Problem: p is very small (true value ~0.00084), so maybe all Xi's are smaller than π.

Example

est_naive = function(n)
{
    x = rnorm(n)
    p_hat = mean(x > pi)
    mu_hat = sum(x * (x > pi)) / sum(x > pi)
    c(p_hat, mu_hat)
}
set.seed(123)
est_naive(n = 100)

## [1]   0 NaN

Motivation

IS attempts to resolve the previous issues:

• It does not discard or waste any sample points

• Instead, it assigns different weights to each point

• By properly choosing the proposal distribution, it is able to more effectively generate sample points around the "important region" A

Basic Idea

The idea of IS is in fact quite simple. It is based on a straightforward identity:

μ=∫f(x)p(x)dx=∫f(x)p(x)q(x)q(x)dx=Eq(f(X)p(X)q(X)),

where q(x) is another density function that is positive on Rp, and Eq(⋅) denotes the expectation for X∼q(x).

Accordingly, the IS estimate for μ is

ˆμq=1MM∑i=1f(Xi)p(Xi)q(Xi),Xi∼q(x), i=1,…,M.

Theorem ^[2]

Suppose that q(x)>0 whenever f(x)p(x)≠0. Then Eq(ˆμq)=μ, and Varq(ˆμq)=σ2q/n, where

σ2q=∫Q[f(x)p(x)]2q(x)dx−μ2=∫Q[f(x)p(x)−μq(x)]2q(x)dx,

and Q={x:q(x)>0}.

Example

Back to the previous example, note that

p=∫+∞−∞I(x>π)ϕ(x)dx=∫+∞πI(x>π)ϕ(x)q(x)q(x)dx=∫+∞πϕ(x)q(x)q(x)dx,q(x)={0,x≤πe−(x−π),x>π,

where q(x) is a shifted exponential distribution. Similarly,

μ=p−1∫+∞πx⋅ϕ(x)q(x)q(x)dx.

Example

est_is = function(n)
{
    x = rexp(n) + pi
    ratio = exp(dnorm(x, log = TRUE) + x - pi)
    p_hat = mean(ratio)
    mu_hat = mean(x * ratio) / p_hat
    c(p_hat, mu_hat)
}
set.seed(123)
est_is(n = 100)

## [1] 0.0007478989 3.4296517837

Optimal q(x) ^[2]

It can be proved that the optimal proposal distribution q∗(x) is given by q∗(x)=|f(x)|p(x)/Ep(|f(X)|).

Proof: for any density function q(x) such that q(x)>0 when f(x)p(x)≠0,

μ2+σ2q∗=∫[f(x)p(x)]2q∗(x)dx=[Ep(|f(X)|)]2=[Eq(|f(X)|p(X)q(X))]2≤Eq([f(X)p(X)]2[q(X)]2)=∫[f(x)p(x)]2q(x)dx=μ2+σ2q.

Optimal q(x)

This means that to approximate ∫f(x)p(x)dx, IS can be better than the simple Monte Carlo estimator!

Self-Normalized IS

Suppose that we can only compute pu(x)∝p(x) and qu(x)∝q(x), the self-normalized IS estimate is given by

˜μ=∑Mi=1f(Xi)w(Xi)∑Mi=1w(Xi),

where w(x)=pu(x)/qu(x) and Xi∼q(x).

Under mild conditions, ˜μ is a consistent estimator of μ, but in general ˜μ is no longer unbiased.

Recap: Inverse Transform Algorithm

• U∼Unif(0,1)

• g=F−1

• X=g(U)⇒X∼f(x)

• f(x)=F′(x) is the density function

Random Variable Transformation

• Continuous random variable X∼pX(x)

• pX(x) density function

• g:R→R a monotone function

• Define Y=g(X), then its density function is given by pY(y)=pX(g−1(y))|ddyg−1(y)|

• Can extend to multivariate case

Random Vector Transformation

• Continuous random vector X∈Rd, X∼pX(x)

• pX(x) density function

• T:Rd→Rd a diffeomorphism (a smooth mapping with smooth inverse)

• Define Y=T(X), then its density function is given by pY(y)=pX(T−1(y))|det(∇(T−1)(y))|=pX(x)|det(∇T(x))|−1,x=T−1(y)

• ∇T (cf. ∇T−1) is the Jacobian matrix of T (cf. T−1)

Recall: Box-Muller Transform

1. (U1,U2)∼Unif([0,1]2)

2. Let Z1=√−2log(U1)cos(2πU2) and Z2=√−2log(U1)sin(2πU2)

3. Then Z1 and Z2 are two independent N(0,1) random variables

Recall: Box-Muller Transform

• Transformation mapping T(U1U2)=(√−2log(U1)cos(2πU2)√−2log(U1)sin(2πU2))

• Jacobian matrix ∇T(U1U2)=(−cos(θ)U1L−2πLsin(θ)−sin(θ)U1L2πLcos(θ))L=√−2log(U1),θ=2πU2

• Determinant det∇T=−2πU−11

Recall: Box-Muller Transform

• Clearly, pU(u1,u2)=1, so pZ(z1,z2)=pU(u1,u2)|det∇T|−1=(2π)−1u1

• To express u1 using z1 and z2, note that Z21+Z22=−2log(U1)

• Therefore, pZ(z1,z2)=12πe−12(z21+z22), which implies that Z1 and Z2 follow independent N(0,1)

Ideas from Inverse Transform

• Fix a source distribution pZ(z), e.g. N(0,I)

• Given a target distribution pX(x)

• Can we learn a mapping T such that X=T(Z)∼pX(x) if Z∼pZ(z)?

Measure Transport Sampler

• If we can obtain such a mapping, sampling would be easy:

  • Simulate Z1,…,Zn∼N(0,I)

  • Set X1=T(Z1),…,Xn=T(Zn)

  • Then Xiiid∼pX(x)

• The key is to find such a mapping T, also called a transport map

Transport Map

For practical use, the transport map T needs to have some "nice" properties:

• T should be invertible and differentiable

• T−1 and det(∇T(x)) should be easy to compute

• T should be flexible enough to characterize sophisticated nonlinear mappings

Modeling Transport Maps

• Polynomials ^[3] (not good enough)

• Normalizing flows ^[4] (tools from the deep learning community)

Training

See [3] for details.

References