Statistical Inference

Department of Computer Science, University of Manchester

Manchester Centre for AI Fundamentals

COMP 64101 – Reasoning and Learning under Uncertainty

(U of Manchester) Statistical Inference 1 / 46

Topics: Variational Inference (VI)

What is VI? (Murphy, 2023, § 10.1)

Gradient-based VI (Murphy, 2023, § 10.2)
Blackbox VI (Murphy, 2023, § 10.2.3)
Coordinate Ascent VI (Murphy, 2023, § 10.3)
Variational Bayes (Murphy, 2023, § 10.3.3)
More VI Goodies (Murphy, 2023, § 10.4 – 10.7)

(U of Manchester) Statistical Inference 2 / 46

Topics: Monte Carlo (MC) Methods

MC integration (Murphy, 2023, § 11.2)

Importance Sampling (Murphy, 2023, § 11.5)
Rao-Blackwellization (Murphy, 2023, § 11.6.2)

(U of Manchester) Statistical Inference 3 / 46

Topics: Markov Chain Monte Carlo (MCMC)

Metropolis-Hastings Algorithm (Murphy, 2023, § 12.2)

Gibbs Sampling (Murphy, 2023, § 12.3)
Hamiltonian MC (Murphy, 2023, § 12.5.3 – 12.5.4)
MCMC Convergence (Murphy, 2023, § 12.6)

(U of Manchester) Statistical Inference 4 / 46

What is VI?

Variational Inference (VI) reduces posterior inference to optimization

Consider a model with unknown (latent) variables z, known variables
x, and fixed parameters θ

(U of Manchester) Statistical Inference 5 / 46

What is VI?

Variational Inference (VI) reduces posterior inference to optimization

Consider a model with unknown (latent) variables z, known variables
x, and fixed parameters θ
Prior: pθ (z)
Likelihood: pθ (x | z)
(Unnormalized) Joint: pθ (x, z) = pθ (x | z)pθ (z)
Posterior:
pθ (x, z)
pθ (z | x) =
pθ (x)

Computing the normalisation constant

Z
pθ (x) = pθ (x, z)dz

is intractable.
Need to approximate the posterior!
(U of Manchester) Statistical Inference 5 / 46
What is VI?

We approximate the posterior by the distribution q(z) that minimizes

the loss

q ⋆ = arg min DKL (q(z)∥pθ (z | x)), (1)

q∈Q

where DKL (·∥·) is the Kullback-Leibler(KL) divergence between

the two distributions.
The KL divergence is defined as
Z
p(x)
DKL (p∥q) = p(x) log dx
q(x)

In practice, we pick a parametric family Q; call ψ its variational

parameters

(U of Manchester) Statistical Inference 5 / 46

What is VI?

(1) can be rewritten as

ψ ⋆ = arg min DKL (qψ (z)∥pθ (z | x))

ψ
  
pθ (x | z)pθ (z)
= arg min Eqψ (z) log qψ (z) − log
ψ pθ (x)
= arg min Eqψ (z) [log qψ (z) − log pθ (x | z) − log pθ (z)] + log pθ (x)
ψ | {z }
=:L(θ,ψ|x)
R
The term log pθ (x) = log pθ (x, z)dz is generally intractable to
compute.
However, it does not depend on ψ so we can drop it.

(U of Manchester) Statistical Inference 5 / 46

ELBO
So, we want to minimize the variational free energy

L(θ, ψ | x) := Eqψ (z) [log qψ (z) − log pθ (x, z)]

Or equivalently, maximize the evidence lower bound (ELBO)

L (θ, ψ | x) := −L(θ, ψ | x) = Eqψ (z) [log pθ (x, z) − log qψ (z)]

L (θ, ψ | x) ≤ log pθ (x), the evidence

(U of Manchester) Statistical Inference 6 / 46

ELBO First Interpretation

We can write the ELBO as follows

L (θ, ψ | x) = Eqψ (z) [log pθ (x, z)] + H(qψ (z)),

| {z } | {z } | {z }
ELBO expected log joint entropy

where the R(differential) entropy of a distribution is defined as

H(p) = − [log p(x)]p(x)dx

Second term encourages the posterior to be maximum entropy; first

encourages it to be a joint MAP configuration

(U of Manchester) Statistical Inference 7 / 46

ELBO Second Interpretation

We can also write the ELBO as follows

L (θ, ψ | x) = Eqψ (z) [log pθ (x | z) + log pθ (z) − log qψ (z)]

