[PPT] - Understanding MCMC Marcel Lthi, University of Basel Slides based PowerPoint Presentation

SLIDE 1

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Understanding MCMC

Marcel Lüthi, University of Basel

Slides based on presentation by Sandro Schönborn

1

SLIDE 2

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The big picture

Markov chain Equilibrium distribution Distribution 𝑞(𝑦) Metropolis Hastings Algorithm induces converges to samples from is If Markov Chain is a- periodic and irreducable it… … which satisfies detailed balance condition for p(x) … an aperiodic and irreducable

SLIDE 3

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Understanding Markov Chains

3

SLIDE 4

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Markov Chain

Sequence of random variables 𝑌𝑗 𝑗=1

𝑂 , 𝑌𝑗 ∈ 𝑇 with joint distribution

𝑄 𝑌1, 𝑌2, … , 𝑌𝑂 = 𝑄 𝑌1 ෑ

𝑗=2 𝑂

𝑄(𝑌𝑗|𝑌𝑗−1)

Simplifications: (for our analysis)
Discrete state space: 𝑇 = {1, 2, … , 𝐿}
Homogeneous Chain: 𝑄 𝑌𝑗 = 𝑚 𝑌𝑗−1 = 𝑛 = 𝑈𝑚𝑛

Initial distribution Transition probability State space

4

Automatically true if we use computers (e.g. 32 bit floats)

1 2 3 1/3 1/2 1/6 1 1

SLIDE 5

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Example: Markov Chain

Simple weather model: dry (D) or rainy (R) hour
Condition in next hour? 𝑌𝑢+1
State space 𝑇 = {𝐸, 𝑆}
Stochastic: 𝑄(𝑌𝑢+1|𝑌𝑢)
Depends only on current condition 𝑌𝑢
Draw samples from chain:
Initial: 𝑌0 = 𝐸
Evolution: 𝑄 𝑌𝑢+1 𝑌𝑢
Long-term Behavior
Does it converge? Average probability of rain?
Dynamics? How quickly will it converge?

DDDDDDDDRRRRRRRRRRRDDDDDDDDDDD DDDDDDDDDDDDDDDDDDDDDDDDDDDDDD DDDDDDDDDRDD...

5

D R 0.05 0.95 0.2 0.8

SLIDE 6

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Discrete Homogeneous Markov Chain

Formally linear algebra:

Distribution (vector):

𝑄 𝑌𝑗 : 𝒒𝒋 = 𝑄(𝑌𝑗 = 1) ⋮ 𝑄(𝑌𝑗 = 𝐿)

Transition probability (transition matrix):

𝑄 𝑌𝑗 𝑌𝑗−1 : 𝑈 = 𝑄 1 ← 1 ⋯ 𝑄 1 ← 𝐿 ⋮ ⋱ ⋮ 𝑄 𝐿 ← 1 ⋯ 𝑄 𝐿 ← 𝐿 𝑈

𝑚𝑛 = 𝑄 𝑚 ← 𝑛 = 𝑄 𝑌𝑗 = 𝑚 𝑌𝑗−1 = 𝑛

6

SLIDE 7

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Evolution of the Initial Distribution

Evolution of 𝑄 𝑌1 → 𝑄(𝑌2):

𝑄 𝑌2 = 𝑚 = ෍

𝑛∈𝑇

𝑄 𝑚 ← 𝑛 𝑄 𝑌1 = 𝑛 𝒒2 = 𝑈𝒒1

Evolution of 𝑜 steps:

𝒒𝑜+1 = 𝑈𝑜𝒒1

Is there a stable distribution 𝒒∗? (steady-state)

𝒒∗ = 𝑈𝒒∗

A stable distribution is an eigenvector of 𝑈 with eigenvalue 𝜇 = 1

7

SLIDE 8

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Steady-State Distribution: 𝒒∗

It exists:
𝑈 subject to normalization constraint: left eigenvector to eigenvalue 1

෍

𝑚

𝑈

𝑚𝑛 = 1

⇔ 1 … 1 𝑈 = 1 … 1

T has eigenvalue 𝜇 = 1 (left-/right eigenvalues are the same)
Steady-state distribution as corresponding right eigenvector

