Markov Chain Monte Carlo (MCMC) Inference Seung-Hoon Na Chonbuk - - PowerPoint PPT Presentation

Markov Chain Monte Carlo (MCMC) Inference

Seung-Hoon Na Chonbuk National University

Monte Carlo Approximation

• Generate some (unweighted) samples from the

posterior:

• Use these samples to compute any quantity of

interest

– Posterior marginal: 𝑞(𝑦1|𝐸) – Posterior predictive: 𝑞 𝑧 𝐸 – …

Sampling from standard distributions

• Using the cdf

– Based on the inverse probability transform

Sampling from standard distributions: Inverse CDF

Sampling from standard distributions: Inverse CDF

• Example: Exponential distribution

cdf

Sampling from a Gaussian: Box-Muller Method

• Sample 𝑨1, 𝑨2 ∈ (−1,1) uniformly
• Discard (𝑨1, 𝑨2) that do not inside the unit-circle,

so satisfying  𝑨1

2 + 𝑨2 2 ≤ 1

• 𝑞 𝑨 =

1 𝜌 𝐽(𝑨 𝑗𝑜𝑑𝑗𝑒𝑓 𝑑𝑗𝑠𝑑𝑚𝑓)

• Define

Sampling from a Gaussian: Box-Muller Method

• Pola

lar coo

• ordinat

dinate상에서 sampl pling ng

• r을 normal distribution에 비례하도록 샘플  Exponential 분포
• 𝜄를 uniform하게 샘플  주어진 r에 대해 x,y 를 uniform하게 분포

Sampling from a Gaussian: Box-Muller Method

• 𝑌 ~ 𝑂 0,1

𝑍 ~ 𝑂 0,1

• 𝑞 𝑦, 𝑧 =

1 2𝜌 exp − 𝑦2+𝑧2 2

= 1

2𝜌 exp − 𝑠2 2

• 𝑠2~𝐹𝑦𝑞

1 2

• 𝐹𝑦𝑞 𝜇 =

− log(𝑉 0,1 ) 𝜇

• 𝑠 ∼

−2 log(𝑉 0,1 )

https://theclevermachine.wordpress.com/2012/09/11/sampling-from-the-normal- distribution-using-the-box-muller-transform/

𝑄 𝑉 ≤ 1 − exp(−0.5𝑠2) =𝑄 𝑠2 ≤ −2𝑚𝑝𝑕𝑉 = 𝑄 𝑠 ≤ −2𝑚𝑝𝑕𝑉

• 𝑨 = 𝑠2
• 𝑄 𝑨 = 𝑄 𝑠

𝑒𝑠 𝑒𝑨 = 0.5 𝑄 𝑠

• 𝑄 𝑠 = 2𝑄 𝑨 = 𝑓𝑦𝑞 −

𝑠2 2

Sampling from a Gaussian: Box-Muller Method

• 1. Draw 𝑣1, 𝑣2 ∼ 𝑉 0,1
• 2. Transform to polar rep.

– 𝑠 = −2log(𝑣1) 𝜄 = 2𝜌𝑣2

• 3. Transform to Catesian rep.

– 𝑦 = 𝑠 cos 𝜄 𝑧 = 𝑠 cos 𝜄

Rejection sampling

• when the inverse cdf method cannot be used
• Create a proposal distribution 𝑟(𝑦) which

satisfies

• 𝑁𝑟(𝑦) provides an upper envelope for ෤

𝑞(𝑦)

• Sample 𝑦 ∼ 𝑟(𝑦)
• Sample 𝑣 ∼ 𝑉(0,1)
• If 𝑣 >

෤ 𝑞(𝑦) 𝑁𝑟(𝑦) , reject the sample

• Otherwise accept it

Rejection sampling

Rejection sampling: Proof

• Why rejection sampling works?
• Let
• The cdf of the accepted point:

𝑞(𝑦)

Rejection sampling

