Using emcee with Gaussian Priors

Question

I'm trying to use a gaussian prior using emcee and can't seem to fully figure it out. Basically I want to replace

def lnprior(theta):
     a, b, c = theta
     if 1.0 < a < 2.0 and 1.0 < b < 2.0 and 1.0 < c < 2.0:
        return 0.0
     return -np.inf

with something that would sample 'a' from a gaussian distribution with a mu and sigma. How would I do that? Like this?

def lnprior(theta):
     a, b, c = theta
     if 1.0 < b < 2.0 and 1.0 < c < 2.0:
        return 0.0
     if 0<a<20:
         mu=10
         sigma=1
         s=np.random.normal(mu, sigma)
         return s
     return -np.inf

That doesn't seem right though?

Answer 1

The previous answer by isinwe is the correct answer, I will only try to explain why it is correct.

Prior role

In the question there is an example of uniform priors (really similar to the example in the docs), which returns 0 if the values of theta are inside some constraints, and minus infinity otherwise.

This is correct because this prior will be used to calculate a posterior probability:

However, as this probabilities tend to be really small numbers, it is better to avoid to much rounding errors to work with their logarithm, which means:

---

Single variable prior

Therefore, a uniform prior like:

Will always have a constant if the condition is fulfilled and zero otherwise. As here only the proportionality is relevant, the normalization constant can be neglected. Thus, when the logarithm is taken, when the condition is fulfilled the logarithm will be zero, and it will be minus infinity otherwise.

---

Multiple priors

In the case of multiple priors, they are multiplied, which once the logarithm is taken becomes a sum. That is, in the case of uniform priors like the example, unless both conditions are fulfilled at the same time, the logarithm of the prior will be zero, -inf+0=-inf.

In the case of more complex combination of priors, we need to fall back to the proper interpretation of the priors, the sum. Therefore, in the case at hand, the prior must return the sum of each of the three priors logarithm, which is exactly what is done in isinwe's answer in an efficient way, that avoids evaluating the gaussian if the contribution of the uniform priors is already -inf.

As a general rule, it is a good idea to check first the uniform priors and return -inf if the condition is not fulfilled, and if the condition is fulfilled, evaluate all the other more complicated priors and return their sum (because the contribution of the uniform priors can be approximated to zero).

Answer 2

The following approach seems to work for me

def lnprior(theta):
    a, b, c = theta
    #flat priors on b, c
    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