Sampling Methods Sara Wade University of Cambridge
Charles University 8-19 April 2013, Prague
Sara Wade
Sampling Methods
1 / 30
Numerical Integration We want to compute G = Ep [g(X)] =
R
g(x)p(x)dx.
Solution 1: approximate integral by ˆ= G
T X
g(xt )p(xt )∆x,
t=1
where xt lie on an equidistant grid.
Problem: the number of grid points needed grows exponentially with the dimension of X . Sara Wade
Sampling Methods
2 / 30
Monte Carlo Solution 2: approximate integral by T X ˆ= 1 g(xt ), G T
iid
where X t ∼ p
t=1
Sara Wade
Sampling Methods
3 / 30
Monte Carlo
Properties: ˆ = 1 PT E[g(X t )] = G. unbiased: E[G] t=1 T ˆ convergence: G → Ep [g(X)] as T → ∞ (a.s.) under mild conditions (SLLN). ˆ = 1 V(g(X)) ← doesn’t depend on the dimension of variance: V(G) T x!
Sara Wade
Sampling Methods
4 / 30
Monte Carlo
Properties: ˆ = 1 PT E[g(X t )] = G. unbiased: E[G] t=1 T ˆ convergence: G → Ep [g(X)] as T → ∞ (a.s.) under mild conditions (SLLN). ˆ = 1 V(g(X)) ← doesn’t depend on the dimension of variance: V(G) T x! But, how do we generate random samples from p(x)?
Sara Wade
Sampling Methods
4 / 30
Probability integral transformation
Solution 1: Probability integral transformation: If X ∼ F , then U = F (X) ∼ Unif(0, 1). → X = F −1 (U ) ∼ F. Note: F −1 is the generalized inverse F −1 (u) = inf x {x : F (x) ≥ u}. We can obtain a sample x from F by sample U ∼ Unif(0, 1), set x = F −1 (u). Problems: need to compute F −1 , which can be expensive or unavailable.
Sara Wade
Sampling Methods
5 / 30
Probability integral transformation
0.4 0.0
0.2
p(x)
0.6
0.8
Ex. truncated normal p(x) ∝ N(x|µ, σ 2 )1(x ≥ 0) sample U ∼ Unif(Φ(−µ/σ), 1), set x = Φ−1 (u)σ + µ.
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●● ● ● ●● ●● ●● ● ● ● ● ●● ●● ● ● ● ● ●● ● ●● ● ● ● ●●
−3
−2
−1
0
1
2
3
x
Sara Wade
Sampling Methods
6 / 30
Rejection sampling Let q(x) be a density (that can be easily sampled from) such that p(x) ≤
q(x) α
∀ x,
for some constant α ∈ (0, 1). Solution 2: rejection sampling algorithm: obtain a sample x from p(x) by 1. sample x∗ from q(x) ∗
) 2. accept x = x∗ with probability α(x∗ ) = min(1, p(x q(x∗ ) α), otherwise go to step 1.
Problems: need to find q(x) that can be easily sampled and such that the number of rejections is small.
Sara Wade
Sampling Methods
7 / 30
Rejection sampling Ex. truncated normal p(x) ∝ N(x|µ, σ 2 )1(x ≥ 0);
q(x) = N(x|µ, σ 2 ),
0.4 0.0
0.2
p(x)
0.6
0.8
1. sample x∗ from N(µ, σ 2 ), 2. set x = x∗ if x∗ ≥ 0, otherwise go to step 1. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●● ● ● ●● ●● ●● ● ● ● ● ●● ●● ● ● ● ● ●● ● ●● ● ● ● ●●
xx x xx x xxxxxxxxxx xxxxxxxxxxxxxxxxx x xxx
−3
−2
−1
0
1
2
3
x
Sara Wade
Sampling Methods
8 / 30
Markov chain Monte Carlo A third solution is to generate a Markov chain {X t }t≥1 with an appropriately defined transition probability such that for large enough t, X t is approximately sampled from p(x). Markov chain (first order): over the state space X , is a sequence of random variables X1 , X2 , . . . such that π(xt+1 |x1 , . . . , xt ) = π(xt+1 |xt ). Homogeneous: a Markov chain is homogeneous if the transition probabilities do not depend on t, π(xt+1 |xt ) ← does not depend on t.
Sara Wade
Sampling Methods
9 / 30
Markov chain Monte Carlo
Stationary distribution: p(x) is stationary with respect to a homogeneous Markov chain if X π(x|xt )p(xt ). p(x) = xt
A sufficient condition to ensure that p(x) is a stationary distribution is the detailed balance condition: p(x)π(x0 |x) = p(x0 )π(x|x0 ).
Sara Wade
Sampling Methods
10 / 30
Markov chain Monte Carlo If {X t }t≥1 is a homogeneous Markov chain such that state space is the support of p(x), irreducible and aperiodic, has stationary distribution p(x), d
from ergodic theory, X t → X, where X has density p(x) and T X 1 g(xt ) → Ep [g(X)] a.s. T − T0 t=T0
Thus, our Markov chain produces (dependent) samples from p can be used to approximate G. → We generate a large Markov chain and discard samples to obtain nearly independent samples from p(x). Sara Wade
Sampling Methods
11 / 30
Metropolis-Hastings How to construct such a chain? Such a chain can be obtained by introducing a proposal density q(x|xt ) and choose a starting value x0 . At each iteration: generate x∗ from q(x|xt ), with acceptance probability p(x∗ ) q(xt |x∗ ) ∗ t α(x |x ) = min 1, , p(xt ) q(x∗ |xt ) set xt+1 = x∗ , otherwise set xt+1 = xt . The transition probability: π(xt+1 |xt ) = q(xt+1 |xt )α(xt+1 |xt ) + r(xt )δxt (xt+1 ), R where r(xt ) = q(x|xt )α(x|xt )dx, satisfies the detailed balance condition with respect to p(x). For large enough T0 , (xt )t≥T0 are approximate samples from p(x). Sara Wade
Sampling Methods
12 / 30
Independent Metropolis Hasting The proposal density q(x) is independent of xt and such that p(x) ≤
q(x) , α
for some constant α ∈ (0, 1). The acceptance probability reduces to p(x∗ ) q(xt ) ∗ t α(x |x ) = min 1, . p(xt ) q(x∗ ) This algorithm is similar to rejection sampling but produces dependent samples, but it can be shows that on average independent MH accepts more often.
Sara Wade
Sampling Methods
13 / 30
Independent MH Ex. truncated normal p(x) ∝ N(x|µ, σ 2 )1(x ≥ 0);
q(x) = N(x|µ, σ 2 ),
0.4 0.0
0.2
p(x)
0.6
0.8
sample x∗ from N(µ, σ 2 ), set xt+1 = x∗ if x∗ ≥ 0, otherwise set xt+1 = xt . ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ●● ●● ● ● ●● ● ● ●● ●●● ●
x x
−3
−2
x xx xxxxxxx xx xxx● x ● −1
0
●
●
1
2
●
3
x
Sara Wade
Sampling Methods
14 / 30
Random walk The candidate value is obtained by perturbing the previous value: x∗ = xt + , where is independent of xt with symmetric distribution centered at 0; for ex. ∼ N(0, σ 2 ). The acceptance probability reduces to p(x∗ ) . α(x |x ) = min 1, p(xt ) ∗
t
The scale of the proposal should be as large as possible without incurring high rejection rates.
Sara Wade
Sampling Methods
15 / 30
Random walk Ex. bivariate normal p(x) = N(x|µ, Σ), q(x|xt ) = N(x1 |xt1 , σq2 )N(x2 |xt2 , σq2 ),
sample x∗ from N(x1 |xt1 , σq2 )N(x2 |xt2 , σq2 ) , sample U ∼ Unif(0, 1) and set ∗ x if u < min(1, N(x∗ |µ, Σ)/N(xt |µ, Σ)) t+1 x = . t x otherwise Sara Wade
Sampling Methods
16 / 30
Gibbs sampling
Assume x = (x1 , . . . , xd ). The Gibbs sampling algorithm iteratively samples xi given x−i = (x1 , . . . , xi−1 , xi+1 , . . . , xn ). Let p(xi |x−i ) be the full conditional density. The chain is initialized at x0 and the Gibbs sampling algorithm procedes by sample xt+1 from p(x1 |xt2 , . . . , xtd ), 1 t t sample xt+1 from p(x2 |xt+1 2 1 , x3 , . . . , xd ), .. . t+1 sample xt+1 from p(xd |xt+1 1 , . . . , xd−1 ), d
The Gibbs sampling procedure is a special case of MH where the proposal is always accepted.
Sara Wade
Sampling Methods
17 / 30
Gibbs sampling Ex. bivariate normal p(x) = N(x|µ, Σ): 2 ), sample xt+1 from N(x1 |β0,1 + β1,1 xt2 , σ1|2 1 2 sample xt+1 from N(x2 |β0,2 + β1,2 xt+1 2 1 , σ2|1 ).
Strong correlations can slow down Gibbs sampling. If groups of variables can be updated jointly, we can improve the mixing; this is known as blocked Gibbs sampling. Sara Wade
Sampling Methods
18 / 30
Convergence of MCMC Samples need to be discarded so that (xt )Tt=1 are independent samples from p(x): burn in: to overcome poor initialization, we may need to throw away the first T0 samples. thinning: to reduce dependency, we may need to keep only every k th sample. It is difficult to monitor convergence, but we can determine lack of convergence via trace plots, autocorrelations, and comparing chains with different initializations.
Sara Wade
Sampling Methods
19 / 30
Summary of MCMC
Advantages: general, applicable to many problems, easy to implement, theoretical guarantees as T → ∞. Disadvantages: can be slow and expensive to compute, difficult to assess convergence.
Sara Wade
Sampling Methods
20 / 30
MCMC in Bayesian inference In Bayesian inference, a closed form for the posterior distribution is often unavailable. MCMC can be used to obtain samples (θt )Tt=1 from the posterior p(θ|x), and T 1X E[g(θ)|x] ≈ g(θt ). T t=1
For example, the posterior distribution can be approximated by P (A|x) ≈
T 1X δθt (A), T t=1
for any measurable set A ⊆ Θ, and the predictive density at x∗ is approximately T 1X p(x∗ |x) ≈ p(x∗ |θt , x). T t=1
Sara Wade
Sampling Methods
21 / 30
Probit regression Model: ind
Yi |xi , w ∼ Bern(p(xi )), p(xi ) = p(Yi = 1|xi , w) = Φ(w0 + w1 xi ). Prior: w ∼ N(µ0 , Σ0 ). Latent model: Introduce a latent variable Y˜i where ind Y˜i |xi , w ∼ N(w0 + w1 xi , 1),
yi |˜ yi , xi , w = 1(˜ yi ≥ 0) Note: if we integrate out Y˜i , we recover the original model. Sara Wade
Sampling Methods
22 / 30
Probit regression - Gibbs sampler
Gibbs sampling algorithm: Initialize w, y ˜ ˆ sample w|˜ y, x ∼ N(w, ˆ Σ), where ˆ = (Σ−1 + X 0 X)−1 , Σ 0 ˆ −1 µ0 + X 0 y). w ˆ = Σ(Σ 0 sample Y˜i |w, xi , yi independently with density 1 N(˜ yi |w0 + w1 xi , 1)(1(˜ yi ≥ 0))yi (1(˜ yi < 0))1−yi . Z
Sara Wade
Sampling Methods
23 / 30
← Trace plot of first 2000 samples for w0 (left) and w1 (right)
0.0
−0.6
0.1
−0.4
0.2
−0.2
0.3
0.0
0.4
0.2
0.5
0.4
0.6
Probit regression
0
500
1000
1500
2000
0
500
1000
1500
2000
1.0 0.5 −0.5
0.0
← Autocorrelation for w0 (left) and w1 (right)
−1.0
−0.5
0.0
Autocorrelation
1.0 0.5
Iterations
−1.0
Autocorrelation
Iterations
0
5
10
15
20 Lag
Sara Wade
25
30
35
0
5
10
15
20
25
30
35
Lag
Sampling Methods
24 / 30
0.6 0.4 0.0
0.2
p(y=1|x)
0.8
1.0
Probit regression
xx xxxxx x x xxxxxxxxxxxxxx xxx x
xx xx xx
x x
0
2
x
x
x
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●● ● ●●●● ●●●●● ●●●●● ● ● ● ●●● ● ●● ●●●●●●● ● ●●●● ● ● ●●● ● ● ●●●●●● ● ●●●●●● ● ●● ● ●●●● ●●● ● ●●● ● ● ● ●● ● ●●
−4
−2
4
x
Posterior mean of p(Y∗ = 1|x∗ , w) at a grid of new values x∗ with 95% pointwise credible intervals (and the true curve in red). The n = 100 data points are represented as crosses for yi = 0 and circles for yi = 1. Sara Wade
Sampling Methods
25 / 30
Logit regression
Model: ind
Yi |xi , w ∼ Bern(p(xi )), exp(w0 + w1 xi ) . p(xi ) = p(Yi = 1|xi , w) = 1 + exp(w0 + w1 xi ) Prior: w ∼ N(µ0 , Σ0 ).
Sara Wade
Sampling Methods
26 / 30
Logit regression - MH random walk MH random walk algorithm: Initialize w sample w∗ |wt with proposal density ˆ −1 )−1 V ), q(w|wt ) = N(w|wt , V (Σ−1 0 +C where Cˆ is variance-covariance matrix of the MLEs and V is a diagonal matrix of tuning parameters. accept wt+1 = w∗ with probability ! n ∗ + w ∗ x ) yi t + wt x ) Y exp(w 1 + exp(w i i 0 1 0 1 α(x∗ |xt ) = min 1, exp(w0t + w1t xi ) 1 + exp(w0∗ + w1∗ xi ) i=1
otherwise set wt+1 = wt . Implemented in MCMClogit function of the R package MCMCpack.
Sara Wade
Sampling Methods
27 / 30
Logit regression Trace of x
0.7
0.5
0.8
Trace of (Intercept)
−0.5
0.3
0.4
0.0
0.5
0.6
← Trace plot of first 2000 samples for w0 (left) and w1 (right)
500
1000
1500
2000
0
1000
(Intercept)
x
1500
2000
0.5 −0.5
0.0
← Autocorrelation for w0 (left) and w1 (right)
−1.0
−0.5
0.0
Autocorrelation
0.5
1.0
Iterations
−1.0
Autocorrelation
500
Iterations
1.0
0
0
5
10
15
20 Lag
Sara Wade
25
30
35
0
5
10
15
20
25
30
35
Lag
Sampling Methods
28 / 30
0.6 0.4 0.0
0.2
p(y=1|x)
0.8
1.0
Logit regression
xxxx xxx xxxxxxxxxxxx xxxx
x xx x x xxxxxx
x
0
2
x x
x x
●●●●●●●●●●●●●● ● ●●●●●●●●●●●● ● ● ● ●●●●●● ● ●●● ● ●●●● ● ●●●●●●●●● ● ●●● ●● ● ● ● ● ●●●●●●●● ● ● ●●●●● ● ● ● ●●●● ● ● ●●● ● ●●● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ●●●●●● ● ●●● ●●●● ●●●
−4
−2
4
x
Posterior mean of p(Y∗ = 1|x∗ , w) at a grid of new values x∗ with 95% pointwise credible intervals (and the true curve in red). The n = 100 data points are represented as crosses for yi = 0 and circles for yi = 1. Sara Wade
Sampling Methods
29 / 30
References
Chib, S. and Greenberg, E. (1995). Understanding the Metropolis-Hasting algorithm. The American Statistician, 49, 4:327-335. Casella, G. and George, E. (1992). Explaining Gibbs sampler. The American Statistician, 46, 3:167-174.
Sara Wade
Sampling Methods
30 / 30
