Evaluating Posterior Distributions by Selectively Breeding Prior Samples

We may say more about the rate of convergence of $\widehat{\Pi}_{1}(A)$ to $\Pi_{1}(A)$.

Proof: Defining $N_{A}$ and $D$ as in the proof of Proposition 1, the delta method gives an asymptotic variance:

Fill in each piece, using the fact that an indicator function is equal to its own square:

Put these together, and pull out the common factor of $n$:

Since $\Pi_{1}(A)=\Pi_{0}(fA)/\Pi_{0}(f)$, after some algebra,

as was claimed. $\Box$

Remarks, or sanity checks: (i) $1\geq 1-2\Pi_{1}(A)\geq-1$, and $\Pi_{0}(f^{2})\geq\Pi_{0}(f^{2}A)$, so the numerator is always $\geq 0$. (ii) As $A\rightarrow\emptyset$ or $A\rightarrow\Theta$, the variance $\rightarrow 0$. (iii) As $\Pi_{1}(A)\rightarrow 1$, the variance tends to $\Pi_{0}(f^{2}A^{c})/\Pi_{0}^{2}(f)$, while as $\Pi_{1}(A)\rightarrow 0$, the numerator approaches $\Pi_{0}(f^{2}A)/\Pi_{0}^{2}(f)$. (iv) An oracle which could sample directly from $\Pi_{1}$ would have an asymptotic variance proportional to $\Pi_{1}(A)-\Pi_{1}^{2}(A)$, so the algorithm is not as efficient as the oracle.

Proof: Since $\Pi_{0}(f^{2}A)(1-2\Pi_{1}(A))\leq\Pi_{0}(f^{2})$, the previous proposition gives us an asymptotic variance of at most $2\Pi_{0}(f^{2})/\Pi_{0}^{2}(f)$ for any set. Now, the Rényi divergence of order $\alpha\geq 0$ is defined (for $\alpha\neq 0,1,\infty$) to be (van Erven and Harremoës, 2014)

Pick $\alpha=2$, and suppose without loss of generality that $P$ and $Q$ have densities $p$ and $q$ with respect to a common measure $M$, such as $(P+Q)/2$, so $\frac{dP}{dQ}(\theta)=p(\theta)/q(\theta)$. Then

Set $P=\Pi_{1}$ and $Q=\Pi_{0}$. The density ratio is

so

Undoing the logarithm proves the result. $\Box$

The preceding results deal with the behavior as the number of Monte Carlo samples $n$ grows; to sum up (Corollary 3), with a large number of draws from the prior, the accuracy of the approximation is controlled by $\Pi_{0}(f^{2})/\Pi_{0}^{2}(f)=1+\mathrm{Var}\left[f\right]/\Pi_{0}^{2}(f)$. There is also a second asymptotic regime of interest, which is when the observation $x$ becomes highly informative. It would be somewhat perverse if the approximation broke down just when the observation became most useful.

These fears can be calmed, at least in some situations, by asymptotic expansions. Suppose that $\Theta=\mathbb{R}^{d}$ and $\Pi_{0}$ has a density $\pi_{0}$ with respect to Lebesgue measure, so $\Pi_{0}(g)=\int{g(\theta)\pi_{0}(\theta)d\theta}$. In particular, write

and

where $t$ gauges, in some relevant sense, how informative the observation $x$ is, with the normalized negative log-likelihood $\ell(\theta)$ converging to a constant-magnitude function as $t\rightarrow\infty$. ($t$ may represent the number of measurements taken, the length of a time series, the precision of a measuring instrument, etc.). Assume that $\ell(\theta)$ has a unique interior minimum at the MLE $\hat{\theta}$. We may then appeal to Laplace approximation to evaluate these integrals: by Eq. 2.3 in Tierney et al. (1989):

where $\Sigma$ is the inverse Hessian matrix of $\ell$ at $\hat{\theta}$. The precise expression for $A_{1}$ in terms of the derivatives of $\pi_{1}$ and $\ell$ is known (Wojdylo, 2006; Nemes, 2013) but not, for the present purposes, illuminating. Appealing to Laplace approximation a second time,

Thus

