Posterior distribution and MCMC [duplicate]

Question

I have read something like 6 articles on Markov Chain Monte carlo methods, there are a couple of basic points I can't seem to wrap my head around.

1. How can you "draw samples from the posterior distribution" without first knowing the properties of said distribution?

2. Again, how can you determine which parameter estimate "fits your data better" without first knowing your posterior distribution?

3. If you already know the properties of your posterior distribution (as is indicated by 1) and 2)), then what's the point of using this method in the first place?

This just seems like circular reasoning to me.

Answer 1

If this was not a clear conflict of interest, I would suggest you invest more time on the topic of MCMC algorithm and read a whole book rather than a few (6?) articles that can only provide a partial perspective.

How can you "draw samples from the posterior distribution" without first knowing the properties of said distribution?

MCMC is based on the assumption that the product$$\pi(\theta)f(x^\text{obs}|\theta)$$can be numerically computed (hence is known) for a given $\theta$, where $x^\text{obs}$ denotes the observation, $\pi(\cdot)$ the prior, and $f(x^\text{obs}|\theta)$ the likelihood. This does not imply an in-depth knowledge about this function of $\theta$. Still, from a mathematical perspective the posterior density is completely and entirely determined by $$\pi(\theta|x^\text{obs})=\dfrac{\pi(\theta)f(x^\text{obs}|\theta)}{\int_ \Theta \pi(\theta)f(x^\text{obs}|\theta)\,\text{d}\theta}\tag{1}$$Thus, it is not particularly surprising that simulation methods can be found using solely the input of the product $$\pi(\theta)\times f(x^\text{obs}|\theta)$$ The amazing feature of Monte Carlo methods is that some methods like Markov chain Monte Carlo (MCMC) algorithms do not formally require anything further than this computation of the product, when compared with accept-reject algorithms for instance, which calls for an upper bound. A related software like Stan operates on this input and still delivers high end performances with tools like NUTS and HMC, including numerical differentiation.

A side comment written later in the light of some of the other answers is that the normalising constant$$\mathfrak{Z}=\int_ \Theta \pi(\theta)f(x^\text{obs}|\theta)\,\text{d}\theta$$is not particularly useful for conducting Bayesian inference in that, were I to "know" its exact numerical value in addition to the function in the numerator of (1), $\mathfrak{Z}=3.17232\,10^{-23}$ say, I would not have made any progress towards finding Bayes estimates or credible regions. (The only exception when this constant matters is in conducting Bayesian model comparison.)

When teaching about MCMC algorithms, my analogy is that in a videogame we have a complete map (the posterior) and a moving player that can only illuminate a portion of the map at once. Visualising the entire map and spotting the highest regions is possible with enough attempts (and a perfect remembrance of things past!). A local and primitive knowledge of the posterior density (up to a constant) is therefore sufficient to learn about the distribution.

Again, how can you determine which parameter estimate "fits your data better" without first knowing your posterior distribution?

Again, the distribution is known in a mathematical or numerical sense. The Bayes parameter estimates provided by MCMC, if needed, are based on the same principle as most simulation methods, the law of large numbers. More generally, Monte Carlo based (Bayesian) inference replaces the exact posterior distribution with an empirical version. Hence, once more, a numerical approach to the posterior, one value at a time, is sufficient to build a convergent representation of the associated estimator. The only restriction is the available computing time, i.e., the number of terms one can call in the law of large numbers approximation.

If you already know the properties of your posterior distribution (as is indicated by 1) and 2)), then what's the point of using this method in the first place?

It is the very paradox of (1) that this is a perfectly well-defined mathematical object such that most integrals related with (1) including its denominator may be out of reach from analytical and numerical methods. Exploiting the stochastic nature of the object by simulation methods (Monte Carlo integration) is a natural and manageable alternative that has proven immensely helpful.

Connected X validated questions:

• Confusion related to MCMC technique
• What are Monte Carlo simulations?
• Is Markov chain based sampling the “best” for Monte Carlo sampling? Are there alternative schemes available?
• MCMC; Can we be sure that we have a ''pure'' and ''large enough'' sample from the posterior? How can it work if we are not?
• How would you explain Markov Chain Monte Carlo (MCMC) to a layperson?
• How to do MC integration from Gibbs sampling of posterior?

Answer 2

How can you "draw samples from the posterior distribution" without first knowing the properties of said distribution?

In Bayesian analysis we usually know that the posterior distribution is proportional to some known function (the likelihood multiplied by the prior) but we don't know the constant of integration that would give us the actual posterior density:

$$\pi( \theta | \mathbb{x} ) = \frac{\overbrace{L_\mathbb{x}(\theta) \pi(\theta)}^{\text{Known}}}{\underbrace{\int L_\mathbb{x}(\theta) \pi(\theta) d\theta}_{\text{Unknown}}} \overset{\theta}{\propto} \overbrace{L_\mathbb{x}(\theta) \pi(\theta)}^{\text{Known}}.$$

So we actually do know one major property of the distribution; that it is proportional to a particular known function. Now, in the context of MCMC analysis, a Markov chain takes in a starting value $\theta_{(0)}$ and produces a series of values $\theta_{(1)}, \theta_{(2)}, \theta_{(3)}, ...$ for this parameter.

