Oh ok, thanks for sharing! Always looking for more covariance functions to add.
Regarding your 2nd question, yes that’s right. This is true for all the stationary covariance functions. For example you can do:
with model:
el = pm.HalfCauchy(“el”, beta=[5, 5], shape=(2,))
cov = pm.gp.cov.ExpQuad(2, el)