Logp questions, synthetic dataset to evaluate modeling

[ricardoV94]

@LeiAsplanatoLorenzo

Also, the reported logp while running find_MAP differs from the one reported from compile_logp()(MAP_point).

This happens because find_MAP ignores the jacobian. You should get the same if you do compile_logp(jacobian=False)(map)

@LeiAsplanatoLorenzo

If the argument “initval” is passed within the model, then both the tuning steps AND the draws are initialized in those values, right?

Yes, although the default sampler adds an amount of jittering around these values. You can choose another init method that does not.

@LeiAsplanatoLorenzo

In this case, the model has a Dirichlet prior. When this happens, the initial_values report a dirichlet_simplex__, and that is what is needed in order to run ```model.compile_logp’‘’. How can I go from the variable’s concentration values to the simplex that is needed? (for example, in any other of the chains’ values, as asked in 1)

Hmm we should probably add a utility to do that. For now you would need to do something manual like

import pymc as pm

with pm.Model() as m:
    x = pm.Uniform("x")
    y = pm.Dirichlet("y", [x, x])
    
transformed_rvs = []
for free_rv in m.free_RVs:
    transform = m.rvs_to_transforms.get(free_rv) 
    if transform is None:
        transformed_rvs.append(free_rv)
    else:
        transformed_rv = transform.forward(free_rv, *free_rv.owner.inputs)
        transformed_rvs.append(transformed_rv)
        
fn = m.compile_fn(inputs=m.free_RVs, outs=transformed_rvs)
res = fn({"x": .75, "y": [0.3, 0.7]})
res_dict = {
    "x_interval__": res[0],
    "y_simplex__": res[1],
}
print(res_dict)
# {'x_interval__': array(1.09861229), 'y_simplex__': array([-0.42364893])}

# Going the other way around is easier:
# unobserved_value_vars includes the original variables and deterministics, so you have
# to pick the outputs manually. Or you could look into the source code and adapt to return only what you need.
fn_inv = m.compile_fn(outs=m.unobserved_value_vars)
res = fn_inv(res_dict)
res_dict = {
    "x": res[2],
    "y": res[3]
}
print(res_dict)
# {'x': array(0.75), 'y': array([0.3, 0.7])}

Alternatively if you don’t care about transforms after you fit, you could recreate your model with all variables (transform=None), so that things like the simplex__ don’t show up. Then you will be able to use the traces to evaluate the logp from model.compile_logp() out of the box.