# Define constraint on the priors

Hello. My first intention of writing this question was to ask about defining constraint on the priors. To explain that I made a very simple model. But I got very strange result. My model simply is `y=sqrt(a-b)` and my data is defined by number 3 plus some noise. My priors are uniform distribution from 1 to 5 for both a and b. Then I expect error, because for example a=3 and b=4 makes the square root argument negative. Surprisingly I do not get any error!!! More surprisingly the results also are reasonable!!! Can someone explain me why I do not get any error, when I must have error running this code?

``````import pymc3 as pm
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
np.random.seed(5)
#%%our data
data = 3 + np.random.normal(0,1,5)
#%% defining bayesian model
with pm.Model() as model:
a = pm.Uniform('a',lower=1,upper=5)
b = pm.Uniform('b',lower=1,upper=5)
sd_lik = pm.HalfNormal('sd_lik', sd=2.0)#prior
y=pm.math.sqrt(a-b)
obs = pm.Normal('obs', mu=y, sd=sd_lik, observed=data)
step1 = pm.Metropolis([a,b])
step2 = pm.NUTS(sd_lik)
trace = pm.sample(draws=5000,tune=1000,step=[step1,step2],njobs=1,start={'a':4.5,'b':2})
#%%
ppc = pm.sample_ppc(trace, model=model, samples=2000)
plt.hist(ppc['obs'],30)
#%%
sns.set(style="white", palette="deep", color_codes=True,font_scale=4)
pm.traceplot(trace, [a,b],figsize=(34,25));
pm.plots.autocorrplot(trace,figsize=(34,25));
pm.plot_posterior(trace, [sd_lik,a,b],figsize=(34,25))
#%%
a_pos=trace.get_values(a)
b_pos=trace.get_values(b)
g = sns.jointplot(a_pos, b_pos, kind="kde", size=16, space=0)
g.ax_joint.set_xlabel("a")
g.ax_joint.set_ylabel("b")``````

You still get valid samples, as proposal like `a=3, b=4` leads to nan or inf logp and rejected. See for example this blog post: https://darrenjw.wordpress.com/2012/06/04/metropolis-hastings-mcmc-when-the-proposal-and-target-have-differing-support/, and this blog post: https://umbertopicchini.wordpress.com/2017/12/18/tips-for-coding-a-metropolis-hastings-sampler/
The MCMC trace would, however, suffer a bit from inefficiency.

Thanks! I see, very interesting. However I have more complicated model that I must define a constraint on priors. Imagine something similar like this example, but if I get a=2 and b=3 my model will generate error. Then I cannot let these values enter to the function, for example here “sqrt”. I remember from pymc2 I need to use “potential”. I did that, but it still does not work for that model, what I did is like adding this line to above mentioned example. Wondering if I did anything wrong?

#%% defining bayesian model
with pm.Model() as model:
a = pm.Uniform(‘a’,lower=1,upper=5)
b = pm.Uniform(‘b’,lower=1,upper=5)
sd_lik = pm.HalfNormal(‘sd_lik’, sd=2.0)#prior
y=pm.math.sqrt(a-b)
const1 = pm.Potential(‘const’, pm.math.switch(pm.math.gt(a-b, 0),0,-np.inf))#a>b
obs = pm.Normal(‘obs’, mu=y, sd=sd_lik, observed=data)
step1 = pm.Metropolis([a,b])
step2 = pm.NUTS(sd_lik)
trace = pm.sample(draws=5000,tune=1000,step=[step1,step2],njobs=1,start={‘a’:2.5,‘b’:2})

Your constraint seems correct to me. I can sample locally without any problem - not sure what do you mean by it does not work.
You might find this related discussion helpful: Sampling uniformly in a triangular support

1 Like

Great! Thank you for confirming it. I will check my code one more time, and then if I see the error again, I will update here. I really appreciate your help.

So this is my model, and deterministic variable “D” cannot be negative. I first copy the model and then report D which is negative.

