Perhaps the simplest numerical approximation of the marginal likelihood is to
draw samples of a model's parameters from their respective prior distributions.
This turns the intractable integral into a sum of the samples' likelihoods.
Because the prior weight of each sample is one in this case, the marginal
likelihood can be approximated by simply calculating the average likelihood of
the prior samples.
Alternatively, if we have a sample of the parameters from the posterior
distribution—like one obtained from a “standard” Bayesian phylogenetic
analysis via MCMC—we can again use summation to approximate the integral.
In this case, the weight of each sample is the ratio of the prior density to
the posterior density.
As a result, the sum simplifies to the harmonic mean (HM) of the likelihoods
from the posterior sample [1]}.
Both of these techniques can be thought of as importance-sampling integral
approximations.
Whereas both provide unbiased estimates of the marginal likelihood in theory,
they can suffer from very large Monte Carlo error due to
the fact that the prior and posterior are often very divergent,
with the latter usually much more peaked than the former due
to the strong influence of the likelihood.
A finite sample from the prior will often yield an underestimate of the
marginal likelihood, because the region of parameter space with high likelihood
is likely to be missed.
In comparison, a finite sample from the posterior will almost always lead to an
overestimate [2]}, [3]}, [4]}, because it will contain too
few samples outside of the region of high likelihood, where the prior weight
“penalizes” the average likelihood.
However, [5]} showed that for trees with 3–6 tips and relatively
simple models, the average likelihood of a very large sample from the prior
(30–50 billion samples) can yield accurate estimates of the marginal
likelihood.
