pythonmcmcemcee

emcee walkers burn in but then remain the same

I'm having an issue using emcee. Its a simple enough 3 parameter fit but occasionally (only has occurred in two scenarios so far despite much use) my walkers burn in just fine but then do not move (see figure below). The acceptance fraction reported is 0.

Has anyone else encountered this issue before? I have tried varying my initial conditions and number of walkers and iterations etc. This piece of code has been running well on very similar data sets. Its not a challenging parameter space and it seems unlikely that the walker would be getting "stuck".

Any ideas? I'm stumped... my walkers are lazy it seems...

enter image description here

Sample code below (and sample data file). This code + data file fail when I run it. 

`
import numpy as n
import math
import pylab as py    
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
from scipy import ndimage
import pyfits
from scipy import stats
import emcee
import corner
import scipy.optimize as op
import matplotlib.pyplot as pl
from matplotlib.ticker import MaxNLocator

def sersic(x, B,r_s,m):
      return B * n.exp(-1.0 * (1.9992*m - 0.3271) * ( (x/r_s)**(1.0/m) - 1.))

def lnprior(theta):
      B,r_s,m, lnf = theta
      if  0.0 < B < 500.0 and 0.5 < m < 10. and r_s > 0. and -10.0 < lnf < 1.0: 
        return 0.0
      return -n.inf

def lnlike(theta, x, y, yerr): #"least squares"
      B,r_s,m, lnf = theta
      model = sersic(x,B, r_s, m)
      inv_sigma2 = 1.0/(yerr**2 + model**2*n.exp(2*lnf))
      return -0.5*(n.sum((y-model)**2*inv_sigma2 - n.log(inv_sigma2)))

def lnprob(theta, x, y, yerr):#kills based on priors
    lp = lnprior(theta)
    if not n.isfinite(lp):
        return -n.inf
    return lp + lnlike(theta, x, y, yerr)


profile=open("testprofile.dat",'r')  #read in the data file    
profilelines=profile.readlines()
x=n.empty(len(profilelines))
y=n.empty(len(profilelines))
yerr=n.empty(len(profilelines))

for i,line in enumerate(profilelines):
    col=line.split()
    x[i]=col[0]
    y[i]=col[1]
    yerr[i]=col[2]

# Find the maximum likelihood value.
chi2 = lambda *args: -2 * lnlike(*args)
result = op.minimize(chi2, [50,2.0,0.5,0.5], args=(x, y, yerr))
B_ml, rs_ml,m_ml, lnf_ml = result["x"]
print("""Maximum likelihood result:
          B   = {0}
          r_s = {1}
          m   = {2}
    """.format(B_ml, rs_ml,m_ml))

# Set up the sampler.
ndim, nwalkers = 4, 4000
pos = [result["x"] + 1e-4*n.random.randn(ndim) for i in     range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
# Clear and run the production chain.
print("Running MCMC...")
Niter = 2000 #2000
sampler.run_mcmc(pos, Niter, rstate0=n.random.get_state())
print("Done.")
# Print out the mean acceptance fraction. 
af = sampler.acceptance_fraction
print "Mean acceptance fraction:", n.mean(af)

# Plot sampler chain
pl.clf()
fig, axes = pl.subplots(3, 1, sharex=True, figsize=(8, 9))
axes[0].plot(sampler.chain[:, :, 0].T, color="k", alpha=0.4)
axes[0].yaxis.set_major_locator(MaxNLocator(5))
axes[0].set_ylabel("$B$")

axes[1].plot(sampler.chain[:, :, 1].T, color="k", alpha=0.4)
axes[1].yaxis.set_major_locator(MaxNLocator(5))
axes[1].set_ylabel("$r_s$")

axes[2].plot(n.exp(sampler.chain[:, :, 2]).T, color="k", alpha=0.4)
axes[2].yaxis.set_major_locator(MaxNLocator(5))
axes[2].set_xlabel("step number")

fig.tight_layout(h_pad=0.0)
fig.savefig("line-time_test.png")

# plot MCMC fit
burnin = 100
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim-1))

B_mcmc, r_s_mcmc, m_mcmc = map(lambda v: (v[0]),
                         zip(*n.percentile(samples, [50],
                                            axis=0)))
print("""MCMC  result:
         B   = {0}
         r_s = {1}
         m   = {2}
    """.format(B_mcmc, r_s_mcmc, m_mcmc))

pl.close()

# Make the triangle plot.
burnin = 50
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim-1))

fig = corner.corner(samples, labels=["$B$", "$r_s$", "$m$"])
fig.savefig("line-triangle_test.png")
Solution

  • Here's a better result. I made the random initial samples not so close to the maximum likelihood value and run the chain for a lot more steps with fewer walkers/chains. Notice that I'm plotting the m parameter and not its exponential, as you did.

    The mean acceptance fraction is ~0.48, and it took about 1 min to run in my laptop. You can of course add more steps and get an even better result.

    enter image description here

    import numpy as n
import emcee
import corner
import scipy.optimize as op
import matplotlib.pyplot as pl
from matplotlib.ticker import MaxNLocator


def sersic(x, B, r_s, m):
    return B * n.exp(
        -1.0 * (1.9992 * m - 0.3271) * ((x / r_s)**(1.0 / m) - 1.))


def lnprior(theta):
    B, r_s, m, lnf = theta
    if 0.0 < B < 500.0 and 0.5 < m < 10. and r_s > 0. and -10.0 < lnf < 1.0:
        return 0.0
    return -n.inf


def lnlike(theta, x, y, yerr):  # "least squares"
    B, r_s, m, lnf = theta
    model = sersic(x, B, r_s, m)
    inv_sigma2 = 1.0 / (yerr**2 + model**2 * n.exp(2 * lnf))
    return -0.5 * (n.sum((y - model)**2 * inv_sigma2 - n.log(inv_sigma2)))


def lnprob(theta, x, y, yerr):  # kills based on priors
    lp = lnprior(theta)
    if not n.isfinite(lp):
        return -n.inf
    return lp + lnlike(theta, x, y, yerr)


profile = open("testprofile.dat", 'r')  # read in the data file
profilelines = profile.readlines()
x = n.empty(len(profilelines))
y = n.empty(len(profilelines))
yerr = n.empty(len(profilelines))

for i, line in enumerate(profilelines):
    col = line.split()
    x[i] = col[0]
    y[i] = col[1]
    yerr[i] = col[2]

# Find the maximum likelihood value.
chi2 = lambda *args: -2 * lnlike(*args)
result = op.minimize(chi2, [50, 2.0, 0.5, 0.5], args=(x, y, yerr))
B_ml, rs_ml, m_ml, lnf_ml = result["x"]
print("""Maximum likelihood result:
          B   = {0}
          r_s = {1}
          m   = {2}
          lnf = {3}
    """.format(B_ml, rs_ml, m_ml, lnf_ml))

# Set up the sampler.
ndim, nwalkers = 4, 10
pos = [result["x"] + 1e-1 * n.random.randn(ndim) for i in range(nwalkers)]

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
# Clear and run the production chain.
print("Running MCMC...")
Niter = 50000
sampler.run_mcmc(pos, Niter, rstate0=n.random.get_state())
print("Done.")
# Print out the mean acceptance fraction.
af = sampler.acceptance_fraction
print("Mean acceptance fraction:", n.mean(af))

# Plot sampler chain
pl.clf()
fig, axes = pl.subplots(3, 1, sharex=True, figsize=(8, 9))
axes[0].plot(sampler.chain[:, :, 0].T, color="k", alpha=0.4)
axes[0].yaxis.set_major_locator(MaxNLocator(5))
axes[0].set_ylabel("$B$")

axes[1].plot(sampler.chain[:, :, 1].T, color="k", alpha=0.4)
axes[1].yaxis.set_major_locator(MaxNLocator(5))
axes[1].set_ylabel("$r_s$")

# axes[2].plot(n.exp(sampler.chain[:, :, 2]).T, color="k", alpha=0.4)
axes[2].plot(sampler.chain[:, :, 2].T, color="k", alpha=0.4)
axes[2].yaxis.set_major_locator(MaxNLocator(5))
axes[2].set_ylabel("$m$")
axes[2].set_xlabel("step number")

fig.tight_layout(h_pad=0.0)
fig.savefig("line-time_test.png")

# plot MCMC fit
burnin = 10000
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim - 1))

B_mcmc, r_s_mcmc, m_mcmc = map(
    lambda v: (v[0]), zip(*n.percentile(samples, [50], axis=0)))
print("""MCMC  result:
         B   = {0}
         r_s = {1}
         m   = {2}
    """.format(B_mcmc, r_s_mcmc, m_mcmc))

pl.close()

# Make the triangle plot.
burnin = 50
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim - 1))

fig = corner.corner(samples, labels=["$B$", "$r_s$", "$m$"])
fig.savefig("line-triangle_test.png")