Square decision region model specification - initial evaluation failed

Questions
6

Good call. I tried a couple of things to no avail, following examples in the Martin text.

The first was ordered like you suggest, I was a bit less familiar with this:

breakpoint_model = pm.Model()

with breakpoint_model:
    # Data
    loc_x = pm.Data("plate_x", right_full["plate_x"])
    loc_z = pm.Data("plate_z", right_full["plate_z"])
    outcome = pm.Data("outcome", right_full["isStrike"])
    
    horiz_bounds = pm.Normal(
        "horiz_bounds",
        mu=[-1.0, 1.0],
        sd=1.0,
        shape=2,
        transform=pm.distributions.transforms.ordered        
    )
    vert_bounds = pm.Lognormal(
        "vert_bounds",
        mu=[1.5, 3.5],
        sigma=1.5,
        shape=2,
        transform=pm.distributions.transforms.ordered        
    )
    
    db = pm.Deterministic(
        "decision",
        tt.switch( 
            (loc_x > horiz_bounds[0])  & (loc_x < horiz_bounds[1]) & (loc_z > vert_bounds[0]) & (loc_z < vert_bounds[1]),
            0.97,
            0.03
        )
    )
    
    yl = pm.Bernoulli('strike', p=db, observed=outcome)
    trace = pm.sample(500, init="advi", n_init=10000)

Which gives the NaN error in initialization.

FloatingPointError: NaN occurred in optimization. 
The current approximation of RV `vert_bounds_ordered__`.ravel()[0] is NaN.

The other try was to use a pm.Potential following one of the examples in the text:

breakpoint_model = pm.Model()

with breakpoint_model:
    # Data
    loc_x = pm.Data("plate_x", right_full["plate_x"])
    loc_z = pm.Data("plate_z", right_full["plate_z"])
    outcome = pm.Data("outcome", right_full["isStrike"])
    

    left_side = pm.Normal("left_boundary", -1.0, 1)
    right_side = pm.Normal("right_boundary", 0.8, 1)
    top_side = pm.Lognormal("top_boundary", 3.5, 1.5 )
    bottom_side = pm.Lognormal("bottom_boundary", 1.5, 2 ) # can't be negative

    db = pm.Deterministic(
        "decision",
        tt.switch( 
            (loc_x > left_side)  & (loc_x < right_side) & (loc_z > bottom_side)  & (loc_z < top_side),
            0.97,
            0.03
        )
    )

    # Force boundaries to be ordered
    order_vert = pm.Potential('order_vert',
                                  tt.switch(top_side - bottom_side < 0,
                                  -np.inf, 0)
                                 )
    order_horiz = pm.Potential('order_horiz',
                                  tt.switch(right_side - left_side < 0,
                                  -np.inf, 0)
                                 )
    
    
    yl = pm.Bernoulli('strike', p=db, observed=outcome)
    trace = pm.sample(500, init="advi", n_init=10000)

Which samples but still has divergences. I’ve only gotten a little time to play with it so I’ll do some more digging to see what I can figure out (no need to burn brainspace before I spin my wheels on it for a while, unless you’re super bored).

Also just for reference, an image of the data

image