Using the compile_logp function to evaluate the log-posterior of a point estimate

Questions v5
1

I am trying to compute the log-posterior of a hierarchical model for a point estimate of the model parameters, which is placed in a dictionary.

I have tried using model.compile_logp(point_estimate_dict) with point_estimate_dict including the estimates of the free random variables of the model, but it returns the following error:

TypeError: Unknown input or state: θ_mu. The function has 8 named inputs (θ_mu_log__, ...

Which of the model variables should be used in the point estimate dictionary? Should I include both the parameters and hyperparameters? Why does the compile_logp function require to add _log_ after my parameter names?

Your help will be appreciated.

Thank you in advance!

2

Should be model.compile_logp()(poont_estimate_dict)

The first call simply compiles and returns a function (you should store this if you intend to do multiple evaluations). The second actually evaluates it.

You should include all the unobserved variables.

Because those variables are transformed on the log scale, so your input should actually be the log of the variable. The name is to avoid surprises. You can disable a variable transformation by passing transform=None when defining it in the Model.

Finally, you can call model.initial_point() to have an example dict of all the variable names and values/shapes that are required for the model logp. model.compile_logp()(model initial_point()), should always work.

3

Thank you, your response has been quite helpful.

Now, I can calculate the log posterior by setting transform=None. Does this option impact the inference?

Also, I can give you some more specifics of what I want to do. My idea is to derive a point estimate (such as the mean of the posterior or MAP) from the posterior trace. This leads to a dictionary with keys being the model.free_RVs.

In case I don’t use transform=None, how would I know the transformed variable names in order to create the point estimate dictionary? Also, would I have to calculate manually their point estimate by transforming the trace accordingly?

4

Depends on what you want to do with the model. NUTS and find_MAP will struggle with transform=None because of the hard parameter constraints when transformations are not used.

And you want to evaluate the model logp at this mean? Or you just want the mean? If you want to use it as an input, you will have to transform it. We should probably add an option to disable transforms in compile_logp, to make this easier.

You can obtain the value variable that corresponds to each rv by doing model.rvs_to_values[rv]. You can then check the name of this variable.

import pymc as pm

with pm.Model() as m:
    x = pm.Normal("x")
    y = pm.HalfNormal("y")
    z = pm.Normal("z", x, y, observed=0)

value_vars = [m.rvs_to_values[rv] for rv in m.free_RVs]
value_var_names = [v.name for v in value_vars]
value_var_transforms = [getattr(v.tag, "transform", None) for v in value_vars]

print(value_var_names)
print(value_var_transforms)

# ['x', 'y_log__']
# [None, <aeppl.transforms.LogTransform object at 0x7f78de2f01c0>]