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).