| {z }
ELBO
= Eqψ (z) [log pθ (x | z)] − DKL (qψ (z)∥pθ (z))
| {z } | {z }
expected log likelihood KL from posterior to prior

The KL term acts like a regularizer, preventing the posterior from

diverging too much from the prior

(U of Manchester) Statistical Inference 8 / 46

Choosing the Form of the Variational Posterior

Fixed-form VI: pick a convenient functional form (e.g. multivariate

Gaussian) and optimize ELBO using gradient-based methods

Free-form VI: make the mean field assumption that the posterior
factorizes as qψ (z) = Jj=1 qj (zj )
Q

qj (zj ) ≡ qψj (zj ) is the posterior over the j-th group of variables
No need to specify the functional form for each qj
Optimal distributional form derived by maximizing the ELBO wrt each
group of variational parameters one at a time, in a coordinate ascent
manner

(U of Manchester) Statistical Inference 9 / 46

Multivariate Normal Example

We want to compute the posterior distribution of the mean z of a

multivariate Gaussian distribution.
Assume we know the covariance matrix of the Gaussian distribution.
The likelihood is given as
N
Y
p(D | z) = N (xn | z, Σ)
n=1

The prior is also Gaussian and follows as

p(z|m̆, V̆) = N (z | m̆, V̆)

We want to compute the posterior distribution p(z|D).

Exact posterior can be computed analytically: Multivariate Normal
conjugate prior.

(U of Manchester) Statistical Inference 10 / 46

Multivariate Normal Example

If q uses a full covariance matrix, it matches the exact posterior;

intractable in high dimensions
If q uses a diagonal covariance matrix (mean field approximation):
the approximation is over confident
Mode-seeking nature of minimizing DKL (q∥p)
If q uses a rank-1 plus diagonal approximation: better approximation;
can be computed efficiently
(U of Manchester) Statistical Inference 11 / 46
Parameter Estimation

If the model parameters θ are not known, can we estimate them?

Maximizing the likelihood of the dataset D = {xn }N
n=1 in the presence
N
of latent variables {zn }n=1 ,
N
Y N
Y
p({xn }, {zn } | θ) = p(xn , zn | θ) = p(zn | θ)p(xn | zn , θ)
n=1 n=1

Since the latent variables zn are hidden, we must marginalize them

out to get the log marginal likelihood
 Z 
log p(xn | θ) = log p(xn | zn , θ)p(zn | θ)dzn
| {z }
intractable: normalization constant of the exact posterior

(U of Manchester) Statistical Inference 12 / 46

Parameter Estimation

ELBO is a lower bound on this, L (θ, ψn | xn ) ≤ log p(xn | θ)

Optimize the model parameters by maximizing
N
X
L (θ, {ψn } | D) := L (θ, ψn | xn ) ≤ log p(D | θ)
n=1

Variational EM algorithm (Murphy, 2023, § [Link]): alternate

between
Maximizing the ELBO wrt the variational parameters {ψn } in the E
step, to give us qψn (zn )
Maximizing the ELBO (using the new ψn ) wrt the model parameters θ
in the M step
Faster alternative: Stochastic VI (Murphy, 2023, § 10.1.4)

(U of Manchester) Statistical Inference 12 / 46

Gradient-Based VI

Choose convenient form for qψ (z), e.g. Gaussian for continuous z, or

product of categoricals for discrete z
Optimize the ELBO using gradient based methods
Stochastic Gradient Descend (SGD): Suppose we have a way of
computing an unbiased estimate gt of the gradient of the ELBO, i.e.

E(gt ) = ∇θ L (θ, ψ | x)|θt

Then, we can use it in a gradient descent procedure

θt+1 = θt − ηt gt ,

where ηt is the learning rate or step size (Murphy, 2023, § [Link])

Same for the gradient wrt ψ

(U of Manchester) Statistical Inference 13 / 46

Gradient wrt Generative Parameters ∇θ L (θ, ψ | x)

Easy to compute: push gradients inside the expectation, and use a

single Monte Carlo sample

∇θ L (θ, ψ | x) = ∇θ Eqψ (z) [log pθ (x, z) − log qψ (z)]

= Eqψ (z) [∇θ ( log pθ (x, z) − log qψ (z) )]
≈ ∇θ log pθ (x, zs ), zs ∼ qψ (z)

This is an unbiased estimate of the gradient, so can be used with SGD

(U of Manchester) Statistical Inference 14 / 46

Gradient wrt Inference Parameters ∇ψ L (θ, ψ | x)

We want an estimate of the following gradient

