I am learning bayesian inference and gaussian process. I am trying to use pymc to do a simple example which involves a partial differential equation in the form of (sorry I don’t know how to type latex here):

frac{partial}{partial r}(A(r)theta(r)partial B(r)/partial r)=C(r)

These are 1-D profiles with A(r), B(r) C(r), r known profiles. I want to infer the theta(r) with bayesian inference. Unlike the example presented in the pymc which is a linear relation, I don’t know how to deal with the derivative in the above equation.

I tried the following code but it doesn’t allow me to perform the derivative. I am wondering how to write the derivatives in the gaussian process with pymc? Or some alternative tools that allows me to do it?

with pm.Model() as model:
    # GP prior for theta
    length_scale = pm.Gamma("length_scale", alpha=2, beta=1)
    cov_func = pm.gp.cov.ExpQuad(1, ls=length_scale)
    theta_gp = pm.gp.Latent(cov_func=cov_func)
    theta = theta_gp.prior("chi", X=rho[:, None])

    # Compute derivatives
    dB_dr = np.gradient(B, r)
    A_theta_dB_dr = A * theta * dB_dr
    derivative = np.gradient(A_theta_dB_dr, r)

    # Likelihood
    S_obs = pm.Normal("S_obs", mu=derivative, sigma=0.1, observed=C)

    # Inference
    trace = pm.sample(1000, return_inferencedata=True)