A Conceptual Introduction to Markov Chain Monte Carlo Methods
Joshua S. Speagle



Combinatorial relations on skew Schur and skew stable Grothendieck polynomials
    [
    December 30, 2023
==============================================================================





Markov Chain Monte Carlo (MCMC) methods have become a cornerstone of
many modern scientific analyses by providing a straightforward
approach to numerically estimate uncertainties in the parameters of a model
using a sequence of random samples.
This article provides a basic introduction to MCMC methods
by establishing a strong conceptual understanding of
what problems MCMC methods are trying to solve,
why we want to use them, and
how they work in theory and in practice.
To develop these concepts, I outline the foundations of Bayesian inference,
discuss how posterior distributions are used in practice,
explore basic approaches to estimate posterior-based quantities, and
derive their link to Monte Carlo sampling and MCMC.
Using a simple toy problem, I then demonstrate how these concepts
can be used to understand the benefits and drawbacks of various MCMC approaches.
Exercises designed to highlight various concepts are also included throughout the
article.








§ INTRODUCTION
 

Scientific analyses generally rest on making inferences about
underlying physical models from various sources of observational data.
Over the last few decades, the quality and quantity of these
data have increased substantially as they become faster and cheaper to 
collect and store. At the same time, the same technology that has
made it possible to collect vast amounts of data has also led to
a substantial increase in the computational power and resources
available to analyze them. 

Together, these changes have made it possible to explore
increasingly complex models using methods that
can exploit these computational resources. This has led to
a dramatic rise in the number of published works that rely on
Monte Carlo methods, which use
a combination of numerical simulation and random number generation
to explore these models.

One particularly popular subset of Monte Carlo methods is known as
Markov Chain Monte Carlo (MCMC). MCMC methods are appealing
because they provide a straightforward, intuitive way to both simulate values
from an unknown distribution and use those simulated values
to perform subsequent analyses.
This allows them to be applicable in a wide variety of domains.

Owing to its widespread use, various overviews of MCMC methods are
common both in peer-reviewed and non-peer-reviewed sources.
In general, these tend to
fall into two groups: articles focused on various statistical underpinnings
of MCMC methods and articles focused on implementation and practical usage.
Readers interested in reading more details on either topic are encouraged
to see <cit.> and <cit.> along with
associated references therein.

This article instead provides an overview of MCMC methods focused
instead on building up a strong conceptual understanding of
the what, why, and how of MCMC based on statistical intuition.
In particular, it tries to systematically answer the following questions:

    
  * What problems are MCMC methods trying to solve?
    
  * Why are we interested in using them?
    
  * How do they work in theory and in practice?


When answering these questions, this article generally assumes that
the reader is somewhat familiar with
the basics of Bayesian inference in theory (e.g., the role of priors)
and in practice (e.g., deriving posteriors), 
basic statistics (e.g., expectation values), 
and basic numerical methods (e.g., Riemann sums).
No advanced statistical knowledge is required. For more details
on these topics, please see <cit.> and <cit.>
along with associated references therein.

The outline of the article is as follows. 
In <ref>, I provide a brief review of Bayesian inference
and posterior distributions. 
In <ref>, I discuss what posteriors are used for in practice,
focusing on integration and marginalization.
In <ref>, I outline a basic scheme to approximate
these posterior integrals using discrete grids.
In <ref>, I illustrate how Monte Carlo methods
emerge as a natural extension of grid-based approaches.
In <ref>, I discuss how MCMC methods fit within
the broader scope of possible approaches and their benefits and drawbacks.
In <ref>, I explore the general challenges MCMC methods face.
In <ref>, I examine how these concepts come together
in practice using a simple example.
I conclude in <ref>.



§ BAYESIAN INFERENCE
 

In many scientific applications, we have access to some data 
that we want to use to make inferences about the world around us.
Most often, we want to interpret these data in light of an underlying
model M that can make predictions about the data we expect
to see as a function of some parameters _M of that
particular model.

We can combine these pieces together to estimate the
probability P(|_M, M) 
that we would actually see that data  we have
collected conditioned on (i.e. assuming)
a specific choice of parameters _M from our model M.
In other words, assuming our model M is right and
the parameters _M describe the data, what is the
likelihood P(|_M, M) of the parameters _M
based on the observed data ? Assuming different values
of _M will give different likelihoods, telling us
which parameter choices appear to best describe the data we
observe.

In Bayesian inference, we are interested in
inferring the flipped quantity, P(_M|, M).
This describes the probability that the 
underlying parameters are actually _M 
given our data  and assuming a particular model M.
By using factoring of probability, we can relate this new
probability P(_M|, M) to the likelihood
P(|_M, M) described above as

    P(_M | , M) P( | M) 
        = P(_M,  | M)
        = P( | _M, M) P(_M | M)

where P(_M,  | M) represents the joint
probability of having an underlying set of parameters _M
that describe the data and observing the particular set of data 
we have already collected.

Rearranging this equality into a more convenient form 
gives us Bayes' Theorem:

    P(_M | , M) 
        = P( | _M, M) P(_M | M)/P( | M)

This equation now describes exactly how our two probabilities relate to each other.

P(_M | M) is often referred to as the prior. This describes
the probability of having a particular set of values _M
for our given model M before conditioning on our data.
Because this is independent of the data, this term
is often interpreted as representing 
our “prior beliefs” about what _M should
be based on previous measurements, physical concerns, 
and other known factors. In practice, this has the effect of essentially
“augmenting” the data with other information.

The denominator

    P( | M) = ∫ P( | _M, M) P(_M | M) _M

is known as the evidence or marginal likelihood
for our model M marginalized (i.e. integrated)
over all possible parameter values _M.
This broadly tries to quantify
how well our model M explains the data  after averaging over
all possible values _M of the true underlying parameters.
In other words, if the observations predicted by our model look similar
to the data , then M is a good model.
Models where this is true more often also tend to be favored 
over models that give excellent agreement occasionally but 
disagree most of the time. Since in most instances we take
 as a given, this often ends up being a constant.

Finally, P(_M | , M) represents our posterior.
This quantifies our belief in _M after combining our
prior intuition P(_M|M) with
current observations P(|_M,M) 
and normalizing by the overall evidence P(|M).
The posterior will be some compromise between the prior and 
the likelihood, with the exact combination depending
on the strength and properties of the prior and the quality
of the data used to derive the likelihood.
A schematic illustration is shown in 
 <ref>.



Throughout the rest of the paper I will write 
these four terms (likelihood, prior, evidence, posterior)
using shorthand notation such that

    () 
        ≡()()/∫()() ≡()()/

where () ≡ P(_M | , M) is the posterior,
() ≡ P( | _M, M) is the likelihood,
() ≡ P(_M | M) is the prior, and the constant
≡ P( | M) is the evidence. I have suppressed the
model M and data  notation
for convenience here since in most cases the data and model are considered
fixed, but will re-introduce them as necessary.

Before moving on, I would like to close by emphasizing that
the interpretation of any result is
only as good as the models and priors that underlie them.
Trying to explore the implications of any particular model
using, for instance, some of the methods described in this article
is fundamentally a secondary concern behind constructing 
a reasonable model with well-motivated priors in the first place. 
I strongly encourage readers to keep this idea in mind
throughout the remainder of this work.



 §.§ Exercise: Noisy Mean
 



  §.§.§ Setup


Consider the case where we have temperature monitoring stations
located across a city. Each station i takes a noisy measurement
T̂_i of the temperature on any given day with some
measurement noise σ_i. We will assume our measurements T̂_i
follow a Normal (i.e. Gaussian) distribution
with mean T and standard deviation σ_i such that

    T̂_i ∼Tσ_i

This translates into a probability of

    P(T̂_i|T,σ_i) ≡Tσ_i
        = 1/√(2πσ_i^2)exp[-1/2(T̂_i - T)^2/σ_i^2]

for each observation and

    P({T̂_i }_i=1^n | T, {σ_i }_i=1^n)
        = ∏_i=1^n P(T̂_i|T,σ_i)

for a collection of n observations.

Let's assume we have five independent noisy measurements
of the temperature (in Celsius) from several monitoring stations

    T̂_1 = 26.3,  T̂_2 = 30.2,  T̂_3 = 29.4, 
        T̂_4 = 30.1,  T̂_5 = 29.8

with corresponding uncertainties

    σ_1 = 1.7,  σ_2 = 1.8,  σ_3 = 1.2, 
        σ_4 = 0.5,  σ_5 = 1.3


Looking at historical data, we find that the typical underlying
temperature T during similar days is roughly Normally-distributed
with a mean T_ prior=25 and variation σ_ prior = 1.5:

    T ∼T_ prior=25σ_ prior=1.5




  §.§.§ Problem


Using these assumptions, compute:

    
  * the prior (T),
    
  * the likelihood (T), and
    
  * the posterior (T)

given our observed data {T̂_i } and errors {σ_i }
over a range of temperatures T. How do the three terms differ?
Does the prior look like a good assumption? Why or why not?



§ WHAT ARE POSTERIORS GOOD FOR?
 

Above, I described how Bayes' Theorem is able to combine
our prior beliefs and the observed data into a new posterior estimate 
() ∝() ().
This, however, is only half of the problem. Once we have
the posterior, we need to then use it to
make inferences about the world around us. 
In general, the ways in which we want to use posteriors
fall into a few broad categories:

    
  * Making educated guesses:
    make a reasonable guess 
    at what the underlying model parameters are.
    
  * Quantifying uncertainty:
    provide constraints on 
    the range of possible model parameter values.
    
  * Generating predictions:
    marginalize over uncertainties in the underlying model parameters 
    to predict observables or other variables that depend on the
    model parameters.
    
  * Comparing models:
    use the evidences from different models
    to determine which models are more favorable.


In order to accomplish these goals,
we are often more interested in trying to use the
posterior to estimate various constraints on the parameters
 themselves or other quantities f() that
might be based on them. This often depends on marginalizing
over the uncertainties characterized by our posterior (via the likelihood
and prior). The evidence , for instance, is
again just the integral of the likelihood and the prior 
over all possible parameters:

    = ∫() () ≡∫()

where () ≡() ()
is the unnormalized posterior.

Likewise, if we are investigating the behavior of a subset of
“interesting” parameters _ int from 
= {_ int, _ nuis},
we want to marginalize over the behavior of the remaining 
“nuisance” parameters _ nuis to see how
they can impact _ int. This process is pretty
straightforward if the entire posterior over  is known:

    (_ int)
        = ∫(_ int, _ nuis)  _ nuis
        = ∫() _ nuis


Other quantities can generally be derived from the expectation value
of various parameter-dependent functions f() with respect to the posterior:

    f()≡∫ f() () /∫()  
        = ∫ f() () /∫()  
        = ∫ f() ()

since ∫()  = 1 by definition 
and () ∝().
This represents a weighted average of f(),
where at each value  we weight the resulting
f() based on to the 
chance we believe that value is correct.

Taken together, we see that in almost all cases we are
more interested in computing integrals over the posterior
rather than knowing the posterior itself.
To put this another way, the posterior is rarely ever useful on its
own; it mainly becomes useful by integrating over it.


























This distinction between estimating expectations and other integrals
over the posterior versus estimating the posterior 
in-and-of-itself is a key element of Bayesian inference.
This distinction is hugely important when it comes
to actually performing inference in practice, since it
is often the case that we can get an excellent estimate of
f() even if we have an
extremely poor estimate of () 
or ().

More details are provided below to further illustrate how
the particular categories described above translate into
particular integrals over the (unnormalized) posterior.
An example is shown in  <ref>.





 §.§ Making Educated Guesses
 

One of the core tenets of Bayesian inference is that
we don't know the true model M_* or its true underlying parameters
_* that characterize the data we observe: the model M
we have is almost always a simplification of what is actually going on.
If we assume that our current model M is correct, however,
we can try to use our posterior ()
to propose a point estimate 
that we think is a pretty good guess for the true value _*.

What exactly counts as “good”? This depends on exactly what we care
about. In general, we can quantify “goodness” by asking the opposite question:
how badly are we penalized if our estimate ≠_*
is wrong? This is often encapsulated through the use of
a loss function L(|_*)
that penalizes us when our point estimate 
differs from _*. An example of a common loss function
is L(|_*) = | - _*|^2
(i.e. squared loss), where an incorrect guess is penalized based on
the square of the magnitude of the separation between the
guess  and the true value _*.

Unfortunately, we don't know what the actual value of
_* is to evaluate the true loss. We can, however,
do the next best thing and compute the expected loss
averaged over all possible values of _* based on
our posterior:

    L_()
        ≡L(|)
        = ∫ L(|) ()

A reasonable choice for  is then
the value that minimizes this expected loss in place of the
actual (unknown) loss:

    ≡_'[L_(')]

