Rrerku

I have a likelihood function $mathcalL(d | theta)$ for the probability of my data $d$ given some model parameters $theta in mathbfR^N$, which I would like to estimate. Assuming flat priors on the parameters, the likelihood is proportional to the posterior probability. I use an MCMC method to sample this probability.

Looking at the resulting converged chain, I find that the maximum likelihood parameters are not consistent with the posterior distributions. For example, the marginalized posterior probability distribution for one of the parameters might be $theta_0 sim N(mu=0, sigma^2=1)$, while the value of $theta_0$ at the maximum likelihood point is $theta_0^ML approx 4$, essentially being almost the maximum value of $theta_0$ traversed by the MCMC sampler.

This is an illustrative example, not my actual results. The real distributions are far more complicated, but some of the ML parameters have similarly unlikely p-values in their respective posterior distributions. Note that some of my parameters are bounded (e.g. $0 leq theta_1 leq 1$); within the bounds, the priors are always uniform.

My questions are:

1. Is such a deviation a problem per se? Obviously I do not expect the ML parameters to exactly coincide which the maxima of each of their marginalized posterior distributions, but intuitively it feels like they should also not be found deep in the tails. Does this deviation automatically invalidate my results?




2. Whether this is necessarily problematic or not, could it be symptomatic of specific pathologies at some stage of the data analysis? For example, is it possible to make any general statement about whether such a deviation could be induced by an improperly converged chain, an incorrect model, or excessively tight bounds on the parameters?

With flat priors, the posterior is identical to the likelihood up to a constant. Thus

1. MLE (estimated with an optimizer) should be identical to the MAP (maximum a posteriori value = multivariate mode of the posterior, estimated with MCMC). If you don't get the same value, you have a problem with your sampler or optimiser.




2. For complex models, it is very common the marginal modes are different from the MAP. This happens, for example, if correlations between parameters are nonlinear. This is perfectly fine, but marginal modes should therefore not be interpreted as the points of highest posterior density, and not be compared to the MLE.




3. In your specific case, however, I suspect that the posterior runs against the prior boundary. In this case, the posterior will be strongly asymmetric, and it doesn't make sense to interpret it in terms of mean, sd. There is no principle problem with this situation, but in practice it often hints towards model misspecification, or poorly chosen priors.

Some possible generic explanations for this perceived discrepancy, assuming of course there is no issue with code or likelihood definition or MCMC implementation or number of MCMC iterations or convergence of the likelihood maximiser (thanks, Jacob Socolar):

1. in large dimensions $N$, the posterior does not concentrate on the
   maximum but something of a distance of order $sqrtN$ from the
   mode, meaning that the largest values of the likelihood function
   encountered by an MCMC sampler are often quite below the value of
   the likelihood at its maximum. For instance, if the posterior is $theta|mathbf xsimmathcal N_N(0,I_n)$, $theta$ is at least at a distance $N-2sqrt2N$ from the mode, $0$.




2. While the MAP and the MLE are indeed confounded under a flat prior, the
   marginal densities of the different parameters of the model may have (marginal) modes
   that are far away from the corresponding MLEs (i.e., MAPs).




3. The MAP is a position
   in the parameter space where the posterior density is highest but
   this does not convey any indication of posterior weight or volume
   for neighbourhoods of the MAP. A very thin spike carries no posterior weight. This is also the reason why MCMC exploration of a posterior may face difficulties in identifying the posterior mode.




4. The fact that most parameters are bounded may lead to some
   components of the MAP=MLE occurring at a boundary.

See, e.g., Druihlet and Marin (2007) for arguments on the un-Bayesian nature of MAP estimators.

As an example of point 1 above, here is a short R code

N=100
T=1e4
lik=dis=rep(0,T)
mu=rmvnorm(1,mean=rep(0,N))
xobs=rmvnorm(1,mean=rep(0,N))
lik[1]=dmvnorm(xobs,mu,log=TRUE)
dis[1]=(xobs-mu)%*%t(xobs-mu)
for (t in 2:T)
 prop=rmvnorm(1,mean=mu,sigma=diag(1/N,N))
 proike=dmvnorm(xobs,prop,log=TRUE)
 if (log(runif(1))<proike-lik[t-1])
 mu=prop;lik[t]=proike
 elselik[t]=lik[t-1]
 dis[t]=(xobs-mu)%*%t(xobs-mu)

which mimics a random-walk Metropolis-Hastings sequence in dimension N=100. The value of the log-likelihood at the MAP is -91.89, but the visited likelihoods never come close:

> range(lik)
[1] -183.9515 -126.6924

which is explained by the fact that the sequence never comes near the observation:

> range(dis)
[1] 69.59714 184.11525