Restricting the value of a Deterministic variable

Questions
1

Is it possible to restrict the value of a Deterministic variable? For example, my model looks like this:

with pm.Model() as myModel:
    sig = pm.HalfNormal('sig', sd=2)
    kappa = pm.Normal('kappa', mu=1, sd=1)
    rhoCp = pm.Normal('rhoCp', mu=4, sd=1)
    alpha = pm.Determinstic('alpha', kappa/rhoCp)
    SE = ComplicatedModel(alpha, data)    # Squared error between model output and data
    pm.Potential('res', -0.5*SE/sig**2)

If alpha is greater than or equal to 0.5, the ComplicatedModel() is unstable and returns NaNs. So is there a way for me to restrict alpha<0.5 without causing the sampler too much grief? I suppose I could use theano ifelse or theano maximum like the following:

from theano import tensor as tt
from theano.ifelse import ifelse
with pm.Model() as myModel:
    sig = pm.HalfNormal('sig', sd=2)
    kappa = pm.Normal('kappa', mu=1, sd=1)
    rhoCp = pm.Normal('rhoCp', mu=4, sd=1)
    a = kappa/rhoCp
    alpha = pm.Determinstic('alpha', ifelse(tt.lt(a, 0.5), a, 0.5))
    # OR
    alpha = pm.Determinstic('alpha', tt.maximum(kappa/alpha, 0.5))

However, is the above likely to mess with the nice geometry of the parameter space and make sampling using NUTS difficult?

2

You can add a potential to your model: pm.Potential('restrict_term', tt.switch(tt.lt(a, 0.5), 0., -np.inf))
This means that when a >= .5, the model will produce -inf log_prob, and the sample will be rejected.

2 Likes
3

Works great. Thanks!

4

On further inspection, it looks like some instances where a>=0.5 are still getting through to the SE = ComplicatedModel(alpha, data) line. Iā€™m wondering at what point is the sample meant to be rejected? For example:

from theano import tensor as tt
with pm.Model() as myModel:
    sig = pm.HalfNormal('sig', sd=2)
    kappa = pm.Normal('kappa', mu=1, sd=1)
    rhoCp = pm.Normal('rhoCp', mu=4, sd=1)
    alpha = kappa/rhoCp
    pm.Potential('restrict_term', tt.switch(tt.lt(alpha, 0.5), 0., -np.inf))
    # Is the sample meant to be rejected at the above line preventing the 
    # non-conforming alpha value from being passed to the line below?
    SE = ComplicatedModel(alpha, data) 
    pm.Potential('res', -0.5*SE/sig**2)

I have included an alpha_printed = tt.printing.Print('alpha')(alpha) line inside the ComplicatedModel(alpha, data) function and found that an alpha value great than 0.5 is indeed finding its way into the function.

5

There is a post about using random variables as bounds for prior distributions (see here) which is also helpful. Using the same method, early experiments indicate that the following works:

from theano import tensor as tt
alphaMax = tt.as_tensor_variable(0.5)
with pm.Model() as myModel:
    sig = pm.HalfNormal('sig', sd=2)
    kappa = pm.Normal('kappa', mu=1, sd=1)
    rhoCp_ = pm.Normal('rhoCp_', mu=4, sd=1)
    rhoCp = pm.Determinstic('rhoCp', tt.maximum(rhoCp, kappa/alphaMax))
    alpha = pm.Determinstic('alpha', kappa/rhoCp)
    SE = ComplicatedModel(alpha, data)    # Squared error between model output and data
    pm.Potential('res', -0.5*SE/sig**2)
6

I find that quite surprising! Sampling with NUTS usually accept-reject all RVs at the same time (ie it is not sequential).
I would evaluate at the point where you got a >= .5, and check the logp:

testpoint = model.test_point
# find the slice where you got a>=.5, and replace the value below
testpoint['a'] = ...
testpoint['sig'] = ...
...
# FYI, there might be better way to do this 
print(model.logp(testpoint))