Can I make out-of-sample predictions using posterior means? Or how to use scipy.optimize with pymc predictions?

If you take the mean of the posterior then optimize you will get the wrong answer due to Jensen’s Inequality. You want \mathbb{E}[f(x)], but you are computing f(\mathbb{E}[x]). You can consider putting the optimization into the model itself, as in here. You would need to add gradients for the optimizer using the implicit value theorem for any non-trivial problem.

If the model is essentially linear though, I would be looking for an approximate closed-form solution rather than running an optimizer.

1 Like