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!
Are you seeing both nans and infs while sampling? Did you set test points explicitly using the
testval argument when you tried
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:
The output of model.logp(model.test_point) is (-inf) and the output of
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
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]
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()
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: Yp.tag.test_value Out: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
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 = 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)
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
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
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
@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
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?