• How efficient is this method?
• 𝑄 𝑏𝑑𝑑𝑓𝑞𝑢 = ׬ 𝑟 𝑦 ׬ 𝐽 𝑦, 𝑣 ∈ 𝑇 𝑒𝑣 𝑒𝑦

= න ෤ 𝑞 𝑦 𝑁𝑟 𝑦 𝑟 𝑦 𝑒𝑦 = 1 𝑁 න ෤ 𝑞 𝑦 𝑒𝑦  We need to choose M as small as possible while still satisfying

න

1

𝐽 𝑣 ≤ 𝑧 = 𝑧

෤ 𝑞 𝑦 𝑁𝑟 𝑦

Rejection sampling: Example

• Suppose we want to sample from a Gamma
• When 𝛽 is an integer, i. e. , 𝛽 = 𝑙, we can use:
• But, for non-integer 𝛽, we cannot use this trick

& instead use rejection sampling

Rejection sampling: Example

• Use as a proposal, where
• To obtain 𝑁 as small as possible, check the ratio

𝑞(𝑦)/𝑟 𝑦 :

• This ratio attains its maximum when

𝑞(𝑦) 𝑟(𝑦) ≤ 𝑁 𝑞(𝑦) 𝑁𝑟(𝑦) ≤ 1

Rejection sampling: Example

• Proposal

Rejection Sampling: Application to Bayesian Statistics

• Suppose we want to draw (unweighted) samples

from the posterior

• Use rejection sampling

– Target distribution – Proposal: – M:

• Accepting probability:

Adaptive rejection sampling

• Upper bound the log density with a piecewise

linear function

Importance Sampling

• MC methods for approximating integrals of the

form:

• The idea: Draw samples 𝒚 in regions which have

high probability, 𝑞(𝒚), but also where |𝑔(𝒚)| is large

• Define 𝑔 𝒚 = 𝐽 𝒚 ∈ 𝐹
• Sample from a proposal

than to sample from itself

Importance Sampling

• Samples from any proposal 𝑟(𝒚) to estimate

the integral:

• How should we choose the proposal?

– Minimize the variance of the estimate መ 𝐽

importance weights.

Importance Sampling

• By Jensen’s inequality, we have 𝐹 𝑣2 𝒚

≥ 𝐹 𝑣 𝒚

2

• Setting 𝑣 𝒚 =

𝑞 𝒚 𝑔 𝒚 𝑟 𝒚

= 𝑥 𝒚 |𝑔 𝒚 |, we have the lower bound:

𝑣 𝒚 =

𝑞 𝒚 𝑔 𝒚 𝑟 𝒚

가 상수이면 equality가 성립 𝑟 𝒚 ∝ 𝑞 𝒚 𝑔 𝒚

Importance Sampling: Handling unnormalized distributions

• When only unnormalized target distribution and

proposals are available without 𝑎𝑞, 𝑎𝑟?

• Use the same set of samples to evaluate 𝑎𝑟/ 𝑎𝑟:

unnormalized importance weight

Ancestral sampling for PGM

• Ancestral sampling

– Sample the root nodes, – then sample their children, – then sample their children, etc. – This is okay when we have no evidence

Ancestral sampling

Ancestral sampling: Example

Rejection sampling for PGM

• Now, suppose that we have some evidence with

interest in conditional queries :

• Rejection sampling (local sampling)

– Perform ancestral sampling, – but as soon as we sample a value that is inconsistent with an observed value, reject the whole sample and start again

• However, rejection sampling is very inefficient

(requiring so many samples) & cannot be applied for real-valued evidences

Importance Sampling for DGM: Likelihood weighting

• Likelihood weighting

– Sample unobserved variables as before, conditional on their parents; But don’t sample

• bserved variables; instead we just use their
• bserved values.

– This is equivalent to using a proposal of the form: – The corresponding importance weight:

the set of observed nodes

Likelihood weighting

Likelihood weighting

Sampling importance resampling (SIR)

• Draw unweighted samples by first using

importance sampling

• Sample with replacement where the

