4

I used a sampling method to fit a model with three parameters to data, by supplying the likelihood function and priors. (I'm using JAGS but I think this applies to any method). I obtain triplets of the parameter values, for each sample (N=10000). For each sample I also have the posterior and likelihood.

I want to get the joint MAP estimate. I know how to get the marginal MAPs for each of the three individual parameters. However I am not sure how to find the joint - i.e. the triplet that maximises the posterior -- which is quite different in my case.

I could, I suppose, just find the single sample with the largest posterior, and take those parameters. But I thought there would be a way of using the samples together to find this maximum point. Maybe I would have to fit a hypersurface somehow and find where its gradient is zero? Or some form of interpolation. I assume this is a common use case but could not find specifics online (perhaps I don't know what it is called). Are there methods and software to do this automatically eg in Matlab or R?

Robin Ryder
  • 1,787
  • 1
  • 11
  • 16
Sanjay Manohar
  • 906
  • 9
  • 16

1 Answers1

3

You are right that just taking the 3 marginal MAPs would be incorrect.

You can take a kernel density estimate of your posterior sample, and then find the argmax of that function. The function kde of R package ks handles up to 6 dimensions. Beware that the accuracy of kernel density estimates decreases rapidly when the dimension increases: you might need a larger posterior sample, and should at least check the sensitivity of your method (e.g. by repeating it several times and checking how the result varies).

As you suggest, another option is to take the sample with highest posterior density. As $N\to\infty$, both methods will agree.

A word of caution: the MAP is not necessarily the best estimator. If you need a point estimate, the posterior mean (or sometimes median) is usually preferred, since the MAP is not associated with a loss function; see e.g. this answer for details.

Robin Ryder
  • 1,787
  • 1
  • 11
  • 16
  • Thanks - this is helpful - I am surprised there is no 'out-of-the-box' solution that automatically does this for e.g. JAGS or MCMC packages. I will use `kde`. But after reading the answers to the question you mentioned, I was still in the dark as to how you can get a posterior mean when model parameters are multivariate - how can I minimise $\int{(\theta-\hat{\theta})^2 p(\theta | X) d\theta }$ when $\theta$ is a vector? I assume you can't just add up $(\theta_1-\hat{\theta_1})^2 + (\theta_2-\hat{\theta_2})^2$ to get the cost, as they're not necessarily normalised. – Sanjay Manohar Sep 26 '18 at 00:21
  • 1
    Or is this related to what @Xi'an refers to in his (difficult) answer as the "choice of the dominating measure"? – Sanjay Manohar Sep 26 '18 at 00:29