Dealing with multiple closely related models

Questions version agnostic
9

This function would work as-in. Inside a with pm.Model() scope, it would automatically add all the declared RVs in the active branch to the model as expected.

While I don’t recommend it, if you really wanted to lean into all the automagic stuff that PyMC does, you could return nothing from the function and access the created RVs using the model, e.g:

with pm.Model() as m:
    hierarchical_normal('mu', dims=(1,), center=True)
    alpha = m['mu'] # it's saved as mu because that's the name in the pm.Normal

This wouldn’t work if center = False, because alpha = mu_alpha + sigma_alpha * alpha_std isn’t registered (We have pm.Deterministic to let users register arbitrary computation to a model).

Obviously it works fine if it is used as you originally intended: alpha = hierarchical_normal('mu', dims=(1,), center=False)