CODEWITHSUNDEEP
pythoncurve-fittingemcee
Sebastiano1991
Sebastiano1991

Reputation: 897

emcee samples initial position

I'm trying to fit a inverse gamma distribution to a histogram. I can easily pick out a set of good fit parameters by playing around with its two parameters (c.f. orange line in Fig1).

However, when I try to estimate the posteriors using emcee I simply sample the initial starting points of the chain.

import numpy as np 
import corner 
import emcee 
from scipy.special import gamma as gamma_function
import matplotlib.pyplot as plt

def _inverse_gamma_distribution(x0, alpha, beta):
    return beta**alpha/gamma_function(alpha)* (1/x0)**(alpha + 1)* np.exp(-beta/x0)


flux =          np.array([0.,            0.,         0.08753462, 0.48757902, 0.59385076, 0.32146012, 0.20280991, 0.06542532, 0.01888836, 0.00369042, 0.00133481,  0.,            0.,         0.10504154, 0.44777665])
bin_center =    np.array([0.06898463,    0.12137053, 0.21353749, 0.37569469, 0.66099163, 1.16293883, 2.04605725, 3.59980263, 6.33343911, 11.1429583, 19.60475495, 0.06898463,    0.12137053, 0.21353749, 0.37569469])
error =         np.array([0.,            0.,         0.03914667, 0.06965415, 0.0579539,  0.03214601, 0.01924986, 0.00824282, 0.00333902, 0.0011127,  0.0005045,   0.,            0.,         0.04288303, 0.0667506 ])


def lnprior(theta):
    alpha, beta         = theta
    if 1.0 < alpha < 2.0 and 1.0 < beta < 2.0:
        return 0.0
    return -np.inf

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

def lnlike(theta, x, y, yerr):
    alpha, beta         = theta
    model               = np.array(_inverse_gamma_distribution(x, alpha, beta))
    return -0.5*np.sum((y - model)**2/yerr**2)

p0                      = [1.5, 1.5]

ndim, nwalkers          = 2, 100
pos                     = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)]

sampler                 = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(bin_center, flux, error))
sampler.run_mcmc(pos, 1000)
samples                 = sampler.chain[:, 500:, :].reshape((-1, ndim))

print("mean posterior: ", samples.T[0].mean(), samples.T[1].mean()) 
print("std posterior: ",  samples.T[0].std(),  samples.T[1].std())

fig                     = corner.corner(samples)

fig, ax                 = plt.subplots()
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip')
ax.set_ylim([0.0001, 10])
x                       = np.linspace(0.1, 20, 1000)

ax.errorbar(bin_center, flux, yerr=error, fmt=".")
ax.plot(x, _inverse_gamma_distribution(x, *p0))
for alpha, beta in samples[np.random.randint(len(samples), size=100)]:
    ax.plot(x, _inverse_gamma_distribution(x, alpha, beta), "-", color="grey", alpha=0.1)
plt.show()

The posterior mean and standard deviation correspond to the values I set in the pos = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)] (i.e. [p0 +- 0.5]). Changing the number of steps and or the number of burned samples doesn't help.

orange initial position walkers stuck at position

If I increase the number of walkers, it just samples the input position better - they never seem to go anywhere.

More walkers

There is this post link, but I don't think it applies here - the inverse gamma distribution is arguably a bit annoying but the curve fitting problem should be superbly defined..?

-------------------EDIT--------------------

I know understand the issue: removing the data points where the flux is zero causes the walkers to start moving and I get a result. I find this behavior a bit weird, so to update the question: why do two single outlier prevent emcee from moving?

Upvotes: 1

Views: 1558

Answers (1)

Sebastiano1991
Sebastiano1991

Reputation: 897

The reason the walkers never started to move lies in the error array:

error = np.array([0., 0., 0.03914667, 0.06965415, ... ])

The first two values of error are 0.. Thus the likelihood function -0.5*np.sum((y - model)**2/yerr**2) is -np.inf for any value of (alpha, beta) and the walkers never start moving. Removing these values or setting the errors to any nonzero value let's walker loose and the fit quickly converges as it should!

Upvotes: 3

Related Questions