probability that we pick 𝒚𝑡 is 𝑥𝑡

Sampling importance resampling (SIR)

• Application: Bayesian inference

– Goal: draw samples from the posterior – Unnormalized posterior: – Proposal: – Normalized weights: – Then, we use SIR to sample from

Typically S’ << S

Particle Filtering

• Simulation based, algorithm for recursive

Bayesian inference

– Sequential importance sampling  resampling

Markov Chain Monte Carlo (MCMC)

• 1) Construct a Markov chain on the state space X

– whose stationary distribution is the target density 𝑞∗(𝒚) of interest

• 2) Perform a random walk on the state space

– in such a way that the fraction of time we spend in each state x is proportional to 𝑞∗(𝒚)

• 3) By drawing (correlated!) samples

𝒚0, 𝒚1, 𝒚2, . . . , from the chain, perform Monte Carlo integration wrt p∗

Markov Chain Monte Carlo (MCMC) vs. Variational inference

• Variational inference

– (1) for small to medium problems, it is usually faster; – (2) it is deterministic; – (3) is it easy to determine when to stop; – (4) it often provides a lower bound on the log liklihood.

Markov Chain Monte Carlo (MCMC) vs. Variational inference

• MCMC

– (1) it is often easier to implement; – (2) it is applicable to a broader range of models, such as models whose size or structure changes depending on the values of certain variables (e.g., as happens in matching problems), or models without nice conjugate priors; – (3) sampling can be faster than variational methods when applied to really huge models or datasets.2

Gibbs Sampling

• Sample each variable in turn, conditioned on the

values of all the other variables in the distribution

• For example, if we have D = 3 variables
• Need to derive full conditional for variable I

Gibbs Sampling: Ising model

• Full conditional

Gibbs Sampling: Ising model

• Combine an Ising prior with a local evidence term 𝜔𝑢

Gibbs Sampling: Ising model

• Ising prior with 𝑋

𝑗𝑘 = 𝐾 = 1

– Gaussian noise model with 𝜏 = 2

Gibbs Sampling: Ising model

Gaussian Mixture Model (GMM)

• Likelihood function
• Factored conjugate prior

Gaussian Mixture Model (GMM): Variational EM

• Standard VB approximation to the posterior:
• Mean field approximation
• VBEM results in the optimal form of 𝑟 𝒜, 𝜾 :

[Ref] Gaussian Models

• Marginals and conditionals of a Gaussian model
• 𝑞 𝒚 ∼ 𝑂(𝝂, 𝚻)
• Marginals:
• Posterior

conditionals:

[Ref] Gaussian Models

• Linear Gaussian systems
• The posterior:
• The normalization constant:

[Ref] Gaussian Models: Posterior distribution of 𝝂

• The likelihood wrt 𝝂:
• The prior:
• The posterior:

[Ref] Gaussian Models: Posterior distribution of 𝚻

• The likelihood form of Σ
• The conjugate prior: the inverse Wishart distribution

The posterior:

Gaussian Mixture Model (GMM): Gibbs Sampling

• Full joint distribution:

Gaussian Mixture Model (GMM): Gibbs Sampling

• The full conditionals:

Gaussian Mixture Model (GMM): Gibbs Sampling

Label switching problem

• Unidentifiability

– The parameters of the model 𝜾, and the indicator functions 𝒜, are unidentifiable – We can arbitrarily permute the hidden labels without affecting the likelihood

• Monte Carlo average of the samples for clusters:

– Samples for cluster 1 may be used for samples for cluster 2

• Labeling switching problem

– If we could average over all modes, we would find 𝐹[𝝂𝑙|𝐸] is the same for all 𝑙

Collapsed Gibbs sampling

• Analytically integrate out some of the unknown

quantities, and just sample the rest.

• Suppose we sample 𝒜 and integrate out 𝜾
• Thus the 𝜾 parameters do not participate in the

Markov chain;

• Consequently we can draw conditionally

independent samples  This will have much lower variance than samples drawn from the joint state space