∇ψ L (θ, ψ | x) = ∇ψ Eqψ (z) [log pθ (x, z) − log qψ (z)]

= ∇ψ Eqψ (z) [fθ,ψ (z)] ,

where we have defined fθ,ψ (z) = log pθ (x, z) − log qψ (z).

If we were to exchange the order of the integration and the
differentiation, it can be shown that

∇ψ L (θ, ψ | x) = ∇ψ Eqψ (z) [fθ,ψ (z)]

Z
= fθ,ψ (z)∇ψ qψ (z)dz + Eqψ (z) [∇ψ fθ,ψ (z)]
z

The first term is not guaranteed to be an expectation.

We can not always have an analytical solution for ∇ψ qθ,ψ (z).

(U of Manchester) Statistical Inference 15 / 46

Two alternatives to compute ∇ψ L (θ, ψ | x)

the reparametrization trick Reparametrized VI (RVI)

REINFORCE estimator Blackbox VI (BBVI)

(U of Manchester) Statistical Inference 16 / 46

Reparametrized VI

Key idea: rewrite z ∼ qψ (z) as differentiable (and invertible)

transformation h of ϵ ∼ p(ϵ), does not depend on ψ, i.e.
z = h(ψ, x, ϵ)

(U of Manchester) Statistical Inference 17 / 46

Reparametrized VI

Key idea: rewrite z ∼ qψ (z) as differentiable (and invertible)

transformation h of ϵ ∼ p(ϵ), does not depend on ψ, i.e.
z = h(ψ, x, ϵ)
E.g. z ∼ N (µ, diag(σ)) ⇐⇒ z = µ + ϵ ⊙ σ, ϵ ∼ N (0, I)

(U of Manchester) Statistical Inference 17 / 46

Reparametrized VI

Key idea: rewrite z ∼ qψ (z) as differentiable (and invertible)

transformation h of ϵ ∼ p(ϵ), does not depend on ψ, i.e.
z = h(ψ, x, ϵ)
E.g. z ∼ N (µ, diag(σ)) ⇐⇒ z = µ + ϵ ⊙ σ, ϵ ∼ N (0, I)

Then, Eqψ (z) [fθ,ψ (z)] = Ep(ϵ) [fθ,ψ (z)], with z = h(ψ, x, ϵ).

(U of Manchester) Statistical Inference 17 / 46

Gradient wrt Inference Parameters ∇ψ L (θ, ψ | x)

Hence

∇ψ L (θ, ψ | x) = ∇ψ Eqψ (z) [fθ,ψ (z)]

= ∇ψ Ep(ϵ) [fθ,ψ (z)] = Ep(ϵ) [∇ψ fθ,ψ (z)]
≈ Ep(ϵ) [∇ψ f (h(ψ, x, ϵs ))], ϵs ∼ p(ϵ)

This is an unbiased estimate of the gradient, so can be used with SGD

(U of Manchester) Statistical Inference 17 / 46

Gradient wrt Inference Parameters ∇ψ L (θ, ψ | x)

Since we are now working with ϵ, need to use the change of variables
formula
 
∂z
log qψ (z) = log p(ϵ) − log det
∂ϵ

∂z/∂ϵ is the Jacobian, i.e. the matrix of all its first-order partial
derivatives
We design transformation z = h(ψ, x, ϵ) s.t. this Jacobian is tractable
to compute

Examples: (Murphy, 2023, § 10.2.1)

(U of Manchester) Statistical Inference 17 / 46

Blackbox VI (BBVI)
Suppose we can evaluate ℓ(ψ, z) := log pθ (x, z) − log qψ (z) pointwise,
but may not be able to take gradients of it

To estimate the gradient of the ELBO, we use the score function

estimator. We have

L (θ, ψ | x) = Eqψ (z) [ℓ(ψ, z)] = Eqψ (z) [fθ,ψ (z)]

We compute ∇ψ L (θ, ψ | x) as

∇ψ L (θ, ψ | x) = ∇ψ Eqψ (z) [fθ,ψ (z)]

Z
= ∇ψ fθ,ψ (z)qψ (z)dz

= Eqψ (z) [fθ,ψ (z)∇ψ log qψ (z)]

(U of Manchester) Statistical Inference 18 / 46
Blackbox VI (BBVI)

Compute a Monte Carlo approximation to this

S
V
1X
∇ψ L (θ, ψ | x) = ℓ(ψ, zs )∇ψ log qψ (zs ) ,
ψt S
s=1 ψt
z1 , . . . , zS ∼ qψ (z)

