Can MCMC iterations after burn in be used for density estimation?

Question

After burn-in, can we directly use the MCMC iterations for density estimation, such as by plotting a histogram, or kernel density estimation? My concern is that the MCMC iterations are not necessarily independent, although they are at most identically distributed.

What if we further apply thinning to the MCMC iterations? My concern is that the MCMC iterations are at most uncorrelated, and not yet independent.

The ground I learned for using an empirical distribution function as an estimation of the true distribution function is based on Glivenko–Cantelli theorem, where the empirical distribution function is calculated based on an iid sample. I seemed to see some grounds (asymptotic results?) for using histograms, or kernel density estimates as density estimations, but I can't recall them.

Answer 1

You can - and people do - estimate densities from MCMC sampling.

One thing to keep in mind is that while histograms and KDEs are convenient, at least in simple cases (such as Gibbs sampling), much more efficient estimates of density may be available.

If we consider Gibbs sampling in particular, the conditional density you're sampling from can be used in place of the sample value itself in producing an averaged estimate of the density. The result tends to be quite smooth.

The approach is discussed in

Gelfand and Smith (1990), "Sampling-Based Approaches to Calculating Marginal Densities"
Journal of the American Statistical Association, Vol. 85, No. 410, pp. 398-409

(though Geyer cautions that if the sampler dependence is high enough it doesn't always reduce the variance and gives conditions for it to do so)

This approach is also discussed, for example, in Robert, C. P. and Casella, G. (1999) Monte Carlo Statistical Methods.

You don't need independence, you're actually computing an average. If you want to compute a standard error of a density estimate (or a cdf), then you have to account for the dependence.

The same notion applies to other expectations, of course, and so it can be used to improve estimates of many other kinds of average.

Answer 2

Resume

You can directly use the MCMC iterations for anything because the average value of your observable will asymptotically approach the true value (because you are after the burn-in).

However, bear in mind that the variance of this average is influenced by the correlations between samples. This means that if the samples are correlated, as is common in MCMC, storing every measurement will not bring any real advantage.

In theory, you should measure after N steps, where N is of the order of the autocorrelation time of the observable you are measuring.

Detailed explanation

Let's define some notation to formally answer your question. Let $x_t$ be the state of your MCMC simulation at time $t$, assumed much higher than the burn-in time. Let $f$ be the observable you want to measure.

For example, $x_t \in \mathbb{R}$, and $f=f_a(x)$: "1 if $x\in[a,a+\Delta]$, 0 else". $x_t$ is naturally being drawn from a distribution $P(x)$, which you do using MCMC.

In any sampling, you will always need to compute an average of an observable $f$, which you do using an estimator:

$$F = \frac{1}{N}\sum_{i=1}^N f(x_i)$$

We see that the average value of this estimator $\langle F\rangle$ (in respect to $P(x)$) is

$$\langle F \rangle = \frac{1}{N}\sum_{i=1}^N \langle f(x_i)\rangle = \langle f(x)\rangle$$

which is what you want to obtain.

The main concern is that when you compute the variance of this estimator, $\langle F^2 \rangle - \langle F \rangle^2$, you will obtain terms of the form

$$\sum_{i=1}^N\sum_{j=1}^N \langle f(x_i)f(x_j)\rangle$$

which do not cancel out if $x_t$ are correlated samples. Moreover, because you can write $j=i+\Delta$, you can write the above double sum as sum of the autocorrelation function of $f$, $R(\Delta)$

So, to recap:

• If computationally it does not cost anything to store every measure, you can do it, but bear in mind that the variance can not be computed using the usual formula.

• If it is computationally expensive to measure at each step of your MCMC, you have to find a way to estimate the cumulative of the autocorrelation time $\tau$ and perform measurements only every $\tau$. In this case, the measurements are independent and thus you can use the usual formula of the variance.