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 '14 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. – Glen_b Mar 30 '14 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 '14 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. – Glen_b Mar 30 '14 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 '14 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 – DeltaIV Apr 16 '17 at 7:34
• @DeltaIV, yes. My point here was that thinning or not, the relevant time-scale is still the autocorrelation time. – Jorge Leitao Apr 16 '17 at 18:26

protected by Glen_bApr 16 '17 at 8:26

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?