Prior-predictive value from fast-growth simulations: Error analysis and bias estimation

Variants of fluctuation theorems recently discovered in the statistical mechanics of nonequilibrium processes may be used for the efficient determination of high-dimensional integrals as typically occurring in Bayesian data analysis. In particular for multimodal distributions, Monte Carlo procedures not relying on perfect equilibration are advantageous. We provide a comprehensive statistical error analysis for the determination of the prior-predictive value (the evidence) in a Bayes problem, building on a variant of the Jarzynski equation. Special care is devoted to the characterization of the bias intrinsic to the method and statistical errors arising from exponential averages. We also discuss the determination of averages over multimodal posterior distributions with the help of a consequence of the Crooks relation. All our findings are verified by extensive numerical simulations of two model systems with bimodal likelihoods.

Figure 1

Verifying that our error analysis is consistent with the result (50) obtained in [11, 28] based on the central limit theorem for the random variable $e^R$. The total number of Markov chains, and consequently of $R$ values, is $N=6\times {10}^7$.

Figure 2

As the number of $R$ values, $N$, gets larger, the confidence interval (28) gets smaller. Nevertheless, for the data set examined here, $D\left(N\right)$ stays within ${\widehat{D}}_{+}\left(N\right)$ and ${\widehat{D}}_{-}\left(N\right)$. The value of $D\left(N\right)$ follows from (27) using the exact value of $p\left(d\right|\mathfrak{M})$ given by (47); the estimated confidence limits ${\widehat{D}}_{\pm }\left(N\right)$ are derived from (41).

Figure 3

Estimation of the log-evidence for the bimodal Gaussian model (44) using the Jarzynski estimator (14) for an increasing number $N$ of $R$ values. The thick line is the analytic result (47), and the symbols use the Jarzynski estimator. The error bars are given by the estimates of the bias, ${\widehat{D}}_{\pm }\left(N\right)={\widehat{B}}_{\pm }(N,N)={\widehat{\alpha }}_{\pm }(N,N)$, from (41).

Figure 4

Verification for the likelihood distribution (53) involving two Cauchy distributions that our error analysis is consistent with the result (50) obtained in [11, 28] based on the central limit theorem for the random variable $e^R$. The total number of Markov chains, and consequently of $R$ values, is $N={10}^8$.

Figure 5

The value of $D\left(N\right)$ defined in (27) together with its confidence interval given by ${\widehat{D}}_{\pm }\left(N\right)$ from (41) as a function of increasing number $N$ of considered samples for $R$ for the likelihood distribution in (53). The determination of $D\left(N\right)$ involves the exact value of $p\left(d\right|\mathfrak{M})$, which follows from using the cumulative Cauchy distribution in (52).

Figure 6

Estimation of the log-evidence for the likelihood distribution (53), which involves two Cauchy distributions, using the Jarzynski estimator (14) for an increasing number $N$ of $R$ values. The thick line is the analytic result obtained from (52); the symbols use the Jarzynski estimator. The error bars are given by the confidence limits ${\widehat{D}}_{\pm }\left(N\right)$ from (41).

Figure 7

Histogram for the average of $x_{\parallel }$ along each Markov chain in our standard Monte Carlo program. Despite their different weights, the peaks of the bimodal posterior distribution are sampled almost equally.

Figure 8

Histogram for the end location $x_{\parallel }\left(T\right)$ of each trajectory in our fast-growth program. Although the Markov chains are nonstationary, they still get trapped around the maxima of the posterior distribution.

Figure 9

Absolute error in the fast-growth estimate of ${\langle x_{\parallel }\rangle }_{\mathrm{post}}$ as the number of trajectories $N$ is varied. The rightmost data point corresponds to (59).