Multiple (uncertain) function observations of the same Gaussian process

That looks fantastic!! I think it would be awesome if you could submit this as a PR! I think you can put this into pymc3/distributions/multivariate.py. Yes, everything looks perfect in your notebook. Githubs a better place to go over code, so we can go over this more there. I’ve been wanting to add this for a while, but haven’t gotten around to it. With a MatrixNormal distribution, it should be pretty straightforward to add a few more GP implementations, maybe called Marginal2D and Latent2D, including the coregional model. Super awesome!