Rao-Blackwellisation

Collapsed Gibbs sampling

• Theorem 24.2.1 (Rao-Blackwell)

– Let 𝒜 and 𝜾 be dependent random variables, and 𝑔(𝒜, 𝜾) be some scalar function. Then

 The variance of the estimate created by analytically integrating out 𝜾 will always be lower (or rather, will never be higher) than the variance of a direct MC estimate.

Collapsed Gibbs sampling

After integrating out the parameters.

Collapsed Gibbs: GMM

• Analytically integrate out the model parameters

𝝂𝑙, 𝚻𝑙 and 𝝆, and just sample the indicators 𝒜

– Once we integrate out 𝝆, all the 𝑨𝑗 nodes become inter-dependent. – once we integrate out 𝜾𝑙, all the 𝑦𝑗 nodes become inter-dependent

• Full conditionals:

Collapsed Gibbs: GMM

• Suppose a symmetric prior of the form 𝝆 =

𝐸𝑗𝑠(𝜷) and 𝛽𝑙 = 𝛽/𝐿

• Integrating out 𝝆, using Dirichlet-multinormial

formula:

෍

𝑙=1 𝐿

𝑂𝑙,−𝑗 = 𝑂 − 1

[ref] Dirichlet-multinormial

• The marginal likelihood for the Dirichlet-

multinoulli model:

Collapsed Gibbs: GMM

All the data assigned to cluster 𝑙 except for 𝒚𝑗

To compute , we remove 𝒚𝑗’s statistics from its current cluster (namely 𝑨𝑗), and then evaluate 𝒚𝑗 under each cluster’s posterior

• predictive. Once we have picked a new cluster, we add 𝒚𝑗’s statistics to this

new cluster

Collapsed Gibbs: GMM

• Collapsed Gibbs sampler for a mixture model

Collapsed Gibbs vs. Vanilla Gibbs

a mixture of K = 4 2-d Gaussians applied to N = 300 data points 20 different random initializations

Collapsed Gibbs vs. Vanilla Gibbs

logprob averaged over 100 different random initializations. Solid line: the median thick dashed: the 0.25 and 0.75 quantiles, thin dashed: the 0.05 and 0.95 quintiles.

Metropolis Hastings algorithm

• At each step, propose to move from the

current state 𝒚 to a new state 𝒚′ with probability 𝑟(𝒚′|𝒚)

– 𝑟: the proposal distribution – E.g.)

Random walk Metropolis algorithm independence sampler

Metropolis Hastings algorithm

• The acceptance probability if the proposal is

symmetric,

• The acceptance probability if the proposal is

asymmetric

– Need the Hastings correction

Given target distribution

Metropolis Hastings algorithm

• When evaluating α, we only need to know

unnormalized density

Metropolis Hastings algorithm

Metropolis Hastings algorithm

• Gibbs sampling is a special case of MH, using

the proposal:

• Then, the acceptance probability is 100%

Metropolis Hastings: Example

An example of the Metropolis Hastings algorithm for sampling from a mixture of two 1D Gaussians

MH sampling results using a Gaussian proposal with the variance 𝑤 = 1

Metropolis Hastings: Example

An example of the Metropolis Hastings algorithm for sampling from a mixture of two 1D Gaussians

MH sampling results using a Gaussian proposal with the variance 𝑤 = 500

Metropolis Hastings: Example

An example of the Metropolis Hastings algorithm for sampling from a mixture of two 1D Gaussians

MH sampling results using a Gaussian proposal with the variance 𝑤 = 8

Metropolis Hastings: Gaussian Proposals

• 1) an independence proposal
• 2) a random walk proposal

MH for binary logistic regression

Metropolis Hastings: Example

• MH for binary logical regression

Joint posterior of the parameters

Metropolis Hastings: Example

• MH for binary logical regression

– Initialize the chain at the mode, computed using IRLS – Use the random walk Metropolis sampler

Metropolis Hastings: Example

• MH for binary logical regression