Hi all,
If you believe that the variables are well set and your data is reliable, what can you do if you are getting nans or infs from model.logp (model.test point) (not allowing you to use NUTS), and some variables getting the chain stuck when using Metropolis or HamiltonianMC?
Thank you very much!
Get nan or inf from model.logp (model.test point) is an attestation of incorrectly configured model?
Are you seeing both nans and infs while sampling? Did you set test points explicitly using the testval argument when you tried model.test_point?
pm.Normal("x", mu=0.0, sd=1.0, testval=10)
It’s hard to say without looking at your code, can you paste a snippet?
Hi @bwengals, here it goes the model:
https://pastebin.com/WBRsVj2q
The output of model.logp(model.test_point) is (-inf) and the output of
model.bijection.rmap(model.dlogp_array(model.dict_to_array(model.test_point))) is:
https://pastebin.com/QXM3DYMg
Thank you!
I’m not sure if this is it, but it looks like @junpenglao recently (~a week ago) made some changes to master which involve a fix for test points being set on transformed variables. See for instance: https://github.com/pymc-devs/pymc3/pull/2328
Could you try downloading master and then running your model?
I don’t think this is related to the transformed test points.
@marcio Could you also provide the data (or some simulation data)? It would be easier to diagnose the problem.
Are you sure you are pasting the model.bijection.rmap(model.dlogp_array(model.dict_to_array(model.test_point))) in https://pastebin.com/QXM3DYMg2? To me it looks like the model.test_point.
1 Like
Yes, I am pasting the output of model.bijection.rmap
loss_shifted40 = [155, 267, 105, 223, 183, 189, 262, 178, 272, 187]
fem_cases = [2, 2, 2, 3, 1, 0, 1, 0, 0, 0]
mz_freq = [0.05, 0.047, 0.045, 0.063, 0.06, 0.062, 0.053, 0.057, 0.061, 0.058]
what about N and pp_mean?
pp_mean: 0.0595759774518
N: 3392425
A few thoughts:
1, In your model, the observed part you are doing:
with model:
...
mz = pm.Uniform("mz", lower=0, upper=1, observed=mz_freq)
M = pm.Binomial("M", n=Yp, p=mz, observed=loss_shifte
However, the node mz is not doing anything, since there is no prior on mz, in effect, this is equivalent to mz = mz_freq
2, Since we got -inf from model.logp (model.logp(model.test_point)), we want to identify which RV node is causing the problem. @aseyboldt suggest doing model.bijection.rmap(model.dlogp_array(model.dict_to_array(model.test_point))) to get the logp of each RV but actually I am not sure that’s correct. I am doing below instead:
# for RV in model.basic_RVs:
# print(RV.name, RV.logp(model.test_point))
model.check_test_point()
which gives:
lower_beta_interval__ -1.3862943611198906
upper_beta_interval__ -1.3862943611198906
beta_interval__ -13.862943611198906
beta_gamma_interval__ -1.3862943611198906
alpha_gamma_interval__ -1.3862943611198906
gamma_log__ 0.14556084518264534
I_interval__ -13.862943611198906
S_interval__ -13.862943611198908
po_interval__ -1.3862943611198906
Y -18.29351833496805
Yp -3.574558647108001
Yo -26.885429018279982
M -inf
So the problem is the node M = pm.Binomial("M", n=Yp, p=mz_freq, observed=loss_shifted40), and looking closer at the test_value of the input parameter:
In[1]: Yp.tag.test_value
Out[1]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
The n of the Binomial is smaller than the observed, which causing the -inf in the log-likelihood. An easy temporally fix is to provide a valid test value to Yp:
Yp = pm.Poisson("Yp", mu=mu_Yp, shape=len(fem_cases), testval=loss_shifted40+1000)
which seems to resolve the problem.
So below is the full code I ran, hope you find it helpful.
import pymc3 as pm
import numpy as np
loss_shifted40 = np.array([155, 267, 105, 223, 183, 189, 262, 178, 272, 187])
fem_cases = np.array([2, 2, 2, 3, 1, 0, 1, 0, 0, 0])
mz_freq = np.array([0.05, 0.047, 0.045, 0.063, 0.06, 0.062, 0.053, 0.057, 0.061, 0.058])
pp_mean = 0.0595759774518
N = 3392425
with pm.Model() as model:
pp = pp_mean
lower_beta = pm.Uniform("lower_beta", lower=0, upper=1)
upper_beta = pm.Uniform("upper_beta", lower=lower_beta, upper=lower_beta+3)
beta = pm.Uniform("beta", lower=lower_beta, upper=upper_beta, shape=len(fem_cases))
beta_gamma = pm.Uniform("beta_gamma", lower=0, upper=120)
alpha_gamma = pm.Uniform("alpha_gamma", lower=(2/3)*beta_gamma/7, upper=(4/3)*beta_gamma/7)
gamma = pm.Gamma("gamma", alpha=alpha_gamma, beta=beta_gamma)
I = pm.Uniform("I", lower=0, upper=fem_cases.sum(), shape=len(fem_cases))
R = pm.Deterministic("R", (gamma*I))
S = pm.Uniform("S", lower=0, upper=N-(I+R), shape=len(fem_cases))
inc = pm.Deterministic("inc", beta*I*(S)/(S+I+R)-R)
po = pm.Uniform("po", lower=0, upper=1)
Y = pm.Poisson('Y', mu=inc, shape=len(fem_cases))
mu_Yp = pm.Deterministic("mu_Yp", var=Y*pp)
Yp = pm.Poisson("Yp", mu=mu_Yp, shape=len(fem_cases), testval=loss_shifted40+1000)
Yo = pm.Binomial('Yo', n=Y, p=po, observed=fem_cases)
M = pm.Binomial("M", n=Yp, p=mz_freq, observed=loss_shifted40)
trace = pm.sample(1e4)
2 Likes
with model:
trace = pm.sample(1e4, njobs=4)
pm.traceplot(trace[7000:]);
Hmm S seems a bit problematic… The next step I would suggest to take is reparameterize the model and try to avoid using bounded (ie. Uniform) variable for the prior.
One more question, it seems obvious, but just checking: does Yp.tag.test_value outputs the start value of this variable, right?
Yes, it can be defined by providing testval in the variable declares (e.g., pm.Normal(..., testval=...), or by the start karg in the calling of sampling (e.g., pm.sample(..., start=start,...)) where start is a dict with the same format as model.test_point.
1 Like
for RV in model.basic_RVs:
print(RV.name, RV.logp(model.test_point))
Using the procedure you taught me, I got this:
lower_beta_interval__ -1.3862943611198906
upper_beta_interval__ -1.3862943611198906
beta_interval__ -40.20253647247681
beta_gamma_interval__ -4.795859992452562
alpha_gamma_interval__ -1.3862943611198904
gamma_log__ -2.300013578603738
I_interval__ -40.20253647247681
S_interval__ -40.20253647247685
po_interval__ -2.407945608651872
Y -160.51534229164463
Yp -119.61861672892917
Yo -4612.920761723422
M -inf
and
Yp.tag.test_value:
array([608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608,
608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608, 608,
608, 608, 608])
Instead of using mz_freq, I used mz_freq.mean(), what solved the problem, even without setting a testval. Thank you very much my friend!
You are welcome! I am glad you find it helpful 
2 Likes
@junpenglao, still in this model, I would like to store the sample, but got ValueError: This backend does not support sampler stats. Why would this happen? Thank you in advance!
What backend are you specified in pm.sample?
db = pm.backends.Text(‘trace’)
trace = pm.sample(5e4, njobs=4, step=[pm.Metropolis(vars=[Y, Yp]), pm.NUTS()], tune=1000, progressbar=True, trace=db)
I see, well I have never work with backend before (and we are thinking of removing it actually), if you want to save it in a file I think there are better ways to do so. @colcarroll?