Question concerning MCMC

[ar-valdez]

Dear @prappleizer ,

Thanks for sharing the explanations linking the concept and the code to perform MCMC evaluations. I have slightly changed your code:

https://github.com/prappleizer/prappleizer.github.io/blob/master/Tutorials/MCMC/MCMC_Tutorial_Solution.ipynb

to my own needs. Therein couple situations appears: 1) The lnprior(theta) routine assumed Uniform distributions for the entries of theta, how can I edit to other prior distributions? Is it relevant this? The priori can be any distribution and the posteriori will always have the same distribution???

2) The routine lnprob(theta, x, y, yerr) returned to my case NaN values, that is because lnlike(theta, x, y, yerr) returned Zero.... How do you suggest to avoid this situation??? The initial guess for theta should be "far" from the real value....

Thats all, thanks for the help

[prappleizer]

Hi, sorry for the late reply -- no one has used this repo directly to ask a question before.

To (roughly) answer your questions: yes, I had a uniform prior on values of parameters in theta in that example, but you can definitely implement whatever prior distributions you want! The posterior should still be roughly gaussian (ideally) for most problems, remember you're trying to sample to get an estimate of the spread around a value, even if your prior has a different shape and is more or less constraining.

Loosely, to implement non-uniform priors in emcee, all you have to do is transform your parameters to the desired distribution within your lnprior function. For example, to adopt a gaussian prior on some parameter a while keeping b and c uniform, you could do something like

def lnprior(theta):
    a, b, c = theta
    if not 1.0 < b < 2.0 and 1.0 < c < 2.0:
        return -np.inf
    #gaussian prior on a
    mu = 10
    sigma = 1
    return np.log(1.0/(np.sqrt(2*np.pi)*sigma))-0.5*(a-mu)**2/sigma**2

The parameters b and c don't need to return any particular value (I usually return 0) because that value is a constant and our lnlike is proportional either way (a constant offset doesn't change our sampling results). But a has different prior likelihoods now (in the context of ln(post) = ln(prior) + ln(like) ) depending on its selected value, which is now captured by the value returned by lnprior.

As to the second question, I'll start with the final part --- no, the initial guess for theta should not be far from the max-likelihood location in parameter space. emcee is a sampler, so while it is useful for fitting, that's only true when you already have reasonable estimates for the parameters. Most fitting routines using emcee will first run your optimizer of choice (least squares, powell minimization, chi-squared, etc) to find the optimum parameters. emcee then effectively explores your parameter space around that point (spreading out considerably) to give you an estimate of, e.g., the parameter uncertainties. Now, you can be somewhat off, because the walkers should move towards the highest likelihood regions, but it will take a very long time if you are very far from the max likelihood region (and walkers might even get stuck trying to get there).

Also, lnlike returning zero shouldn't trigger a NaN -- the function lnprob just returns lnlike + lnprior (but with a check first to make sure lnprior is not -inf), and returns -inf if that's the case. Probably you are getting NaNs because you are attempting to take the log of a negative number, which you'll have to look into in your likelihood function. Note there is more than one function you can write down to parametrize your likelihood.

Hope that helps (or, if you're long gone, someone else who stumbles on thi).

[ar-valdez]

Thank you very much @prappleizer . This was really useful and helpful