1
$\begingroup$

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

enter image description here

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()
$\endgroup$
4
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ You are absolutely correct, it is the log-likelihood (actually it is the pseudo-likelihood in my previous question). Two questions if you don't mind: 1. could you expand on the association with simulated annealing?. 2. Is there a general approach for this flattening in multi-dimensional models? I find that the performance of the MCMC is greatly affected by it, in the sense that if I flatten it too much I get a flat parameter distribution, and if I don't flatten it enough the MCMC gets stuck in a tiny region forever. $\endgroup$ – Gabriel Sep 7 '18 at 15:15
  • 2
    $\begingroup$ Actually the right term is tempering, which is related to simulated annealing. I have updated the answer; hopefully the blog post by Darren Wilkinson that I linked to, and the references therein, will be helpful. $\endgroup$ – Robin Ryder Sep 7 '18 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.