# Normalize likelihood for better MCMC performance?

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


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