# Find joint maximum of a sampled density

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?

As you suggest, another option is to take the sample with highest posterior density. As $N\to\infty$, both methods will agree.
• 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 0:21