# Custom likelihood from a python function

**URL:** <https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637>\
**Category:** v5\
**Tags:** modeling\
**Created:** [March 17, 2023, 10:32am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637 "2023-03-17T10:32:33Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![geop](https://avatars.discourse-cdn.com/v4/letter/g/58f4c7/32.png) [@geop](https://discourse.pymc.io/u/geop)\
**Post date:** [March 17, 2023, 10:32am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/1 "2023-03-17T10:32:33Z")

</div>

Hello,

I have generated some data **D** (blue dots) from the non-linear model (red) that takes in a set of model parameters **M** =(m1,m2,…), where I want to estimate the posterior for, for example m1 and m2: P(m1,m2|D).

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/7/7a4004dd4b019a54f989305be26ef139dc779978.png)

How should the likelihood function be defined without a simple expression?

Should the following code in principle be valid? (the “bg” refers to another library for the model I use):

```auto
K_data = ... # The "K-noisy" data in the above plot
with pm.Model() as model_cem:
    phi = pm.Data('phi', phi_data, mutable=True)

    # Define priors
    m1 = pm.Normal("Kf", mu=0, sigma=1)

    # Standard deviation - Note HalfNormal!
    s = pm.HalfNormal("sigma", sigma=1)

    def testfunc(phi,m1): # This function calls on the model used to generate K_data
      K_d_cem,MU_d_cem = bg.rockphysics.contact_cement(K_qz, MU_qz, phi, phi_c=phic, Cn=Cn, Kc=K_qz, Gc=MU_qz, scheme=1)
      return bg.rockphysics.fluidsub.vels(K_d_cem,MU_d_cem,K_qz,RHO_qz,x1,RHO_b,phi)[-1]
    
    # Define the likelihood
    likelihood = pm.Normal("y", mu=testfunc(phi,m1), sigma=s, observed=K_data)

    step = pm.NUTS()

    trace = pm.sample(1000, tune=500, init=None, step=step, cores=2)
    

```

Thanks in advance,  
Kenneth

---

<div class="post-metadata">

**Author:** ![BrynWalton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/brynwalton/32/3921_2.png) [@BrynWalton](https://discourse.pymc.io/u/BrynWalton)\
**Post date:** [March 17, 2023, 11:43am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/2 "2023-03-17T11:43:41Z")

</div>

Hi, at a glance, I think i have had to do something similar to this. I found the “custon likelihood” tutorial in the docs very useful. Not sure how easy it would be to get the gradient of testfunc, so if you want to just get things running, you may need to use a metropolis hastings sampling step rather than NUTS. Hope this helps

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [March 18, 2023, 12:36am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/3 "2023-03-18T00:36:27Z")

</div>

[Here](https://www.pymc.io/projects/examples/en/latest/case_studies/blackbox_external_likelihood_numpy.html) is the custom pytensor `Op` documentation that @BrynWalton mentioned.

It might be possible that you don’t need to do this. It depends entirely on what is happening inside of `contact_cement` and `fluidsub.vels`. If they only do vectorized numpy-style operations (that is, no control flow, no loops, no exotic linear algebra), you can use the function out-of-the-box by just passing symbolic tensors as inputs. If there is some optional control flow (perhaps related to the `scheme` argumnet?) you could rip out only the bits you need and re-implement it in native pytensor. Making an `Op` should be your last resort.

---

<div class="post-metadata">

**Author:** ![BaptisteDE](https://avatars.discourse-cdn.com/v4/letter/b/3d9bf3/32.png) [@BaptisteDE](https://discourse.pymc.io/u/BaptisteDE)\
**Post date:** [May 16, 2024, 9:26am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/4 "2024-05-16T09:26:31Z")

</div>

> [@BrynWalton](#):
>
> custon

Hello !  
I have a similar problem, I would like to use the result of an FMU (Fuctional Mockup Unit) calculation as the mean of my likelyhood function. I already wrapped the FMU in a python function, but it doesn’t accept Pytensors.

I think this is similar to this discussion:

> [@Using PyMC3 with external functions](https://discourse.pymc.io/t/using-pymc3-with-external-functions/6403):
>
> I am trying to use PyMC3 to replace emcee as my MCMC sampler of choice. However, I depend on external functions to draw the model curves, which I import from existing packages that are not compatible with Theano’s style. How can I still use PyMC3 as a sampler? For example: #::: Option 1: How the external fct currently looks def external\_fct\_1(x, xmid, depth): y = np.zeros\_like(x) y[(x\>xmid-0.1) & (x\<xmid+0.1)] = -depth return y #::: Option 2: Rewriting the external fct in Theano…

I tried to find the tutorials listed here :  
_“Specifically, the tutorial on using a [blackbox likelihood function](https://www.pymc.io/projects/examples/en/latest/case_studies/blackbox_external_likelihood_numpy.html), the tutorial on [custom distributions](https://www.pymc.io/projects/examples/en/latest/howto/custom_distribution.html), the tutorial on [wrapping jax functions into PyTensor Ops](https://www.pymc.io/projects/examples/en/latest/case_studies/wrapping_jax_function.html).”_  
from [Likelihood Approximations with Neural Networks in PyMC - PyMC Labs](https://www.pymc-labs.com/blog-posts/likelihood-approximations-through-neural-networks/)

The one you mention @jessegrabowski is among them, but all the links seem to be broken.

Do you have any idea where I should look ?

Thank you !

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [May 17, 2024, 9:42pm UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/5 "2024-05-17T21:42:51Z")

</div>

Please find the correct links here:

- [How to wrap a JAX function for use in PyMC](https://www.pymc.io/projects/examples/en/latest/howto/wrapping_jax_function.html)
- [Using a “black box” likelihood function](https://www.pymc.io/projects/examples/en/latest/howto/blackbox_external_likelihood_numpy.html)
- [Defining a Custom Distribution in PyMC3](https://www.pymc.io/projects/examples/en/2022.01.0/pymc3_howto/custom_distribution.html)

---

<div class="post-metadata">

**Author:** ![BaptisteDE](https://avatars.discourse-cdn.com/v4/letter/b/3d9bf3/32.png) [@BaptisteDE](https://discourse.pymc.io/u/BaptisteDE)\
**Post date:** [May 27, 2024, 8:28am UTC](https://discourse.pymc.io/t/custom-likelihood-from-a-python-function/11637/6 "2024-05-27T08:28:50Z")

</div>

Thank you !

[Using a “black box” likelihood function](https://www.pymc.io/projects/examples/en/latest/howto/blackbox_external_likelihood_numpy.html) is particularly usefull.