PyMC4: custom likelihood

Hello pymc community,

Using PyMC4, I am trying to use a custom likelihood wrapped in a python function call (it computes a misfit between simulated data from a physics PDE and observations). It also has several variables. I have not been able to figure out the syntax to do that yet. Below is a minimalistic example of the syntax I have pieced together. The code runs BUT:

  • my output idata has no likelihood section, only posterior and sample_stats
    image

  • in the as_op, I need to add an additional variable for the name, otherwise the code returns an error
    image

What should I be doing instead to use my custom likelihood, with pymc 4, and multiple variables? A code snippet would be deeply appreciated as I am still new to pymc.

Thanks in advance

import pymc as pm
from aesara.compile.ops import as_op
import numpy as np
import aesara.tensor as at

basic_model = pm.Model()

def custom_log_likelihood(x,y):
    # do something, for simplicity here let just put a quadratic so that the code run
    #return log_likelihood
    return -0.5 * (x**2 + y**2)  # placeholder log-likelihood

with basic_model:

    # Priors for unknown model parameters
    x0p = pm.Normal("x0p", mu=0, sigma=10)
    y0p = pm.Normal("y0p", mu=-30, sigma=10)
    
    # I somehow need an additional variable here for the name
    @as_op(itypes=[at.dscalar, at.dscalar,at.dscalar], otypes=[at.dscalar])
    def logp(name, x0p, y0p):
        return custom_log_likelihood_quadratic(x0p, y0p)

    likelihood = pm.DensityDist("likelihood",x0p,y0p,logp=logp)
        
    start = pm.find_MAP()
    # draw 1000 posterior samples
    idata = pm.sample(step=pm.Metropolis(),cores=4, start=start) # 1, tune=1, start=start, 

I think you might be confusing the pm.DensityDist with pm.Potential:
The DensityDist is designed for situations where you want to “manually draw the PDF” of a prior distribution. It has nothing to do with the likelihood.

It sounds like you want to use a custom logp/custom likelihood which you can do using pm.Potential.
For that you can make the necessary computations with pure Aesara Ops if you want them to be differentiable, or you can implement a custom Aesara Op with or without a grad() implementation. Or you can use the as_op decorator like you’re doing already.

So try to replace the DensityDist with this:

pm.Potential("likelihood", logp("Dummy name", x0p, y0p))