Switch function has little constraint in sampling

Hi there :)
As the title says, I’m looking for a way to build a model where the variables are in same ranges with different values, however the result sampling suggest the sampling points are away from the switch constraints.
So, what do you think about this?

    with pm.Model() as model:
        rs = pm.Normal('rs', 800, 50)
        rp = pm.Normal('rp', 800, 50)
        re = pm.Normal('re', 800, 50)
        pm.Potential('const1', tt.switch(tt.gt(rp, rs), 0., -1E10))
        pm.Potential('const2', tt.switch(tt.gt(re, rp), 0., -1E10))
        rhop = pm.Uniform('rhop', 1, 5)
        rhoe = pm.Uniform('rhoe', 1, 5)
        pm.Potential('const3', tt.switch(tt.gt(rhop, rhoe), 0., -1E10))
        aa = pm.Deterministic('aa', np.sqrt(-np.log(rhoe / rhop)))
        trace = pm.sample(1000)
    data = az.from_pymc3(trace)
    pm.summary(data, hdi_prob=0.95)
         mean      sd  hdi_2.5%  hdi_97.5%  ...  mcse_sd  ess_bulk  ess_tail  r_hat
rs    793.549  48.173   698.735    888.375  ...    9.795      13.0     283.0   1.13
rp    803.784  49.538   710.618    904.680  ...    2.813     157.0     294.0   1.02
re    805.033  41.822   727.193    890.052  ...    2.876     105.0     188.0   1.00
rhop    2.821   1.165     1.000      4.724  ...    0.154      32.0     126.0   1.06
rhoe    3.064   1.203     1.165      4.995  ...    0.140      38.0     109.0   1.05
aa        NaN     NaN     0.027        NaN  ...      NaN       NaN       NaN    NaN

What exactly do you expect to see from this model that you are not currently seeing?

Thank you for your quickly answer. I want to transfer aa to another function. however the prior sampling points show Nan. do you have any idea about it?

Potentials don’t play a role in prior or posterior predictive sampling since they are only part of the logp graph, not the random one.

So if I want to transfer aa to a pm.Simulator module, how to add the constraints properly to make sure the aa is always float not nan?

You can try to parametrize your model in a way that it respects the constraints naturally. For instance, instead of modelling rp and re directly, you can model the shift from the previous variable.

with pm.Model():
  rs = pm.Normal('rs', 800, 50)
  rp_shift = pm.LogNormal('rp_shift', 0, 1)
  re_shift = pm.Lognormal('re_shift', 0, 1)
  rp = rs + rp_shift
  re = rs + rp_shift + re_shift

The priors will be trickier to reason about, but it isn’t too bad. You can wrap rp and re in pm.Deterministic if you need them in the trace.

Thanks! I will try it.