Hello everyone, I have a few concerns about using PyMC4 for my MCMC needs, and I would like to outline them here in hopes of getting some guidance.

## Brief summary of the model

I have a model that gives me a function u'(t | \theta), which I have to solve using an ODE solver to obtain u(t | \theta), which I will then have to integrate and compute \int_0^T u(t | \theta) dt, where T also depends on the parameters.

### ODE solvers

In my search I came across sunode, a Python package that allows me to solve an IVP for ODEs. I havenâ€™t tried it out, but its usage looks somewhat straightforward and Iâ€™m currently learning how to apply it to my model.

### Integration routines

What concerns me is that there seems to be no integration routine, or at least one that is obvious to me that is compatible with PyMC4.

I say this because I came across this post on StackOverflow, named â€śCustom Theano Op to do numerical integrationâ€ť, which suggest a particularly complicated solution based on *Theano*, which was deprecated in favor of *Aesara*. The hyperlink that I use to go to the repository of the latter actually leads me to *PyTensor*, which I imagine is the replacement for *Aesara*.

Posts such as â€śCustom theano Op to do numerical integrationâ€ť, " PyMC3 : Using a parameter as a limit of integral" and â€śUsing Pymc3 to do forecasting and numerical integrationâ€ť (which actually leads to the first forum post) all work based on *Theano*.

### Lambert function

Additionally, in order to compute one of the quantities in my model, I need to have the Lambert W function in PyStan.

This function is available in SciPy, however, I canâ€™t use it directly in PyMC. I imagine the problem is the same as any other external function that I want to use, and in order to have it work correctly I will most likely face the same challenges as the people trying to implement the integration routines in the posts I have mentioned in the previous section.

### Questions

My questions are then as follows:

- Has anybody used this ODE solver successfully in PyMC4?
- What integration routine should I use in PyMC4?
- How can I use the Lambert W function in PyMC4?

If both questions 2. and 3. do not have an answer, is it recommended to downgrade to a previous version of PyMC? If so, which one? And which integration routine and Lambert W functions would I use?

If you are interested in the model itself, I will share it with you here in further detail.

Feel free to skip this section if you donâ€™t want to hear about it.

First, I must solve the ODE

\left.u^{\prime}(t | \theta) = \frac{e^{-\lambda u^2}}{\lambda u\left(u^{-2}-2 \lambda\right)-u^{-3}}\left[\frac{3}{2}\Omega_m(1+t)^2+2 \Omega_r(1+t)^3\right)\right],

where \Omega_m and \Omega_r are the two parameters in this model, and \lambda is derived from the parameters by relying on the main branch of the Lambert function

\lambda = \frac{1}{2} + W_0\left( - \frac{\Omega_m + \Omega_r}{2e^{1/2}} \right) .

The final step is to compute

\int_0^T u(t | \theta) dt,

where T is a function of the parameters \Omega_m and \Omega_r with a rather complicated expression which I will omit here for brevity, as I believe it only matters that it depends on the parameters.

This integral can then be directly compared with the observations.

I would also like to share with you the reason why Iâ€™m using PyMC4: because Stan has a bug, which I found, that doesnâ€™t have an ETA on its solution.

Apparently, using an integration routine on top of an ODE solver when applied to this system somehow gives a segmentation fault.

If you are interested, I opened a bug report in the Stan math repository: Segmentation fault on making use of `ode_rk45` with `integrate_1d` Â· Issue #2848 Â· stan-dev/math Â· GitHub

If something isnâ€™t clear, do let me know.

Thank you in advance!