# Reliability of Mode from an MCMC sample

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?

• I thik John Kruschke speaking about algorithm for mode estimation which is based on kernel density estimation. – Andrey Kolyadin Feb 1 '17 at 12:06
• Unless I'm mistaken being new to this area of statistics, JAGS outputs a set of samples from the posterior distribution rather than a probability density function so not sure kernal density estimation comes into it. Thanks for the link though. – Morgan Ball Feb 1 '17 at 12:45
• I think this is probably more to do with how you obtain a mode from a large sample of a continuous variable where there might not be more than one of any particular value and so you have to group (or smooth) the sample. – Morgan Ball Feb 1 '17 at 12:49
• You can obtain mode as value with maximum density on kernel density estimation. (At least this is what I do, and if I'm not mistaken J. Kruschke use same approach in his examples) – Andrey Kolyadin Feb 1 '17 at 12:58

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))


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.