where  indicates the value (argument) of '
that minimizes the expected loss L_(').

While this strategy can work for any arbitrary loss function, solving
for  often requires using numerical methods and repeated
integration over (). However, analytic solutions
do exist for particular loss functions. For example,
it is straightforward to show 
(and an insightful exercise for the interested reader) 
that the optimal point estimate 
under squared loss is simply the mean.



 §.§ Quantifying Uncertainty
 

In many cases we are not just interested in computing
a prediction  for _*, but also constraining
a region () of possible values 
within which _* might lie with some amount of certainty.
In other words, can we construct a region _X
such that we believe there is an X% chance that it contains
_*?

There are many possible definitions
for this credible region. One common definition
is the region above some posterior threshold _X
where X% of the posterior is contained, i.e. where

    ∫_ ∈ _X() 
        = X/100

given

    _X
        ≡{ : () ≥_X }


In other words, we want to integrate our posterior over all 
where the value () > _X is greater than some threshold
_X, where _X is set so 
that this integral encompasses X% of the full posterior.
Common choices for X include 68% and 95% (i.e.
“1-sigma” and “2-sigma” credible intervals).

In the special case where our (marginalized) posterior is 1-D, credible intervals 
are often defined using percentiles rather than thresholds, where
the location x_p of the pth percentile is defined as

    ∫_-∞^x_p(x)  x = p/100

We can use these to define a credible region [x_ low, x_ high]
containing Y% of the data by taking x_ low = x_(1-Y)/2 and
x_ high = x_(1+Y)/2. While this leads to asymmetric thresholds
and does not generalize to higher dimensions, it has the benefit of
always encompassing the median value x_50 and having equal tail
probabilities (i.e. (1-Y)/2% of the posterior on each side).

In general, when referring to “credible intervals” throughout the text
the percentile definition should be assumed unless explicitly stated otherwise.



 §.§ Making Predictions
 

In addition to trying to estimate the underlying parameters of our
model, we often also want to make predictions of other observables or variables
that depend on our model parameters.
If we think we know the underlying true model parameters _*, then
this process is straightforward. Given that we only have access to the
posterior distribution () over possible
values _* could take, however, to predict what will happen we
will need to marginalize over this uncertainty.

We can quantify this intuition using the posterior predictive
P(|), which represents the probability
of seeing some new data  based on our existing data :

    P(|) 
        ≡∫ P(|) P(|) ≡∫() () 
        = ()

In other words, for hypothetical data , 
we want to compute the expected value of the likelihood
() over all possible
values of  based on the current posterior ().



 §.§ Comparing Models
 

One final point of interest in many Bayesian analyses is
trying to investigate whether the data particularly favors
any of the model(s) we are assuming in our analysis.
Our choice of priors or the particular way we parameterize the
data can lead to substantial differences in the way we
might want to interpret our results.

We can compare two models by computing the Bayes factor:

    ^1_2≡P(M_ 1|)/P(M_ 2|)
        = P(|M_ 1)P(M_ 1)/P(|M_ 2)P(M_ 2)≡_ 1/_ 2_ 1/_ 2

where _M is again the evidence for model M 
and _M is our prior belief that M 
is correct relative to the competing model.
Taken together, the Bayes factor  tells us
how much a particular model is favored over another given the
observed data, marginalizing over all possible values of the underlying
model parameters _M, and our previous relative confidence in the model.

Again, note that computing _M requires computing the
integral ∫() of the
unnormalized posterior () over . 
Combined with the other examples outlined in this section,
it is clear that many common use cases in Bayesian analysis rely on
computing integrals over the (possibly unnormalized) posterior.



 §.§ Exercise: Noisy Mean Revisited
 



  §.§.§ Setup


Let's return to our temperature posterior (T) from
<ref>. We want to use this result to derive interesting
estimates and constraints on the possible underlying temperature T.



  §.§.§ Point Estimates


The mean can be defined as the
point estimate  that minimizes the expected loss
L_() under squared loss:

    L_ mean(|_*) = | - _*|^2

The median can be defined as the point estimate that minimizes
L_() under absolute loss:

    L_ med(|_*) = | - _*|

And the mode can be defined as the point estimate that
minimizes L_() under “catastrophic” loss:

    L_ mode(|_*) = -δ(|-_*|)

where δ(·) is the Dirac delta function defined such that

    ∫ f(x) δ(x-a)  x = f(a)


Given these expressions for the mean, median,
and mode, estimate the corresponding temperature point estimate
T_ mean, T_ med, and T_ mode from our corresponding
posterior. Feel free to experiment with various analytic and numerical
methods to perform these calculations.

We might expect that the historical data we used for our priors
might not hold as well today if there have been some long-term changes
in the average temperature.
For instance, we expect that the average temperature has increased
over time, and so we might not want to penalize hotter temperatures
T ≥ T_ prior as much as cooler ones T < T_ prior.
We can encode this information in an asymmetric loss function such as

    L(T̂|T_*) = 
        
        |T̂ - T_*|^3     T < T_ prior
    
        |T̂ - T_*|     T ≥ T_ prior

What is the optimal point estimate T_ asym that minimizes
the expected loss in this case?



  §.§.§ Credible Intervals


Next, let's try to quantify the uncertainty. Given the posterior
(T), compute the 50%, 80%, and 95% credible intervals
using posterior thresholds _X. Next, compute these
credible intervals using percentiles. Are there are differences between
the credible intervals computed from the two methods? Why or why not?



  §.§.§ Posterior Predictive


To propagate our uncertainties into the next observations,
compute the posterior predictive P(T̂_6|{T̂_1, …, T̂_5 }) 
over a range of possible temperature measurements T̂_6 
for the next observations given the previous five
{T̂_1, …, T̂_5 } assuming
an uncertainty of σ_6=0, σ_6=0.5, and σ_6=2.



  §.§.§ Model Comparison


Finally, we want to investigate whether our prior appears to be a good
assumption. Using numerical methods, compute the
evidence  for our default prior with mean T_ prior = 25
and standard deviation σ_ prior = 1.5. Then compare this to the
evidence estimated based on an alternative prior
where we assume the temperature has risen by roughly five
degrees with mean T_ prior = 30 but with a corresponding larger uncertainty
σ_ prior = 3. Is one model particularly favored over the other?



§ APPROXIMATING POSTERIOR INTEGRALS WITH GRIDS
 

I now want to investigate methods for estimating posterior integrals.
While in some cases (e.g., conjugate priors) 
these can be computed analytically, this is not true in general.
To properly estimate quantities such as those outlined in
<ref> therefore requires the use of numerical methods
(highlighted in the previous exercises).

To start, I will first focus on the case where our integral over
 is 1-D. In that case, we can approximate it using standard
numerical techniques such as a Riemann sum over a discrete grid
of points:

    f() 
        = ∫ f() () ≈∑_i=1^n f(_i) 
        (_i) Δ_i

where

    Δ_i = _j+1 - _j

is simply the spacing between the set of j=1,…,n+1 points
on the underlying grid and

    _i = _j+1 + _j/2

is just defined to be the mid-point between _j and
_j+1.[Choosing _i to be one of the end-points 
gives consistent behavior (see <ref>) as 
the number of grid points n →∞ but generally leads
to larger biases for finite n.] As shown in 
 <ref>, 
this approach is akin to trying to approximate the
integral using a discrete set of n rectangles with heights of
f(_i) (_i) and widths of Δ_i.



This idea can be generalized to higher dimensions. In that case,
instead of breaking up the integral into n 1-D segments, 
we instead can decompose it into a set of n N-D cuboids.
The contribution of each of these pieces is then proportional
to the product of the “height” f(_i) (_i)
and the volume

    Δ_i = ∏_j=1^dΔΘ_i,j

where ΔΘ_i,j is the width of the ith cuboid in the jth dimension.
See  <ref>
for a visual representation of this procedure.

Substituting () = ()/
into the expectation value and replacing any integrals 
with their grid-based approximations then gives:

    f()
        = ∫ f() () /∫() 
        = ∫ f() () /∫() ≈∑_i=1^n f(_i) 
        (_i) Δ_i/∑_i=1^n(_i) Δ_i

Note the denominator is now an estimate for the evidence:

    = ∫() ≈∑_i=1^n(_i) Δ_j

This substitution of the unnormalized posterior
()
for the posterior () is a crucial part of computing 
expectation values in practice since we can compute 
() = () ()
directly without knowing .



 §.§ The Curse of Dimensionality
 

While this approach is straightforward, it has one immediate
and severe drawback: the total number of grid points increases
exponentially as the number of dimensions increases. 
For example, assuming we have roughly k ≥ 2
grid points in each dimensions, the total number of points n in our
grid scales as

    n ∼∏_j=1^d k = k^d

This means that even in the absolute best case where k=2,
we have 2^d scaling.

This awful scaling is often referred to as the curse of dimensionality.
This exponential dependence turns out to be a generic feature 
of high-dimensional distributions
(i.e. posteriors of models with larger numbers of parameters)
that I will return to later in <ref>.



 §.§ Effective Sample Size
 

Apart from this exponential scaling of dimensionality, there is
a more subtle drawback to using grids. Since we do not
know the shape of the distribution ahead of time, the contribution of
each portion of the grid (i.e. each N-D cuboid)
can be highly uneven depending on the structure of the grid. 
In other words, the effectiveness of this approach not
only depends on the number of grid points n but 
also where they are allocated. 
If we do not specify our grid points well,
we can end up with many points located in regions where
() and/or f()()
is relatively small. This then implies that their respective sums 
will be dominated by a small number of points with 
much larger relative “weights”. Ideally, we would
want to increase the resolution of the grid
in regions where the posterior is large and decrease it elsewhere
to mitigate this effect. 

Note that our use of the term “weights” in the preceding paragraph 
is quite deliberate. Looking back at our original approximation,
the form of equation (<ref>) 
is quite similar to one which might be used to compute a 
weighted sample mean of f(). 
In that case, where we have n observations { f_1, …, f_n }
with corresponding weights { w_1, …, w_n }, the
weighted mean is simply:

    f̂_ mean≡∑_i=1^n w_i f_i/∑_i=1^n w_i

Indeed, if we define

    f_i ≡ f(_i),
           
        w_i ≡(_i) Δ_i

then the connection between the weighted sample mean
in equation (<ref>) and the expectation value from our grid
in equation (<ref>) becomes explicit:

    f()≈∑_i=1^n f(_i) 
        (_i) Δ_i/∑_i=1^n(_i) Δ_i≡∑_i=1^n w_i f_i/∑_i=1^n w_i


Thinking about our grid as a set of n samples
also allows us to consider an associated effective sample size (ESS)
n_ eff≤ n. The ESS encapsulates the idea that not all of our
samples contribute the same amount of information:
if we have n samples that are very similar to each other, we
expect to have a substantially worse estimate
than if we have n samples that are quite different.
This is because the information in correlated samples are 
at least partially redundant with one another, with the amount of
redundancy increasing with the strength of the correlation: while
two independent samples provide completely unique information about
the distribution and no information about each other, 
two correlated samples instead provide some information about
each other at the expense of the underlying distribution.



Returning to grids, this correspondence means that we can in theory
come up with an estimate of the expectation value
f() that is at least as good
as the one we might currently have using a smaller number
n_ eff≤ n of grid points if we were able to allocate them
more efficiently. This distinction matters because errors on 
our estimate of the expectation value generally scale as a function
of n_ eff rather than n. For instance, the error on the
mean typically goes as ∝ n_ eff^-1/2 rather than
∝ n^-1/2.

We can quantify the ideas behind the ESS as discussed above
by introducing a formal definition following <cit.>:

    n_ eff≡(∑_i=1^n w_i)^2/∑_i=1^n w_i^2

In line with our intuition, the best case under this definition is
one where all the weights are equal (w_i = w): 

    n_ eff^ best
        = (∑_i=1^n w_i)^2/∑_i=1^n w_i^2
        = (nw)^2/∑_i=1^n w^2
        = n^2w^2/nw^2 = n

Likewise, the worst case is one where
all the weight is concentrated around a single sample 
(w_i = w for i=j and w_i=0 otherwise):

    n_ eff^ worst
        = (∑_i=1^n w_i)^2/∑_i=1^n w_i^2
        = (w)^2/w^2
        = 1

This former situation (with n_ eff^ best)
would be the case where each of the elements of our grid
all have roughly the same contribution to the integral, while the
latter (with n_ eff^ worst) would be where the entire
integral is essentially contained in just one of our n N-D cuboid regions.
An illustration of this behavior is shown in
 <ref>.



 §.§ Convergence and Consistency
 

Now that I have outlined the relationship between the structure
of our grid and the ESS, I want to examine
two final issues: convergence and consistency.
Convergence is the idea that, while our estimates using n
samples (grid points) might be noisy, it approaches some fiducial value
as n →∞:

    lim_n →∞∑_i=1^n f(_i) 
        (_i) Δ_i/∑_i=1^n(_i) Δ_i
        = C

Consistency is subsequently the idea that the value we converge to
is the true value we are interested in estimating:

    lim_n →∞∑_i=1^n f(_i) 
        (_i) Δ_i/∑_i=1^n(_i) Δ_i
        = f()




It is straightforward to show
that if the expectation value is well-defined (i.e. it exists)
and the grid covers the entire domain of  (i.e. spans
the smallest and largest possible values in every dimension)
then using a grid is a consistent way to estimate 
the expectation value.
This should make intuitive sense: provided our grid is
expansive enough in  so that we're not “missing”
any region of parameter space, we should be able to estimate
f() to arbitrary precision
by simply increasing the resolution in Δ.

Unfortunately, we do not know beforehand what range of 
values of  our grid should span. While parameters
can range over (-∞, +∞), grids rely on finite-volume elements
and so we have to choose some finite sub-space to grid up.
So while grids may give estimates that converge to some
value over the range spanned by the grid points,
there is always a possibility that a significant portion of the posterior
lies outside that range. In these cases,
grids are not guaranteed to be consistent estimators
of f(). An illustration of this issue is
shown in  <ref>.
This fundamental problem is not shared by Monte Carlo
methods, which I will cover in <ref>.



 §.§ Exercise: Grids over a 2-D Gaussian
 



  §.§.§ Setup


Consider an unnormalized posterior 
well-approximated by a 2-D Gaussian (Normal)
distribution centered on (μ_x,μ_y) 
with standard deviations (σ_x, σ_y):

    (x,y) 
        = exp{-1/2[(x-μ_x)^2/σ_x^2
        + (y-μ_y)^2/σ_y^2]}

Assume that we expect to find our posterior has a mean of 0 and
a standard deviation of 1.
In reality, however, our posterior actually has means (μ_x,μ_y)=(-0.3,0.8) and
standard deviations (σ_x^2,σ_y^2)=(2,0.5), mimicking the common
case where our prior expectations and posterior inferences somewhat disagree.



  §.§.§ Grid-based Estimation


We want to use a 2-D grid to estimate various forms of posterior integrals.
Starting with an evenly-spaced 5 × 5 grid from [-2, 2], compute:

    
  * the evidence ,
    
  * the means x
    and y,
    
  * the 68% credible intervals (or closest approximation) 
    [x_ low, x_ high] and [y_ low, y_ high],
    
  * and the effective sample size n_ eff.

How accurate are each of these quantities with the values we might
expect? What does n_ eff/n tell us about how efficiently
we have allocated our grid points?



  §.§.§ Convergence


Repeat the above exercise using an evenly-spaced grid of
20 × 20 points and 100 × 100 points. Comment
on any differences. How much has the overall accuracy improved? Do
the estimates appear convergent?



  §.§.§ Consistency


Next, expand the bounds of the grid to be from [-5, 5] and
perform the same exercise as above. Do the answers change substantially?
If so, what does this tell us about the consistency of our previous
estimates? Adjust the density and bounds of the grid until the 
answers appear both convergent and consistent.
Remember that we do not know the exact shape of the posterior 
ahead of time. What does this imply
about general concerns when applying grids in practice?



  §.§.§ Effective Sample Size


Finally, explore whether there is a straightforward scheme to adjust the
locations of the x and y grid points to maximize the effective
sample size based on the definition outlined in <ref>.
If so, can you explain why it works? If not, why not?
Compared to equivalent evenly-spaced grids, how much can adaptively adjusting
the grid spacing improve n_ eff and the overall accuracy of our estimates?



§ FROM GRIDS TO MONTE CARLO METHODS
 



 §.§ Connecting Grid Points and Samples



Earlier, I outlined how we can relate estimating
f() using a grid of n points to
an equivalent estimate using a set of n samples { f_1, …, f_n}
and a series of associated weights { w_1, …, w_n }.
The main result is that there is an intimate connection between the
structure of the posterior and the grid to the relative amplitude
of the weights w_i ≡(_i)Δ_i
for each point f_i ≡ f(_i). Adjusting the resolution
of the grid then affects these weights, with a more uniform distribution
of weights leading to a larger ESS which can improve our estimate.

The fact that decreasing the spacing (making grid denser)
also decreases the weights makes sense: we have more points 
located in that region, so each point should in general
get less relative weight when computing f(). 
Likewise, if we have the same spacing
but change the relative shape of the posterior, 
the weight of that point when estimating f()
should also change accordingly.

I now want to extend this basic relationship further.
In theory, adaptively increasing the resolution of our grid allows
us more control over the volume elements Δ_i used to
derive our weights. If we knew the shape of our posterior sufficiently
well, for large n 
we should in theory be able to adjust Δ_i such that
the weights w_i = (_i)Δ_i
are uniform to some amount of desired precision.
By inspection, this should happen when

    Δ_i ∝1/(_i)

for all i.

Taking this reasoning to its conceptual limit, as n →∞
we can imagine estimating the posterior using a larger and larger number of 
grid points whose spacing Δ changes as a function of .
Using this, we can now define the density of points ()
based on the varying resolution Δ() of our
infinitely-fine grid as a function of :

    () ∝1/Δ ()

This result suggests that, in the continuum limit where n →∞, 
the structure of our infinite-resolution grid is equivalent to a new
continuous distribution (). An illustration
of this concept is shown in  <ref>. 
Using (), we can 
then rewrite our original expectation value as

    f()≡∫ f() () /∫() 
        = ∫ f() ()/()() /∫()/()() 
        = f() 
        ()/()/()/()

For reasons that will soon become clear, I will refer to 
() as the proposal distribution.



At this point, this may mostly seem like a mathematical trick:
all I have done is rewrite our original single 
expectation value with respect to the (unnormalized)
posterior () in terms of two
expectation values with respect to the proposal distribution
(). 
This substitution, however, actually allows us to fully realize 
the connection between grid points and samples.

Earlier, I showed that the estimate for the expectation
value from grid points is exactly analogous to the estimate
we would derive assuming the grid points were random samples
{ f_1, …, f_n } with associated weights
{ w_1, …, w_n }. Once we have defined our
expectation with respect to (), however,
this statement can become exact assuming we can explicitly generate
samples from ().

Let's quickly review what this means. Initially,
we looked at trying to estimate f() over a grid
with n points. In the limit of infinite resolution, however, our grid becomes
equivalent to some distribution (). Using (),
we can then rewrite our original expression in terms of two expectations,
f() ()/()
and ()/(), over
() instead of ().
This helps us because we can in theory estimate these final expressions
explicitly using a series of n randomly generated samples from
(). Due to the randomness inherent in this approach, 
this is commonly referred to as a Monte Carlo approach for estimating
f() due to historical connections with
randomness and gambling.

On the face of it, this should come across as a surprising claim. When we
compute an integral of a function f() on a bounded grid, we
know that there is some error in our approximation having to do with the
discretization of the grid. This error is entirely deterministic:
given a number of grid points n and an a particular discretization density
() ∝ 1/Δ(), we will get the
same result (and error) for f() every time.

By contrast, drawing n samples {_1,…,_n}
from () is an inherently random (i.e. stochastic) 
process that seems to look nothing like a
grid of points. And because these points are inherently random,
the actual deviation between our estimate and the true value of
f() will also be random.
The “error” from random samples then tells us something about
how much we expect our estimate can differ over many possible
realizations of our random process given a particular number of samples
n generated from (). The fact that we can derive
roughly equivalent estimates from these these very different
approaches as we adjust n and () lies at the heart of the
connection between grid points and samples.

There are three primary benefits from moving from an adaptively-spaced grid
to a continuous distribution ().
First, a grid will always have some minimum 
resolution Δ_i that makes it difficult
to get our weights to be roughly uniform, limiting our maximum
ESS in practice. By contrast, we can in theory get
() to more closely match the posterior (),
giving a larger ESS at fixed n.

Second, because we are now working with distributions
rather than a finite number of grid points, we are no longer limited
to some finite volume when estimating expectations. Since distributions can
range over (-∞, +∞), we can guarantee ()
will provide sufficient coverage over all possible  
values that our posterior () could be defined over.
This means that some of the theoretical issues raised in <ref>
associated with applying grids to posteriors 
that range over (-∞, +∞) no longer apply. 
Monte Carlo methods therefore can serve as a consistent
estimator for a wider range of possible posterior expectations
than grid-based methods, making them substantially more flexible.

Finally, the minimum number of grid points always scales exponentially
with dimensionality (see <ref>), regardless of how many 
parameters we are interested in marginalizing over. Since 
Monte Carlo methods do not rely on these, they can take full
advantage of marginalizing over parameters 
when estimating expectations f(). 
They are therefore less susceptible to this effect 
(although see <ref>).



 §.§ Importance Sampling
 




As I have tried to emphasize previously, the core tenet of this article is that
we do not know what () looks like beforehand.
This means we do not know what grid structure will
provide an optimal estimate (i.e. maximum ESS) for
f(, let alone how this should
behave as () in the continuum limit. This gives
us ample motivation to choose ()
in such a way to make generating samples from it easy and straightforward.

Assuming we have chosen such a (), we can 
subsequently generate a series of n samples from it. Assuming
these samples have weights q_i associated with them and defining

    f(_i) ≡ f_i,   (_i)/(_i) 
        ≡w̃(_i) ≡w̃_i

our original expression reduces to

    f() 
        = f() 
        w̃()/w̃()≈∑_i=1^n f_i w̃_i q_i/∑_i=1^nw̃_i q_i

If we further assume that we have chosen () so that
we can simulate samples that are 
independently and identically distributed (iid)
(i.e. each sample has the same probability distribution 
as the others and all the samples are mutually independent),
then the corresponding sample weights immediately reduce to 
q_i = 1/n and our result becomes

    f()≈n^-1∑_i=1^n f_i w̃_i/n^-1∑_i=1^nw̃_i

As with the previous case using grids (<ref>),
the denominator of this
expression is again a direct approximation for the evidence

    = ∫() ≈ n^-1∑_i=1^nw̃_i


This gives a straightforward recipe for estimating our original expectation
value:

    
  * Draw n iid samples {_1, …, _n }
    from ().
    
  * Compute their corresponding weights w̃_i =
    (_i)/(_i).
    
  * Estimate f()
    by computing w̃() and
    f()w̃() using the
    weighted sample means.

Since this process just involves “reweighting” the samples
based on w̃_i, these weights are often referred to as
importance weights and the method as Importance Sampling.
A schematic illustration of Importance Sampling
is highlighted in  <ref>.



We can interpret the importance weights as ways to correct for how
“far off” our original guess () 
is from the truth ().
If the posterior density is higher at position
_i relative to the proposal density, 
then we were less likely to generate
a sample at that position compared to what we would have
seen if we had drawn samples directly from the posterior.
As a result, we should increase its corresponding weight
to account for this expected deficit of samples at a given position.
If the posterior density is lower relative to the proposal density, 
then the alternative is true and we want to lower the weight of the
corresponding sample to account for the expected excess of samples
at a given position.



 §.§ Examples of Sampling Strategies
 

Importance Sampling serves as a useful first step for understanding how
the weights {w̃_1, …, w̃_n } for the corresponding set 
of n samples are related to different Monte Carlo sampling strategies.

As an example, one common approach is to generate samples uniformly within some 
cuboid with volume V. The proposal distribution for this
will then be

    ^ unif() =
        
        1/V      in cuboid
    
        0     otherwise

The corresponding importance weights subsequently will
just be proportional to the posterior at a given position:

    w̃_i^ unif
        = (_i)/^ unif(_i)
        = V (_i)
        ∝(_i)


Another possible approach would be if we instead take our proposal
to be our prior:

    ^ prior() = ()

This seems like a well-motivated choice: the prior characterizes
our knowledge before looking at the data, so it should serve
as a useful first guess and encompass the range of all possibilities.
Under this assumption, we now find our weights will be 
equal to the likelihood () at each position:

    w_i^ prior 
        = (_i)/^ prior(_i)
        = (_i) (_i)/(_i)
        = (_i)


Finally, notice that the optimal sampling strategy is to
assume that we can take our proposal to be identical to our
posterior:

    ^ post() = ()

The corresponding weights will then just be constant and
equal to the evidence :

    w_i^ post 
        = (_i)/^ post(_i)
        = (_i)/(_i)
        =


As expected, this final result guarantees the maximum
possible ESS of n_ eff = n. Getting
() to be as “close” as possible
to () therefore becomes a crucial part of
analyses when trying to use Importance Sampling 
to estimate expectation values. It is this result in particular 
that motivates the use of Markov Chain Monte Carlo (MCMC) methods
discussed from <ref> onward:
if we can somehow generate samples directly
from () or something close to it, 
then we can achieve an optimal estimate of our
corresponding expectation values.



 §.§ Exercise: Importance Sampling over a 2-D Gaussian
 




  §.§.§ Setup


Let's return to our exercise from <ref>, in which
our unnormalized posterior is well-approximated by a 2-D Gaussian (Normal)
distribution:

    (x,y) 
        = exp{-1/2[(x-μ_x)^2/σ_x^2
        + (y-μ_y)^2/σ_y^2]}

where (μ_x,μ_y)=(-0.3,0.8) and (σ_x^2,σ_y^2)=(2,0.5).



  §.§.§ Importance Sampling


We want to use Importance Sampling to approximate various posterior
integrals from this distribution.
We will start by choosing our proposal distribution (x,y)
to be a 2-D Gaussian with a mean of 0 and standard deviation of 1:

    (x,y) = (μ_x,μ_y)=(0,0)(σ_x,σ_y)=(1,1)


Using n=25 iid random samples drawn 
from the proposal distribution, compute an estimate for:

    
  * the evidence ,
    
  * the means x
    and y,
    
  * the 68% credible intervals (or closest approximation) 
    [x_ low, x_ high] and [y_ low, y_ high],
    
  * and the effective sample size n_ eff.

How accurate are each of these quantities with the values we might
expect? What does n_ eff/n tell us about how well
our proposal (x,y) traces the underlying posterior
(x,y)?



  §.§.§ Uncertainty


Repeat the above exercise m=100 times to get an
estimate for how much our estimates of each quantity can vary.
Is the variation in line with what might be expected given 
the typical effective sample size? Why or why not?



  §.§.§ Convergence


Now repeat the above exercise using n=100, n=1000, and n=10000 
points rather than n=25 points and comment on any differences.
How much has the overall accuracy improved? Do
the estimates appear convergent and consistent as n_ eff increases?
How much do the errors on quantities shrink as a function
of n and/or n_ eff? Is this behavior expected? Why or why not?



  §.§.§ Consistency


Next, let's expand our proposal distribution to instead have
(σ_x,σ_y)=(2,2) to get more coverage in the “tails” of the
posterior. Perform the same exercise as above 
with n={100,1000,10000} iid random samples. 
Do the answers change substantially? Why or why not?

While in theory we can choose (x,y) ≈(x,y) 
so that n_ eff≈ n, we do not know the exact shape of the posterior 
ahead of time. Given that (x,y) may differ from
our initial expectations, what does this exercise imply about general concerns
applying Importance Sampling in practice?



§ MARKOV CHAIN MONTE CARLO
 

Now that we see how the weights relate to various Monte Carlo sampling strategies
(e.g., generating samples from the prior), I will now outline
the idea behind Markov Chain Monte Carlo (MCMC). In brief,
MCMC methods try to generate samples in such a way that the importance
weights {w̃_1, …, w̃_n } associated with each sample
are constant. Based on the results from <ref>, this
means MCMC seeks to generate samples proportional to
the posterior () in order to arrive 
at an optimal estimate for our expectation value.

MCMC accomplishes this by creating a chain of 
(correlated) parameter values 
{_1 →…→_n }
over n iterations such that the number of iterations m(_i)
spent in any particular region δ__i 
centered on _i is proportional to the posterior 
density (_i) contained within that region.
In other words, the “density” of samples generated from MCMC

    ρ() ≡m()/n

at position  integrated over δ_ is approximately

    ∫_∈δ_() ≈∫_∈δ_ρ() ≈ n^-1∑_j=1^n_j ∈δ_

where · is the indicator function which evaluates to
1 if the inside condition is true and 0 otherwise. We can therefore
approximate the density by simply
adding up the number of samples within δ_
and normalizing by the total number of samples n.
A schematic illustration of this concept is shown in
 <ref>.



While this will just be approximately true for any finite n,
as the number of samples n →∞ this procedure generally
guarantees that ρ() →() 
everywhere.[Discussing the details of exactly when/where this
condition holds in theory and in practice
is beyond the scope of this paper but can be found in
other references such as <cit.> and <cit.>.]
In theory then, once we have a reasonable enough
approximation for ρ(), we can also use the 
samples {_1 →…→_n }
generated from ρ()
to get an estimate for the evidence using the 
same substitution trick introduced in <ref>:

    = ∫()/ρ()ρ() ≡()/ρ()ρ≈ n^-1∑_i=1^n(_i)/ρ(_i)

This is just the average of the ratio 
between (_i) and ρ(_i)
over all n samples.

Finally, since our MCMC procedure gives us a series of n
samples from the posterior, our expectation value simply reduces to

    f()≈n^-1∑_i=1^n f_i w̃_i/n^-1∑_i=1^nw̃_i
        = n^-1∑_i=1^n f_i/n^-1∑_i=1^n 1
        = n^-1∑_i=1^n f_i

This is just the sample mean of
the corresponding { f_1, …, f_n } values over
our set of n samples. 

I wish to take a moment here to highlight two features of the above results
related to common misconceptions surrounding MCMC methods. First, there
is a widespread belief that because MCMC methods generate a chain of 
samples whose behavior follows the posterior, we do not
have any ability to use them to estimate normalizing constants such
as the evidence . As shown above, this is not true at all: not only
can we do this using ρ(), but the estimate we derive
is actually a consistent one (although it will converge slowly;
see <ref>).

The second misconception is that the primary goal of MCMC is to
“approximate” or “explore” the posterior. In other words,
to estimate ρ(). However, as shown above,
the ability of MCMC methods to estimate ρ()
is really only useful for estimating the evidence . In fact, by tracing
its heritage from Importance Sampling-based methods, we see its primary purpose
is actually to estimate expectation values (i.e.
integrals over the posterior). I have explicitly tried to avoid
introducing any mention of “approximating the posterior” up to this point
in order to avoid this misconception, but will spend some time discussing
this point in more detail in <ref>.

To summarize, the idea behind MCMC is to simulate a series of
values {_1 →…→_n }
in a way that their density ρ() after a given amount
of time follows the underlying posterior (). We can
then estimate the posterior within any particular region δ_
by simply counting up how many samples we simulate there and normalizing by
the total number of samples n we generated. Because we are also
simulating values directly from the posterior, any expectation values
also reduce to simple sample averages. This procedure is incredibly
intuitive and part of the reason MCMC methods have 
become so widely adopted.



 §.§ Generating Samples with the Metropolis-Hastings Algorithm
 

There is a vast literature on various approaches to generating samples
(see, e.g., cites). Since this article focuses on
building up a conceptual understanding of MCMC methods, 
exploring how the majority of these methods behave both in theory and
in practice is beyond the scope of this paper.

Instead of an overview, I aim to clarify the basics of how these methods operate.
The central idea is that we want a way to generate new samples 
_i →_i+1 such that the
distribution of the final samples ρ() as n →∞
(1) is stationary (i.e. it converges to something) and 
(2) is equal to the (). These are essentially analogs
to the convergence and consistency constraints 
discussed in <ref>.

We can satisfy the first condition by invoking detailed balance.
This is the idea that probability is conserved when moving
from one position to another (i.e. the process is reversible). More formally,
this just reduces to factoring of probability:

    P(_i+1|_i) P(_i) 
        = P(_i+1, _i)
        = P(_i|_i+1) P(_i+1)

where P(_i+1|_i) is the probability of moving
from _i →_i+1 and P(_i|_i+1) 
is the probability of the reverse move from 
_i+1→_i. 
Rearranging then gives the following constraint:

    P(_i+1|_i)/P(_i|_i+1) 
        = P(_i+1)/P(_i)
        = (_i+1)/(_i)

where the final equality comes from the fact that the distribution
we are trying to generate samples from is the posterior ().

We now need to implement a procedure that enables us to actually
move to new positions by computing this probability. 
We can do this by breaking each move into two steps. First,
we want to propose a new position 
_i →_i+1' based on a
proposal distribution (_i+1'|_i)
similar in nature to the () used in to 
Importance Sampling (<ref>). 
Then we will either decide to accept the new position
(_i+1=_i+1') or reject the new position
(_i+1=_i) with some transition probability
T(_i+1'|_i).
Combining these terms together then gives us the probability
of moving to a new position:

    P(_i+1|_i) 
        ≡(_i+1|_i) T(_i+1|i)


As with Importance Sampling, we can choose (_i+1'|_i)
so that it is straightforward to propose new samples
_i+1' by numerical simulation. 
We then need to determine the transition probability
T(_i+1'|_i) of whether we should accept 
or reject _i+1'. Substituting into
our expression for detailed balance, we find that our form for the
transition probability must satisfy the following constraint:

    T(_i+1|_i)/T(_i|_i+1) 
        = (_i+1)/(_i)(_i|_i+1)/(_i+1|_i)

It is straightforward to show that the 
Metropolis criterion <cit.>

    T(_i+1|_i)
        ≡min[1, (_i+1)/(_i)(_i|_i+1)/(_i+1|_i)]

satisfies this constraint.

Generating samples following this approach can be done using the
Metropolis-Hastings (MH) Algorithm <cit.>:

    
  * Propose a new position _i →_i+1'
    by generating a sample from the proposal distribution
    (_i+1'|_i).
    
  * Compute the transition probability
    T(_i+1'|_i)
    = min[1, (_i+1')/(_i)(_i|_i+1')/(_i+1'|_i)].
    
  * Generate a random number u_i+1 from [0, 1].
    
  * If u_i+1≤ T(_i+1'|_i), accept
    the move and set _i+1 = _i+1'. 
    If u_i+1 > T(_i+1'|_i), reject the move
    and set _i+1 = _i.
    
  * Increment i = i + 1 and repeat this process.

See  <ref> for a schematic illustration of this process.



Because algorithms like the MH algorithm generate a chain of states where the
next proposed position only depends on the current position rather than
any of its past positions (i.e. it “forgets” the past),
they are known as Markov processes.
Combining these two terms with the Monte Carlo nature of simulating new
positions is what gives Markov Chain Monte Carlo (MCMC) its namesake.

An issue with generating a chain of samples in practice is the
fact that our chain only has finite length and a starting position _0. 
If our chain were infinitely long, we would expect it
to visit every possible position in parameter space, rendering
the exact starting position is unimportant. However, since in practice we terminate 
sampling after only n iterations, starting from a location _0 that has
an extremely low probability means an inordinate fraction
of our n samples will occupy this low-probability region, possibly biasing
our final results. Since we have limited knowledge beforehand about where
_0 is relative to our posterior, 
in practice we generally want to remove the initial
chain of states once we are confident our chain has begun sampling from
higher-probability regions. Discussing various approaches
for identifying and removing samples from this burn-in period 
is beyond the scope of this article; for additional information,
please see <cit.>, <cit.>, 
and <cit.> along with references therein.



 §.§ Effective Sample Size and Auto-Correlation Time
 

At this point, MCMC seems like it should be the optimal method for
any situation: by simulating samples directly from the (unknown) posterior,
we can achieve an optimal estimate for any expectation values we wish
to evaluate. In practice, however, this does not hold true.
MCMC values rely on specific algorithmic procedures such as 
the MH algorithm to generate samples,
whose limiting behavior reduces to a chain of samples
{_1 →…→_n }
whose distribution follows the posterior. Any given sample _i,
however, is more likely than not to be 
correlated with both the previous sample in the 
sequence _i-1 and the subsequent sample in the sequence
_i+1.

This occurs for two reasons. First, new positions _i drawn
from (_i|_i-1) by construction tend to depend
on the current position _i-1. This means that the position
we propose at iteration i+1 from will 
be correlated with the position at iteration i, 
which itself will be correlated with the position at
iteration i-1, etc. 

Second, even if we set
('|)=(') so that all
of our proposed positions are uncorrelated, 
our transition probability T('|) still ensures 
that we will eventually reject the new position 
so that _i+1=_i.
Since samples at exactly the same position are maximally correlated,
this ensures that samples from our chain will “on average” have non-zero
correlations. Note that having low acceptance fractions 
(i.e. the fraction of proposals that are accepted rather than rejected)
will lead to a larger fraction of the chain containing 
these perfectly correlated samples, increasing the overall correlation.

As mentioned in <ref>, correlated samples
provide less information about the underlying distribution they
are sampled from since their behavior doesn't just depend on the
underlying distribution but also the neighboring samples in the
sequence. Samples that are more highly correlated then should 
lead to a reduced ESS.

This intuition can be quantified by introducing
the auto-covariance C(t) 
for some integer lag t. Assuming that we have an infinitely long chain
{_1→…},
the auto-covariance C(t) is:

    C(t) ≡(_i - )
        · (_i+t - )i
        = lim_n→∞1/n∑_i=1^n (_i - ) · (_i+t - )

where · is the dot product. In other words, 
we want to know the covariance between _i at some iteration
i and _i+t at some other iteration i+t, averaged
over all all possible pairs of samples (_i, _i+t)
in our infinitely long chain. Note that the amplitude |C(t)|
will be maximized at |C(t=0)|,
where the two samples being compared are identical,
and minimized with |C(t)|=0 when 
_i and _i+t are completely independent
from each other.

Using the auto-covariance, we can define the corresponding
auto-correlation A(t) as

    A(t) ≡C(t)/C(0)

This now measures the average degree of correlation between 
samples separated by an integer lag t. In the case where t=0,
both samples are identical and A(t=0) = 1. In the case
where the samples are uncorrelated over lag t,
A(t) = 0.

The overall auto-correlation time
for our chain is just the auto-correlation A(t) summed over
all non-zero lags (t ≠ 0):

    τ≡∑_t=-∞^∞ A(t) - 1
        = 2 ∑_t=1^∞ A(t)

where the -1 comes from the fact that
the auto-correlation with no lag is just A(t=0) = 1
(i.e. each sample perfectly correlates with itself) and
the substitution arises from the fact that A(t) = A(-t) by symmetry.
If τ = 0, then it takes
no time at all for samples to become uncorrelated and the
samples can be assumed to be iid. If τ > 0, then
it takes on average τ additional iterations for
samples to become uncorrelated. 
An illustration of this process is shown in 
 <ref>.



Incorporating the auto-correlation time leads directly
to a modified definition for the ESS:

    n_ eff' ≡n_ eff/1 + τ

In practice, we cannot precisely compute τ since
we do not have an infinite number of samples and do not know
(). Therefore we often need to
generate an estimate τ̂ of the auto-correlation time
using the existing set of n samples we have. While discussing
various approaches taken to derive τ̂ is beyond the
scope of this work, please see <cit.> 
for additional details.

The fact that MCMC methods are subject to non-negative
auto-correlation times (τ≥ 0) but have optimal importance
weights w̃_i = 1 give an ESS of

    n_ eff,MCMC' = n_ eff,MCMC/1 + τ
        = n/1 + τ≤ n

This means that there is no guarantee
that MCMC is always the optimal choice to achieve
the largest ESS. In particular,
Importance Sampling methods, which can generate
fully iid samples with no auto-correlation time (τ = 0)
but non-optimal importance weights w̃_i, instead have
an ESS of

    n_ eff,IS' = n_ eff,IS/1 + τ
        = n_ eff,IS 
        = (∑_i=1^nw̃_i)^2/∑_i=1^nw̃_i^2≤ n

which can be greater than n_ eff, MCMC' at fixed n.

Given the results above, it should now be clear
that the central motivating concern of MCMC methods
is whether they can generate a chain of samples
with an auto-correlation time small enough to outperform
Importance Sampling. Whether or not
this is true will depend on the posterior, the approach used to generate
the chain of samples (see <ref> and <ref>)
and the proposal distribution () 
used for Importance Sampling (see <ref>).



 §.§ Exercise: MCMC over a 2-D Gaussian
 



  §.§.§ Setup


Let's again return to our examples from 
<ref> and <ref>, in which
our unnormalized posterior is well-approximated by a 2-D Gaussian (Normal)
distribution:

    (x,y) 
        = exp{-1/2[(x-μ_x)^2/σ_x^2
        + (y-μ_y)^2/σ_y^2]}

where (μ_x,μ_y)=(-0.3,0.8) and (σ_x^2,σ_y^2)=(2,0.5).

We want to use MCMC to approximate various posterior
integrals from this distribution.
We will start by choosing our proposal distribution (x',y'|x,y)
to be a 2-D Gaussian with a mean of 0 and standard deviation of 1:

    (x',y'|x,y) = (μ_x,μ_y)=(x,y)(σ_x,σ_y)=(1,1)




  §.§.§ Parameter Estimation


Using the above proposal, generate n=1000 samples following the MH algorithm starting
from the position (x_0,y_0)=(0,0). Using these samples,
compute an estimate of the means x
and y as well as the corresponding 68% credible
intervals (or closest approximation)
[x_ low, x_ high] and [y_ low, y_ high].
How accurate are each of these quantities compared with the values we might
expect?



  §.§.§ Evidence Estimation


Next, use a set of 10 × 10 bins from
x=[-5, 5] and y=[-5, 5] to construct an estimate ρ(x,y) 
from the resulting set of samples. Using this estimate
for the density, compute an estimate of the evidence .
How accurate is our approximation? Does it substantially
change if we adjust the number and/or size of the bins?



  §.§.§ Auto-Correlation Time and Effective Sample Size


Use numerical methods to compute an estimate of the
auto-correlation time τ and the corresponding effective sample size n_ eff.
How efficient is our sampling (n_ eff/n) compared to the default
Importance Sampling approach from the exercise in <ref>?
Does this mirror what we'd expect given the acceptance fraction of our proposals?
What do these quantities this tell us about how well 
our proposal (x,y) matches the structure of the underlying posterior
(x,y)?



  §.§.§ Uncertainties


Repeat the above exercises m=30 times to get an
estimate for how much our estimates of each quantity can vary.
Is the variation in line with what might be expected given 
the typical effective sample size?



  §.§.§ Consistency and Convergence


Now repeat the above exercise using n=2500 and n=10000 samples
points and comment on any differences.
How much has the overall accuracy improved? Do
the estimates appear convergent and consistent as n_ eff increases? 
How much do the errors on quantities shrink as a function 
of n and/or n_ eff? Is this similar or different
from the observed dependence from the Importance Sampling exercise
in <ref>?



  §.§.§ Sampling Efficiency


Next, adjust the (σ_x, σ_y) of the proposal distribution
to try and improve n_ eff at fixed n. How close is the final
ratio σ_x/σ_y of our proposal to that of the underlying
posterior? Are there any additional scaling differences between the rough
size of our proposal (x',y'|x,y) relative to the
underlying posterior (x,y)?
Given that (x,y) may differ from the
structure assumed when picking (x',y'|x,y),
can you think of any possible scheme to try and adjust our proposal
using an existing set of samples?



  §.§.§ Burn-In


Finally, adjust the starting position to be at
(x_0,y_0)=(10,10) instead of (0,0) and generate a new chain of 
samples. Plot the x and y positions of the chain over time.
Are there any obvious signs of the burn-in period? How many samples
roughly should be assigned to burn-in and subsequently removed from
our chain? Are there any possible heuristics that might help to identify
the initial burn-in period?



§ SAMPLING THE POSTERIOR WITH MCMC
 

The approach by which MCMC methods are able to generate a chain of
samples immediately gives a mental image of our chain “exploring” the
posterior. While it is true that the density of samples from the chain
ρ() →() as n →∞, 
the primary purpose of MCMC is estimating expectation values
f().
Although this might seem like a subtle difference,
this distinction is actually crucial for understanding how MCMC
algorithms (should) behave in practice. We discuss this in more
detail below.



 §.§ Approximating the Posterior
 

Although algorithms such as MH (<ref>) are constructed to
ensure the density of the chain of samples ρ()
generated by MCMC converges to the posterior ()
as n →∞, this does not necessarily translate
into an efficient method to approximate the posterior in practice. In
other words, n might need to be extremely large for this constraint to hold.
So how many samples do we need to ensure ρ() is a good approximation
to ()?

To start, we first need to define some metric for what a “good”
approximation is. A reasonable one might be that we would like
to know the posterior within some region δ_
to within some precision ϵ so that

    | 1/n∑_i=1^n_i ∈δ_ 
        - ∫_δ_() |
        ≡ |p̂(δ_) - p(δ_)| < ϵ

where p(δ_) is the total probability contained within
δ_ and p̂(δ_) is the fraction of
the MCMC chain of samples contained within the same region. While it might
seem strange to only estimate this for one region, I will shortly
generalize this to encompass the
entire[Technically the procedure outlined in this section
only works for finite volumes. The basic intuition, however, holds even
when parameters are unbounded although proving those results is beyond
the scope of this work.] posterior.

In the ideal case where our samples are iid and drawn from (),
our samples each have a probability p(δ_) 
of being within δ_.
The probability that p̂(δ_) = m/n then follows the
binomial distribution:

    P(p̂(δ_) = m/n) 
        = nm[p(δ_)]^m
        [1 - p(δ_) ]^n-m

In other words, our samples end up inside δ_ a total of m times
with probability p(δ_) and outside δ_ a
total of n-m times with probability 1 - p(δ_). The
additional binomial coefficient nm for “n choose m” 
accounts for all possible unique cases where m samples can
end up within δ_ out of our total sample size of n.

This distribution has a mean of p(δ_), so for
any finite n we expect p̂(δ_) to be an
unbiased estimator of p(δ_):

    p̂(δ_) - p(δ_) 
        = p(δ_) - p(δ_) = 0

The variance, however, depends on the sample size:

    |p̂(δ_) - p(δ_)|^2
        = p(δ_) [1 - p(δ_)]/n


In practice, we can expect there to be some non-zero auto-correlation
time τ > 0. This will increase the number of MCMC samples
we will need to generate to be confident that our estimate
p̂(δ_) is well-behaved. Inserting a factor
of 1+τ and substituting our expectation value from above
into our accuracy constraint then gives a
rough constraint for the number of samples n we would require
as a function of ϵ:

    n ≳p(δ_) [1 - p(δ_)]/ϵ^2/(1+τ)∼p̂(δ_) 
        [1 - p̂(δ_)]/ϵ^2× (1+τ̂)

The final substitution of p(δ_) and τ with their
noisy estimates p̂(δ_) and τ̂
arises from the fact that
in practice we don't know p(δ_) or 
τ (both of which require full knowledge of the posterior).
We are therefore forced to rely on estimators
derived from our set of n samples.

Let's now examine this result more closely. As expected,
the total number of samples is proportional
to 1 + τ̂: if it takes longer to generate independent
samples, then we need more samples to be confident we have characterized
the posterior well in a given region. 
We also see that n ∝ϵ^-2, so that
if we want to reduce the error by a factor of x we need
to increase our sample size by a factor of x^2.

The behavior in the numerator is more interesting. Note
that p̂(δ_) [1 - p̂(δ_)]
is maximized for p̂(δ_) = 0.5, and so the largest
sample size needed is when we have split our posterior directly in half.
In all other cases the sample size needed will be smaller because
there will be more samples outside or inside the region of interest
whose information we can leverage. The exact value of p̂(δ_)
of course depends on both the posterior ()
and the target region δ_: the sample size needed to approximate
the posterior to some ϵ 
near the peak of the distribution (the small region where () is large)
will likely be different than
the sample size needed to accurately estimate the tails of the
distribution (the large region where () is small).

While the above argument holds if we are looking to estimate the
posterior in just one region, “converging to the posterior”
implies that we want ρ()
to become a good approximation to () everywhere.
We can enforce this new requirement by splitting our posterior into m
different sub-regions {δ__1, …, δ__m} 
and requiring that each sub-region is well constrained:

    |p̂(δ__1) - p(δ__1)| < ϵ_1
            …    
        |p̂(δ__m) - p(δ__m)| < ϵ_m

Substituting in the expected errors on each of these
constraints then gives us an approximate limit
on the number of samples n_j that we need
to estimate the posterior in each region δ__j:

    n_j ≳p̂(δ__j) 
        [1 - p̂(δ__j)]/ϵ_j^2× (1+τ̂)

The total number of samples we need is then simply:

    n ≳∑_j=1^m n_j


This approach of dividing up our posterior into sub-regions
is conceptually similar to the grid-based approaches
described in <ref>. As such,
it is also subject to the same drawbacks:
we expect the number of regions m to increase
exponentially with the number of dimensions d.
For instance, if we just wanted to divide our posterior up
into m orthants we would end up with m=2^d regions:
2 in 1-D (left-right), 4 in 2-D (upper-left, lower-left,
upper-right, lower-right), 8 in 3-D, etc.

This effect implies that we should in general expect
the number of samples required to ensure ρ() is a good
approximation to () for some specified accuracy ϵ
to scale as

    n ≳ k^d

where k is a constant that depends on the accuracy requirements.
This puts approximating the full posterior firmly in the 
“curse of dimensionality” regime (see <ref>).[
A direct corollary of this result is that, while the
evidence estimates from MCMC are consistent, 
the rate of convergence to the underlying value 
will proceed exponentially more slowly
as d increases.]

While many practitioners talk about MCMC being 
an efficient method to “approximate the posterior”, in practice it is rarely
used to approximate () directly.
As discussed in <ref> and shown in
 <ref>, almost all quantities 
that are reported in the literature do not rely on approximations to
the full d-dimensional posterior, but rather approximations to marginalized
distributions that are almost always restricted to no more 
than k≲3 parameters at a time.
The act of marginalizing over the remaining d-k parameters
helps to counteract the curse of dimensionality illustrated here.
While it is technically fair to say that MCMC can “explore” 
the marginalized k-D posteriors for certain limited sets of parameters, 
this type of language can often lead to more misconceptions than insights.



 §.§ Posterior Volume
 

The basic consequences outlined in <ref> 
are more general than the specific case where we imagine dividing
up the posterior into orthants or other regions.
Fundamentally, computing any expectation
over the posterior f() requires
integrating over the entire domain of our parameters 
. We therefore want to understand how the volume
of this domain behaves (i.e. how many parameter combinations there are).
Once we have a grasp on how this behaves, we can then starting trying
to quantify how this will impact our estimates.




To start, let's consider the
d-dimensional hyper-cube (the d-cube)
with side length ℓ in all d dimensions. Its volume scales as

    V(ℓ) = ∏_i=1^dℓ = ℓ^d

The differential volume element between ℓ and
ℓ + ℓ is

    V(ℓ) = (d ×ℓ^d-1) × (ℓ) ∝ℓ^d-1


This exponential scaling with dimensionality means that volume becomes
increasingly concentrated in thin shells located in regions located
progressively further away from the center of the d-cube. 
As an example, consider the length-scale

    ℓ_50 = 2^-1/dℓ

that divides the d-cube into two equal-sized regions
with 50% of the volume contained interior to ℓ_50 and
50% of the volume exterior to ℓ_50. 
In 1-D, this gives ℓ_50/ℓ = 0.5 as we'd
expect. In 2-D, this gives ℓ_50/ℓ≈ 0.7. In 3-D,
ℓ_50/ℓ≈ 0.8. In 7-D, ℓ_50/ℓ≈ 0.9.
By the time we get to 15-D, we have ℓ_50 /ℓ≈ 0.95,
which means that 50% of the volume is located in the last 5% of the 
length-scale near the boundary of the d-cube.
While the constants may change when considering other shapes (e.g., spheres),
in general this exponential scaling as a function of d is a generic feature
of higher-dimensional volumes. In other words, increasing the number
of parameters leads to an exponential increase in the number of
available parameter combinations that we have to explore.



In addition to affecting the long-term behavior of MCMC, 
this exponential increase in volume also directly
impacts how MCMC methods operate. To see why this is the case, we
need look no further than the transition probability used in the
MH algorithm discussed in <ref>:

    T(_i+1|_i)
        ≡min[1, (_i+1)/(_i)(_i|_i+1)/(_i+1|_i)]

The non-trivial portion of this expression cleanly splits into two terms.
The first is dependent on the volume and is related to
how we proposed our next position from ('|).
The second is dependent on the density and is related to
how the posterior density changes between the two positions.

In practice, our transition probability can be interpreted as
a basic corrective approach: after proposing
a new position from some nearby volume, we then try to “correct”
for differences between our proposal and the underlying posterior
by only accepting these moves sometimes based on changes in the underlying density.
In high dimensions, this basic “tug of war” between the volume (proposal) and
the density (posterior) can break down as the vast majority of an object's 
volume becomes concentrated near the outer edges.[Alternative
methods such as Hamiltonian Monte Carlo <cit.> can get around
this problem by smoothly incorporating changes in the density and volume.]
For instance, in the case where our
proposal ('|) is a cube with side-length ℓ
centered on , this leads to a median length-scale of
ℓ_50 = 2^-1/dℓ, which increases rapidly from 0.5 ℓ
to ≈ℓ as the dimensionality increases. The same logic
also applies to other proposal distributions (see <ref>).
This focus on positions either far away or with very similar separation
length-scales as ℓ_50→ℓ 
means that many choices of ('|)
have a tendency to “overshoot”, proposing new positions with much
smaller posterior densities compared to the current position. These
new positions are then almost always rejected, leading to extremely
low acceptance fractions and correspondingly
long auto-correlation times. An example of this effect is
illustrated in  <ref>.

One of the main ways to counteract this behavior is to
adjust the size/shape of the proposal ('|)
so that the fraction of proposed positions that are accepted remains
sufficiently high. This helps to ensure the posterior density 
() does not change too drastically when proposing 
positions new positions, leading to lower overall auto-correlation times.
Details of how to implement these schemes in practice are beyond 
the scope of this article; please see
citation for additional details.



 §.§ Posterior Mass and Typical Sets
 

Above, I described how the behavior of volume in high dimensions
can impact the performance of our MCMC MH sampling algorithm,
possibly leading to inefficient proposals and 
low acceptance fractions. Let's assume that we have resolved this
problem and have an efficient way of generating our chain
of samples. We now have a secondary question: where are these
samples located?

From our discussion in <ref>, we know that
the highest density of samples ρ()
will be located where the posterior density ()
is also correspondingly high. However, this region δ_
might only correspond to a small portion of the posterior. 
Indeed, given there is exponentially more volume as the dimensionality
increases, it is almost guaranteed that models with many parameters
 will have the vast majority of the posterior located
outside the region of highest density. 

A consequence of this is that
the majority of samples in our chain will be located
away from the peak density. As a result, our chain
spends the majority of its time
generating samples in these regions.
This has a huge impact in the way our chain is expected to behave:
while the highest concentration
of samples will be located in the regions of highest
posterior density, the largest amount of samples
will actually be located in the regions of highest posterior
mass (i.e. density times volume). 
Since this implies that a “typical” sample 
(picked at random) will most likely be located in this
region of high posterior mass,
this region is also commonly referred to
as the typical set.

To make this argument a little easier to conceptualize,
let's imagine that we have a 3-parameter model = (x, y, z) 
and (x,y,z) is spherically symmetric.
While we could imagine trying to
integrate over (x,y,z) directly in terms of
x  y  z, it is almost always
easier to instead integrate over such a distribution
in “shells” with differential volume
V(r) = 4 π r^2  r
as a function of radius
r = √(x^2 + y^2 + z^2). This allows us
to rewrite the 3-D integral over (x,y,z) as a 1-D integral over r:

    ∫(x,y,z)  x  y  z
        = ∫(r) 4 π r^2  r
        ≡∫'(r)  r

where '(r) ≡ 4 π r^2 (r) is now
the 1-D density as a function of r. This “boosts” the contribution
as a function of r by the differential volume element of the
shell associated with (r), and implies that the
the posterior should have some sort of shell-like structure (i.e.
'(r) is maximized for r > 0).

Although not all posterior densities can be expected to
be spherically-symmetric in this way, in general we can
rewrite the d-D integral over  as a 1-D volume integral
over V defined by some unknown iso-posterior 
contours[Indeed, alternative Monte Carlo methods such as 
Nested Sampling <cit.> or 
Bridge/Path Sampling <cit.> actually are designed to 
evaluate this type of volume integral explicitly.]

    ∫() 
        = ∫(V)  V

As outlined in <ref>, we generically
expect the size of each volume element to go as
V ∼ r^d-1 r where r is the
distance from the peak of posterior. So the basic
intuition we get from the simple spherically-symmetric case
still applies and we expect

    ∫(V)  V 
        ∼∫(r) r^d-1 r
        = ∫'(r)  r


As before, the differential volume element of the
shell associated with (r) “boosts” its overall
contribution as a function of r. This boost also becomes
exponentially stronger as d increase. 
For even moderately-sized d,
we therefore expect the posterior mass
to be mostly contained in a thin shell located at a radius r'
with some width Δ r'.
See  <ref> for an illustration of this effect
based on the toy problem presented in <ref>.



This result has two immediate implications. First, the majority
of our samples are not located where the posterior
density is maximized. This is the result of an exponentially
increasing number of parameter combinations, which allow a
small handful of excellent fits to the data to
be easily overwhelmed by a substantially larger number
of mediocre fits. MCMC methods are therefore generally
inefficient at locating and/or characterizing the region of
peak posterior density.

Second, as d increases we generally would expect the radius of the
shell containing the bulk of the posterior mass to increase, moving further
and further away from the peak density due to the exponentially increasing
available volume. Since the majority of our samples are located in this
region, our chain will spend the vast majority 
of time generating samples from this shell.

This allows us to now outline exactly why it is challenging
to propose samples efficiently in high dimensions:

    
  * To make sure our acceptance fractions remain reasonable, 
    we need to ensure our proposed positions mostly lie within 
    this shell of posterior mass.
    
  * However, obtaining an independent sample requires 
    being able to (in theory) propose any position within this shell.
    
  * This means that our auto-correlation time
    will principally be set by how long it takes to “wander around”
    the shell, which will be a function of its overall size r',
    its width Δ r', and the number of dimensions d.




§ APPLICATION TO A SIMPLE TOY PROBLEM
 

I now consider a concrete, detailed example
to illustrate how all the concepts discussed in <ref> and
<ref> come together in practice. 
Throughout this section, I will outline a number of analytic
results and utilize several different MCMC sampling strategies
to generate chains of samples. I strongly encourage interested readers
to implement their own versions of the methods outlined here,
which can be used to reproduce the numerical results from this section
in their entirety.



 §.§ Toy Problem
 

In this toy problem, we will take our (unnormalized) 
posterior to be a d-dimensional 
Gaussian (Normal) distribution with a mean
of μ = 0 and a standard deviation of σ
in all dimensions:

    () 
        = exp[-1/2||^2/σ^2]

where ||^2 = ∑_i=1^dΘ_i^2 is the squared magnitude
of the position vector.

Based on the results from <ref>, 
we can better understand the properties of this distribution by
rewriting the posterior density
in terms of the “radius”
r ≡ |  | = √(∑_i=1^dΘ_i^2)
away from the center:

    (r) = exp[-r^2/2σ^2]

The corresponding volume contained within a given radius r is then

    V(r) ∝ r^d

The corresponding posterior mass is '(r)
is then defined via

    (V)  V(r)
        ∝ e^-r^2/2σ^2 r^d-1 r ≡'(r)  r

Note that this is closely related
to the chi-square distribution.

The typical radius r_ peak 
where the posterior mass peaks (i.e. is maximized)
and a sample is most likely to be located
can be derived by setting '(r)/ r = 0.
Solving this gives

    r_ peak = √(d-1)σ

In other words, while in 1-D a typical sample
is most likely to be located at the peak of the 
distribution with r_ peak = 0, in higher dimensions
this changes quite drastically. While r_ peak = 1σ
in 2-D, it is 2σ in 5-D, 3σ in 10-D, and 5σ
in 26-D. This is a direct consequence of the
huge amount of volume at larger radii in high dimensions:
although a sample at r=5σ has a posterior density (r)
orders of magnitude worse than a sample at r=0, the enormous number of
parameter combinations (volume) available at r=5σ
more than makes up for it.

In general, we expect the posterior mass to comprise
a “Gaussian shell” centered at some radius

    r_ mean≡r'
        = ∫_0^∞ r '(r)  r
        = √(2)Γ(d+1/2)/Γ(d/2)σ≈√(d)σ

with a standard deviation of

    Δ r_ mean≡√((r - r_ mean)^2')
        = σ√(d - 2 
        (Γ(d+1/2)/Γ(d/2))^2)≈σ/√(2)

where Γ(d) is the Gamma function and the
approximations are taken for large d.
See  <ref> for an
illustration of this behavior.



 §.§ MCMC with Gaussian Proposals
 

Let us now consider a chain of samples
{_1 →…→_n }.
The distance between two samples
_m and _m+t 
separated by some lag t will be

    | - '|
        = √(∑_i=1^d (Θ_m,i - Θ_m+t,i)^2)

Assuming that the lag t ≫τ is
substantially larger than the auto-correlation time τ,
we can assume each sample is approximately iid distributed following
our Gaussian posterior. This then gives an expected separation of

    Δ r_ sep≡√(|_m - _m+t|^2)
        = √(∑_i=1^d(Θ_m,i - Θ_m+t,i)^2)
        = √(2d)σ≈√(2) r_ mean


We can in theory propose samples in such a way so that the separation
|_i+1 - _i|
between a proposed position _i+1 
and the current position _i follows the
ideal separation of √(2)r_ mean derived above
by using a simple Gaussian proposal distribution:

    (_i+1|_i)
        ∝exp[-1/2|_i+1 - _i|^2/2σ^2]

While this proposal has the same shape as the posterior,
it is centered on _i rather than 0. Using our intuition for
how volume behaves based on <ref>, we can
conclude that the majority of samples proposed from 
this choice of ('|) will probably have little overlap
with the posterior. 



Indeed, numerical simulation suggests
the typical fraction of positions that will be accepted given the above proposal
roughly scales as

    ⟨ f_ acc(d) ⟩≡exp[ln T(_i+1|_i),]
        ∼exp[-d/4 - 1/2]

which decreases exponentially as the dimensionality increases,
similar to  <ref>. Likewise, we find
the auto-correlation time roughly scales as

    ⟨τ(d) ⟩≡exp[lnτ,]
        ∼exp[d/4 + 7/4]

This exponential dependence
arises because the overlap between the typical Gaussian proposal
('|) and the underlying posterior 
() essentially reduces to the small volume
where two thin shells overlap. Since the radii of the shells
goes as ∝√(d) while the widths remain roughly constant,
the “fractional size” of the shell (and the corresponding overlap)
ends up decreasing exponentially.

To counteract this effect, we need to adjust the σ of our
proposal distribution by some factor γ:

    _γ(_i+1|_i)
        ∝exp[-1/2|_i+1 - _i|^2/(γσ)^2]

where our previous proposal assumes γ = √(2).
If we want to ensure our typical acceptance fraction will remain
roughly constant as a function of dimension d, γ needs to scale as

    ⟨ f_ acc(γ(d)) ⟩≈ C
          ⇒  γ(d) ∝1/√(d)

which inversely tracks the expected radius r_ mean
of the typical set. We find that taking

    γ = δ/√(d)

leads to a typical acceptance fraction of

    ⟨ f_ acc(δ/√(d)) ⟩≈exp[-(δ^2/4)^2 
        - δ/2]

as d becomes large with a typical auto-correlation time of

    ⟨τ(δ/√(d)) ⟩≈ 3 d

for reasonable choices of δ.
This linear dependence is a substantial improvement over
our earlier exponential scaling.



  §.§.§ Numerical Tests
 

To confirm these results, I sample from this d-dimensional
Gaussian posterior (assuming σ=1 for simplicity) using
two MH MCMC algorithms for n=20,000 iterations
based on these proposal distributions.
The first proposes new points 
assuming γ=√(2). The second
assumes γ=2.5/√(d) in order
to maintain a roughly constant acceptance fraction of 25%. 
As shown in  <ref>, the chains
behave as expected given our theoretical predictions as a function
of dimensionality, with the constant proposal quickly becoming stuck
while the adaptive proposal continues sampling normally. While
the auto-correlation time τ increases in both cases, 
the increase in the latter case (where it is driven by decreasing
size/scale of the proposal distribution)
is much more manageable than the former (where it is driven by
the exponentially decreasing acceptance fraction).



 §.§ MCMC with Ensemble Proposals
 

One drawback to the Gaussian proposals explored above is that
we have to specify the structure of the distribution ahead of time.
In this specific case, we assumed that:

    
  * the width of the posterior in each dimension (parameter)
    was constant such that σ_1 = σ_2 = … = σ_n = σ and
    
  * the parameters were entirely uncorrelated with each other such that
    the correlation coefficient ρ_ij = 0 between any two dimensions
    i and j.


In general, there is no good reason to assume that either of these are true.
This means we have to also estimate the entire set of
d(d+1)/2 free parameters that determine the overall covariance
structure of our unknown posterior distribution.
Trying to adjust the covariance structure in order
to improve our sampling efficiency and decrease the auto-correlation
time (see <ref> and <ref>)
becomes one of the most difficult parts of running MCMC algorithms in practice.

While there are schemes to perform these adjustments during
an extended burn-in period (see, e.g., ), there
is significant appeal in methods that can “auto-tune” without
much additional input from the user. One class of such approaches
are known as ensemble or particle methods.
These methods attempt to use many m chains running
simultaneously (i.e. in parallel) to improve the performance 
of any individual chain.

We explore three variations of ensemble methods here
that attempt to exploit m ≳ d(d+1)/2 chains running
simultaneously:

    
  * using the ensemble of particles to condition
    a Gaussian proposal distribution,
    
  * using trajectories from multiple particles along with
    Gaussian “jitter”, and
    
  * using affine-invariant transformations of
    trajectories from multiple particles.

A schematic illustration of these 
methods is shown in  <ref>.

As we might expect, an immediate drawback of these methods is
they rely on having enough particles to characterize
the overall structure of the space (i.e. the curse of dimensionality).
While this limits their utility when sampling from high-dimensional spaces, 
they can be attractive options in moderate-dimensional spaces 
(d ≲ 25) where a few hundred particles
are often sufficient to ensure reasonable performance.





  §.§.§ Gaussian Proposal
 

The first approach is simply a modified
Gaussian proposal: at any iteration i for any chain
j, we propose a new position _i+1^j
based on the current position _i^j
using a Gaussian proposal

    _γ^j(_i+1^j|_i^j)
        ∝exp[-1/2
        (_i+1^j - _i^j)^ T
        (γ^2_i^j)^-1
        (_i+1^j - _i^j)]

where T is the transpose operator and

    _i^j =  Cov[{_i^1, …, 
        _i^j-1, _i^j+1, …, _i^m}]

is the empirical covariance matrix estimated from the
current positions of the m chains excluding chain j.
We repeat this process for each of the m chains in turn.

In other words, at each iteration i we want to update all
m chains. We do so by updating each chain j in turn based on what
the other chains are currently doing. Assuming the current position of
each chain is distributed following the underlying posterior
(), it is straightforward to show that
_i^j is a reasonable approximation to the unknown
covariance structure of our posterior. In addition,
because we exclude j when computing _i^j, this proposal
is symmetric going from _i^j→_i+1^j
and from _i+1^j→_i^j.
This means that we satisfy detailed balance and do not
have to incorporate any proposal-dependent factors when
computing the transition probability.



  §.§.§ Ensemble Trajectories with a Gaussian Proposal
 


The approach taken in <ref>
solves the problem of trying to tune
the covariance of our initial Gaussian proposal. However, it still
assumes that a Gaussian proposal is the optimal solution.
A more general approach is one that does not rely on assuming a
proposal explicitly, but rather only relies on the distribution of
the remaining particles. 

One such approach used in the literature is
Differential Evolution MCMC <cit.>.
The main idea behind DE-MCMC is to rely on the relative positions
of the chains at a given iteration i when making new proposals. We first
randomly select two other particles k and l where
^j_i ≠^k_i ≠^l_i. We then propose
a new position based on the vector distance between the other
two particles ^k_i - ^l_i with some scaling γ
along with some additional “jitter” ϵ:

    _i+1^j 
        = _i^j + γ× (^k_i - ^l_i + ϵ)


In the case where the behavior of chains k and l are approximately
independent of each other and assuming the
underlying posterior distribution ()
is Gaussian with some unknown
mean  and covariance 
(and “standard deviation” ^1/2), it is straightforward to
show that the distribution of ^k_i - ^l_i will then follow

    ^k_i - _l ∼0(2)^1/2

Typically, the jitter ϵ is chosen to also be Gaussian distributed
with covariance _ϵ such that

    ϵ∼0_ϵ^1/2

In general, _ϵ is mostly
used to try and avoid issues caused by finite particle
sampling: since the number of unique trajectories (ignoring symmetry) is

    n_ traj 
        = m-12 
        = (m-1)!/2!(m-3)!
        = (m-1)(m-2)/2

if m is sufficiently small the DE-MCMC
procedure can only explore a small number of possible trajectories
at any given time, leading to extremely inefficient sampling.

Combined, this implies that the
proposed position has a distribution of

    _i+1^j∼_i^jγ× (2+_ϵ)^1/2

This shows that the 3-particle DE-MCMC procedure can generate
new positions in a manner analogous to the ensemble Gaussian proposal
we first discussed.



  §.§.§ Affine-Invariant Transformations of Ensemble Trajectories



Another approach used in the literature 
<cit.>
is the Affine-Invariant “stretch move” from 
<cit.>. This 
uses only one additional particle _i^k rather than two:

    _i+1^j 
        = _i^k + γ× (^j_i - ^k_i)

In place of the jitter term ϵ from DE-MCMC,
the stretch move instead injects some amount of 
randomness by allowing γ to vary.
By sampling γ from some probabilty distribution 
g(γ), we allow the proposals
to explore various “stretches” of the direction vector.
As shown in <cit.>, if this function is chosen such that

    g(γ^-1) = γ× g(γ)

then this proposal is symmetric. Typically, g(γ) is chosen
to be

    g(γ|a) = 
        γ^-1/2    a^-1≤γ≤ a 
    
        0     otherwise

where a=2 is often taken as a typical value.
Note that when γ = 1, this move leaves 
_i+1^j = _i^j unchanged.

Compared to DEMCMC, the stretch move appears to have one clear advantage:
it doesn't have any reliance on some “jitter” term ϵ that reintroduces
scale-dependence into the proposal. That makes the proposal invariant to
affine transformations and only
sensitive to a single parameter a, which governs the range of
scales the stretch factor γ is allowed to explore.

This lack of jitter, however, is not substantially advantageous in practice. 
As noted in <ref>,
ϵ is really designed to avoid possible
degeneracies due to the limited number of available trajectories. In that
case we had (m-1)(m-2)/2 ∼ m^2/2 possible trajectories; here, however, we only
have m (since ^j_i is always included). This is a much
smaller number of possible trajectories at a given m, 
making this particular proposal more susceptible to that particular effect.



In addition, because this proposal involves adjusting γ and therefore
the length of the trajectory itself, we need to consider how changing γ
affects the total volume of the sphere centered on ^j_i
with radius ^k_i-^j_i. As discussed in <ref>,
the differential volume increases as r^d-1. Therefore, increasing
or decreasing γ substantially adjusts the differential volume
in our proposal. This involves introducing a steep boost/penalty into our
transition probability, which now becomes:

    T(_i+1^j|_i^j, γ)
        = min[1, γ^d-1(^j_i+1)/(_i^j)]

This heavily favors proposals with γ > 1 (outwards) and heavily
disfavors proposals with γ < 1 as d increases to account for the
exponentially increasing volume at larger radii.

Finally, while this stretch move actually generates proposals 
in the right overall direction, it is not efficient at generating
samples within the bulk of the posterior mass
as the dimensionality increases.
As discussed in <ref>, given the typical
position of _i^j, the typical length-scale of
the proposed positions needs to shrink by
∝ 1/√(d) in order to guarantee our new sample
remains within the bulk of the posterior mass.
However, the form for g(γ|a) specified above 
instead ensures that γ will always be between 1/a and a.
Even if we attempt to account for this effect by letting 
a(d) → 1 as d →∞ in order to
target a constant acceptance fraction and ensure more overlap, 
the asymmetry of our proposal and the γ^d-1 term
in the transition probability systematically biases our
proposed and accepted positions compared with the ideal distribution.
This subsequently leads to lager auto-correlation times, mostly
counteracting any expected gains.



  §.§.§ Numerical Tests
 

To confirm these results, I sample from this d-dimensional
Gaussian posterior (assuming σ=1 for simplicity) 
using each of these ensemble MH MCMC algorithms with for n=1500
iterations with m=100 chains.
In the first case, I propose a new position for chain j 
at iteration i using a Gaussian distribution with
a covariance γ^2 _i^j computed over the remaining ensemble of
k ≠ j chains, where the scale factor γ=2.5/√(d) is
chosen to target a constant acceptance fraction of roughly 25%. 
In the second case, I propose new positions using the DE-MCMC algorithm
with a scale factor of γ=1.7/√(d) 
and additional Gaussian jitter with covariance _ϵ = _i^j / 5
derived from the remaining chains in the ensemble, again targeting
an acceptance fraction of roughly 25%.
In the third case, I propose new positions using the affine-invariant stretch
move assuming the typical form for g(γ|a) with a=2.[Allowing
a(d) to vary as a function of dimensionality to target a roughly constant
acceptance fraction gives similar results.]

As shown in  <ref>, the chains
behave as expected given our theoretical predictions as a function
of dimensionality. Similar to the adaptive Gaussian case,
the first two approaches continue 
sampling efficiently even as d increases. The affine-invariant
stretch move, however, experiences exponentially-decreasing efficiency
and struggles to sample the posterior effectively.



 §.§ Additional Comments
 

Before concluding, I wish to emphasize that the toy problem explored in this
section should only be interpreted as a tool to build intuition surrounding
how certain methods are expected to behave in a controlled environment.
While the behavior as a function of dimensionality
helps to illustrate common issues, in practice 
the performance of any method will depend on the specific problem, tuning parameters,
the time spent on tuning, and many other possible factors. Since it is always possible
to find problems for which any particular method
will perform well or poorly, I encourage users to try out a variety of
approaches to find the ones that work best for their problems.



§ CONCLUSION
 

Bayesian statistical methods have become increasingly prevalent in
modern scientific analysis as models have become more complex.
Exploring the inferences we can draw from these models often
requires the use of numerical techniques, the most popular of
which is known as Markov Chain Monte Carlo (MCMC).

In this article, I provide a conceptual introduction
to MCMC that seeks to highlight the what, why,
and how of the overall approach. I first give
an overview of Bayesian inference and discuss what
types of problems Bayesian inference generally is trying to solve, 
showing that most quantities we are interested in computing
require integrating over the posterior density.
I then outline approaches to computing these integrals using
grid-based approaches, and illustrate how adaptively changing
the resolution of the grid naturally transitions into 
the use of Monte Carlo methods. I illustrate how
different sampling strategies affect the overall efficiency
in order to motivate why we use MCMC methods.
I then discuss various details related to how MCMC methods 
work and examine their expected overall behavior based on
simple arguments derived from how volume and posterior density behave
as the number of parameters increases. 
Finally, I highlight the impact this conceptual understanding
has in practice by comparing the performance of various MCMC methods 
on a simple toy problem.

I hope that the material in this article, 
along with the exercises and applications,
serve as a useful resource that helps
build up intuition for how MCMC and other
Monte Carlo methods work.
This intuition should be helpful
when making decisions over when to apply MCMC methods to your own
problems over possible alternatives,
developing novel proposals and sampling strategies, 
and characterizing what issues you might
expect to encounter when doing so.



§ ACKNOWLEDGEMENTS


JSS is grateful to Rebecca Bleich for continuing to tolerate his
(over-)enthusiasm for sampling during their time together. 
He would also like to thank a number of people for
helping to provide much-needed feedback during earlier stages
of this work, including Catherine Zucker, Dom Pesce, Greg Green,
Kaisey Mandel, Joel Leja, David Hogg, Theron Carmichael, and Jane Huang.
He would also like to thank Ana Bonaca, Charlie Conroy, Ben Cook,
Daniel Eisenstein, Doug Finkbeiner, Boryana Hadzhiyska, 
Will Handley, Locke Patton, and Ioana Zelko
for helpful conversations surrounding the material.

JSS also wishes to thank Kaisey Mandel and
the Institute of Astronomy at the University
of Cambridge, Hans-Walter Rix and
the Galaxies and Cosmology Department at
the Max Planck Institute for Astronomy, and Renée
Hloz̆ek, Bryan Gaensler, and the Dunlap Institute
for Astronomy and Astrophysics at the University of Toronto
for their kindness and hospitality while hosting him over
the period where a portion of this work was being completed.

JSS acknowledges financial support from the National Science Foundation
Graduate Research Fellowship Program (Grant No. 1650114)
and the Harvard Data Science Initiative.


aasjournal