This is an unbiased estimate of the gradient, so can be used with SGD

(U of Manchester) Statistical Inference 19 / 46

Coordinate Ascent VI (CAVI)

Mean field approximation in VI: assume that all latent variables are
independent, qψ (z) = Jj=1 qj (zj ) ≡ Jj=1 qψj (zj )
Q Q

J is the number of hidden variables

ψj are the variational parameters for the j-th distribution

From previous slides, the ELBO becomes

Z J
X
L (θ, ψ | x) = qψ (z) log pθ (x, z)dz + H(qj )
j=1

(U of Manchester) Statistical Inference 20 / 46

Coordinate Ascent VI (CAVI)

From previous slides, the ELBO becomes

Z J
X
L (θ, ψ | x) = qψ (z) log pθ (x, z)dz + H(qj )
j=1

First term decomposes according to Markov properties of the model

This allows us to use a coordinate ascent optimization scheme to

estimate each ψj
Optimize L (θ, ψ | x) wrt each qj , one at a time, keeping others fixed
Convergence is guaranteed since the bound is concave wrt each of the
factors qj
Functional form of the qj ’s need not be specified in advance:
determined by the form of the log joint

CAVI algorithm: (Murphy, 2023, Algorithm 10.4)

(U of Manchester) Statistical Inference 21 / 46
Variational Bayes

Bayesian modeling: treat the parameters θ as latent variables

Goal: approximate parameter posterior p(θ | D) ∝ p(θ)p(D | θ)
No latent variables
Q except for the shared global parameters,
p(θ, D) = p(θ) N n=1 p(xn | θ)
Fit the variational posterior by maximizing the ELBO

L (θ, ψ | D) = Eqψ (θ) [log p(θ, D)] + H(qψθ (θ))

θ

QJ
Assume the variational posterior factorizes as qψθ (θ) = j=1 qψθj (θj )
Update each ψθj via CAVI (Murphy, 2023, Algorithm 10.4)
Examples: (Murphy, 2023, § 10.3.4)
Of possible interest: Variational Bayes EM (Murphy, 2023, § 10.3.5 –
10.3.6)

(U of Manchester) Statistical Inference 22 / 46

More VI Goodies

Improve tightness of ELBO lower bound – reducing KL of our

posterior approximation – if we use more flexible posterior families
Optimizing within more flexible families may be slower, and can incur
statistical error if the sample size is low (Murphy, 2023, § 10.4)

(U of Manchester) Statistical Inference 23 / 46

More VI Goodies

Improve tightness of ELBO lower bound – reducing KL of our

posterior approximation – if we use more flexible posterior families
Optimizing within more flexible families may be slower, and can incur
statistical error if the sample size is low (Murphy, 2023, § 10.4)

Improve quality of posterior approximation: optimize q wrt a bound

that is a tighter approximation to the log marginal likelihood
compared to ELBO (Murphy, 2023, § 10.5)

(U of Manchester) Statistical Inference 23 / 46

More VI Goodies

Problem with lower bound maximization (standard VI): minimizing

DKL (q∥p), which induces zero-forcing behavior
This means that q(z | x) tends to be too compact (over-confident),
to avoid the situation in which q(z | x) > 0 but p(z | x) = 0,
Would incur infinite KL penalty
Zero-forcing can be desirable for multi-modal posteriors (e.g. mixture
models); not so reasonable for unimodal posteriors
Avoid this problem: minimize DKL (p∥q), which is zero-avoiding
Result in broad posteriors, which avoids overconfidence
Expectation propagation: local approximation to DKL (p∥q) (Murphy,
2023, § 10.7)

(U of Manchester) Statistical Inference 23 / 46

Monte Carlo Methods

Monte Carlo methods are a stochastic approach to solving numerical

integration problems

(U of Manchester) Statistical Inference 24 / 46

Monte Carlo Integration

Integration of interest
Z
E[φ(x)] = φ(x)π(x)dx
Rn

φ : Rn → Rm
π(x) is the target distribution, often a posterior

Approximate it via Monte Carlo integration

Ns
1 X
E[φ(x)] ≈ φ(xn ), x1 , . . . , xNs ∼ π(x)
Ns
n=1

(U of Manchester) Statistical Inference 25 / 46

Monte Carlo Integration

Approximate it via Monte Carlo integration

Ns
1 X
E[φ(x)] ≈ φ(xn ), x1 , . . . , xNs ∼ π(x)
Ns
n=1

The function is only evaluated in places where there is non-negligible

