How to extend an existing model

If I have a model such as:

y = np.random.randn(10)

with pm.Model() as model_1:
    noise = pm.Gamma('noise', alpha=2, beta=1)
    y_observed = pm.Normal('y_observed', mu=0, sigma=noise, observed=y)

and after exploring this model I want to change the model, to:

with pm.Model() as model_2:
    noise = pm.Gamma('noise', alpha=2, beta=1)
    y_observed = pm.Exponential('y_observed', mu=noise, observed=y)

but I don’t want to rewrite the whole model again (which in reality would be larger, and 99 lines may be the same and one line changed or added). How do I do this? Doing:

with model_1 as model_2:
    y_observed = pm.Exponential('y_observed', mu=noise, observed=y)

results in an error because y_observed already exists. I just want to update its contents. I believe I can add variable by putting the model definition into a custom class and instantiating and building on that, but I don’t know how to update a model by changing variables. This includes Deterministics as well.

You can create the def function that will choose either of two models depending on the argument (e.g. you put the if-statement inside that function)

def F(parameter):
  with pm.Model() as model:
    ...
  if (parameter):
    with model:
      ...
  else:
    with model:
      ...
  ...
  return model

Thank you, this sounds like a good way to specify different configurations of a model without duplicating code. But is there any way to replace parts of a model (not data) that already exist?

That’s a really interesting question! I like using a “functional” pattern, though you could do this with graph rewriting in Theano if you have a really specific need. Here’s an example of my approach:

def hierarchical_normal(name, model):
    with model:
        mean = pm.Normal(f'{name}_mean', 0, 1)
        sd = pm.HalfNormal(f'{name}_sd', 1)
    return model

def make_priors(model):
    model = hierarchical_normal('a', model)
    model = hierarchical_normal('b', model)
    return model

def model_version_one():
    model = pm.Model()
    model = make_priors(model)
    with model:
        c = pm.Normal('c', model['a'] + model['b'], observed=data)
    return model
1 Like