The Markov chain has a stationary distribution which is the distribution that preserves itself if you run it through the chain. Under certain broad assumptions (e.g., the chain is irreducible, aperiodic), the stationary distribution will also be the limiting distribution of the Markov chain, so that regardless of how you choose the starting value, this will be the distribution that the outputs converge towards as you run the chain longer and longer. It turns out that it is possible to design a Markov chain with a stationary distribution equal to the posterior distribution, even though we don't know exactly what that distribution is. That is, it is possible to design a Markov chain that has $\pi( \theta | \mathbb{x} )$ as its stationary limiting distribution, even if all we know is that $\pi( \theta | \mathbb{x} ) \propto L_\mathbb{x}(\theta) \pi(\theta)$. There are various ways to design this kind of Markov chain, and these various designs constitute available MCMC algorithms for generating values from the posterior distribution.

Once we have designed an MCMC method like this, we know that we can feed in any arbitrary starting value $\theta_{(0)}$ and the distribution of the outputs will converge to the posterior distribution (since this is the stationary limiting distribution of the chain). So we can draw (non-independent) samples from the posterior distribution by starting with an arbitrary starting value, feeding it into the MCMC algorithm, waiting for the chain to converge close to its stationary distribution, and then taking the subsequent outputs as our draws.

This usually involves generating $\theta_{(1)}, \theta_{(2)}, \theta_{(3)}, ..., \theta_{(M)}$ for some large value of $M$, and discarding $B < M$ "burn-in" iterations to allow the convergence to occur, leaving us with draws $\theta_{(B+1)}, \theta_{(B+2)}, \theta_{(B+3)}, ..., \theta_{(M)} \sim \pi( \theta | \mathbb{x} )$ (approximately).

If you already know the properties of your posterior distribution ... then what's the point of using this method in the first place?

Use of the MCMC simulation allows us to go from a state where we know that the posterior distribution is proportional to some given function (the likelihood multiplied by the prior) to actually simulating from this distribution. From these simulations we can estimate the constant of integration for the posterior distribution, and then we have a good estimate of the actual distribution. We can also use these simulations to estimate other aspects of the posterior distribution, such as its moments.

Now, bear in mind that MCMC is not the only way we can do this. Another method would be to use some other method of numerical integration to try to find the constant-of-integration for the posterior distribution. MCMC goes directly to simulation of the values, rather than attempting to estimate the constant-of-integration, so it is a popular method.

Answer 3

Your confusion is understandable. Surely, if you already know $p(\theta|X)$, why would you need to draw samples of $\theta$ under this distribution? The answer is usually that the distribution is multivariate, and you want to marginalize over some dimensions of $\theta$ but not others. So for instance, $\theta$ might be a vector of 10 parameters, and you're interested in the marginal distribution $p(\theta_1|X)=\int p(\theta|X)d\theta_{2:10}$. The integrals required to do this marginalization are often very hard to compute exactly. They may be analytically intractable, and (deterministic) numerical integration is often cumbersome in high dimensions.

This is where MCMC can help. As long as you know $p(\theta|X)$ up to a constant of multiplication, you can generate samples of $\theta$ that follow this distribution. Then, given a sufficient number of such samples, you can simply look at the distribution of sampled values of $\theta_1$ (e.g. by making a histogram), and those samples will approximate the desired marginal distribution. Compared to numerical integration methods, MCMC is more efficient because it spends more time exploring parts of the distribution where more of the probability mass is concentrated. Also, many MCMC algorithms (such as the classic Metropolis Hastings algorithm) only require that you know the target distribution up to a constant of proportionality, which is helpful if you don't know the normalization constant required to make the distribution proper (which is very often the case, because to compute that constant itself often requires computing a multivariate integral just as complex as the one you're interested in).

Edit: it occurred to me that this perhaps doesn't fully answer your first question. The answer to this is that MCMC only requires that you can calculate the posterior probability (density) of a certain parameter value (up to a constant of proportionality). So all you need is a function where, if you put a parameter value in, it gives you its probability under the target distribution (or a value proportional to that probability). That is the sense in which the target distribution must be 'known'. But you don't need to know anything else about it. You can be blissfully ignorant about the mean & covariance of the distribution, or about the little squiggles and bumps that it has here or there, or any number of other things (although some of those things can be helpful to know in order to make MCMC run more smoothly).

Answer 4

Just one example to address part (1).

Sometimes you can evaluate the posterior up to a partition function only.

For example, you know that $p(x)= \frac{1}{z}f(x)$, but $z$ is unknown.

The metropolis hasting algorithm:

-Initialize $x_0$

-Choose some distribution $q$

Repeat:

-Sample $y$ from $q(x_{i-1})$

-Accept $y$ if $p(y)$ is large (essentially) via an "acceptance rule"

-if accepted set $x_i=y$

But at each step we don't know $p(y)$, we only know $f(y)$ because $z$ is unknown. However, The acceptance rule can be written (essentially) as a ratio of $p(x_{i-1})$ and $p(y)$ so $z$ cancels.

The final output of the sampling then provides $p(x), z $ included, but you never had to compute (or know) $z$.