"log_likelihood" not found in InferenceData

Questions
1

The arviz documentation suggests that PyMC3 should be storing a log-likelihood with the InferenceData when “log_likelihood = True” is one of the `idata_kwargs.

step_func = pm.DEMetropolisZ(vars = params_list)
trace = pm.sample(
    draws = NUM_DRAWS, step = step_func, cores = NUM_CORES,
    chains = NUM_CORES, return_inferencedata = True,
    idata_kwargs = {'log_likelihood': True})

However, the output trace (InferenceData`) does not contain a “log_likelihood” group.

>>> trace
Inference data with groups:
	> posterior
	> sample_stats

Am I missing something? Thanks!

2

So, thankfully, there is a way to:

These are sufficient for now, but better support for InferenceData might be worth pursuing. From the PyMC3 trunk, it’s apparent that:

3

which versions of pymc and arviz are you using?

Yes, when converting to InferenceData the log likelihood data is automatically computed and included in the resulting InferenceData if possible. I think there should be a warning printed if that fails, but I might be wrong. You can see that this works in most of the example notebooks, for example: https://pymc-examples.readthedocs.io/en/latest/diagnostics_and_criticism/model_comparison.html

That being the default also means that using idata_kwargs={"log_likelihood": True} is equivalent to not using any arguments.

If either of these options has worked, then it seems like the auto conversion should work too, provided you have relatively new versions of both pymc and arviz. If so, could you share a minimal example in which it fails?

Here you are wrong, which is not strange given the actual change and large updates happening to the codebase as version 4.x approaches.

The to_inference_data you mention is a method of the converter class, not a function. The function that is actually called from the sample method in the 2nd link you provided is a few lines below that. And it does have keyword arguments. However, those are only used for >=4 versions, which you’ll only have if you have installed from github as it has not been released yet. You can also see how in the methods above there is code similar to the answers you shared above that automatically store pointwise log likelihood values.

If you are using pymc3 installed from pip or from conda, the conversion is then being done with an ArviZ function instead of with this one, which is why you originally found that info in the ArviZ docs. Both are similar and both store log likelihood by default, but there are some important differences.

1 Like
4

I am using ArviZ version 0.11.4, PyMC3 version 3.11.4.

Here is a minimal example in which the log_likelihood is not saved in the InferenceData (i.e., not in trace).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm
import theano.tensor as tt
from matplotlib import pyplot

NUM_CORES = 4
NUM_DRAWS = 25000

class BlackBoxLikelihood(tt.Op):
    '''
    Specify what type of object will be passed and returned to the Op when it is
    called. In our case we will be passing it a vector of values (the parameters
    that define our model) and returning a single "scalar" value (the
    log-likelihood)
    '''
    itypes = [tt.dvector] # Expects a vector of parameter values when called
    otypes = [tt.dscalar] # Outputs a single scalar value (the log likelihood)

    def __init__(self, model, observed, x, loglik = None):
        '''
        Initialise the Op with various things that our log-likelihood function
        requires. Below are the things that are needed in this particular
        example.

        Parameters
        ----------
        model : Callable
            An arbitrary "black box" function that takes two arguments: the
            model parameters ("params") and the forcing data ("x")
        observed : numpy.ndarray
            The "observed" data that our log-likelihood function takes in
        x:
            The forcing data (input drivers) that our model requires
        loglik : Callable or None
            The log-likelihood (or whatever) function we've defined
        '''
        self.model = model
        self.observed = observed
        self.x = x
        # Optionally overwrite the Gaussian log-likelihood function
        if loglik is not None:
            self.loglik = loglik

    def loglik(self, params, x, observed):
        predicted = self.model(params, x)
        sigma = params[-1]
        return -0.5 * np.log(2 * np.pi * sigma**2) - (0.5 / sigma**2) *\
            np.nansum((predicted - observed)**2)

    def perform(self, node, inputs, outputs):
        # The method that is used when calling the Op
        (params,) = inputs
        logl = self.loglik(params, self.x, self.observed)
        outputs[0][0] = np.array(logl) # Output the log-likelihood


def main():
    def model(params, x, n = 1000):
        'Our black-box model; really, a simulator of Gaussian random vars'
        return np.random.normal(*params, n)

    # Observed data
    data = np.random.normal(loc = 0, scale = 1, size = 1000)
    loglik = BlackBoxLikelihood(model, data, None)

    basic_model = pm.Model()
    with basic_model:
        # (Stochstic) Priors for unknown model parameters
        a = pm.Normal('a', mu = 0, sd = 5)
        b = pm.HalfNormal('b', sd = 1)

        # Convert model parameters to a tensor vector
        params = tt.as_tensor_variable([a, b])
        pm.Potential('likelihood', loglik(params))

        trace = pm.sample(
            draws = NUM_DRAWS, step = pm.DEMetropolisZ(vars = [a, b]),
            cores = NUM_CORES, chains = NUM_CORES,
            return_inferencedata = True, start = {'a': 1, 'b': 10},
            idata_kwargs = {'log_likelihood': True})

Output trace looks like the following, even if idata_kwargs is omitted:

>>> trace
Inference data with groups:
	> posterior
	> sample_stats
5

Log likelihood values can only be stored automatically for those variables that have been defined with the observed kwarg. You are using pm.Potential which can be used for blackbox likelihood terms but can also be used for any other goal with nothing to do with observed values such as breaking symmetry or imposing extra prior constraints on the parameters.

This is also the reason why you don’t see any warning. Potential terms are not interpreted as observed variables (not they can be by design) and so the issue is not that auto conversion fails but instead that it’s not triggered at all.

If you wish to have pointwise log likelihood data retrieved automatically, you should instead use pm.DensityDist for example, which might require to reformat your log likelihood function. After all, it needs to be somewhat compatible to a distribution and take an observed kwarg. Using a random variable as observed - #5 by OriolAbril covers multiple working ways of calling DensityDist. In your case you do have “proper” observed data, so anything that allows for DensityDist(..., observed=data) should work and trigger the converter to store log likelihood values automatically.

1 Like
6

Worked! I am on PyMC v5.3, and I guess the default was to not save it…

3 Likes