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

Thanks in advance,

Kenneth