The upshot of this is calculation²²2Marginally further simplified, when $d=1$, by the fact that then $A_{2}=A_{1}/2$; whether this is true more generally I do not know. is that as $t\rightarrow\infty$, that is, as the data becomes increasingly informative, the ratio $\Pi_{0}(f^{2})/\Pi_{0}^{2}(f)$ tends to a limit proportional to the prior density at the MLE, with vanishing corrections. Heuristically, as $t\rightarrow\infty$, both $\Pi_{0}(f^{2})$ and $\Pi_{0}(f)$ become integrals dominated by their slowest-decaying integrand. The rate of exponential decay is twice as fast for $f^{2}$ as for $f$, so $\Pi_{0}(f^{2})/\Pi_{0}^{2}(f)$ tends to a $t$-independent limit.³³3Bengtsson et al. (2008) analyzes a related problem, in the context of particle filtering, and concludes that increasing $t$, as from a growing number of observations, leads to “collapse” of the filter, in the sense that $\Pi_{0}(\max_{i}{f(\theta_{i})/D})\rightarrow 1$, unless $n$ grows exponentially in some (fractional) power of $t$. They do not, however, analyze the quality of the approximation to the posterior distribution, which is my concern here. The limit they consider is one where the likelihood is extremely sharply peaked over a region of very small prior probability mass, so it is indeed unsurprising that if the number of samples is small, only a few fall within it, and many samples would be needed to accurately sketch a very rapidly-changing posterior. The real source of a difficulty, if there is one, is $\pi_{0}(\hat{\theta})$ being small.

It’s worth remarking that (in the present notation) $\log{\pi_{1}(\hat{\theta})/\pi_{0}(\hat{\theta})}$ was proposed by Haldane (1954, pp. 483, 486) as a measure of the intensity of natural selection on a population.

It would be nice if there were finite-sample concentration results to go along with these asymptotics. One approach to this would be to use results developed for weighting training samples to cope with distributional shift. The most nearly applicable result is that of Cortes et al. (2010, Theorem 7 and Corollary 1), which runs (in the present notation) as follows. If the samples $\theta_{1},\ldots\theta_{n}$ are drawn iidly from $Q$, and $dP/dQ=w$, then over any function class $H$ bounded between 0 and 1,

where $G_{H}$ is the growth function of $H$, the number of ways a set of $n$ points can be dichotomized by sets from $H$ (Anthony and Bartlett, 1999). As in this result, we sample from one distribution and want to compute integrals under another probability measure, so we would like to set $Q=\Pi_{0}$ and $P=\Pi_{1}$. (This would make the denominator on the left hand side $\exp{\left\{D_{2}(\Pi_{1}\|\Pi_{0})/2\right\}}$.) Standing in the way is the annoyance that the conversion weights $w=f(\theta)/\Pi_{0}(f)$ are unavailable, because instead of $\Pi_{0}(f)$ we only have a Monte Carlo approximation $D$. Introduce a new random measure $R$, defined through $dR/d\Pi_{0}=f/D$. ($R$ is not, generally, a properly normalized probability measure.) The theorem of Cortes et al. then bounds the deviation of $\widehat{\Pi}_{1}(h)$ from $R(h)$. Since $d\Pi_{1}/dR=D/\Pi_{0}(f)$, a finite-sample deviation inequality for $D$ will yield a bound on the relative deviation of $\widehat{\Pi}_{1}(h)$ from $\Pi_{1}(h)$. Controlling the convergence of $D$ to $\Pi_{0}(f)$ will need assumptions on the distribution (under the prior) of likelihoods, which seem too model-specific to pursue here.

For any measurable $A\subseteq\Theta$, $\widehat{\Pi}_{1}(A)=\Pi_{1}(A)+O_{P}(1/\sqrt{n})$. Thus Algorithm 1 provides a rapidly-converging way of estimating posterior probabilities. The error of approximation is controlled by $D_{2}(\Pi_{1}\|\Pi_{0})$, or more concretely by the ratio $\Pi_{0}(f^{2})/\Pi_{0}^{2}(f)=1+\mathrm{Var}\left[f\right]/\Pi_{0}^{2}(f)$. This is to say, the more (prior) variance there is in the likelihood, the worse this will work. (The dimensionality of $\Theta$ is not, however, directly relevant.) Algorithms 2 and 3 simply introduce additional approximation error, which however can be made as small as desired either by letting $c\rightarrow\infty$ or $m\rightarrow\infty$.

These algorithms are all closely related to particle filters (Doucet et al., 2001), which however usually aim at approximating a sequence of probability measures which themselves evolve randomly, and where convergence results often rely on the ergodicity of the underlying process. Here, however, we have only one (possibly high-dimensional) observation, and only one posterior distribution to match.

Artificial intelligence employs many population-based optimization techniques to explore fitness landscapes, by copying more successful solutions in proportion to their fitness, most notably genetic algorithms (Holland, 1975; Mitchell, 1996). All of these, however, have a vital role for operations which introduce new members of the population, corresponding to mutation, cross-over, sex, etc., which are necessarily lacking in Bayesian updating⁴⁴4From a strictly Bayesian point of view, innovation is merely an invitation to be Dutch Booked, and so “irrational” and/or “incoherent” (Gelman and Shalizi, 2013)., and avoided in the algorithms set forth above.

