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