pythonpython-2.7statisticsemcee

Define a custom prior for each parameter in emcee


I have a function with three parameters a,b and c and I want to define different priors for each of these parameters. I am using the emcee package.

I started with the simple uniform (non-informative) prior:

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

I would like to have for each parameter a different prior. For instance for a I would like to have a Normal(mu,sigma) prior, while for b an uniform and for c a Jeffreys prior (1/c). Up to now I come out with the following:

def lnprior(theta):
    a, b, c = theta

    mu = 0.5 # mean of the Normal prior
    sigma = 0.1 # standard deviation of the Normal prior

if not (1.0 < b < 2.0): # the bound on the uniform
    return -np.inf
if c < 0.0:             # the bound on the Jeffreys
    return -np.inf
return .... # total log-prior to be determined

As far as I understood in log-scale I have to add together all the probabilities to define the total one (the return value of lnprior). So let's start with the Normal on a:

log_Pr(a) = np.log( 1.0 / (np.sqrt(2*np.pi)*sigma) ) - 0.5*(a - mu)**2/sigma**2;

then the prior on c:

log_Pr(c) = -log(c).

Thus the total log-prior should be: Pr(a)+Pr(c). My question, is this approach correct?

Thanks


Solution

  • Try the following one:

    def lnprior(theta):
        a, b, c = theta
        #flat priors on b, c
        if not 1.0 < b < 2.0 and c > 0:
            return -np.inf
        #gaussian prior on a and c
        mu = 0.5
        sigma = 0.1
        ### your prior is gaussian * (1/c), take natural log is the following:
        return np.log(1.0/(np.sqrt(2*np.pi)*sigma))-0.5*(a-mu)**2/sigma**2 - np.log(c)