I ran into this as well when trying to track down the first bug, I was hoping it was just something wrong with my system. I think there’s a glitch that sets the shape of the observed data as the number of classes in the categorical (rather than inferring it from the length of p).
Could you test that:
There’s no error if the length observed data (still with missing values) is exactly equal to the length of p, and;
If (1) works successfully, that if you sample from data with fewer observations than the number of classes, the largest class you get in your samples is the length of the data, not the length of p?
If I’m right I think I can track the bug down pretty quick.