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.
Thank you for your comments, it was a great discussion.

Oh I see what you mean. Yes it will be enter into the model as the values are valid proposal. But it will get rejected as the logp evaluated on these invalid points would be inf.
Usually it is fine, it might compromise the efficiency a bit but in general, the sample would most likely stay in the valid region of the posterior space.

1 Like

I see, yes it affects the efficiency because especially my model is very computationally expensive.
Thanks!