probability
Not need to uniformly cover the entire space

Accuracy is independent of the dimensionality of x; only depends on

the number of samples Ns .

The catch: being able to generate the samples xn ∼ π(x).

(U of Manchester) Statistical Inference 26 / 46

Sampling from the Target Distribution

How to sample from π(x)?

One possible method: Rejection Sampling (Murphy, 2023, § 11.4)

Drawbacks in high dimensions (Murphy, 2023, § 11.4.4)

A better one: Importance Sampling (Murphy, 2023, § 11.5)

In general, it is difficult to sample from π(x)

Sample from a proposal distribution q(x)
Adjust for this inaccuracy by associating weights with each sample
We end up with a weighted MC approximation
Ns
X
E[φ(x)] ≈ Wn φ(xn ) (2)
n=1

(U of Manchester) Statistical Inference 27 / 46

Direct Importance Sampling

We assume we can evaluate the normalized target distribution π(x),

but we cannot sample from it

Instead, we sample from the proposal q(x),

Z Z
π(x)
E[φ(x)] = φ(x)π(x)dx = φ(x) q(x)dx
q(x)

The proposal must be non-zero whenever the target is non-zero

the support of q(x) is greater or equal to the support of π(x)

(U of Manchester) Statistical Inference 28 / 46

Direct Importance Sampling

Sample x1 , . . . , xNs ∼ q(x); then,

Ns Ns
1 X π(xn ) 1 X
E[φ(x)] ≈ φ(xn ) = w̃n φ(xn ),
Ns q(xn ) Ns
n=1 n=1

where the importance weights are defined as

π(xn )
w̃n =:
q(xn )

Unbiased estimate of the true mean E[φ(x)]

w̃n π(xn )
Comparing to expression (2), Wn = Ns = Ns ·q(xn ) .

(U of Manchester) Statistical Inference 28 / 46

Self-Normalized Importance Sampling

Say we can only evaluate

R the unnormalized target distribution
γ̃(x) = Z π(x), Z = Rn γ̃(x)dx

E.g. π(x) = p(x | y), γ̃(x) = p(x, y), and Z = p(y)

Self-normalized importance sampling: key idea also approximate

normalization constant Z with importance sampling

Resulting estimate is a ratio of two estimates: biased

However as Ns → ∞, the bias goes to zero, under some weak
assumptions

(U of Manchester) Statistical Inference 29 / 46

Self-Normalized Importance Sampling

Sample x1 , . . . , xNs ∼ q(x); then,

R
n φ(x)γ̃(x)dx
Z
E[φ(x)] = φ(x)π(x)dx = RR
Rn Rn γ̃(x)dx
R h γ̃(x) i
1 PNs
Rn q(x) φ(x) q(x)dx Ns n=1 w̃n φ(xn )
= R h i ≈ 1 PNs ,
γ̃(x) w̃
Rn q(x) q(x)dx N s n=1 n

where w̃n = γ̃(xn )/q(xn )

PNs w̃
From Eq. (2), E[φ(x)] ≈ n=1 Wn φ(xn ), and Wn = PNs n
n′ =1
w̃n′

1 PNs
Approximation to the normalization constant: Z ≈ Ns n=1 w̃n =: Ẑ

(U of Manchester) Statistical Inference 29 / 46

Choosing the Proposal

Performance of importance sampling depends crucially on quality of

proposal distribution q

We need the support of q to be larger than the support of π, but we

do not want it to be “too large”

The support of q should also take into account properties of the

target function φ as well

However, usually the target function φ is unknown or ignored

Try to find a “generally useful” approximation to the target
One way to come up with a good proposal is to learn one, by
optimizing the ELBO

(U of Manchester) Statistical Inference 30 / 46

Controlling Monte Carlo Variance
√
The standard error in a Monte Carlo estimate is O(1/ S)
S ≤ Ns is the number of (independent) samples
It may take many samples to reduce the variance to a sufficiently
small value
Rao and Blackwell come to the rescue!

(U of Manchester) Statistical Inference 31 / 46

Rao-Blackwellization

Consider two rv’s, θ and X ; we want to estimate f¯ = E[f (θ, X )]

Naı̈ve MC approx:
S
1X
fˆMC = f (θs , Xs ),
S
s=1

with

(θ1 , X1 ), . . . , (θS , XS ) ∼ p(θ, X ) i.i.d.

Unbiased, but may have high variance

(U of Manchester) Statistical Inference 32 / 46

Rao-Blackwellization

