3

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?

spnlola
  • 31
  • 1
  • 2

2 Answers2

7

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
isinwe
  • 71
  • 4
5

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:

posterior prod

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:

log posterior

Single variable prior

Therefore, a uniform prior like:

prior example

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).

OriolAbril
  • 7,315
  • 4
  • 29
  • 40
  • 1
    you mean for two parameters(a and b)which we want to be gaussian we should write this one?`np.log(1.0/(np.sqrt(2*np.pi)*sigma))-0.5*(a-mu)**2/sigma**2 + np.log(1.0/(np.sqrt(2*np.pi)*sigma1))-0.5*(b-mu1)**2/sigma1**2` – Ma Y Mar 05 '19 at 22:00
  • 1
    Exactly. In addition, as it is only the proportionality dependance that is needed, you can use only `-0.5*(a-mu)**2/sigma**2-0.5*(b-mu1)**2/sigma1**2` and omit the normalization constants. Both expressions will yield the same result. – OriolAbril Mar 06 '19 at 11:43
  • I have a problem in this regards. I am getting stuck by the way – Ma Y Apr 21 '19 at 21:08
  • What is the problem? I would be happy to clarify anything that is not clearly explained. – OriolAbril Apr 22 '19 at 13:49
  • Thank you. There are several problems, but the important one is: When I choose the mean value and the erro for a parameter. for example `mu_d=0.623` and `sig_d=0.001` and the initial point for `d` equal to `0.527` the final results is exactly what I considered as mean value `0.623` for example `0.6229` – Ma Y Apr 24 '19 at 19:50
  • I mean the contour has Gaussian behavior but the results of the parameter is very close to the mean value I considered in Gaussian equation. – Ma Y Apr 24 '19 at 19:52
  • As explained, the posterior is the likelihood times the prior. In your case, the prior imposes a really tight constraint. Therefore, unless your likelihood is of the same order of magnitude as the prior, you may approximate the posterior with the prior. From your results it looks like this is what is happening. This however, has nothing to do with your implementation nor with emcee. In fact, from this result it seems your implementation of the gaussian prior is correct. – OriolAbril Apr 24 '19 at 22:30