ValueError: Probability function returned NaN when using emcee package

Question

I am using the emcee package to determine the optimal parameters of a measured dataset that should follow a Poisson distribution. The code I use is

def lnL_Poisson(theta,x,y,yerr):

    logA,beta = theta
    A = 10**logA
    model = the Poisson likelihood
    return np.sum(model)

def lnprior(theta):

    logA,beta = theta
    if -5 < logA < 0 and -2 < beta < 4:
        return 0.0
    return -np.inf


def lnprob_Poisson(theta, x, y, yerr):

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

However, when running this code it returns

ValueError                                Traceback (most recent call last)
<ipython-input-81-460e20ecdf72> in <module>
      1 sampler_2 = emcee.EnsembleSampler(nwalkers, ndim, lnprob_Poisson, args=(x, y_obs, dy))
----> 2 tmp = sampler_2.run_mcmc(pos, 500) #Run the sampler 500 times
      3 samples_2 = sampler_2.chain[:, 50:, :].reshape((-1, 2))
      4 fig = corner.corner(samples_2, labels=[r"$\log(A)$", r"$\beta$"],quantiles=[0.16, 0.5, 0.84], show_titles=True,label_kwargs=dict(fontsize=15))

~\Anaconda3\lib\site-packages\emcee\ensemble.py in run_mcmc(self, initial_state, nsteps, **kwargs)
    382 
    383         results = None
--> 384         for results in self.sample(initial_state, iterations=nsteps, **kwargs):
    385             pass
    386 

~\Anaconda3\lib\site-packages\emcee\ensemble.py in sample(self, initial_state, log_prob0, rstate0, blobs0, iterations, tune, skip_initial_state_check, thin_by, thin, store, progress)
    283             state.blobs = blobs0
    284         if state.log_prob is None:
--> 285             state.log_prob, state.blobs = self.compute_log_prob(state.coords)
    286         if np.shape(state.log_prob) != (self.nwalkers,):
    287             raise ValueError("incompatible input dimensions")

~\Anaconda3\lib\site-packages\emcee\ensemble.py in compute_log_prob(self, coords)
    454         # Check for log_prob returning NaN.
    455         if np.any(np.isnan(log_prob)):
--> 456             raise ValueError("Probability function returned NaN")
    457 
    458         return log_prob, blob

ValueError: Probability function returned NaN

The code did work when using a Gaussian log liklihood. I'm guessing it has something to do with a probability being 0 somewhere and then dividing by this value. However, I have no clue how to solve this. Anybody got experience with this?

Answer 1

Check that the x array and y array's dtype. I faced the same problem, then discovered that my x array is float32 while the y array is float64. After changing x to float64, problem solved.