Reliability of Mode from an MCMC sample

Question

In his book Doing Bayesian Data Analysis, John Kruschke states that in using JAGS from R

...the estimate of the mode from an MCMC sample can be rather unstable because the estimate is based on a smoothing algorithm that can be sensitive to random bumps and ripples in the MCMC sample. (Doing Bayesian Data Analysis, page 205, section 8.2.5.1)

While I have a grasp on Metropolis algorithm and exact forms like Gibbs sampling I'm not familiar with the smoothing algorithm that is alluded too and why it would mean the estimate of the mode from the MCMC sample is unstable. Is anyone able to give an intuitive insight into what the smoothing algorithm is doing and why it makes the estimate of the mode unstable?

Answer 1

I don't have the book at hand so I'm not sure what smoothing method Kruschke uses, but for intuition consider this plot of 100 samples from a standard normal, along with Gaussian kernel density estimates using various bandwidths from 0.1 to 1.0. (Briefly, Gaussian KDEs are a sort of smoothed histogram: They estimate density by adding a Gaussian for each data point, with mean at the observed value.)

You can see that even once smoothing creates a unimodal distribution, the mode is generally below the known value of 0.

More, here's a plot of the estimated mode (y-axis) by kernel bandwidth used to estimate density, using the same sample. Hopefully this lends some intuition to how the estimate varies with smoothing parameters.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  1 09:35:51 2017

@author: seaneaster
"""

import numpy as np
from matplotlib import pylab as plt
from sklearn.neighbors import KernelDensity

REAL_MODE = 0
np.random.seed(123)

def estimate_mode(X, bandwidth = 0.75):
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    return u[np.argmax(log_density)]

X = np.random.normal(REAL_MODE, size = 100)[:, np.newaxis] # keeping to standard normal

bandwidths = np.linspace(0.1, 1., num = 8)

plt.figure(0)
plt.hist(X, bins = 100, normed = True, alpha = 0.25)

for bandwidth in bandwidths:
    kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
    u = np.linspace(-3,3,num=1000)[:, np.newaxis]
    log_density = kde.score_samples(u)
    plt.plot(u, np.exp(log_density))

bandwidths = np.linspace(0.1, 3., num = 100)
modes = [estimate_mode(X, bandwidth) for bandwidth in bandwidths]
plt.figure(1)
plt.plot(bandwidths, np.array(modes))

Answer 2

Sean Easter provided a nice answer; here's how it's actually done by the R scripts that come with Kruschke's book. The plotPost() function is defined in the R script named DBDA2E-utilities.R. It displays the estimated mode. Inside the function definition, there are these two lines:

mcmcDensity = density(paramSampleVec)
mo = mcmcDensity$x[which.max(mcmcDensity$y)]

The density() function comes with the base stats package of R, and implements a kernel density filter of the sort Sean Easter described. It has optional arguments for the bandwidth of the smoothing kernel, and for the type of kernel to use. It defaults to a Gaussian kernel, and has some internal magic for finding a nice bandwidth. The density() function returns an object with a component named y that has the smoothed densities at various values x. The second line of code, above, just finds the x value where y is maximum.