1

I have a model with 6 parameters with uniform priors:

parameter1 = pm.Uniform('parameter1',0.01,1)
parameter2 = pm.Uniform('parameter2',0,2)
parameter3 = pm.DiscreteUniform('parameter3',1,50)
parameter4 = pm.Uniform('parameter4',0,1.75)
parameter5 = pm.Uniform('parameter5', 0.005, 0.25)
parameter6 = pm.Uniform('parameter6', 0.005, 0.15)

I have a custom likelihood function that returns log likelihood value:

@pm.potential
def log_l(experiment=experiment,parameter1=parameter1,parameter2=parameter2,parameter3=parameter3,parameter4=parameter4,parameter5=parameter5,parameter6=parameter6):

    if parameter5<parameter4:
        return -np.inf

    parameters=[parameter1, parameter2, parameter3]

    log_l=calculate_probability(parameters, t_m, tol, n_integration, parameter4, parameter5, parameter6, experiment.decon_all[freq,:,:])

    return log_l

Where calculate_probability is a my function that returns log likelihood for this model given parameter values and observed data. For some reason when MCMC samples:

model = pm.MCMC([parameter1,parameter2,parameter3,parameter4,parameter5,parameter6,log_l])
model.sample(100)

and the programme satisfies the if condition (parameter5<parameter4) I get this error:

pymc.Node.ZeroProbability: Potential log_l forbids its parents' current values

I was wondering if anyone knows what I may be doing wrong?

Ale
  • 133
  • 7
  • The question here is: why are you using `pymc` instead of `pymc3`? – David Aug 17 '18 at 01:35
  • I tried using pymc3 first, but during sampling I need all parameter values to be floats or integers so that the calculate_probability function works, and I didn't know how to convert theano objects into floats. Are you suggesting that this type of model would be easier to handle in pymc3? – Ale Aug 17 '18 at 06:02
  • Yes, it should be easier in `pymc3`. Maybe this site can help you https://discourse.pymc.io/ with this: " and I didn't know how to convert theano objects into floats" – David Aug 17 '18 at 06:28

1 Answers1

1

Per the docs:

Potentials have one important attribute, logp, the log of their current probability or probability density value given the values of their parents. The only other additional attribute of interest is parents, a dictionary containing the potential’s parents.

So the parents of your log_l are the arguments of log_l:

In [11]: list(log_l.parents.keys())
Out[11]: ['experiment', 'parameter2', 'parameter3', 'parameter1', 'parameter4', 
          'parameter6', 'parameter5']

Per this answer (my emphasis):

When a random variable is defined as a function of another random variable, PyMC checks that no value of the parent distribution leads to an impossible value for the child distribution.

Making log_l return -np.inf implies certain values of the parent variables are impossible. Hence, PyMC raises a ZeroProbability exception.

So instead of constraining the model by using 

if parameter5<parameter4:
    return -np.inf

you could define parameter5 with

parameter4 = pm.Uniform('parameter4', 0, 1.75)
parameter5 = pm.Uniform('parameter5', parameter4, 0.25)

to ensure that parameter5 > parameter4.

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
  • Great, thanks, this should work! I kept playing around with the code and put in a `if 4<5: return -np.inf` into the log_l function and it actually gives the same error. I was wondering if you would know why pymc is unhappy with the function returning `-np.inf`? – Ale Aug 17 '18 at 06:38
  • @Ale: returning `-np.inf` implies that certain values of the parent variables (namely `parameter4` and `parameter5`) are impossible (i.e., have zero probability). PyMC wants all values of the variables (defined by the domain of their distributions) to be *possible* even if improbable. If you want values to be impossible, they must be built into the definition of the distribution (e.g. `pm.Uniform('parameter5', 0.005, parameter4)`), not after the fact. – unutbu Aug 17 '18 at 11:32
  • Now I understand. Thank you for your help! – Ale Aug 17 '18 at 11:43