The Mauna Loa example uses find_MAP(), but I think since that example was written, a warning was added to pymc3 that "find_MAP should not be used to initialize the NUTS sampler. Simply call pymc3.sample() and it will automatically initialize NUTS in a better way."
Is find_MAP() still the recommended way to fit a gaussian process model? Or should I always be using NUTS instead?
The warning says that you shouldn’t use find_MAP to initialize the chain. It will automatically initialize it smartly if you are using pm.sample with NUTS. But they are two independent inference methods altogether. You can either use find_MAP to infer the parameters of the model or use MCMC methods according to the size of your data and your needs.
Coming to your question, Yes, you can use find_MAP to fit a GP model. In fact, that’s what most of the papers and practical methods tell you to do.
OTOH, PyMC3 provides nice MCMC methods that can also be used on a GP model. Again, it’s upto you what you choose.
The main difference between these methods is that find_MAP may get stuck at a local minima, plateau (flat surfaces) or other non-convex high dimensional optimization problems while NUTS performs very good on high dimensional spaces and provides unbiased estimates of model parameters in a long run (with a downside of being slow on larger datasets).