TP.Marginal implementation or workaround for proper student's t process model

[AndrewFalkowski]

I am looking to implement a student's t process as described in Shah et al.. The TP implementation in PyMC lacks a .Marginal method needed to properly implement regression with student t processes. This was discussed in a 2018 PyMC forum post but hasn't been implemented to my knowledge.

Following the paper, I would expect a model to look something like this, with the noise term applied to the covariance matrix. I'm not entirely sure what the observational model would look like since I am adding in noise at the covariance matrix.

with pm.Model() as tp_model:

    ell = pm.Gamma("ell", alpha=2, beta=1)
    eta = pm.HalfNormal("eta", sigma=5)

    noise = pm.HalfNormal("noise", sigma=1.0)

    k_theta = eta**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=ell)
    # k = kθ + δ
    k = k_theta + pm.gp.cov.WhiteNoise(sigma=noise) # independent noise

    nu = 1 + pm.Gamma("nu", alpha=2, beta=1)
    tp = pm.gp.TP(scale_func=k, nu=nu)
    
    ...

How challenging would this be to implement or are there any workarounds available? The paper provides clear formulas for both marginal likelihood (equations 4-5) and conditional distributions (equation 6).

Thanks in advance!