Moving the RVs out of the scan seems to have done the trick.
I had to bundle the varying-length time-series into a rectangular matrix so that one response RV of shape (N, M) could replace the N response RVs of shape (M - k,).
To mask out the invalid padding values, I set them to be NaNs in the observed vectors and am replacing the NaNs in the result of logp with zero using pt.where; wrapped in a pm.Potential.
Using an indexing mask seems to not work with “ragged” indexes.
@ricardoV94 , I have marked your recommendation as the solution. Thank you, sir!