Suppose we can analytically marginalize out X , provided we know θ,

i.e., we can tractably compute
Z
fθ (θs ) = f (θs , X )p(X | θs )dX = E[f (θ, X ) | θ = θs ]
Rn

PS
Rao-Blackwellized estimator: fˆRB = 1
S s=1 fθ (θs ),
θ1 , . . . , θS ∼ p(θ)

fˆRB is unbiased, and has lower variance than fˆMC

We are now sampling in a reduced dimensional space

(U of Manchester) Statistical Inference 33 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions

(U of Manchester) Statistical Inference 34 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions
Construct a Markov chain on the state space X whose stationary
distribution is the target density π(x) of interest, e.g. posterior.

(U of Manchester) Statistical Inference 34 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions
Construct a Markov chain on the state space X whose stationary
distribution is the target density π(x) of interest, e.g. posterior.
Perform a random walk on the state space so that the fraction of
time we spend in each state x is proportional to π(x)

(U of Manchester) Statistical Inference 34 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions
Construct a Markov chain on the state space X whose stationary
distribution is the target density π(x) of interest, e.g. posterior.
Perform a random walk on the state space so that the fraction of
time we spend in each state x is proportional to π(x)
Initial samples from the chain do not come from the stationary
distribution: discarded

(U of Manchester) Statistical Inference 34 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions
Construct a Markov chain on the state space X whose stationary
distribution is the target density π(x) of interest, e.g. posterior.
Perform a random walk on the state space so that the fraction of
time we spend in each state x is proportional to π(x)
Initial samples from the chain do not come from the stationary
distribution: discarded
Mixing time: time it takes to reach stationarity

(U of Manchester) Statistical Inference 34 / 46

Markov Chain Monte Carlo
Popular method for sampling from high-dimensional distributions
Construct a Markov chain on the state space X whose stationary
distribution is the target density π(x) of interest, e.g. posterior.
Perform a random walk on the state space so that the fraction of
time we spend in each state x is proportional to π(x)
Initial samples from the chain do not come from the stationary
distribution: discarded
Mixing time: time it takes to reach stationarity

(U of Manchester) Statistical Inference 34 / 46

Metropolis-Hastings Algorithm

MH Algorithm: (Murphy, 2023, Algorithm 12.1)

At each step, we propose to move from the current state x to a new

state x′ w.p. q(x′ | x)
q is the proposal distribution
Common proposal distributions: (Murphy, 2023, § 12.3)
Good proposal design depends on the form of target distribution π(x)

The user is free to use any kind of proposal they want, subject to
supp(π) ⊆ ∪x∈X supp(q(· | x))
q is valid if it “covers” the support of the target π: MH flexible method

(U of Manchester) Statistical Inference 35 / 46

Metropolis-Hastings Algorithm

Having proposed a move to x′ , decide whether to accept this

proposal, or to reject it
The long-term fraction of time spent in each state is ∝ π(x)

If the proposal is accepted, the new state is x′ , otherwise same as the

current state x (repeat the sample).

(U of Manchester) Statistical Inference 36 / 46

MH Algorithm: Symmetric Proposal

Proposal is symmetric: q(x′ | x) = q(x | x′ )

′)
 
