Walkers not "walking" in model fit using MCMC

Question

I am having difficulty in performing an MCMC analysis of a model. I believe it relates to the fact I have an incomplete gamma function within the model.

I trying to minimise a Gaussian log-likelihood but it appears the walkers are stuck in their well and not trying to minimize the likelihood function. This is demonstrated in the below picture, where the y axis are the parameters of the model and the x axis is the step number. The image shows how the walkers are not exploring the parameter space. I have added another image to demonstrate what a correct exploration of the parameter space looks like.

Incorrect exploration of parameter space and Correct exploration of parameter space

I have added some code below to demonstrate what I am doing, where x, y and yerr are arrays of about 4000 points. The code works for other models but only breaks on this model so it must be something intrinsic to this that the others do not have. The most obvious change to the other models is the addition of the incomplete gamma function, otherwise it has a very similar functional form to the other models.

The model I am fitting has the form:

def singinhomobremss(freq,S_norm,alpha,p,freq_peak): # Single inhomogeneous free-free emission model
    return S_norm*(p+1)*(np.power((freq/freq_peak),(2.1*(p+1)-alpha))*special.gammainc((p+1),(freq/freq_peak)**(-2.1)))*special.gamma(p+1)

Note I am using the python package emcee (I would post the link but supposedly I do not have enough reputation...). I really do not see why the walkers refuse to 'walk' for this model when they do for the other models. Any help is much appreciated but I understand it is reasonably niche area.

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
import scipy.special as special # For access to the incomplete gamma function.
import emcee
import triangle 
import inspect

def singinhomobremss(freq,S_norm,alpha,p,freq_peak): # Single inhomogeneous free-free emission model
    return S_norm*(p+1)*(np.power((freq/freq_peak),(2.1*(p+1)-alpha))*special.gammainc((p+1),(freq/freq_peak)**(-2.1)))*special.gamma(p+1)

# initial guess for a fit

p0guess = [7, 0.6, -0.5, 3.5]

# Defining log-likelihood function
def lnlike(theta,x,y,yerr):

    S_norm,alpha,p,freq_peak = theta
    model = singinhomobremss(x,S_norm,alpha,p,freq_peak)
    inv_sigma = 1.0/(yerr**2)

    return -0.5*(np.sum((y-model)**2*inv_sigma - np.log(inv_sigma))) 
    # Use the scipy.opt model to find the optimum of this likelihood function

nll = lambda *args: -lnlike(*args)
result = opt.fmin(nll,p0guess,args=(x,y,yerr),full_output='true')
S_norm_ml,alpha_ml,p_ml,freq_peak_ml = result[0]

# Defining priors
def lnprior(theta):
    S_norm,alpha,p,freq_peak = theta
    if S_norm_ml/100. < S_norm < S_norm_ml/100. and alpha_ml/100. < alpha < alpha_ml*100. and p_ml/100. < p < p_ml*100. and freq_peak_ml/100. < freq_peak < freq_peak_ml*100:
        return 0.00 
    return -np.inf      

# Combining this prior with the definition of the likelihood function, the probablity fucntion is:

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

# Now implement emcee

ndim, nwalkers, nsteps = len(inspect.getargspec(singinhomobremss)[0])-1, 200, 2000 

# Initialising the walkers in a Gaussian ball around maximum likelihood result
#pos = [result['x'] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
pos = [result[0] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args = (x,y,yerr))
sampler.run_mcmc(pos, nsteps) # This is the workhorse step.

# Now plotting the walks of the walkers with each step, for each parameter. If they have converged on a good value they should have clumped together.

fig = plt.figure(2,figsize=(10, 10))
fig.clf()
for j in range(ndim):
    ax = fig.add_subplot(ndim,1,j+1)
    ax.plot(np.array([sampler.chain[:,i,j] for i in range(nsteps)]),"k", alpha = 0.3)
    ax.set_ylabel((r'$S_{norm}$',r'$\alpha$',r'$p$',r'$\nu_{peak}$')[j], fontsize = 15)
plt.xlabel('Steps', fontsize = 15)
fig.show()

# To me it looks like the burn in period is well and truly over by 400 steps. So I will exclude those. 
print 'The burnin applied was 400. Make sure the walkers have converged after that many steps.'
samples = sampler.chain[:,400:,:].reshape((-1,ndim))

# Plotting the histograms of the fit.
trifig = triangle.corner(samples, labels = [r'$S_{norm}$',r'$\alpha$',r'$p$',r'$\nu_{peak}$'])

# Finally to get the 1 sigma final uncertainties you do
S_norm_singinhomobremss_mcmc, alpha_singinhomobremss_mcmc, p_singinhomobremss_mcmc, freq_peak_singinhomobremss_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples,[16,50,84], axis = 0))) # Uncertainites based on the 16th, 50th and 84th percentile.

plt.figure()
plt.clf()
plt.errorbar(nu, flux, flux_err, marker = '.', color = 'gray', linestyle='none', label = 'Data',alpha=0.2)  
plt.loglog(nu_peaklong,singinhomobremss(nu_peaklong, *poptsinginhomobremss), 'saddlebrown',label="Best fit inhomogeneous model from least-square")
plt.loglog(nu_peaklong,singinhomobremss(nu_peaklong,S_norm_singinhomobremss_mcmc[0], alpha_singinhomobremss_mcmc[0], p_singinhomobremss_mcmc[0], freq_peak_singinhomobremss_mcmc[0]),color = 'r', linestyle='-', label="Best fit inhomogeneous free-free model.")
plt.title('PKS 1718-649 Epoch - '+filename, fontsize = 15)
minnu = np.array(min(nu))-0.05*np.array(min(nu))
plt.legend(loc='lower center', fontsize=10) # make a legend in the best location
plt.xlabel('Frequency (GHz)', fontsize = 15)
plt.ylabel('Flux Density (Jy)', fontsize = 15)
plt.axis([min(nu)-0.1*min(nu), max(nu)+0.1*max(nu), min(flux)-0.1*min(flux), max(flux)+0.1*max(flux)])
plt.rc('xtick',labelsize=15)
plt.rc('ytick',labelsize=15)
plt.xticks([2,4,6,8,10]) #Setting grid line positions.
plt.yticks([3,4,5])
plt.grid(True)
plt.show()

Answer 1

I have solved the problem! Just writing up the solution here for future reference if anyone else in the future stumbles over a similar situation. The walkers were not walking because the lnprob was -inf. The solution is actually really remedial and has nothing really to do with my coding.

I isolated it down to lnprior because it was outputting -infs. It turns out that the parameter p is negative. Therefore, the condition:

p_ml/100. < p < p_ml*100.

never occurs. This statement should read

p_ml/100. > p > p_ml*100.

Once that edit is made the above code works. Pretty silly oversight!