Multiple Gaussian processes without for loop?

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