𝑈𝒒∗ = 𝒒∗

Does any arbitrary initial distribution evolve to 𝒒∗?
Convergence?
Uniqueness?

8

SLIDE 9

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Equilibrium Distribution: 𝒒∗

Additional requirement for 𝑈: Tn 𝑚𝑛 > 0 for 𝑜 > 𝑂0

The chain is called irreducible and aperiodic (implies ergodic)

All states are connected using at most 𝑂0 steps
Return intervals to a certain state are irregular
Perron-Frobenius theorem for positive matrices:
PF1: 𝜇1 = 1 is a simple eigenvalue with 1d eigenspace (uniqueness)
PF2: 𝜇1 = 1 is dominant, all 𝜇𝑗 < 1, 𝑗 ≠ 1 (convergence)
𝒒∗ is a stable attractor, called equilibrium distribution

𝑈𝒒∗ = 𝒒∗

9

SLIDE 10

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Convergence

Time evolution of arbitrary distribution 𝒒0

𝒒𝑜 = 𝑈𝑜𝒒0

Expand 𝒒0 in Eigen basis of 𝑈:

𝑈𝒇𝑗 = 𝜇𝑗𝒇𝑗, 𝜇𝑗 < 𝜇1 = 1, 𝜇𝑙 ≥ |𝜇𝑙+1| 𝒒0 = ෍

𝑗

𝐿

𝑑𝑗𝒇𝑗 𝑈𝒒0 = ෍

𝑗

𝐿

𝑑𝑗𝜇𝑗𝒇𝑗 𝑈𝑜𝒒0 = ෍

𝑗 𝐿

𝑑𝑗𝜇𝑗

𝑜𝒇𝑗 = 𝑑1𝒇1 + 𝜇2 𝑜𝑑2𝒇2 + 𝜇3 𝑜𝑑3𝒇3 + ⋯

10

SLIDE 11

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Convergence (II)

𝑈𝑜𝒒0 = ෍

𝑗 𝐿

𝑑𝑗𝜇𝑗

𝑜𝒇𝑗 = 𝑑1𝒇𝟐 + 𝜇2 𝑜𝑑2𝒇2 + 𝜇3 𝑜𝑑3𝒇3 + ⋯

≈ 𝒒∗ + 𝜇2

𝑜𝑑2𝒇2

We have convergence:

𝑈𝑜𝒒0

𝑜→∞ 𝒒∗

Rate of convergence:

𝒒𝑜 − 𝒒∗ ≈ 𝜇2

𝑜𝑑2𝒇2

= 𝜇2 𝑜 𝑑2

Normalizations: 𝒇1 = 1 σ𝑗 𝑞𝑗

∗ = 1

𝑑1𝒇𝟐 = 𝒒∗

(𝑜 ≫ 1)

11

SLIDE 12

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Example: Weather Dynamics

Rain forecast for stable versus mixed weather:

𝑋

𝑡 = 0.95

0.2 0.05 0.8 stable 𝑋

𝑛 = 0.85

0.6 0.15 0.4 mixed 𝒒∗ = 0.8 0.2 𝒒∗ = 0.8 0.2 Eigenvalues: 1, 0.75 0.75 Eigenvalues: 1, 0.25 0.25 RDDDDDDDDDDDDDDD RDDDRDDDDDDDD... RRRRDDDDDDDDDDDD DDDDDDDDDDDDD... Rainy now, next hours? Rainy now, next hours?

Long-term average probability of rain: 20% 20%

12

D R

SLIDE 13

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Markov Chain: First Results

Aperiodic and irreducible chains are ergodic:

(every state reachable after > 𝑂 steps, irregular return time)

Convergence towards a unique equilibrium distribution 𝒒∗
Equilibrium distribution 𝒒∗
Eigenvector of 𝑈 with eigenvalue 𝜇 = 1:

𝑈𝒒∗ = 𝒒∗

Rate of convergence:

Exponential decay with second largest eigenvalue ∝ 𝜇2 𝑜

13

Only useful if we can design chain with desired equilibrium distribution!

SLIDE 14

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Detailed Balance

Special property of some Markov chains

