# Custom likelihood from a python function

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).

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):

``````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)

``````

Here 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.