Thank you, @ricardoV94! I suspect the ndim_supp is the culprit here. We had to define ndim_supp=0 even though the RandomVariable is multivariate, because all its parameters are scalars, and setting ndim_supp=1 would result in errors. I know this is kind of a hack and violates some assumptions in PyTensor, but it has worked well for us until now. Do you know if there are ways to get around this?