I understand. Thank you.
And the mode of the marginals is not necessarily the MAP of the joint posterior, right? From what I see, the find_MAP function in PyMC3 gives values that are close enough to the mode of each marginal, but not the same (sometimes quite different actually).
I have seen several examples (e.g., the hierarchical modelling example in the PyMC3 documentation) in which the means of the marginals are taken as reference values for the parameters in the regression model:
In the figure, the black dash line is created by:
bp = varying_intercept_trace[a].mean(axis=0)
mp = varying_intercept_trace[b].mean()
I am then wondering whether the line created by taking the mean of each parameter is a reasonable summary of the model.