Acceptance probability: A = min 1, π(x
π(x)

If x′ is more probable than x, we move there, since π(x′ )/π(x) > 1

If x′ is less probable, we may still move there anyway, depending on
the relative probabilities
Instead of greedily moving to only more probable states, we
occasionally allow “downhill” moves to less probable states
This procedure ensures that the fraction of time we spend in each
state x is equal to π(x) (Murphy, 2023, § 12.2.2)

(U of Manchester) Statistical Inference 37 / 46

MH Algorithm: Asymmetric Proposal
Proposal is asymmetric: q(x′ | x) ̸= q(x | x′ )

(U of Manchester) Statistical Inference 38 / 46

MH Algorithm: Asymmetric Proposal
Proposal is asymmetric: q(x′ | x) ̸= q(x | x′ )
Need Hastings correction
π(x′ )q(x | x′ ) π(x′ )/q(x′ | x)
A = min(1, α), α= =
π(x)q(x′ | x) π(x)/q(x | x′ )

Correction needed: the proposal distribution itself (rather than just the
target distribution) might favor certain states

(U of Manchester) Statistical Inference 38 / 46

MH Algorithm: Asymmetric Proposal
Proposal is asymmetric: q(x′ | x) ̸= q(x | x′ )
Need Hastings correction
π(x′ )q(x | x′ ) π(x′ )/q(x′ | x)
A = min(1, α), α= =
π(x)q(x′ | x) π(x)/q(x | x′ )

Correction needed: the proposal distribution itself (rather than just the
target distribution) might favor certain states
When evaluating α, we only need to know the target density up to a
normalization constant
1 (π̃(x′ )/Z )q(x | x′ ) π̃(x′ )q(x | x′ )
π(x) = π̃(x) =⇒ α = =
Z (π̃(x)/Z )q(x′ | x) π̃(x)q(x′ | x)

We can sample from π even if Z is unknown

(U of Manchester) Statistical Inference 38 / 46

MH Algorithm: Asymmetric Proposal
Proposal is asymmetric: q(x′ | x) ̸= q(x | x′ )
Need Hastings correction
π(x′ )q(x | x′ ) π(x′ )/q(x′ | x)
A = min(1, α), α= =
π(x)q(x′ | x) π(x)/q(x | x′ )

Correction needed: the proposal distribution itself (rather than just the
target distribution) might favor certain states
When evaluating α, we only need to know the target density up to a
normalization constant
1 (π̃(x′ )/Z )q(x | x′ ) π̃(x′ )q(x | x′ )
π(x) = π̃(x) =⇒ α = =
Z (π̃(x)/Z )q(x′ | x) π̃(x)q(x′ | x)

We can sample from π even if Z is unknown

Stochastic MH
Very active line of research still nowadays, see e.g. Bieringer et al. (2023)
(U of Manchester) Statistical Inference 38 / 46
Gibbs Sampling

Problems with MH: (i) need to choose the proposal distribution, (ii)
acceptance rate may be low

Gibbs Sampling: MH method exploiting conditional independence to

automatically create a good proposal; acceptance probability 1
(Murphy, 2023, Eq.’s (12.21), (12.22))

(U of Manchester) Statistical Inference 39 / 46

Gibbs Sampling: Algorithm
(0) (0) (0)
We begin with some initial value x(0) = (x1 , x2 , x3 ) sampled from
a joint distribution p on R3
We want the next sample x(1)
(1)
To do so, for each component xj , we condition on x (1) ’s
(1)
components up to xj−1
(0) (0)
Thereafter condition on x(0) ’s components, starting from xj+1 until x3

(1) (0) (0)

x1 ∼ p(x1 | x2 , x3 )
(1) (1) (0)
x2 ∼ p(x2 | x1 , x3 )
(1) (1) (1)
x3 ∼ p(x3 | x1 , x2 )

This immediately generalizes to RD , D > 3 and x(s) , s > 1

(s)
p(xj | x−j ) is called the full conditional for variable j
(U of Manchester) Statistical Inference 40 / 46
Gibbs Sampling as a Specias Case of MH

Gibbs sampling is a special case of MH where we use a sequence of

proposals of the form

qj (x′ | x) = p(xj′ | x−j )I(x′−j = x−j )

We move to a new state where xj is sampled from its full conditional,

but x−j is left unchanged
The acceptance rate of each proposal is 100%,
′
p xi′ | x−i ′ ′




p (x ′ ) qi (x | x ′ ) p x−i p xi | x−i
α= =
p(x)qi (x ′ | x) p (xi | x−i ) p (x−i ) p xi′ | x−i

p (xi′ | x−i ) p (x−i ) p (xi | x−i )

=
= 1,
p (xi | x−i ) p (x−i ) p xi′ | x−i
′ =x .
where we have use the fact that x−i −i

(U of Manchester) Statistical Inference 41 / 46

Gibbs Sampling as a Specias Case of MH

The fact that the acceptance rate is 100% does not necessarily mean
that Gibbs will converge rapidly
Only updates one coordinate at a time

If we can group together correlated variables, then we can sample

them as a group: help mixing (Murphy, 2023, § 12.3)

(U of Manchester) Statistical Inference 42 / 46

Metropolis Within Gibbs

What should we do if we cannot sample from full conditionals?

(s+1) (s+1) (s)
To sample xj ∼ p(xj | x1:j−1 , xj+1:D ), we
(s)
Propose xj′ ∼ q(xj′ | xj )
Compute acceptance probability Aj = min(1, αj ),
(s+1) (s) (s)
p(x1:j−1 , xj′ , xj+1:D )/q(xj′ | xj )
αj = (s+1) (s) (s) (s)
p(x1:j−1 , xj , xj+1:D )/q(xj | xj′ )

Sample u ∼ U(0, 1), and set

(
(s+1) xj′ if u < Aj
xj = (s)
xj if u ≥ Aj

(U of Manchester) Statistical Inference 43 / 46

Hamiltonian Monte Carlo (HMC)

Many MCMC algorithms perform poorly in high dimensional spaces,

because they rely on a form of random search based on local
perturbations
Hamiltonian Monte Carlo leverages gradient information to guide the
local moves
Important in applications
Algorithm: (Murphy, 2023, § 12.5.3)
Tuning: (Murphy, 2023, § 12.5.4)
Implementation: BlackJAX

(U of Manchester) Statistical Inference 44 / 46

MCMC Convergence

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence

We start MCMC from an arbitrary initial state

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence

We start MCMC from an arbitrary initial state

The samples will be coming from the chain’s stationary distribution π
only when the chain has “forgotten” where it started from

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence

We start MCMC from an arbitrary initial state

The samples will be coming from the chain’s stationary distribution π
only when the chain has “forgotten” where it started from
Mixing time: amount of time it takes to enter the stationary
distribution

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence

We start MCMC from an arbitrary initial state

The samples will be coming from the chain’s stationary distribution π
only when the chain has “forgotten” where it started from
Mixing time: amount of time it takes to enter the stationary
distribution
Samples collected before the chain has reached its stationary
distribution do not come from π, and are usually thrown away

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence

We start MCMC from an arbitrary initial state

The samples will be coming from the chain’s stationary distribution π
only when the chain has “forgotten” where it started from
Mixing time: amount of time it takes to enter the stationary
distribution
Samples collected before the chain has reached its stationary
distribution do not come from π, and are usually thrown away
Burn-in phase: initial period, whose samples will be ignored

(U of Manchester) Statistical Inference 45 / 46

MCMC Convergence Diagnostic

(U of Manchester) Statistical Inference 46 / 46

MCMC Convergence Diagnostic
Assessing if the method has converged:
Run multiple chains (typically 3 or 4) from very different overdispersed
starting points
Trace plot: plot the samples of some quantity of interest
E.g. value of a certain component of the state vector, or some event
such as the value taking on an extreme value

(U of Manchester) Statistical Inference 46 / 46

MCMC Convergence Diagnostic
Assessing if the method has converged:
Run multiple chains (typically 3 or 4) from very different overdispersed
starting points
Trace plot: plot the samples of some quantity of interest
E.g. value of a certain component of the state vector, or some event
such as the value taking on an extreme value
If the chain has mixed, it should have “forgotten” where it started
from
Trace plots should converge to the same distribution, and thus overlap
with each other

(U of Manchester) Statistical Inference 46 / 46

MCMC Convergence Diagnostic
Assessing if the method has converged:
Run multiple chains (typically 3 or 4) from very different overdispersed
starting points
Trace plot: plot the samples of some quantity of interest
E.g. value of a certain component of the state vector, or some event
such as the value taking on an extreme value
If the chain has mixed, it should have “forgotten” where it started
from
Trace plots should converge to the same distribution, and thus overlap
with each other

(U of Manchester) Statistical Inference 46 / 46

MCMC Convergence Diagnostic

MCMC lets us draw samples from a target distribution (assuming it

has converged), but the samples are not independent
We may need to draw a lot of them to get a reliable estimate

How to compute the effective sample size from a set of (possibly

correlated) samples: (Murphy, 2023, § 12.6.3)

How to improve the speed of convergence: (Murphy, 2023, § 12.6.4)

(U of Manchester) Statistical Inference 46 / 46

Probabilistic programming

Plenty of software for Bayesian inference

How do they work?

declarative specification of probabilistic models (likelihood and priors)
inference is performed automatically

Some popular options include PyMC, Stan, TensorFlow Probability,

Pyro, and NumPyro.

Most of these PPs implement an efficient version of HMC known as

No-U-Turn sampler and Automatic Differentiation Variational
Inference (ADVI).

The main difference among them is on the type of automatic

differentiation framework they use, i.e. TensorFlow, PyTorch, Theano,
JAX, etc.
(U of Manchester) Statistical Inference 46 / 46
References

Sebastian Bieringer, Gregor Kasieczka, Maximilian F. Steffen, and Mathias

Trabs. Statistical guarantees for stochastic Metropolis-Hastings.
Available at arXiv:2310.09335, 2023.
Kevin P. Murphy. Probabilistic Machine Learning: Advanced Topics. MIT
Press, 2023. URL [Link]

(U of Manchester) Statistical Inference 1/1