You can do multi-dimensional GPs quite easily: Mean and Covariance Functions — PyMC3 3.11.4 documentation