I believe this is redundant in INLA and is only necessary for Discrete marginal stuff (and thus it is fine to define this as None). Please correct me if I'm wrong.

Yes I think in the refactor WIP PR the base marginalRV class doesn't have this property

In lines 613-619, marginalization_op is defined consistently regardless of the MarginalRV it's built from. For compactness, I think it's fine to set it to None through this else block rather than having a seperate marginalize_constructor line?

marginalization_op = marginalize_constructor(
    inputs=inner_inputs,
    outputs=inner_outputs,
    dims_connections=dependent_rvs_dim_connections,
    dims=dims,
    **marginalize_kwargs,
)

Not within the scope of this PR, but is there a nice (i.e. pytensor-native) way to specify an RNG vector of shape (d,) without np.random? I believe that rng.random(d) crashes when we set d = marginalized_vv.shape[0] because marginalized_vv is a pt tensor, so it's unhappy with d being something symbolic rather than an int. For now I've worked around it by hardcoding d.

If it's a one-liner fix, I'll include that in this PR of course.

pt.random.uniform(size=x0.shape)

Throws the following:

ValueError: Random variables detected in the logp graph: {uniform_rv{"(),()->()"}.out, uniform_rv{"(),()->()"}.out}.
This can happen when mixing variables from different models, or when CustomDist logp or Interval transform functions reference nonlocal variables.

Yes I think in the refactor WIP PR the base marginalRV class doesn't have this property