Thank you for the recommendation. I will try that out and get back on the result.
I am actually trying to replicate the Gaussian Process model in this paper in PyMC3
If you look at page 9 of the paper, it requires summing covariance functions that produces matrices of different sizes. Since we are on this, I have a follow up question. PyMC3 provides the function ExpQuad(input_dim,lengthscales) so that it is easy to model the exponentiated quadratic covariance function. My question is, in higher dimensions of lengthscales how is the covariance computed. Is it
cov(x,x’) = exp(-\sum_{k=1}^p \frac{(x-x’)^2}{2\ell[k]^2})
where p is the number of dimensions of lengthscales