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

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.

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.

• Thanks! Do you mean that, because marginal distributions are expectations wrt the joint distribution, it doesn't matter to use correlated MCMC iterations to estimate marginal distributions? What if using the correlated iterations to estimate the joint distribution? Still ok?
– Tim
Mar 30, 2014 at 18:09
• No that's what I mean. I mean that the estimators we're dealing with are averages of things, and are being used to estimate population quantities that may be construed in turn as expectations of those things. Yes, you can use dependent draws to estimate a joint distribution in the same sense. Mar 30, 2014 at 22:43
• Why can we use the correlated iterations to estimate the joint distribution? I think no, because joint distribution is not expectation of something. Note that in Glivenko–Cantelli theorem, the empirical cdf is calcuated on iid sample.
– Tim
Mar 30, 2014 at 22:45
• For the density, you might consider something like the sample estimate described here for example (and might be considered as the limit of a histogram with increasingly-narrow bins); it's an average, and I believe its expectation is the density. In respect of the cdf you might like to consider whether you can do something with the empirical cdf to make it in the form of an average. Both ideas would seem to work with samples from a joint distribution. Mar 30, 2014 at 23:04

## 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.

• This doesn't answer the particular question, which concerned using the samples from a Markov chain to construct a valid density estimator of the posterior. The point that the standard error of our estimate of a linear functional is higher than a naive estimate based on independence is well appreciated, but OP still wouldn't know on the basis of this answer if it is a good idea to construct a density estimator using (say) kernel-smoothing (which even under iid sampling wouldn't converge at a $\sqrt n$ rate).
– guy
Mar 30, 2014 at 6:51
• Thinning is just a waste of useful data. It doesn't reduce the variance of the estimate. See the comments to this question: stats.stackexchange.com/a/258529/58675 Apr 16, 2017 at 7:34
• @DeltaIV, yes. My point here was that thinning or not, the relevant time-scale is still the autocorrelation time. Apr 16, 2017 at 18:26