Beyond these general points of inspiration, Rubin (1987) proposed what is basically Algorithm 3 as a way of doing multiple imputation, as opposed to sampling from the posterior of the parameters. Cappé et al. (2004) is an even more precise anticipation, though their use of sequential importance sampling means they aren’t just sampling from the prior, with attendant loss of simplicity and analytical tractability.

Best practices in Markov chain Monte Carlo counsel discarding all the samples from an initial burn-in period, running diagnostic tests to check converge to the equilibrium, etc. (Admittedly, some doubt the value of these steps: see Geyer 1992.) None of this is needed or even applicable here: all that’s needed is for $n$ to be large enough that $\widehat{\Pi}_{1}$ and/or $\tilde{\Pi}_{1}$ is a good approximation to the posterior distribution $\Pi_{1}$. Moreover, there is no need to pick or adjust proposal distributions. Indeed, there aren’t really any control settings that need to be tuned, since it’s straightforward that $n$, $c$ and $m$ should all be as large as possible.

Implementing this scheme requires only the ability to draw IID samples from the prior $\Pi_{0}$, and to compute the likelihood function $f$. The former tends to be straightforward, at least when the prior is expressed hierarchically. Actually, examination of the proofs will confirm that we do not even need to calculate $f$; it would be enough to calculate something proportional to the likelihood, $r(x)f(\theta)$, where the factor of proportionality $r(x)>0$ and is strictly constant in $\theta$. That being said, calculating the likelihood can be difficult for complex models, but this is outside the scope here⁵⁵5If the difficulty in calculating $f(\theta;x)$ comes from integrating out latent variables, say $Y\in\Upsilon$, so that $f(\theta;x)=\int_{\Upsilon}{p(x|y;\theta)p(y|\theta)dy}$, there is a straightforward work-around if $p(x|y;\theta)$ is computationally tractable and $p(y|\theta)$ is easy to simulate from. Formally expand the parameter space to $\Theta^{\prime}=\Theta\times\Upsilon$, drawing $\theta_{i}$ from $\Pi_{0}$ and $Y_{i}$ from $p(y|\theta_{i})$. Then use $p(x|y;\theta)$ as the likelihood, and apply the posterior approximation only to sets of the form $A\times\Upsilon$..

One special advantage of this scheme is that it lends itself readily to parallelization, which MCMC does not. Drawing iidly from the prior is “embarrassingly” parallel. Calculating $f(\theta_{i})$ is also embarrassingly parallel, requiring no communication between the different $\theta_{i}$, as is making copies of them. Every step in algorithm 2 is embarrassingly parallel, though with large $c$ we might need a lot of memory in each node. Algorithms 1 and 3 require us to compute $\sum_{i=1}^{n}{f(\theta_{i})}$, which does require communication, but finding the sum of $n$ numbers is among the first parallelized algorithms taught to beginners (Burch, 2009, §2.1), so this might not be much of a drawback.

The great advantage of all three algorithms is that the only distribution they ever need to draw from is the prior, and priors are typically constructed to make that easy. The corresponding disadvantage is that they can spend a lot of time evaluating the likelihood of parameter values which have in fact very low likelihood, and so end up making a small contribution to the posterior. MCMC, by contrast, preferentially spends most of its time in regions of high posterior density.

Three observations may, however, lessen this contrast. First, while the output of MCMC follows the posterior distribution, its proposals do not, and the likelihood must be evaluated once for each proposal. Thus it is not clear that, in general, MCMC has a computational advantage on this score. Second, while I have not explored it in this brief note, it should be possible to short-cut the evaluation of the likelihood of very hopeless parameter points, as in the “go with the winners” of Grassberger (2002)⁶⁶6Grassberger traces the idea of randomly eliminating bad parameter values, and making more copies of the good ones, to Kahn (1956). The latter, in turn, gave this trick the arresting name of “Russian Roulette and splitting”, and attributed it to an unpublished paper by John von Neumann and Stanislaw Ulam., or the “austerity” of Korattikara et al. (2014). Third, the problematic case is where most of the posterior probability mass is concentrated in a region of very low prior mass. Numerical evidence from trying out more complicated situations than the Gaussian prior-and-likelihood suggests that this can be a serious issue, which is why I do not advocate these algorithms, in their present form, as a replacement for MCMC (note that the non-Markov chain “population Monte Carlo” of Cappé et al. (2004) does not suffer from this drawback). It is, however, at least arguable that what goes wrong in these situation isn’t that these algorithms converge slowly, but that the priors are bad, and should be replaced by less misleading ones (Gelman and Shalizi, 2013).