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