2

I am trying to use the emcee python package to draw samples from a distribution and maximize the likelihood that my data came from the sampled parameters.

So for example, I have a parameter N and I trying to find a value for N that maximizes the posterior likelihood. (I'm actually using 3 parameters, but I'm using 1 in this example for simplicity).

I ran: sampler = emcee.EnsembleSampler(100, 3, logL, args=[new_data])

The I chose initial positions p0 for my parameters.

And then I ran:

pos, prob, state = sampler.run_mcmc(p0, 100)  # burn in
sampler.reset()
pox, prob, state = sampler.run_mcmc(pos, 100, rstate0=state)  # sample

It's mostly working, but sometimes the sampler chooses values for N that are nonphysical.

So how do I place restrictions on the range of values that are chosen by the sampler? For example, perhaps I want the sampler to stop trying N as a negative number or to stop N being greater than 100.

I understand that I can change my own likelihood function to make the nonphysical values pay a big penalty and be disfavoured - but I don't want the sampler to be allowed to even choose them in the first place.

I now think that I am supposed to build my likelihood function such that the numbers I don't want the sampler to choose (e.g. negative numbers) are penalised and given a very low likelihood.

I just want someone to confirm this is what I should be doing, in case I am missing a much more simple way to restrict which numbers are chosen in Emcee itself.

user1551817
  • 6,693
  • 22
  • 72
  • 109

1 Answers1

0

You can include a prior that the posterior function can use to determine whether to return a value that indicates the non-physical nature of the parameters, e.g. negative infinity. Example for two parameters:

def logprior(theta):
    # Extract values
    lam_pr, H0_pr = theta
    # If the values are in the domain, taken here to be uniform
    if 0.0 < lam_pr < 1.0 and 20.0 < H0_pr < 120.0:
        return 0.0
    # Otherwise
    return -numpy.inf

def logposterior(theta):
    # Calculate log prior
    lnpr = lnprior(theta)
    # If prior is unphysical
    if not numpy.isfinite(lnpr):
        return -numpy.inf

    # Else, calculate posterior as normal
SiderealTaxi
  • 1
  • 1