You’re just doing a deterministic transformation of variables, so I think you’re overthinking it? As least as written in the PyMC2 code, the function will know both `r`

and `f`

when it’s time to do the integral, so you just need to re-write the function to be “pytensor aware”.

Actually the `cumulative_trapezoid`

function is quite simple. Here is a stripped down version with all the input checks removed:

```
def tupleset(t, i, value):
l = list(t)
l[i] = value
return tuple(l)
def cumulative_trapezoid_simple(y, x=None, dx=1.0, axis=-1):
y = np.asarray(y)
if x is None:
d = dx
else:
x = np.asarray(x)
if x.ndim == 1:
d = np.diff(x)
else:
d = np.diff(x, axis=axis)
nd = len(y.shape)
slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
return res
```

It’s just some fancy slicing (which pytensor supports), differencing (supported), and cumsum (supported). Actually the only thing that won’t work is `len(y.shape)`

, but it’s bizarre they used this anyway (`ndim`

exists, and they know it since it was used for `x`

). So switch that, and switch the `np`

to `pt`

and get the following result:

```
from pytensor import tensor as pt
def cumulative_trapezoid_pt(y, x, dx=1.0, axis=-1):
if x is None:
d = dx
else:
if x.ndim == 1:
d = pt.diff(x)
else:
d = pt.diff(x, axis=axis)
nd = y.ndim
slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
res = pt.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
return res
```

This function you can drop into your PyMC model. Here’s a small numerical check to make sure it works:

```
f = lambda x: 2 * x
x = np.linspace(0, 10, 100)
y = f(x)
# This code compiles the pytensor function, you **will not** need this in your pymc model
# it's included for completeness/transparency
x_pt = pt.dvector('x')
y_pt = pt.dvector('y')
res = cumulative_trapezoid_pt(y_pt, x_pt)
cumint = pytensor.function([y_pt, x_pt], res)
# Here are the actual checks
from scipy import integrate
np.allclose(integrate.cumulative_trapezoid(y, x), cumint(y, x)) # Out: True
np.allclose(cumint(y, x), x[1:] ** 2) # Out: True
```