Re: Highest Posterior Density

Reef Fish wrote:
AB wrote:
how can I find an HPD using a statistical software (e.g. Mathematica)?

By HPD do you mean VALUE of the variable that corresponds to
the highest point of the posterior density in a problem using a
Bayesian method of inference?

You need to be able to supply the Bayesian ingredients of a prior
distribution for your parameter, be able to perform the integration to
obtain your posterior distribution of your PARAMETER

Actually, you do not need to perform the integration as the
normalizing constant does not depend on the parameter you
wish to estimate, and hence does not affect the
maximization.

Most statistical
software should have some method for finding the maximum of a
function -- that part should be easy.

Whether the optimization is easy very much depends on the
dimensionality of the parameter space on which the density
is defined, and on the complexity of the density itself. It
can be an almost impossible task to find a global (as
opposed to local) maximum in a high-dimensional space,
hence the large amount of research devoted to this topic.

There is also a conceptual difficulty with MAP (HPD)
estimation. While Bayes' theorem supplies you with a
measure on the parameter space, in order to perform the
maximization you need a function, i.e. a density. You
therefore need an underlying measure with respect to which
you can define the density of the posterior. The problem is
that different underlying measures give rise to different
estimates. (A similar thing happens with MMSE estimates
too.) Since the most common procedure is simply to drop the
'd\theta' symbol in the density (\theta is the parameter),
i.e. to use 'Lebesgue measure' in the current coordinate
system on the parameter space as the underlying density,
different choices of coordinates, e.g. using \theta^{2}
instead of \theta for a positive parameter, lead to
different estimates. In order to get around this, it is
necessary to choose an underlying density that is invariant
to the choice of coordinates on the parameter space.

The most sensible choice is to use a Riemannian metric to
define the underlying measure, since a metric generates
coherent notions of volume (MAP estimates) and distance
(MMSE estimates). The question is then, how to choose the
metric? If we are just given a density on some space, there
seems to be no good way to do this. However, in the context
of Bayes' theorem, there are good arguments for choosing
the Fisher information metric constructed from the
likelihood, as this is the only choice that does not
introduce extra information. In the case that the prior is
Jeffreys' prior (the volume element of the Fisher
information metric), this leads to maximum likelihood as
the estimation method.

To return to the OP's question: the first thing to know is
that there is no algorithm that will guarantee to find the
global maximum for an arbitrary function. The difficulty of
your problem largely depends first, on the dimensionality
of your space, and second, on the complexity of your
density.

Matlab has an fminsearch function that looks for an
unconstrained minimum of a multivariable function, but may
only return a local minimum. (Obviously finding a maximum
is the same as finding a minimum of the negative of the
function.) The Matlab Optimization toolbox provides
fminunc, which addresses the same problem, as well as other
optimization tools.

Mathematica has the functions NMinimize and FindMinimum,
with similar goals.

There is no real substitute, though, for analysing your
function and trying to find out as much as you can about
it, and then writing your own code. Many optimization
algorithms are not hard to implement (gradient descent,
simulated annealing), and can often provide good results
with all your assumptions known, because you wrote the
code.

illywhacker;

.