Sampler (NUTS) proposes the same value

Questions
1

Hi! I have the following problem and couldn’t find any way to debug it further. When I run the code:

with pm.Model() as model:
    B1 = pm.Lognormal('B1', -2, 2)
    B1_print = tt.printing.Print('B1')(B1)
    jm.B1 = B1_print
    image_model = jm.image()
    image_obs = pm.Normal('image_obs', mu=image_model, sd=1e-6, observed=image_data)
    
    trace = pm.sample(chains=1, njobs=1, step=pm.NUTS())

where jm is a largeish class so I don’t include it here (but it works correctly when I e.g. set jm.B1 = 0.1 and call jm.image.eval()) it alternates between two outputs for B1:

B1 __str__ = 0.15222451156570388
B1 __str__ = nan
B1 __str__ = 0.15222451156570388
B1 __str__ = nan
B1 __str__ = 0.15222451156570388
B1 __str__ = nan

the number changes between code runs, but remains the same during a single run. Any ideas or suggestions what to try?

2

Seems to me the sampler is rejecting all proposal. You should not call the step directly, instead do trace = pm.sample() which use the auto initialization with multiple chain.

Also, are you sure the input to B1 is valid?

3

I also used plain pm.sample without any parameters and got the same result.

And what do you mean by B1 validity? As I say it works when I set jm.B1 = 0.1 and call jm.image().eval().

4

Then what likely happen is that the starting value is at some pathological region, which the sampler could not propose a new statue that get out from the region.

Try with other starting value (in this case using init='advi+adapt_diag' might helps).

5

Well, that seems unlikely, but will try. The B1 used to generate image_data is exactly in the middle of its prior (currently I test the inference on model data).

6

If B1 is a neural network, you can easily have wrong prior that makes sampling impossible because the model is multimodal and ill defined.

7

Nope, nothing of this kind :slight_smile: its an analytical physical model and I expect a unimodal posterior.

8

Hmm, it would be easier to check if you can share the full model/data in a notebook.

9

It turned out that the issue has something to do with theano broadcasting, and I managed to fix it (even through not completely understanding).

10

aplavin, can you please share what the error was in regard to theano broadcasting and how you fixed it. That will allow others to overcome this error if they face it in the future. Thank you

11

As I say, it’s not clear to me what the problem was. I fixed it with manually broadcasting two numpy arrays which were 2nd and 3rd arguments to tt.where. Everything worked correctly before this fix when I set parameters like jm.B1 = 1 and call jm.image().eval(), but somehow failed when using this with pymc. After this fix I managed to sample this model and get correct results.

12

thank you