Help with `logp` of a threshold model

Questions
2

A direct translation of the notation would be something like:

eps = 0.00001  # <== set it to 0 to have the exact notation, but
               #     you might ran into problem sampling
with pm.Model() as m:
    beta = ...
    psi = pm.Normal('psi', X@beta, 1.)
    observed = pm.Bernoulli('y', tt.switch(psi <= 0, eps, 1.-eps), observed=y)

However, a switch breaks the gradient and you might ran into problem sampling with NUTS, which you should replace it with a smoothstep function