Distribution 𝑞 satisfies detailed balance if the total flow of probability between every pair of states is equal, (we have a local equilibrium):

𝑄 𝑚 ← 𝑛 𝑞 𝑛 = 𝑄 𝑛 ← 𝑚 𝑞 𝑚

Detailed balance implies: 𝑞 is the equilibrium distribution

𝑈𝒒 𝑚 = ෍

𝑛

𝑈𝑚𝑛𝑞𝑛 = ෍

𝑛

𝑈𝑛𝑚𝑞𝑚 = 𝑞𝑚

Most MCMC methods construct chains which satisfies detailed balance.

14

SLIDE 15

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The Metropolis-Hastings Algorithm

MCMC to draw samples from an arbitrary distribution

16

SLIDE 16

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Idea of Metropolis Hastings algorithm

Design a Markov Chain, which satisfies the detailed balance condition

𝑈𝑁𝐼 𝑦′ ← 𝑦 𝑄 𝑦 = 𝑈𝑁𝐼 𝑦 ← 𝑦′ 𝑄 𝑦′

Ergodicity ensures that chain converges to this distribution

SLIDE 17

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Attempt 1: A simple algorithm

Initialize with sample 𝒚
Generate next sample, with current sample 𝒚

1. Draw a sample 𝒚′ from 𝑅(𝒚′|𝒚) (“proposal”) 2. Emit current state 𝒚 as sample

It’s a Markov chain
Need to choose Q for every P to satisfy detailed balance

𝑅 𝑦′ ← 𝑦 𝑄 𝑦 = 𝑅 𝑦 ← 𝑦′ 𝑄 𝑦′

SLIDE 18

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

Attempt 2: More general solution

Initialize with sample 𝒚
Generate next sample, with current sample 𝒚

1. Draw a sample 𝒚′ from 𝑅(𝒚′|𝒚) (“proposal”) 2. With probability 𝛽(x, x′) emit 𝒚′ as new sample 3. With probability 1 − 𝛽(x, x′) emit 𝑦 as new sample

It’s a Markov chain
Decouples Q from P through acceptance rule a
How to choose a?

SLIDE 19

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

What is the acceptance function a?

𝑈𝑁𝐼 𝑦′ ← 𝑦 𝑄 𝑦 = 𝑈𝑁𝐼 𝑦 ← 𝑦′ 𝑄 𝑦′ 𝑏 𝑦′ 𝑦 𝑅 𝑦′ 𝑦 𝑄 𝑦 = 𝑏 𝑦 𝑦′ 𝑅 𝑦 𝑦′ 𝑄 𝑦′ Case A: x’ = x

Detailed balance trivially satisfied for every a(x’,x)

Case B: 𝑦′ ≠ 𝑦

We have the following requirement

𝑏 𝑦′ 𝑦 𝑏 𝑦 𝑦′ = 𝑅 𝑦 𝑦′ 𝑄 𝑦′ 𝑅 𝑦′ 𝑦 𝑄 𝑦

SLIDE 20

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

What is the acceptance function a?

Requirement: Choose 𝑏(𝑦′|𝑦) such that 𝑏 𝑦′ 𝑦 𝑏 𝑦 𝑦′ = 𝑅 𝑦 𝑦′ 𝑄 𝑦′ 𝑅 𝑦′ 𝑦 𝑄 𝑦

𝑏 𝑦 𝑦′ is probability distribution 𝑏 𝑦 𝑦′ ≤ 1 and 𝑏 𝑦 𝑦′ ≥ 0
Easy to check that:

𝑏 𝑦′ 𝑦 = min 1, 𝑅 𝑦 𝑦′ 𝑄 𝑦′ 𝑅 𝑦′ 𝑦 𝑄 𝑦 satisfies this property.

SLIDE 21

> DEPARTMENT OF MATHEMATICS AND COMPUTER SCIENCE GRAVIS 2018 | BASEL

The big picture

Markov chain Equilibrium distribution Distribution 𝑞(𝑦) Metropolis Hastings Algorithm induces converges to samples from is If Markov Chain is a- periodic and irreducable it… … which satisfies detailed balance condition for p(x) … an aperiodic and irreducable