I’m trying to specify a model in which I have a roughly square decision region, inside the square is true, outside is false. My approach was to basically glue two invlogit functions together in either dimensions - one as a decision “turn on” and one as a “turn off.” I managed to make this with a theano switch:

```
x = np.linspace(-5,5,2000)
change1 = -2
change2 = 2
turnon_coeff = 5
turnoff_coeff = 5
dist = tt.switch(x< change1 + (change2-change1)/2,
pm.invlogit(turnon_coeff*(x - change1)),
1-pm.invlogit(turnoff_coeff*(x - change2))
).eval()
```

Now, I’d like to model this and infer the change points and the slopes of the invlogits, but I keep getting a failure for the inital evaluation of the model, and can’t seem to locate what’s misspecified.

The model looks like the following:

```
sigmoid_model = pm.Model()
with sigmoid_model:
# Data
loc_x = pm.Data("loc_x", df["loc_x"])
loc_z = pm.Data("loc_z", df["loc_z"])
outcome = pm.Data("outcome", df["isTrue"])
# priors
left_side = pm.Normal("left_boundary", -1, 0.75)
right_side = pm.Normal("right_boundary", 1, 0.75)
top_side = pm.Normal("top_boundary", 3.5, 1.5 )
bottom_side = pm.Lognormal("bottom_boundary", 1.5, 1.5 ) # can't be negative
# My only prior here is that I think it's *fairly* steep and positive, but not anything strong
turnon_coeff_horiz = pm.Lognormal("turnon_coeff_x", 5, 5)
turnoff_coeff_horiz = pm.Lognormal("turnoff_coeff_x", 5, 5)
turnon_coeff_vert = pm.Lognormal("turnon_coeff_z", 5, 4)
turnoff_coeff_vert = pm.Lognormal("turnoff_coeff_z", 5, 4)
horiz_p = pm.Deterministic("horiz_p",
tt.switch(
loc_x < left_side + (right_side - left_side)/2,
pm.invlogit(turnon_coeff_horiz*(loc_x - left_side)),
1-pm.invlogit(turnoff_coeff_horiz*(loc_x - right_side))
)
)
vert_p = pm.Deterministic("vert_p",
tt.switch(
loc_z < bottom_side + (top_side - bottom_side)/2,
pm.invlogit(turnon_coeff_horiz*(loc_z - bottom_side)),
1-pm.invlogit(turnoff_coeff_horiz*(loc_z - top_side))
)
)
# Making the probability the sum of the probabilities, however clipping at 0.99
# to cover random judgement mistakes inside/outside the decision region
p_add = pm.Deterministic("p_add" , pm.math.clip(horiz_p + vert_p, 0.01, 0.99))
yl = pm.Bernoulli('isTrue', p=p_add, observed=outcome)
trace = pm.sample(1000, init="advi", n_init=10000)
```

The error it gives me is the following:

```
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'left_boundary': array(-1.), 'right_boundary': array(1.), 'top_boundary': array(3.5), 'bottom_boundary_log__': array(1.5), 'turnon_coeff_x_log__': array(5.), 'turnoff_coeff_x_log__': array(5.), 'turnon_coeff_z_log__': array(5.), 'turnoff_coeff_z_log__': array(5.)}
Initial evaluation results:
left_boundary -0.63
right_boundary -0.63
top_boundary -1.32
bottom_boundary_log__ -1.32
turnon_coeff_x_log__ -2.53
turnoff_coeff_x_log__ -2.53
turnon_coeff_z_log__ -2.31
turnoff_coeff_z_log__ -2.31
isTrue -inf
Name: Log-probability of test_point, dtype: float64
```

Any suggestions here? Something blatently obviously wrong with the model? (Or is there a less silly way to go about this problem?)