Normalize likelihood for better MCMC performance?

Question

I'm using the emcee package to sample the distribution of a single parameter, using a uniform prior and 8 chains. In this toy example, my likelihood is defined numerically since my actual model is rather complicated to reproduce.

The original shape of this likelihood is shown below (top left plot) and consists of negative values for the entire parameter range. If I try to apply the MCMC sampler using this likelihood shape, the chains do mix (top middle plot) but they are stuck in a mode and do not explore the full parameter space (top right plot).

If instead I normalize the likelihood (left bottom plot) and apply the same exact sampling, the full parameter is explored (middle bottom plot) and the shape of the likelihood is properly recovered (right bottom plot).

I have two questions:

1. Why does this work? Shouldn't the MCMC sampler explore the posterior independently of the likelihood's range?
2. I can perform this normalization because this is a single-parameter toy example and I can explore the likelihood previously, but: how do I generalize this procedure for a $P$-parameters model?

---

Python code

import numpy as np
from scipy.interpolate import interp1d
import emcee
from scipy import stats
import matplotlib.pyplot as plt


def main():
    # Data that defines the shape of my likelihood.
    y = -np.array([
        5715.75, 5592.3, 5548.33, 5638.97, 5586.43, 5703.21, 5660.6,
        5714.96, 5637.59, 5599.72, 5631.14, 5684.31, 5586.08, 5617.43,
        5629.58, 5530.08, 5540.53, 5475.53, 5505.21, 5500.96, 5500.58,
        5474.65, 5462.45, 5443.82, 5441.77, 5463.53, 5512.18, 5395.85,
        5389.87, 5432.94, 5366.31, 5284.45, 5176.52, 5221.89, 5182.52,
        5084.92, 5084.3, 4972.78, 4968.32, 4818.19, 4789.56, 4872.02,
        4809.45, 4855.06, 4806.77, 4717.93, 4741.29, 4822.45, 4760.51,
        4698.31, 4744.1, 4797.08, 4777.43, 4785.02, 4687.61, 4820.73,
        4753.5, 4777.99, 4812.5, 4856.53, 4859.69, 4905.37, 4838.71,
        5058.49, 5053.58, 5057., 5159.58, 5155.03, 5079.21, 5228.57,
        5257.26, 5409.64, 5505.87, 5511.82, 5471.4, 5478.47, 5530.9,
        5578.88, 5705.87, 5633.66, 5740.72, 5760.05, 5801.39, 5808.52,
        5803.22, 5832.76, 5867.51, 5837.56, 5923.97, 5933.75, 5945.04,
        5932.16, 5909.68, 5951.29, 5958.6, 5958.07, 5970.75, 5931.93,
        5947.53, 5956.36])
    x = np.linspace(0., .6, 100)

    # Define the likelihood as functions.
    # Original
    lnlike_orig = interp1d(x, y)
    # Normalized
    norm_y = y + abs(min(y))
    norm_y = norm_y / norm_y.max()
    lnlike_norm = interp1d(x, norm_y)

    # MCMC sampler data.
    nwalkers, ndim, nsteps, nburn = 8, 1, 2000, 500
    p0 = [.3 + 1e-4 * np.random.randn(ndim) for i in range(nwalkers)]

    # Define output plot size.
    plt.style.use('seaborn-darkgrid')
    fig = plt.figure(figsize=(10, 10))

    for i, lnlike in enumerate([lnlike_orig, lnlike_norm]):
        # Run MCMC sampler
        sampler = emcee.EnsembleSampler(
            nwalkers, ndim, lnprob, args=[lnlike])
        sampler.run_mcmc(p0, nsteps)
        # Make plots.
        flat_samples = sampler.chain[:, nburn:, :].reshape(-1, ndim).T
        makePlot(
            i, x, nsteps - nburn, sampler.chain[:, nburn:, :], flat_samples,
            lnlike)

    # Store plot file.
    fig.tight_layout()
    plt.savefig("test.png", dpi=300)


def lnprior(x):
    if .0 < x < .6:
        return 0.0
    return -np.inf


def lnprob(x, lnlike):
    lp = lnprior(x)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(x)


def makePlot(j, x, N, samples, flat_samples, lnlike):
    plt.subplot(int("33" + str(j * 3 + 1)))
    titles = ["Original", "Normalized"]
    plt.title("{} likelihood".format(titles[j]))
    plt.plot(x, lnlike(x))
    plt.subplot(int("33" + str(j * 3 + 2)))
    for c in samples:
        plt.plot(range(N), c, lw=.5)
    plt.subplot(int("33" + str(j * 3 + 3)))
    for c in samples:
        hist, _ = np.histogram(c, bins=25)
        n, edges, _ = plt.hist(
            c, histtype='step', bins=25,
            weights=np.ones_like(c) / hist.max())

    # KDE for flat samples.
    x_kde = np.linspace(edges[0], edges[-1], 100)
    kernel_cl = stats.gaussian_kde(flat_samples)  # , bw_method=bw)
    kde = np.reshape(kernel_cl(x_kde).T, x_kde.shape)
    plt.plot(x_kde, kde / max(kde), color='k', lw=1.5)


if __name__ == '__main__':
    main()

Answer 1

In the top plot, what you call likelihood is presumably the log-likelihood. As such, the output you get is not surprising: the values around 0.297 are 800 log-likelihood points above values around 0.1 for example; we would not expect the chain to visit states with such low likelihood. However, the likelihood looks multimodal, so different initializations might lead to finding different modes. If you run the chain for long enough, it will find several modes, but with a basic algorithm this will take a long time.

What you have called the normalization step is therefore not really a normalization: you are multiplying the log-likelihood by a constant, hence taking the likelihood to a certain power. This will have the effect of flattening the function; notice the scales are vastly different on the y-axis in the left-hand plots. Since the likelihood is less peaked, you explore the whole space.

This is called tempering or annealing. The power to which you put the likelihood $L$ is the inverse temperature $1/T$: you are thus exploring $L^{1/T}$. If $T\to \infty$, $L^{1/T}$ becomes flat and is easy to explore. As you decrease $T$, $L^{1/T}$ becomes more peaked at its mode(s).

There are several schemes to handle the temperature: you can decrease it gradually, which is the basic idea behind simulated annealing; you can have chains run in parallel at different temperatures and interact, as in parallel tempering... Darren Wilkinson has written a nice introduction to parallel tempering.

For a multi-dimensional model, you can still take a power of the likelihood to explore it better: the ideas are exactly the same.