#%% defining bayesian model
with pm.Model() as model:
CGK1 = pm.Uniform(‘First component GK’,lower=0.97911392*.5e10-0.203312385e10,upper=0.979113925e10-0.20331238*.5e10)#prior14478996000
CGK2 = pm.Uniform(‘Second component GK’,lower=-0.203312385e10-0.979113925e10,upper=-0.20331238*.5e10-0.97911392*.5e10)#prior=11721092000 Km=CGK10.97911392+CGK2(-0.20331238)
Km=pm.Deterministic(‘Bulk modulus’,CGK1*(0.97911392)+CGK2*(-0.20331238))
Gm=pm.Deterministic(‘Shear modulus’,CGK1*(-0.20331238)+CGK2*(-0.97911392))
CDW1 = pm.Uniform(‘First component DW’,lower=.5e-10/(1.5192534852562095e-10)0.70710678+.015/0.0737650507479504640.70710678,upper=20e-10/(1.5192534852562095e-10)0.70710678+0.99/0.0737650507479504640.70710678)
CDW2 = pm.Uniform(‘Second component DW’,lower=-20e-10/(1.5192534852562095e-10)0.70710678+.015/0.0737650507479504640.70710678,upper=-0.5e-10/(1.5192534852562095e-10)0.70710678+0.99/0.0737650507479504640.70710678)
D=pm.Deterministic(‘D’,(CDW1-CDW2)1.5192534852562095e-10/(20.70710678))
W=pm.Deterministic(‘W’,(CDW1+CDW2)0.073765050747950464/(20.70710678))
sd_lik = pm.HalfNormal(‘sd_lik’, sd=2.0)#prior
const1 = pm.Potential(‘const’, pm.math.switch(pm.math.gt(D, 0),0,-np.inf))#D>0
acc = beamacc(Gm,Km,D,W)
obs = pm.Normal(‘obs’, mu=acc, sd=sd_lik, observed=data)
step1 = pm.Metropolis([CGK1,CGK2,CDW1,CDW2])
step2 = pm.NUTS(sd_lik)
db = pm.backends.Text(“savedchain”)
trace = pm.sample(draws=1100,tune=1000,trace=db,step=[step1,step2],njobs=1,start={‘First component GK’:0.9791139214478996000-0.2033123811721092000,‘Second component GK’:-0.2033123814478996000-0.9791139211721092000,‘First component DW’:4.6412e-10/(1.5192534852562095e-10)0.70710678+0.42/0.0737650507479504640.70710678,‘Second component DW’:-4.6412e-10/(1.5192534852562095e-10)0.70710678+0.42/0.0737650507479504640.70710678,‘sd_lik’:0.07})

and the reason that error happened is this input value in my model for D: -1.1760875837940131e-11
apparently potential do not work. However I got an idea from you previous message that maybe I should define NAN output for my model in these cases.
Need to add that `beamacc(Gm,Km,D,W)` is an external model, and It was checked and it works within pymc.

As a conclusion, I think pymc3 let the model takes the forbidden values, but then assign the zero probability to its results. I want to forbid the model to take these value from beginning. Based on your comment, I added an “if” statement to my model to give me “nan” output when its input is not what I want, and so far the chain is running.
Thanks!

Glad that you found a solution. I am still a bit surprise that in your more complex code it doesnt work.

1 Like

The code now works. With this trick that my core model (the model for mu) gives NAN as output when inputs are forbidden values. I think based on your comment, pymc will ignore them. However wondering if you confirm my previous comment? the way that potential enters to the model? Am I correct?
Thanks!

In theory it should always rejected if a<b, as the potential makes the logp goes -inf.
In the simple model it seems to work as `np.testing.assert_array_less(trace['b'], trace['a'])` always pass.

I think it will be rejected, because of the same reason you just mentioned, however the value executed by the model and then they would be rejected. In my model I wrote a line of code that if not good values entered, see a warning, and I see that sometimes forbidden values are used by my model.
So I conclude these values enters to the model, then they would be rejected. I do not like it though.