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 posteriorbased 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.
['Statistics - Other Statistics', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Physics - Data Analysis', 'Statistics and Probability']
A Conceptual Introduction to Markov Chain Monte Carlo Methods
{'AConceptual Introduction to Markov Chain Monte Carlo Methods': 'Joshua S. Speagle  \nCenter for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA02138, USA  \njspeagle@cfa.harvard.edu', 'Abstract': '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.', '1 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.  \nTogether, 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.  \nOne 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.  \nOwing 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 Brooks et al. (2011) and Hogg & ForemanMackey (2018) along with associated references therein.  \nThis 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:  \n- 1. What problems are MCMC methods trying to solve?\n- 2. Why are we interested in using them?\n- 3. How do they work in theory and in practice?  \nWhen 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 Gelman et al. (2013) and Blitzstein &Hwang (2014) along with associated references therein.  \nThe outline of the article is as follows. In § 2, I provide a brief review of Bayesian inference and posterior distributions. In § 3, I discuss what posteriors are used for in practice, focusing on integration and marginalization. In § 4, I outline a basic scheme to approximate these posterior integrals using discrete grids. In § 5, I illustrate how Monte Carlo methods emerge as a natural extension of grid-based approaches. In § 6, I discuss how MCMC methods fit within the broader scope of possible approaches and their benefits  \nand drawbacks. In § 7, I explore the general challenges MCMC methods face. In § 8, I examine how these concepts come together in practice using a simple example. I conclude in § 9.', '2 Bayesian Inference': "In many scientific applications, we have access to some data D 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.  \nWecancombinethese pieces together to estimate the probability P ( D | Θ M , M ) that we would actually see that data D 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 ( D | Θ M , M ) of the parameters Θ M based on the observed data D ? Assuming different values of Θ M will give different likelihoods, telling us which parameter choices appear to best describe the data we observe.  \nIn Bayesian inference, we are interested in inferring the flipped quantity, P ( Θ M | D , M ) . This describes the probability that the underlying parameters are actually Θ M given our data D and assuming a particular model M . By using factoring of probability, we can relate this new probability P ( Θ M | D , M ) to the likelihood P ( D | Θ M , M ) described above as  \nP ( Θ M | D , M ) P ( D | M ) = P ( Θ M , D | M ) = P ( D | Θ M , M ) P ( Θ M | M ) (1)  \nwhere P ( Θ M , D | M ) represents the joint probability of having an underlying set of parameters Θ M that describe the data and observing the particular set of data D we have already collected.  \nRearranging this equality into a more convenient form gives us Bayes' Theorem :  \nP ( Θ M | D , M ) = P ( D | Θ M , M ) P ( Θ M | M ) P ( D | M ) (2)  \nThis equation now describes exactly how our two probabilities relate to each other.  \nP ( Θ 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.  \nThe denominator  \nP ( D | M ) = ∫ P ( D | Θ M , M ) P ( Θ M | M )d Θ M (3)  \nis 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  \nFigure 1: An illustration of Bayes' Theorem. The posterior probability P ( Θ ) (black) of our model parameters Θ is based on a combination of our prior beliefs π ( Θ ) (blue) and the likelihood L ( Θ ) (red), normalized by the overall evidence Z = ∫ π ( Θ ) L ( Θ )d Θ (purple) for our particular model. See § 2 for additional details.  \n<!-- image -->  \nour model M explains the data D 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 D , 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 D as a given, this often ends up being a constant.  \nFinally, P ( Θ M | D , M ) represents our posterior . This quantifies our belief in Θ M after combining our prior intuition P ( Θ M | M ) with current observations P ( D | Θ M , M ) and normalizing by the overall evidence P ( D | 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 Figure 1.  \nThroughout the rest of the paper I will write these four terms (likelihood, prior, evidence, posterior) using shorthand notation such that  \nP ( Θ ) ≡ L ( Θ ) π ( Θ ) ∫ L ( Θ ) π ( Θ )d Θ ≡ L ( Θ ) π ( Θ ) Z (4)  \nwhere P ( Θ ) ≡ P ( Θ M | D , M ) is the posterior, L ( Θ ) ≡ P ( D | Θ M , M ) is the likelihood,  \nπ ( Θ ) ≡ P ( Θ M | M ) is the prior, and the constant Z ≡ P ( D | M ) is the evidence. I have suppressed the model M and data D notation for convenience here since in most cases the data and model are considered fixed, but will re-introduce them as necessary.  \nBefore 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.", 'Setup': "Let's again return to our examples from § 4 and § 5, in which our unnormalized posterior is well-approximated by a 2-D Gaussian (Normal) distribution:  \n˜ P ( x, y ) = exp { -1 2 [ ( x -µ x ) 2 σ 2 x + ( y -µ y ) 2 σ 2 y ]}  \nwhere ( µ x , µ y ) = ( -0 . 3 , 0 . 8) and ( σ 2 x , σ 2 y ) = (2 , 0 . 5) .  \nWe want to use MCMC to approximate various posterior integrals from this distribution. We will start by choosing our proposal distribution Q ( x ' , y ' | x, y ) to be a 2-D Gaussian with a mean of 0 and standard deviation of 1 :  \nQ ( x ' , y ' | x, y ) = N [( µ x , µ y ) = ( x, y ) , ( σ x , σ y ) = (1 , 1)]", 'Problem': 'Using these assumptions, compute:  \n- 1. the prior π ( T ) ,\n- 2. the likelihood L ( T ) , and\n- 3. the posterior P ( T )  \ngiven 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?', '3 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 P ( Θ ) ∝ L ( Θ ) π ( Θ ) . 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:  \n- 1. Making educated guesses : make a reasonable guess at what the underlying model parameters are.\n- 2. Quantifying uncertainty : provide constraints on the range of possible model parameter values.\n- 3. Generating predictions : marginalize over uncertainties in the underlying model parameters to predict observables or other variables that depend on the model parameters.\n- 4. Comparing models : use the evidences from different models to determine which models are more favorable.  \nIn 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 Z , for instance, is again just the integral of the likelihood and the prior over all possible parameters:  \nZ = ∫ L ( Θ ) π ( Θ )d Θ ≡ ∫ ˜ P ( Θ )d Θ (5)  \nwhere ˜ P ( Θ ) ≡ L ( Θ ) π ( Θ ) is the unnormalized posterior.  \nLikewise, 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  \nstraightforward if the entire posterior over Θ is known:  \nP ( Θ int ) = ∫ P ( Θ int , Θ nuis ) d Θ nuis = ∫ P ( Θ )d Θ nuis (6)  \nOther quantities can generally be derived from the expectation value of various parameterdependent functions f ( Θ ) with respect to the posterior:  \nE P [ f ( Θ )] ≡ ∫ f ( Θ ) P ( Θ )d Θ ∫ P ( Θ )d Θ = ∫ f ( Θ ) ˜ P ( Θ )d Θ ∫ ˜ P ( Θ )d Θ = ∫ f ( Θ ) P ( Θ )d Θ (7)  \nsince ∫ P ( Θ )d Θ = 1 by definition and ˜ P ( Θ ) ∝ P ( Θ ) . 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.  \nTaken 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.  \nThis 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 E P [ f ( Θ )] even if we have an extremely poor estimate of P ( Θ ) or ˜ P ( Θ ) .  \nMore 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 Figure 2.", '3.1 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 P ( Θ ) to propose a point estimate ˆ Θ that we think is a pretty good guess for the true value Θ ∗ .  \nglyph[negationslash]  \nWhat 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 Θ ∗ .  \nUnfortunately, 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:  \nL P ( ˆ Θ ) ≡ E P [ L ( ˆ Θ | Θ ) ] = ∫ L ( ˆ Θ | Θ ) P ( Θ )d Θ (8)  \nFigure 2: A 'corner plot' showing an example of how posteriors are used in practice. Each of the top panels shows the 1-D marginalized posterior distribution for each parameter (grey), along with associated median point estimates (red) and 68% credible intervals (blue). Each central panel shows the 10%, 40%, 65%, and 85% credible regions for each 2-D marginalized posterior distribution. See § 3 for additional details.  \n<!-- image -->  \nA reasonable choice for ˆ Θ is then the value that minimizes this expected loss in place of the actual (unknown) loss:  \nˆ Θ ≡ argmin Θ ' [ L P ( Θ ' )] (9)  \nwhere argmin indicates the value (argument) of Θ ' that minimizes the expected loss L P ( Θ ' ) While this strategy can work for any arbitrary loss function, solving for ˆ Θ often requires using numerical methods and repeated integration over P ( Θ ) . 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.  \n.", '3.2 Quantifying Uncertainty': "In many cases we are not just interested in computing a prediction ˆ Θ for Θ ∗ , but also constraining a region C ( Θ ) of possible values within which Θ ∗ might lie with some amount of certainty. In other words, can we construct a region C X such that we believe there is an X % chance that it contains Θ ∗ ?  \nThere are many possible definitions for this credible region . One common definition is the region above some posterior threshold P X where X % of the posterior is contained, i.e. where  \n∫ Θ ∈C X P ( Θ )d Θ = X 100 (10)  \ngiven  \nC X ≡ { Θ : P ( Θ ) ≥ P X } (11)  \nIn other words, we want to integrate our posterior over all Θ where the value P ( Θ ) > P X is greater than some threshold P X , where P 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).  \nIn 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 p th percentile is defined as  \n∫ x p -∞ P ( x )d x = p 100 (12)  \nWe 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).  \nIn general, when referring to 'credible intervals' throughout the text the percentile definition should be assumed unless explicitly stated otherwise.", '3.3 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 P ( Θ ) over possible values Θ ∗ could take, however, to predict what will happen we will need to marginalize over this uncertainty.  \nWe can quantify this intuition using the posterior predictive P ( ˜ D | D ) , which represents the probability of seeing some new data ˜ D based on our existing data D :  \nP ( ˜ D | D ) ≡ ∫ P ( ˜ D | Θ ) P ( Θ | D )d Θ ≡ ∫ ˜ L ( Θ ) P ( Θ )d Θ = E P [ ˜ L ( Θ ) ] (13)  \nIn other words, for hypothetical data ˜ D , we want to compute the expected value of the likelihood ˜ L ( Θ ) over all possible values of Θ based on the current posterior P ( Θ ) .', '3.4 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.  \nWe can compare two models by computing the Bayes factor :  \nR 1 2 ≡ P ( M 1 | D ) P ( M 2 | D ) = P ( D | M 1 ) P ( M 1 ) P ( D | M 2 ) P ( M 2 ) ≡ Z 1 Z 2 π 1 π 2 (14)  \nwhere Z 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 R 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.  \nAgain, note that computing Z M requires computing the integral ∫ ˜ P ( Θ )d Θ of the unnormalized posterior ˜ P ( Θ ) 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.', 'Point Estimates': "The mean can be defined as the point estimate ˆ Θ that minimizes the expected loss L P ( ˆ Θ ) under squared loss :  \nL mean ( ˆ Θ | Θ ∗ ) = | ˆ Θ -Θ ∗ | 2  \nThe median can be defined as the point estimate that minimizes L P ( ˆ Θ ) under absolute loss :  \nL med ( ˆ Θ | Θ ∗ ) = | ˆ Θ -Θ ∗ |  \nAnd the mode can be defined as the point estimate that minimizes L P ( ˆ Θ ) under 'catastrophic' loss :  \nL mode ( ˆ Θ | Θ ∗ ) = -δ ( | ˆ Θ -Θ ∗ | )  \nwhere δ ( · ) is the Dirac delta function defined such that  \n∫ f ( x ) δ ( x -a )d x = f ( a )  \nGiven 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.  \nWe 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  \nL ( ˆ T | T ∗ ) = { | ˆ T -T ∗ | 3 T < T prior | ˆ T -T ∗ | T ≥ T prior  \nWhat 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 P ( T ) , compute the 50%, 80%, and 95% credible intervals using posterior thresholds P 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 Z 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?', '4 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 § 3 therefore requires the use of numerical methods (highlighted in the previous exercises).  \nTo 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:  \nE P [ f ( Θ )] = ∫ f ( Θ ) P ( Θ )d Θ ≈ n ∑ i =1 f ( Θ i ) P ( Θ i )∆ Θ i (15)  \nwhere  \n∆ Θ i = Θ j +1 -Θ j (16)  \nis simply the spacing between the set of j = 1 , . . . , n + 1 points on the underlying grid and  \nΘ i = Θ j +1 + Θ j 2 (17)  \nis just defined to be the mid-point between Θ j and Θ j +1 . 1 As shown in Figure 3, this approach is akin to trying to approximate the integral using a discrete set of n rectangles with heights of f ( Θ i ) P ( Θ i ) and widths of ∆ Θ i .  \nThis 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 ) P ( Θ i ) and the volume  \n∆ Θ i = d ∏ j =1 ∆Θ i,j (18)  \nwhere ∆Θ i,j is the width of the i th cuboid in the j th dimension. See Figure 3 for a visual representation of this procedure.  \nFigure 3: An illustration of how to approximate posterior integrals using a discrete grid of points. We break up the posterior into contiguous regions defined by a position Θ i (e.g., an endpoint or midpoint) with corresponding posterior density P ( Θ i ) and volume ∆ Θ i over a grid with i = 1 , . . . , n elements. Our integral can then be approximated by adding up each of these regions proportional to the posterior mass P ( Θ i ) × ∆ Θ i contained within it. In 1-D (top), these volume elements ∆ Θ i correspond to line segments while in 2-D (middle), these correspond to rectangles. This can be generalized to higher dimensions (bottom), where we instead used N-D cuboids. See § 4 for additional details.  \n<!-- image -->  \nSubstituting P ( Θ ) = ˜ P ( Θ ) / Z into the expectation value and replacing any integrals with their grid-based approximations then gives:  \nE P [ f ( Θ )] = ∫ f ( Θ ) P ( Θ )d Θ ∫ P ( Θ )d Θ = ∫ f ( Θ ) ˜ P ( Θ )d Θ ∫ ˜ P ( Θ )d Θ ≈ ∑ n i =1 f ( Θ i ) ˜ P ( Θ i )∆ Θ i ∑ n i =1 ˜ P ( Θ i )∆ Θ i (19)  \nNote the denominator is now an estimate for the evidence:  \nZ = ∫ ˜ P ( Θ )d Θ ≈ n ∑ i =1 ˜ P ( Θ i )∆ Θ j (20)  \nThis substitution of the unnormalized posterior ˜ P ( Θ ) for the posterior P ( Θ ) is a crucial part of computing expectation values in practice since we can compute ˜ P ( Θ ) = L ( Θ ) π ( Θ ) directly without knowing Z .", '4.1 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  \nn ∼ d ∏ j =1 k = k d (21)  \nThis means that even in the absolute best case where k = 2 , we have 2 d scaling.  \nThis 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 § 7.', '4.2 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 ˜ P ( Θ ) and/or f ( Θ ) ˜ P ( Θ ) 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.  \nNote that our use of the term 'weights' in the preceding paragraph is quite deliberate. Looking back at our original approximation, the form of equation (19) is quite similar to one which might be used to compute a weighted sample mean of f ( Θ ) . In that case,  \nwhere we have n observations { f 1 , . . . , f n } with corresponding weights { w 1 , . . . , w n } , the weighted mean is simply:  \nIndeed, if we define  \nˆ f mean ≡ ∑ n i =1 w i f i ∑ n i =1 w i (22)  \nf i ≡ f ( Θ i ) , w i ≡ ˜ P ( Θ i )∆ Θ i (23)  \nthen the connection between the weighted sample mean in equation (22) and the expectation value from our grid in equation (19) becomes explicit:  \nE P [ f ( Θ )] ≈ ∑ n i =1 f ( Θ i ) ˜ P ( Θ i )∆ Θ i ∑ n i =1 ˜ P ( Θ i )∆ Θ i ≡ ∑ n i =1 w i f i ∑ n i =1 w i (24)  \nThinking 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.  \nReturning to grids, this correspondence means that we can in theory come up with an estimate of the expectation value E P [ 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 -1 / 2 eff rather than ∝ n -1 / 2 .  \nWe can quantify the ideas behind the ESS as discussed above by introducing a formal definition following Kish (1965):  \nn eff ≡ ( ∑ n i =1 w i ) 2 ∑ n i =1 w 2 i (25)  \nIn line with our intuition, the best case under this definition is one where all the weights are equal ( w i = w ):  \nn best eff = ( ∑ n i =1 w i ) 2 ∑ n i =1 w 2 i = ( nw ) 2 ∑ n i =1 w 2 = n 2 w 2 nw 2 = n (26)  \nLikewise, 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):  \nn worst eff = ( ∑ n i =1 w i ) 2 ∑ n i =1 w 2 i = ( w ) 2 w 2 = 1 (27)  \nFigure 4: An example of how changing the spacing (volume elements) of the grid can dramatically affect its associated estimate of posterior integrals. On a toy 2-D posterior P ( Θ ) , simply changing the spacing of the associated 2-D 30 × 30 grid dramatically affects the effective sample size (ESS) (see § 4.2). Differences between poor spacing (left), uniform spacing (middle), and optimal spacing (right) leads to an order of magnitude difference in the ESS, as highlighted by the distribution of weights (bottom) associated with the volume elements of each grid. See § 4.2 for additional details.  \n<!-- image -->  \nFigure 5: An illustration of how grid-based estimates can be convergent (i.e. converge to a single value as the number of grid points increases) but not consistent (i.e. the value it converges to is not the correct answer). Our toy 2-D unnormalized posterior ˜ P ( Θ ) has two modes that are well-separated with a total evidence of Z = 200 . If we are not aware of the second mode, we might define a grid region that only encompasses a subset of the entire parameter space (left). While increasing the resolution of the grid within this region allows the estimated Z to converge to an single answer (left to right), this is not equal to the correct answer of Z = 200 because we have neglected the contribution of the other component (right). See § 4.3 for additional details.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nIncreasing resolution  \n<!-- image -->  \nThis former situation (with n best eff ) 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 worst eff ) would be where the entire integral is essentially contained in just one of our n N-Dcuboid regions. An illustration of this behavior is shown in Figure 4.", '4.3 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 →∞ :  \nlim n →∞ ∑ n i =1 f ( Θ i ) ˜ P ( Θ i )∆ Θ i ∑ n i =1 ˜ P ( Θ i )∆ Θ i = C (28)  \nConsistency is subsequently the idea that the value we converge to is the true value we are interested in estimating:  \nlim n →∞ ∑ n i =1 f ( Θ i ) ˜ P ( Θ i )∆ Θ i ∑ n i =1 ˜ P ( Θ i )∆ Θ i = E P [ f ( Θ )] (29)  \nIt 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  \nvalues 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 E P [ f ( Θ )] to arbitrary precision by simply increasing the resolution in ∆ Θ .  \nUnfortunately, 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 E P [ f ( Θ )] . An illustration of this issue is shown in Figure 5. This fundamental problem is not shared by Monte Carlo methods, which I will cover in § 5.", '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:  \n- 1. the evidence Z ,\n- 2. the means E P [ x ] and E P [ y ] ,\n- 3. the 68% credible intervals (or closest approximation) [ x low , x high ] and [ y low , y high ] ,\n- 4. and the effective sample size n eff .  \nHow 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': '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?  \nWhile in theory we can choose Q ( x, y ) ≈ P ( x, y ) so that n eff ≈ n , we do not know the exact shape of the posterior ahead of time. Given that ˜ P ( x, y ) may differ from our initial expectations, what does this exercise imply about general concerns applying Importance Sampling 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 § 4.2. 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?', '5.1 Connecting Grid Points and Samples': "Earlier, I outlined how we can relate estimating E P [ 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 ≡ ˜ P ( Θ 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.  \nThe 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 E P [ f ( Θ )] . Likewise, if we have the same spacing but change the relative shape of the posterior, the weight of that point when estimating E P [ f ( Θ )] should also change accordingly.  \nI 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 = ˜ P ( Θ i )∆ Θ i are uniform to some amount of desired precision. By inspection, this should happen when  \n∆ Θ i ∝ 1 ˜ P ( Θ i ) (30)  \nfor all i .  \nTaking 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  \na function of Θ . Using this, we can now define the density of points Q ( Θ ) based on the varying resolution ∆ Θ ( Θ ) of our infinitely-fine grid as a function of Θ :  \nQ ( Θ ) ∝ 1 ∆ Θ ( Θ ) (31)  \nThis result suggests that, in the continuum limit where n →∞ , the structure of our infiniteresolution grid is equivalent to a new continuous distribution Q ( Θ ) . An illustration of this concept is shown in Figure 6. Using Q ( Θ ) , we can then rewrite our original expectation value as  \nE P [ f ( Θ )] ≡ ∫ f ( Θ ) ˜ P ( Θ )d Θ ∫ ˜ P ( Θ )d Θ = ∫ f ( Θ ) ˜ P ( Θ ) Q ( Θ ) Q ( Θ )d Θ ∫ ˜ P ( Θ ) Q ( Θ ) Q ( Θ )d Θ = E Q [ f ( Θ ) ˜ P ( Θ ) / Q ( Θ ) ] E Q [ ˜ P ( Θ ) / Q ( Θ ) ] (32)  \nFor reasons that will soon become clear, I will refer to Q ( Θ ) as the proposal distribution .  \nAt 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 ˜ P ( Θ ) in terms of two expectation values with respect to the proposal distribution Q ( Θ ) . This substitution, however, actually allows us to fully realize the connection between grid points and samples.  \nEarlier, 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 Q ( Θ ) , however, this statement can become exact assuming we can explicitly generate samples from Q ( Θ ) .  \nLet's quickly review what this means. Initially, we looked at trying to estimate E P [ f ( Θ )] over a grid with n points. In the limit of infinite resolution, however, our grid becomes equivalent to some distribution Q ( Θ ) . Using Q ( Θ ) , we can then rewrite our original expression in terms of two expectations, E Q [ f ( Θ ) ˜ P ( Θ ) / Q ( Θ ) ] and E Q [ ˜ P ( Θ ) / Q ( Θ ) ] , over Q ( Θ ) instead of P ( Θ ) . This helps us because we can in theory estimate these final expressions explicitly using a series of n randomly generated samples from Q ( Θ ) . Due to the randomness inherent in this approach, this is commonly referred to as a Monte Carlo approach for estimating E P [ f ( Θ )] due to historical connections with randomness and gambling.  \nOn 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 Q ( Θ ) ∝ 1 / ∆ Θ ( Θ ) , we will get the same result (and error) for E P [ f ( Θ )] every time.  \nBy contrast, drawing n samples { Θ 1 , . . . , Θ n } from Q ( Θ ) 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 E P [ 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  \nIncreasing resolution  \n<!-- image -->  \nFigure 6: An illustration of the connection between grids and continuous density distributions. As we increase the number of grid points, our estimate of the posterior P ( Θ ) improves (top). Since the spacing between the grid points varies to maximize the effective sample size (see Figure 4 and § 4.2), the differential volume elements ∆ Θ i change depending on our location (middle). As we continue to increase the number of volume elements, the density of grid points at any particular location ρ ( Θ i ) = [∆ Θ i ] -1 behaves like a continuous function Q ( Θ ) whose distribution is similar to P ( Θ ) (bottom). This implies we should be able to use Q ( Θ ) in some way to estimate P ( Θ ) . See § 5 for additional details.  \nQ ( Θ ) . The fact that we can derive roughly equivalent estimates from these these very different approaches as we adjust n and Q ( Θ ) lies at the heart of the connection between grid points and samples.  \nThere are three primary benefits from moving from an adaptively-spaced grid to a continuous distribution Q ( Θ ) . 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 Q ( Θ ) to more closely match the posterior P ( Θ ) , giving a larger ESS at fixed n .  \nSecond, 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 Q ( Θ ) will provide sufficient coverage over all possible Θ values that our posterior P ( Θ ) could be defined over. This means that some of the theoretical issues raised in § 4.3 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.  \nFinally, the minimum number of grid points always scales exponentially with dimensionality (see § 4.1), 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 E P [ f ( Θ )] . They are therefore less susceptible to this effect (although see § 7.2).", '5.2 Importance Sampling': "As I have tried to emphasize previously, the core tenet of this article is that we do not know what P ( Θ ) looks like beforehand . This means we do not know what grid structure will provide an optimal estimate (i.e. maximum ESS) for E P [ f ( Θ ] , let alone how this should behave as Q ( Θ ) in the continuum limit. This gives us ample motivation to choose Q ( Θ ) in such a way to make generating samples from it easy and straightforward.  \nAssuming we have chosen such a Q ( Θ ) , we can subsequently generate a series of n samples from it. Assuming these samples have weights q i associated with them and defining  \nf ( Θ i ) ≡ f i , ˜ P ( Θ i ) / Q ( Θ i ) ≡ ˜ w ( Θ i ) ≡ ˜ w i (33)  \nour original expression reduces to  \nE P [ f ( Θ )] = E Q [ f ( Θ ) ˜ w ( Θ )] E Q [ ˜ w ( Θ )] ≈ ∑ n i =1 f i ˜ w i q i ∑ n i =1 ˜ w i q i (34)  \nIf we further assume that we have chosen Q ( Θ ) 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  \nE P [ f ( Θ )] ≈ n -1 ∑ n i =1 f i ˜ w i n -1 ∑ n i =1 ˜ w i (35)  \nFigure 7: A schematic illustration of Importance Sampling. First, we take a given proposal distribution Q ( Θ ) (left) and generate a set of n iid samples from it (middle left). We then weight each sample based on the corresponding 'importance' ˜ P ( Θ ) / Q ( Θ ) it has at that location (middle right). We then can use these weighted samples to approximate posterior expectations (right). See § 5.2 for additional details.  \n<!-- image -->  \nAs with the previous case using grids ( § 4), the denominator of this expression is again a direct approximation for the evidence  \nZ = ∫ ˜ P ( Θ )d Θ ≈ n -1 n ∑ i =1 ˜ w i (36)  \nThis gives a straightforward recipe for estimating our original expectation value:  \n- 1. Draw n iid samples { Θ 1 , . . . , Θ n } from Q ( Θ ) .\n- 2. Compute their corresponding weights ˜ w i = ˜ P ( Θ i ) / Q ( Θ i ) .\n- 3. Estimate E P [ f ( Θ )] by computing E Q [ ˜ w ( Θ )] and E Q [ f ( Θ ) ˜ w ( Θ )] using the weighted sample means.  \nSince 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 Figure 7.  \nWe can interpret the importance weights as ways to correct for how 'far off' our original guess Q ( Θ ) is from the truth P ( Θ ) . 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  \nthe 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.", '5.3 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.  \nAs an example, one common approach is to generate samples uniformly within some cuboid with volume V . The proposal distribution for this will then be  \nQ unif ( Θ ) = { 1 /V Θ in cuboid 0 otherwise (37)  \nThe corresponding importance weights subsequently will just be proportional to the posterior at a given position:  \n˜ w unif i = ˜ P ( Θ i ) Q unif ( Θ i ) = V ˜ P ( Θ i ) ∝ P ( Θ i ) (38)  \nAnother possible approach would be if we instead take our proposal to be our prior:  \nQ prior ( Θ ) = π ( Θ ) (39)  \nThis 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 L ( Θ ) at each position:  \nw prior i = ˜ P ( Θ i ) Q prior ( Θ i ) = L ( Θ i ) π ( Θ i ) π ( Θ i ) = L ( Θ i ) (40)  \nFinally, notice that the optimal sampling strategy is to assume that we can take our proposal to be identical to our posterior:  \nQ post ( Θ ) = P ( Θ ) (41)  \nThe corresponding weights will then just be constant and equal to the evidence Z :  \nw post i = ˜ P ( Θ i ) Q post ( Θ i ) = ZP ( Θ i ) P ( Θ i ) = Z (42)  \nAs expected, this final result guarantees the maximum possible ESS of n eff = n . Getting Q ( Θ ) to be as 'close' as possible to P ( Θ ) 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 § 6 onward: if we can somehow generate samples directly from P ( Θ ) or something close to it, then we can achieve an optimal estimate of our corresponding expectation values.", 'Importance Sampling': 'We want to use Importance Sampling to approximate various posterior integrals from this distribution. We will start by choosing our proposal distribution Q ( x, y ) to be a 2-D Gaussian with a mean of 0 and standard deviation of 1 :  \nQ ( x, y ) = N [( µ x , µ y ) = (0 , 0) , ( σ x , σ y ) = (1 , 1)]  \nUsing n = 25 iid random samples drawn from the proposal distribution, compute an estimate for:  \n- 1. the evidence Z ,\n- 2. the means E P [ x ] and E P [ y ] ,\n- 3. the 68% credible intervals (or closest approximation) [ x low , x high ] and [ y low , y high ] ,\n- 4. and the effective sample size n eff .  \nHow accurate are each of these quantities with the values we might expect? What does n eff /n tell us about how well our proposal Q ( x, y ) traces the underlying posterior P ( 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?', '6 Markov Chain Monte Carlo': "Nowthat 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 § 5.3, this means MCMC seeks to generate samples proportional to the posterior P ( Θ ) in order to arrive at an optimal estimate for our expectation value.  \nMCMCaccomplishes 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 P ( Θ i ) contained within that region. In other words, the 'density' of samples generated from MCMC  \nρ ( Θ ) ≡ m ( Θ ) n (43)  \nat position Θ integrated over δ Θ is approximately  \n∫ Θ ∈ δ Θ P ( Θ )d Θ ≈ ∫ Θ ∈ δ Θ ρ ( Θ )d Θ ≈ n -1 n ∑ j =1 1 [ Θ j ∈ δ Θ ] (44)  \nwhere 1 [ · ] 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 Figure 8.  \nWhile this will just be approximately true for any finite n , as the number of samples n → ∞ this procedure generally guarantees that ρ ( Θ ) → P ( Θ ) everywhere. 2 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 § 5:  \nZ = ∫ ˜ P ( Θ ) ρ ( Θ ) ρ ( Θ )d Θ ≡ E ρ [ ˜ P ( Θ ) /ρ ( Θ ) ] ≈ n -1 n ∑ i =1 ˜ P ( Θ i ) ρ ( Θ i ) (45)  \nFigure 8: A schematic illustration of Markov Chain Monte Carlo (MCMC). MCMC tries to create a chain of n (correlated) samples { Θ 1 →··· → Θ n } (top) such that the number of samples m in some particular volume δ gives a relative density m/n (middle) comparable to the posterior P ( Θ ) integrated over the same volume (bottom). See § 6 for additional details.  \n<!-- image -->", 'Increasing resolution': "This is just the average of the ratio between ˜ P ( Θ i ) and ρ ( Θ i ) over all n samples.  \nFinally, since our MCMC procedure gives us a series of n samples from the posterior, our expectation value simply reduces to  \nE P [ f ( Θ )] ≈ n -1 ∑ n i =1 f i ˜ w i n -1 ∑ n i =1 ˜ w i = n -1 ∑ n i =1 f i n -1 ∑ n i =1 1 = n -1 n ∑ i =1 f i (46)  \nThis is just the sample mean of the corresponding { f 1 , . . . , f n } values over our set of n samples.  \nI 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 Z . 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 § 7.1).  \nThe 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 Z . 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 § 7.1.  \nTo 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 P ( Θ ) . 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.", '6.1 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.  \nInstead 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 P ( Θ ) . These are essentially analogs to the convergence and consistency constraints discussed in § 4.3.  \nWe 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:  \nP ( Θ i +1 | Θ i ) P ( Θ i ) = P ( Θ i +1 , Θ i ) = P ( Θ i | Θ i +1 ) P ( Θ i +1 ) (47)  \nwhere 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:  \nP ( Θ i +1 | Θ i ) P ( Θ i | Θ i +1 ) = P ( Θ i +1 ) P ( Θ i ) = P ( Θ i +1 ) P ( Θ i ) (48)  \nwhere the final equality comes from the fact that the distribution we are trying to generate samples from is the posterior P ( Θ ) .  \nWe 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 Q ( Θ ' i +1 | Θ i ) similar in nature to the Q ( Θ ) used in to Importance Sampling ( § 5.2). 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:  \nP ( Θ i +1 | Θ i ) ≡ Q ( Θ i +1 | Θ i ) T ( Θ i +1 | Θ i ) (49)  \nAs with Importance Sampling, we can choose Q ( Θ ' 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:  \nT ( Θ i +1 | Θ i ) T ( Θ i | Θ i +1 ) = P ( Θ i +1 ) P ( Θ i ) Q ( Θ i | Θ i +1 ) Q ( Θ i +1 | Θ i ) (50)  \nIt is straightforward to show that the Metropolis criterion Metropolis et al. (1953)  \nT ( Θ i +1 | Θ i ) ≡ min [ 1 , P ( Θ i +1 ) P ( Θ i ) Q ( Θ i | Θ i +1 ) Q ( Θ i +1 | Θ i ) ] (51)  \nsatisfies this constraint.  \nGenerating samples following this approach can be done using the Metropolis-Hastings (MH) Algorithm (Metropolis et al., 1953; Hastings, 1970):  \n- 1. Propose a new position Θ i → Θ ' i +1 by generating a sample from the proposal distribution Q ( Θ ' i +1 | Θ i ) .\n- 2. Compute the transition probability T ( Θ ' i +1 | Θ i ) = min [ 1 , P ( Θ ' i +1 ) P ( Θ i ) Q ( Θ i | Θ ' i +1 ) Q ( Θ ' i +1 | Θ i ) ] .\n- 3. Generate a random number u i +1 from [0 , 1] .  \nFigure 9: A schematic illustration of the Metropolis-Hastings algorithm. At a given iteration i , we have generated a chain of samples { Θ 1 → ··· → Θ i } (white) up to the current position Θ i (red) whose behavior follows the underlying posterior P ( Θ ) (viridis color map). We then propose a new position Θ ' i +1 (yellow) from the proposal distribution (orange shaded region). We then compute the transition probability T ( Θ ' i +1 | Θ i ) (white) based on the posterior Q ( Θ ) and proposal Q ( Θ ' | Θ ) densities. We then generate a random number u i +1 uniformly from 0 to 1. If u i +1 ≤ T ( Θ ' i +1 | Θ i ) , we accept the move and make our next position in the chain Θ i +1 = Θ ' i +1 . If we reject the move, then Θ i +1 = Θ i . See § 6.1 for additional details.  \n<!-- image -->  \n- 4. 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 .\n- 5. Increment i = i +1 and repeat this process.  \nSee Figure 9 for a schematic illustration of this process.  \nBecause 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.  \nAn 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  \nbiasing 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 Gelman & Rubin (1992), Gelman et al. (2013), and Vehtari et al. (2019) along with references therein.", '6.2 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 .  \nThis occurs for two reasons. First, new positions Θ i drawn from Q ( Θ 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.  \nSecond, even if we set Q ( Θ ' | Θ ) = Q ( Θ ' ) 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.  \nAs mentioned in § 4.2, 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.  \nThis 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:  \nC ( t ) ≡ E i [ ( Θ i -¯ Θ ) · ( Θ i + t -¯ Θ ) ] = lim n →∞ 1 n n ∑ i =1 ( Θ i -¯ Θ ) · ( Θ i + t -¯ Θ ) (52)  \nwhere · 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.  \nUsing the auto-covariance, we can define the corresponding auto-correlation A ( t ) as  \nA ( t ) ≡ C ( t ) C (0) (53)  \nThis 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 .  \nglyph[negationslash]  \nThe overall auto-correlation time for our chain is just the auto-correlation A ( t ) summed over all non-zero lags ( t = 0 ):  \nτ ≡ ∞ ∑ t = -∞ A ( t ) -1 = 2 ∞ ∑ t =1 A ( t ) (54)  \nwhere 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 Figure 10.  \nIncorporating the auto-correlation time leads directly to a modified definition for the ESS:  \nn ' eff ≡ n eff 1 + τ (55)  \nIn practice, we cannot precisely compute τ since we do not have an infinite number of samples and do not know P ( Θ ) . 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 Brooks et al. (2011) for additional details.  \nThe 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  \nn ' eff , MCMC = n eff , MCMC 1 + τ = n 1 + τ ≤ n (56)  \nThis 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  \nn ' eff , IS = n eff , IS 1 + τ = n eff , IS = ( ∑ n i =1 ˜ w i ) 2 ∑ n i =1 ˜ w 2 i ≤ n (57)  \nwhich can be greater than n ' eff , MCMC at fixed n .  \nGiven 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 § 6.1 and § 8) and the proposal distribution Q ( Θ ) used for Importance Sampling (see § 5.3).  \nFigure 10: A schematic illustration of the auto-correlation associated with MCMC. MCMC methods generate a chain of samples { Θ 1 →··· → Θ n } (top), but these tend to be strongly correlated on small length scales (top middle). We can quantify the degree of correlation by computing the corresponding auto-correlation A ( t ) over our set of samples and all possible time lags t (bottom middle). This quantity is 1 when t = 0 and drops to 0 as t →±∞ . The overall auto-correlation time τ associated with our chain of samples is then just the integrated auto-correlation over t = 0 . See § 6.2 for additional details.  \n<!-- image -->  \nglyph[negationslash]", '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 E P [ x ] and E P [ 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 Z . 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 § 5? 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 Q ( x, y ) matches the structure of the underlying posterior P ( 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 § 5?', '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 Q ( x ' , y ' | x, y ) relative to the underlying posterior P ( x, y ) ? Given that ˜ P ( x, y ) may differ from the structure assumed when picking Q ( 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?', '7 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 ρ ( Θ ) → P ( Θ ) as n → ∞ , the primary purpose of MCMC is estimating expectation values E P [ 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.", '7.1 Approximating the Posterior': "Although algorithms such as MH ( § 6.1) are constructed to ensure the density of the chain of samples ρ ( Θ ) generated by MCMC converges to the posterior P ( Θ ) 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 P ( Θ ) ?  \nTo 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 δ Θ  \nto within some precision glyph[epsilon1] so that  \n∣ ∣ ∣ ∣ ∣ 1 n n ∑ i =1 1 [ Θ i ∈ δ Θ ] -∫ δ Θ P ( Θ )d Θ ∣ ∣ ∣ ∣ ∣ ≡ | ˆ p ( δ Θ ) -p ( δ Θ ) | < glyph[epsilon1] (58)  \nwhere p ( δ Θ ) is the total probability contained within δ Θ and ˆ p ( δ Θ ) is the fraction of the MCMCchain 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 3 posterior.  \nIn the ideal case where our samples are iid and drawn from P ( Θ ) , our samples each have a probability p ( δ Θ ) of being within δ Θ . The probability that ˆ p ( δ Θ ) = m/n then follows the binomial distribution :  \nP ( ˆ p ( δ Θ ) = m n ) = ( n m ) [ p ( δ Θ )] m [1 -p ( δ Θ )] n -m (59)  \nIn 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 ( n m ) 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 .  \nThis distribution has a mean of p ( δ Θ ) , so for any finite n we expect ˆ p ( δ Θ ) to be an unbiased estimator of p ( δ Θ ) :  \nE [ˆ p ( δ Θ ) -p ( δ Θ )] = p ( δ Θ ) -p ( δ Θ ) = 0 (60)  \nThe variance, however, depends on the sample size:  \nE [ | ˆ p ( δ Θ ) -p ( δ Θ ) | 2 ] = p ( δ Θ ) [1 -p ( δ Θ )] n (61)  \nIn 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 glyph[epsilon1] :  \nn glyph[greaterorsimilar] p ( δ Θ ) [1 -p ( δ Θ )] glyph[epsilon1] 2 / (1 + τ ) ∼ ˆ p ( δ Θ ) [1 -ˆ p ( δ Θ )] glyph[epsilon1] 2 × (1 + ˆ τ ) (62)  \nThe 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.  \nLet'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 ∝ glyph[epsilon1] -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 .  \nThe 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 P ( Θ ) and the target region δ Θ : the sample size needed to approximate the posterior to some glyph[epsilon1] near the peak of the distribution (the small region where P ( Θ ) is large) will likely be different than the sample size needed to accurately estimate the tails of the distribution (the large region where P ( Θ ) is small).  \nWhile 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 P ( Θ ) 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:  \n| ˆ p ( δ Θ 1 ) -p ( δ Θ 1 ) | < glyph[epsilon1] 1 . . . | ˆ p ( δ Θ m ) -p ( δ Θ m ) | < glyph[epsilon1] m (63)  \nSubstituting 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 :  \nn j glyph[greaterorsimilar] ˆ p ( δ Θ j ) [ 1 -ˆ p ( δ Θ j ) ] glyph[epsilon1] 2 j × (1 + ˆ τ ) (64)  \nThe total number of samples we need is then simply:  \nn glyph[greaterorsimilar] m ∑ j =1 n j (65)  \nThis approach of dividing up our posterior into sub-regions is conceptually similar to the grid-based approaches described in § 4. 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 wewouldendupwith 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.  \nThis effect implies that we should in general expect the number of samples required to ensure ρ ( Θ ) is a good approximation to P ( Θ ) for some specified accuracy glyph[epsilon1] to scale as  \nn glyph[greaterorsimilar] k d (66)  \nwhere k is a constant that depends on the accuracy requirements. This puts approximating the full posterior firmly in the 'curse of dimensionality' regime (see § 4.1). 4  \nWhile many practitioners talk about MCMC being an efficient method to 'approximate the posterior', in practice it is rarely used to approximate P ( Θ ) directly. As discussed in § 3 and shown in Figure 2, 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 glyph[lessorsimilar] 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.", '7.2 Posterior Volume': "The basic consequences outlined in § 7.1 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 E P [ 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.  \nTo start, let's consider the d -dimensional hyper-cube (the d -cube) with side length glyph[lscript] in all d dimensions. Its volume scales as  \nV ( glyph[lscript] ) = d ∏ i =1 glyph[lscript] = glyph[lscript] d (67)  \nThe differential volume element between glyph[lscript] and glyph[lscript] +d glyph[lscript] is  \nd V ( glyph[lscript] ) = ( d × glyph[lscript] d -1 ) × (d glyph[lscript] ) ∝ glyph[lscript] d -1 (68)  \nThis 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  \nglyph[lscript] 50 = 2 -1 /d glyph[lscript] (69)  \nthat divides the d -cube into two equal-sized regions with 50% of the volume contained interior to glyph[lscript] 50 and 50% of the volume exterior to glyph[lscript] 50 . In 1-D, this gives glyph[lscript] 50 /glyph[lscript] = 0 . 5 as we'd expect. In 2-D, this gives glyph[lscript] 50 /glyph[lscript] ≈ 0 . 7 . In 3-D, glyph[lscript] 50 /glyph[lscript] ≈ 0 . 8 . In 7-D, glyph[lscript] 50 /glyph[lscript] ≈ 0 . 9 . By the time we get to 15-D, we have glyph[lscript] 50 /glyph[lscript] ≈ 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.  \nIn 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  \nFigure 11: A schematic illustration of how the curse of dimensionality affects MCMC acceptance fractions via posterior volume. At a given position position Θ , the volume increases ∝ r d as a function of distance r away from that position (left). As the dimensionality increases, this implies volume becomes concentrated progressively further out, leading to larger distances between proposed positions Θ ' and the current position Θ . Most of these positions have significantly lower posterior probabilities P ( Θ ' ) compared to the current value P ( Θ ) , leading to an exponential decline in the typical acceptance fraction (and a corresponding increase in the auto-correlation time) as the dimensionality increases (right). Adjusting the size and/or shape of the proposal Q ( Θ ' | Θ ) can help to counteract this behavior. See § 7.2 for additional details.  \n<!-- image -->  \ncase, we need look no further than the transition probability used in the MH algorithm discussed in § 6.1:  \nT ( Θ i +1 | Θ i ) ≡ min [ 1 , P ( Θ i +1 ) P ( Θ i ) Q ( Θ i | Θ i +1 ) Q ( Θ i +1 | Θ i ) ]  \nThe 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 Q ( Θ ' | Θ ) . The second is dependent on the density and is related to how the posterior density changes between the two positions.  \nIn 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. 5 For instance, in the case where our proposal Q ( Θ ' | Θ ) is a cube with side-length glyph[lscript] centered on Θ , this leads to a median length-scale of glyph[lscript] 50 = 2 -1 /d glyph[lscript] , which increases rapidly from 0 . 5 glyph[lscript] to ≈ glyph[lscript] as the dimensionality increases. The same logic also applies to other proposal distributions (see § 8). This focus on positions either far away or with very similar separation length-scales as glyph[lscript] 50 → glyph[lscript] means that many choices of Q ( Θ ' | Θ ) 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 Figure 11.  \nOne of the main ways to counteract this behavior is to adjust the size/shape of the proposal Q ( Θ ' | Θ ) so that the fraction of proposed positions that are accepted remains sufficiently high. This helps to ensure the posterior density P ( Θ ) 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.", '7.3 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?  \nFrom our discussion in § 7.1, we know that the highest density of samples ρ ( Θ ) will be located where the posterior density P ( Θ ) 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  \nthat models with many parameters Θ will have the vast majority of the posterior located outside the region of highest density.  \nAconsequence 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 .  \nTo make this argument a little easier to conceptualize, let's imagine that we have a 3-parameter model Θ = ( x, y, z ) and P ( x, y, z ) is spherically symmetric. While we could imagine trying to integrate over P ( x, y, z ) directly in terms of d x d y d z , it is almost always easier to instead integrate over such a distribution in 'shells' with differential volume d V ( r ) = 4 πr 2 d 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 :  \n∫ P ( x, y, z )d x d y d z = ∫ P ( r )4 πr 2 d r ≡ ∫ P ' ( r )d r (70)  \nwhere P ' ( r ) ≡ 4 πr 2 P ( 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 P ( r ) , and implies that the the posterior should have some sort of shell-like structure (i.e. P ' ( r ) is maximized for r > 0 ).  \nAlthough 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 6  \n∫ P ( Θ )d Θ = ∫ P ( V )d V (71)  \nAs outlined in § 7.2, we generically expect the size of each volume element to go as d V ∼ r d -1 d 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  \n∫ P ( V )d V ∼ ∫ P ( r ) r d -1 d r = ∫ P ' ( r )d r (72)  \nAs before, the differential volume element of the shell associated with P ( 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 Figure 12 for an illustration of this effect based on the toy problem presented in § 8.1.  \nIncreasing dimensionality  \n<!-- image -->  \nFigure 12: A schematic illustration of how the posterior mass behaves as a function of dimensionality using a d -dimensional Gaussian. The top panel shows the posterior density P ( r ) ∝ e -r 2 / 2 (red) plotted as a function of distance r from the maximum posterior density at r = 0 as the number of dimensions d increases (left to right). As expected, this distribution remains constant. The middle panel shows the differential volume element d V ( r ) ∝ r d -1 d r (blue) of the corresponding shell at radius r . This illustrates the exponentially increasing volume contributed by shells further away from the maximum. The bottom panel shows corresponding 'posterior mass' as a function of radius P ' ( r ) ∝ r d -1 P ( r ) ∝ r d -1 e -r 2 / 2 (purple). Due to the increasing amount of volume located further away from the maximum posterior density, we see that the majority of the posterior mass (and therefore of any samples we generate with MCMC) are actually located a shell located far away from the r = 0 . See § 7.3 for additional details.  \nThis 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.  \nSecond, 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 .  \nThis allows us to now outline exactly why it is challenging to propose samples efficiently in high dimensions:  \n- 1. To make sure our acceptance fractions remain reasonable, we need to ensure our proposed positions mostly lie within this shell of posterior mass.\n- 2. However, obtaining an independent sample requires being able to (in theory) propose any position within this shell.\n- 3. 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 .", '8 Application to a Simple Toy Problem': 'I now consider a concrete, detailed example to illustrate how all the concepts discussed in § 6 and § 7 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.', '8.1 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:  \n˜ P ( Θ ) = exp [ -1 2 | Θ | 2 σ 2 ] (73)  \nwhere | Θ | 2 = ∑ d i =1 Θ 2 i is the squared magnitude of the position vector.  \nBased on the results from § 7.3, we can better understand the properties of this distribution by rewriting the posterior density in terms of the 'radius' r ≡ | Θ | = √ ∑ d i =1 Θ 2 i  \naway from the center:  \n˜ P ( r ) = exp [ -r 2 2 σ 2 ] (74)  \nThe corresponding volume contained within a given radius r is then  \nV ( r ) ∝ r d (75)  \nThe corresponding posterior mass is ˜ P ' ( r ) is then defined via  \n˜ P ( V )d V ( r ) ∝ e -r 2 / 2 σ 2 r d -1 d r ≡ ˜ P ' ( r )d r  \nNote that this is closely related to the chi-square distribution .  \nThe 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 d ˜ P ' ( r ) / d r = 0 . Solving this gives √  \nr peak = d -1 σ (76)  \nIn 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 P ( 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.  \nIn general, we expect the posterior mass to comprise a 'Gaussian shell' centered at some radius  \nr mean ≡ E P ' [ r ] = ∫ ∞ 0 r P ' ( r )d r = √ 2 Γ ( d +1 2 ) Γ ( d 2 ) σ ≈ √ dσ (77)  \nwith a standard deviation of  \n∆ r mean ≡ √ E P ' [( r -r mean ) 2 ] = σ √ √ √ √ d -2 ( Γ ( d +1 2 ) Γ ( d 2 ) ) 2 ≈ σ √ 2 (78)  \nwhere Γ( d ) is the Gamma function and the approximations are taken for large d . See Figure 12 for an illustration of this behavior.", '8.2 MCMCwith 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  \n| Θ -Θ ' | = √ √ √ √ d ∑ i =1 (Θ m,i -Θ m + t,i ) 2 (79)  \nAssuming that the lag t glyph[greatermuch] τ 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  \n∆ r sep ≡ √ E P [ | Θ m -Θ m + t | 2 ] = √ √ √ √ d ∑ i =1 E P [(Θ m,i -Θ m + t,i ) 2 ] = √ 2 dσ ≈ √ 2 r mean (80)  \nWe 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:  \nQ ( Θ i +1 | Θ i ) ∝ exp [ -1 2 | Θ i +1 -Θ i | 2 2 σ 2 ] (81)  \nWhile 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 § 7.2, we can conclude that the majority of samples proposed from this choice of Q ( Θ ' | Θ ) will probably have little overlap with the posterior.  \nIndeed, numerical simulation suggests the typical fraction of positions that will be accepted given the above proposal roughly scales as  \n〈 f acc ( d ) 〉 ≡ exp [ E P , Q [ln T ( Θ i +1 | Θ i )]] ∼ exp [ -d 4 -1 2 ] (82)  \nwhich decreases exponentially as the dimensionality increases, similar to Figure 11. Likewise, we find the auto-correlation time roughly scales as  \n〈 τ ( d ) 〉 ≡ exp [ E P , Q [ln τ ]] ∼ exp [ d 4 + 7 4 ] (83)  \nThis exponential dependence arises because the overlap between the typical Gaussian proposal Q ( Θ ' | Θ ) and the underlying posterior P ( Θ ) 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.  \nTo counteract this effect, we need to adjust the σ of our proposal distribution by some factor γ :  \nQ γ ( Θ i +1 | Θ i ) ∝ exp [ -1 2 | Θ i +1 -Θ i | 2 ( γσ ) 2 ] (84)  \nwhere 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  \n〈 f acc ( γ ( d )) 〉 ≈ C ⇒ γ ( d ) ∝ 1 √ d (85)  \nFigure 13: Numerical results showcasing the performance of a simple MH MCMC sampler with Gaussian proposals on our toy problem, a d -dimensional Gaussian with mean µ = 0 and standard deviation σ = 1 in every dimension. The top series of panels show snapshots of a random parameter from the chain as a function of dimensionality (increasing from left to right) assuming an unchanging proposal with constant scale factor γ = √ 2 (blue) and a shrinking proposal with γ = 2 . 5 / √ d designed to target a constant acceptance fraction of ∼ 25% (red). The bottom panels show the corresponding acceptance fractions (left), auto-correlation times (middle), and effective sample sizes (right) from our chains (colored points) as a function of dimensionality. The approximations from § 8.2 are shown as light colored lines. Shrinking the size of the proposal helps to keep samples within the bulk of the posterior mass, substantially reducing the auto-correlation time and increasing the effective sample size. Failing to do so leads to an exponentially decreasing fraction of good proposals and a corresponding exponential increase/decrease in the auto-correlation time/effective sample size. See § 8.2 for additional discussion.  \n<!-- image -->", 'Increasing dimensionality': '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': 'glyph[negationslash]  \nTo 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 C j i 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 C glyph[epsilon1] = C j i / 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 . 7  \nAs shown in Figure 15, 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.', '8.3 MCMCwith 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:  \n- 1. the width of the posterior in each dimension (parameter) was constant such that σ 1 = σ 2 = · · · = σ n = σ and\n- 2. the parameters were entirely uncorrelated with each other such that the correlation coefficient ρ ij = 0 between any two dimensions i and j .  \nIn 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 autocorrelation time (see § 5.3 and § 6.2) becomes one of the most difficult parts of running MCMCalgorithms in practice.  \nWhile there are schemes to perform these adjustments during an extended burn-in period (see, e.g., Brooks et al. 2011), there is significant appeal in methods that can 'autotune' 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.  \nWe explore three variations of ensemble methods here that attempt to exploit m glyph[greaterorsimilar] d ( d +1) / 2 chains running simultaneously:  \n- 1. using the ensemble of particles to condition a Gaussian proposal distribution,\n- 2. using trajectories from multiple particles along with Gaussian 'jitter', and\n- 3. using affine-invariant transformations of trajectories from multiple particles.  \nAschematic illustration of these methods is shown in Figure 14.  \nAs 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 glyph[lessorsimilar] 25 ) where a few hundred particles are often sufficient to ensure reasonable performance.", '8.3.1 Gaussian Proposal': 'The first approach is simply a modified Gaussian proposal: at any iteration i for any chain j , we propose a new position Θ j i +1 based on the current position Θ j i using a Gaussian proposal  \nQ j γ ( Θ j i +1 | Θ j i ) ∝ exp [ -1 2 ( Θ j i +1 -Θ j i ) T ( γ 2 C j i ) -1 ( Θ j i +1 -Θ j i ) ] (89)  \nwhere T is the transpose operator and  \nC j i = Cov [ { Θ 1 i , . . . , Θ j -1 i , Θ j +1 i , . . . , Θ m i } ] (90)  \nis 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.  \nIn 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 P ( Θ ) , it is straightforward to show that C j i is a reasonable approximation to the unknown covariance structure of our posterior. In addition, because we exclude j when computing C j i , this proposal is symmetric going from Θ j i → Θ j i +1 and from Θ j i +1 → Θ j i . This means that we satisfy detailed balance and do not have to incorporate any proposal-dependent factors when computing the transition probability.', 'Ensemble Gaussian': "Figure 14: A schematic illustration of the three ensemble MCMC methods described in § 8.3. The current state of the chain we are interested in updating (red) and the other chains in the ensemble (gray) are shown on the left. In the top panels (ensemble Gaussian; § 8.3.1), we compute the covariance of the other k = j chains (middle) and use a scaled version to subsequently propose a new position. In the middle panels (3-chain shift + jitter; § 8.3.2), we use two additional chains k = l = j to compute a trajectory. We then propose a new position based on this scaled trajectory plus a small amount of 'jitter'. In the bottom panels (2-chain stretch; § 8.3.3), we use only one additional chain k = j to propose a new trajectory. We then propose a random position along a scaled version of this trajectory with the proposal probability varying as a function of scale. See § 8.3 for additional details.  \n<!-- image -->  \nglyph[negationslash]  \nglyph[negationslash]  \nglyph[negationslash]", '8.3.2 Ensemble Trajectories with a Gaussian Proposal': "The approach taken in § 8.3.1 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.  \nglyph[negationslash]  \nglyph[negationslash]  \nOnesuchapproach used in the literature is Differential Evolution MCMC(DE-MCMC; Storn & Price, 1997; Ter Braak, 2006). 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' glyph[epsilon1] :  \nΘ j i +1 = Θ j i + γ × ( Θ k i -Θ l i + glyph[epsilon1] ) (91)  \nIn the case where the behavior of chains k and l are approximately independent of each other and assuming the underlying posterior distribution P ( Θ ) is Gaussian with some unknown mean µ and covariance C (and 'standard deviation' C 1 / 2 ), it is straightforward to show that the distribution of Θ k i -Θ l i will then follow  \nΘ k i -Θ l ∼ N [ 0 , (2 C ) 1 / 2 ] (92)  \nTypically, the jitter glyph[epsilon1] is chosen to also be Gaussian distributed with covariance C glyph[epsilon1] such that  \nglyph[epsilon1] ∼ N [ 0 , C 1 / 2 glyph[epsilon1] ] (93)  \nIn general, C glyph[epsilon1] is mostly used to try and avoid issues caused by finite particle sampling: since the number of unique trajectories (ignoring symmetry) is  \nn traj = ( m -1 2 ) = ( m -1)! 2!( m -3)! = ( m -1)( m -2) 2  \nif 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.  \nCombined, this implies that the proposed position has a distribution of  \nΘ j i +1 ∼ N [ Θ j i , γ × (2 C + C glyph[epsilon1] ) 1 / 2 ] (94)  \nThis shows that the 3-particle DE-MCMC procedure can generate new positions in a manner analogous to the ensemble Gaussian proposal we first discussed.", '8.3.3 Affine-Invariant Transformations of Ensemble Trajectories': "Another approach used in the literature (e.g., emcee ; Foreman-Mackey et al., 2013) is the Affine-Invariant 'stretch move' from Goodman & Weare (2010). This uses only one additional particle Θ k i rather than two:  \nΘ j i +1 = Θ k i + γ × ( Θ j i -Θ k i ) (95)  \nIn place of the jitter term glyph[epsilon1] 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 Goodman & Weare (2010), if this function is chosen such that  \ng ( γ -1 ) = γ × g ( γ ) (96)  \nthen this proposal is symmetric. Typically, g ( γ ) is chosen to be  \ng ( γ | a ) = { γ -1 / 2 a -1 ≤ γ ≤ a 0 otherwise (97)  \nwhere a = 2 is often taken as a typical value. Note that when γ = 1 , this move leaves Θ j i +1 = Θ j i unchanged.  \nCompared to DEMCMC, the stretch move appears to have one clear advantage: it doesn't have any reliance on some 'jitter' term glyph[epsilon1] 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.  \nThis lack of jitter, however, is not substantially advantageous in practice. As noted in § 8.3.2, glyph[epsilon1] 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.  \nIn 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 § 7.2, 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:  \nT ( Θ j i +1 | Θ j i , γ ) = min [ 1 , γ d -1 P ( Θ j i +1 ) P ( Θ j i ) ] (98)  \nThis 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.  \nFinally, 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 § 8.2, given the typical position of Θ j i , 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  \nFigure 15: Numerical results showcasing the performance of several ensemble MH MCMC samplers on our toy problem, a d -dimensional Gaussian with mean µ = 0 and standard deviation σ = 1 in every dimension. The top series of panels show snapshots of a random parameter from the collection of chains (with a few chains highlighted) as a function of dimensionality (increasing from left to right) assuming ensemle Gaussian proposals with γ = 2 . 5 / √ d (blue), 3-chain 'shift and jitter' proposals with γ = 1 . 7 / √ d (red), and 2-chain 'stretch' proposals with γ drawn from the distribution g ( γ | a ) with a = 2 as described in § 8.3.3 (orange). The bottom panels show the corresponding acceptance fractions (left), auto-correlation times (middle), and effective sample sizes (right) from our chains (colored points) as a function of dimensionality. Approximations based on the § 8.2 are shown as light solid colored lines, with dashed lines showing rough fits. The first two methods, which allow the size of the proposal to shrink, are able to propose samples within the bulk of the posterior mass. The last method, which is unable to do so, instead proposes exponentially fewer good positions as the dimensionality increases. See § 8.3 for additional details.  \n<!-- image -->", '8.4 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.', '9 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) .  \nIn 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  \ninference 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.  \nI 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.  \nJSS 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'ee Hlo˘zek, 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.  \nJSS acknowledges financial support from the National Science Foundation Graduate Research Fellowship Program (Grant No. 1650114) and the Harvard Data Science Initiative.", 'References': 'Asmussen, S., & Glynn, P. W. 2011, Statistics & Probability Letters, 81, 1482 , doi: https: //doi.org/10.1016/j.spl.2011.05.004  \nBlitzstein, J., & Hwang, J. 2014, Introduction to Probability, Chapman & Hall/CRC Texts  \nin Statistical Science (CRC Press/Taylor & Francis Group). https://books.google. com/books?id=ZwSlMAEACAAJ  \nBrooks, S., Gelman, A., Jones, G., & Meng, X.-L. 2011, Handbook of Markov Chain Monte Carlo (CRC press)  \nForeman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Pub. of the Astron. Sco. of the Pac., 125, 306, doi: 10.1086/670067  \nGelman, A., Carlin, J., Stern, H., et al. 2013, Bayesian Data Analysis, Third Edition, Chapman & Hall/CRC Texts in Statistical Science (Taylor & Francis). https://books. google.com/books?id=ZXL6AQAAQBAJ  \nGelman, A., & Meng, X.-L. 1998, Statist. Sci., 13, 163, doi: 10.1214/ss/1028905934  \nGelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457, doi: 10.1214/ss/ 1177011136  \nGoodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computer Science, 5, 65, doi: 10.2140/camcos.2010.5.65  \nHastings, W. 1970, Biometrika, 57, 97, doi: 10.1093/biomet/57.1.97  \nHogg, D. W., & Foreman-Mackey, D. 2018, The Astrophys. Journal Supp., 236, 11, doi: 10. 3847/1538-4365/aab76e  \nKish, L. 1965, Survey sampling, Wiley classics library (J. Wiley). https://books. google.com/books?id=xiZmAAAAIAAJ  \nMetropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, Journal of Chem. Phys., 21, 1087, doi: 10.1063/1.1699114  \nNeal, R. M. 2012, arXiv e-prints, arXiv:1206.1901. https://arxiv.org/abs/1206. 1901  \nSkilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395405  \nSkilling, J. 2006, Bayesian Anal., 1, 833, doi: 10.1214/06-BA127  \nStorn, R., & Price, K. 1997, Journal of global optimization, 11, 341  \nTer Braak, C. J. 2006, Statistics and Computing, 16, 239  \nVehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Burkner, P.-C. 2019, arXiv e-prints, arXiv:1903.08008. https://arxiv.org/abs/1903.08008'}