 # Restricting the value of a Deterministic variable

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?

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

Works great. Thanks!

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.

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

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