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?