Hi, bwengals, thank you for your reply. I had previously raised a question about the possibility of implementing an ‘cumtrapz’ within the model and received a definitive answer. I have recently revisited this issue and would like to re-edit this question to make it appear more complete. The formula is as follows:

{U(x)=\frac{s}{f_{1}(x)(f_{4}(x)-f_{5}(x))}\int_{0}^{x}f_{1}(x^{\prime})\Big(f_{2}(x^{\prime})-f_{3}(x^{\prime})\Big)dx^{\prime}}

f1-f5~gp, each f with observed values that have a fixed variance., s~Normal, finally I am looking to obtain the posterior of f.

The v2code can be roughly described as:

```
M = pm.gp.Mean(...)
C = pm.gp.Covariance(...)
# Some observations
pm.gp.observe(M, C, r, y, obs_variance)
# Setting initial values
init_vals = pm.gp.Realization(M, C)(self.x)
gp1 = pm.gp.GPSubmodel(..., M, C, self.x, init_vals=init_vals)
...
# gp1 through gp5 are all created using the gp model mentioned above
@pm.deterministic
def func(gp1, gp2, gp3, gp4, gp5, s, x):
f = gp1.f_eval * (gp2.f_eval - gp3.f_eval)
um = cumtrapz(f, x)
us = s * um / (gp1.f_eval * (gp4.f_eval - gp5.f_eval))
return us
u = func
v = pm.Normal(..., mu=u, tau=..., observed=True, value=...)
# sample
```

I have looked briefly at the PyMC2 documentation on Gaussian processes, which states that after `pm.gp.observe`

, M and C are transformed into M0 and C0, in `gp1 = pm.gp.GPSubmodel(..., M, C, ...)`

, M and C are actually M0 and C0.

So, what do I need to pass to this ‘cumtrapz’ function? If I use `latent`

, should I define its mean and variance as M0 and M0, for example:

```
M0 = pm.gp.mean...
C0 = pm.gp.cov...
gp1 = pm.gp.Latent(M0, C0)
f1 = gp.prior("f1", X=x)
# f1 through f5 are obtained by the same procedure
f = f1 * (f2 - f3)
u_m = cumulative_trapezoid_pt(f, r) / (f4 * (f5 - f6))
u_s = pm.Deterministic("u_s", s * u_m)
u = pm.Normal("u", mu=u_s, sigma=10, observed=observe_value)
pm.sample()
```

Or using `marginal`

, for example:

```
M = pm.gp.Mean(...)
C = pm.gp.Covariance(...)
gp = pm.gp.Marginal(M, C)
y_ = gp.marginal_likelihood("y", X=r, y=y, sigma=10)
f1 = gp.conditional("f1", x)
##
f = f1 * (f2 - f3)
u_m = cumulative_trapezoid_pt(f, x) / (f1 * (f4 - f5))
u_s = pm.Deterministic("u_s", s * u_m)
u = pm.Normal("u", mu=u_s, sigma=10, observed=observe_value)
pm.sample()
```

Forgive my immature ideas, my knowledge in this area is quite limited. If you have any further suggestions or guidance, I would greatly appreciate it. Thank you!