Extracting Maximum A Posteriori (MAP) estimates from MC samples I notice in some academic papers the use of MAP estimators in Bayesian settings. The samples are generated through some MC method, but the MAP calculation from these samples is never elaborated on. 
Are there any standard ways of extracting MAP estimates from continuous samples? I feel like I am missing something - since googling doesn't help either. 
 A: [The following is copied from earlier posts on my blog.]
I have never found MAP estimators very appealing for many reasons, one being indeed that the MAP estimator cannot correctly be expressed as the solution to a minimisation problem. I also find the point-wise nature of the estimator quite a drawback: the estimator is only associated with a local property of the posterior density, not with a global property of the posterior distribution. This is in particular striking when considering the MAP estimates for two different parametrisations. The estimates often are quite different, just due to the Jacobian in the change of parameterisation. For instance, the MAP of the usual normal mean $\mu$ under a flat prior is $x$, for instance x=2, but if one use a logit parameterisation instead
$$
\mu = \log \eta/(1-\eta)
$$
the MAP in $\eta$ can be quite distinct from $1/(1+\exp-x)$, for instance leading to $\mu=3$ when $x=2$… Another bad feature is the difference between the marginal MAP and the joint MAP estimates. This is not to state that the MAP cannot be optimal in any sense, as I suspect it could be admissible as a limit of Bayes estimates (under a sequence of loss functions). 
Here are the details for the normal example. I am using a flat prior on $\mu$ when $x\sim\mathcal{N}(\mu,1)$. The MAP estimator of $\mu$ is then $\hat\mu=x$. If I consider the change of variable $\mu=\text{logit}(\eta)$, the posterior distribution on $\eta$ is
$$
\pi(\eta|x) = \exp[ -(\text{logit}(\eta)-x)^2/2 ] / \sqrt{2\pi} \eta (1-\eta)
$$
and the MAP in $\eta$ is then obtained numerically. For instance, the R code
f=function(x,mea) dnorm(log(x/(1-x)),mean=mea)/(x*(1-x))
g=function(x){ a=optimise(f,int=c(0,1),maximum=TRUE,mea=x)$max;log(a/(1-a))}
plot(seq(0,4,.01),apply(as.matrix(seq(0,4,.01)),1,g),type="l")
abline(a=0,b=1,col="tomato2",lwd=2)

shows the divergence between the MAP estimator \hat\mu and the reverse transform of the MAP estimator $\hat\eta$ of the transform… The second estimator is asymptotically (in $x$) equivalent to $x+1$.
An example I like very much in The Bayesian Choice is Example 4.1.2, when observing $x\sim\text{Cauchy}(\theta,1)$ with a double exponential prior on $\theta\sim\exp\{-|\theta|\}/2$. The MAP is then always $\hat\theta=0$!
The dependence of the MAP estimator on the dominating measure is also studied in a BA paper by Pierre Druihlet and Jean-Michel Marin, who propose a solution that relies on Jeffreys’ prior as the reference measure.

An interesting paper by Burger and Lucka compares MAP and posterior mean, even though I do not share their concern that we should pick between those two estimators (only or at all), since what matters is the posterior distribution and the use one makes of it. I thus disagree there is any kind of a "debate concerning the choice of point estimates”. If Bayesian inference reduces to producing a point estimate, this is a regularisation technique and the Bayesian interpretation is both incidental and superfluous.
Maybe the most interesting result in the paper is that the MAP is expressed as a proper Bayes estimator! I was under the opposite impression, mostly because the folklore (and even The Bayesian Core)  have it that it corresponds to a 0-1 loss function does not hold for continuous parameter spaces and also because it seems to conflict with the results of Druihlet and Marin (BA, 2007), who point out that the MAP ultimately depends on the choice of the dominating measure. (Even though the Lebesgue measure is implicitly chosen as the default.) The authors of this arXived paper start with a distance based on the prior; called the Bregman distance. Which may be the quadratic or the entropy distance depending on the prior. Defining a loss function that is a mix of this Bregman distance and of the quadratic distance
$$
||K(\hat u-u)||^2+2D_\pi(\hat u,u)
$$
produces the MAP as the Bayes estimator. So where did the dominating measure go? In fact, nowhere: both the loss function and the resulting estimator are clearly dependent on the choice of the dominating measure… (The loss depends on the prior but this is not a drawback per se!)
