Inverse problem with Gaussian Processes

Hi!

Gaussian Processes are a well-known way to solve inverse problems of the form L[f]=g where g is a function, my observed data is a set \{g(x_i)\}, and the mapping is of the form L[f](y)=\int K(x,y) f(x)dx. The interesting property is that if f\sim{\cal GP} then L[f]\sim{\cal GP} as well.

This for instance would be one of the starting points for what I would like to do : Phys. Rev. D 105, 036014 (2022) - Reconstructing QCD spectral functions with Gaussian processes
What makes it non-trivial to implement in PyMC are the integrals in Eqs. (6) and (7). Without those, the problem would be just as simple as interpolating g with GPs.

This replaces older methods such as Backus–Gilbert, and has plenty of applications in many different domains.

I was wondering if anybody is aware of any similar application based on PyMC, or maybe some extra packages which could deal with that.