Calculation of joint density from pymc model

Questions v5
1

Hello,
I’m trying to compute the joint log density of a PyMC model with respect to different input data.

As far as I understand, once observed is set in a PyMC model, it’s fixed and cannot be easily changed when evaluating the log density. Is that correct?

The function below does return the correct result, but it’s too slow since it rebuilds and recompiles the model every time. Is there a better or more efficient way to do this?

def cal_log_likelihood(theta, data):
    with pm.Model() as test_model:
        mu = pm.Normal('mu', 0, 1, shape=(2,))
        Sigma = pm.InverseGamma('Sigma', 2, 1, shape=(2,))
        _ = pm.MvNormal('data', mu, pt.diag(Sigma), observed=data, shape=(data.shape[0], 2))
    return test_model.compile_logp()(theta)

for i in tqdm(range(100)):
    cal_log_likelihood(
        {'mu': np.array([1, 1]), 'Sigma_log__': np.log(np.array([2, 3]))},
        np.random.multivariate_normal(mean, cov, size=1)
    )
2

If you use pm.Data for the observed data you can change it between evaluations or you can define the model without making the observed variable observed.

model.compile_logp() gives you the compiled logp function. The inputs have the format of model.initial_point().

If you remove the observed kwarg, the data is part of that initial_point dictionary.

Otherwise, if you use pm.Data, data will be a shared variable that can be modified with data.set_value(new_value) between calls.

1 Like