Monte Carlo Fusion

TL;DR

This paper introduces Monte Carlo Fusion, a novel theoretical framework and methodology for unifying distributed analyses on shared parameters without approximation errors, applicable to various complex data settings.

Contribution

It provides the first general, approximation-error-free approach for combining distributed inferences, with a focus on theoretical foundations and Monte Carlo interpretations.

Findings

First general approach avoiding approximation errors in inference fusion

Theoretical underpinnings of Monte Carlo Fusion methodology

Potential for tailored, efficient parallel implementations

Abstract

This paper proposes a new theory and methodology to tackle the problem of unifying distributed analyses and inferences on shared parameters from multiple sources, into a single coherent inference. This surprisingly challenging problem arises in many settings (for instance, expert elicitation, multi-view learning, distributed 'big data' problems etc.), but to-date the framework and methodology proposed in this paper (Monte Carlo Fusion) is the first general approach which avoids any form of approximation error in obtaining the unified inference. In this paper we focus on the key theoretical underpinnings of this new methodology, and simple (direct) Monte Carlo interpretations of the theory. There is considerable scope to tailor the theory introduced in this paper to particular application settings (such as the big data setting), construct efficient parallelised schemes, understand the approximation and computational efficiencies of other such unification paradigms, and explore new theoretical and methodological directions.

Tables

Table 1. Table 1: System running times in seconds for simulating 10,000 realisations.

Algorithms	CMC	Algorithm 1	Algorithm 2
Running time	0.05	0.36	4.05

Equations

$$f\left(\mathbf{x}\right)\propto f_1\left(\mathbf{x}\right)\cdots f_C\left(\mathbf{x}\right),$$

$$g({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y})$$

$$h^{bm}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y})$$

$$\mathrm{d}{\mathbf{X}}_t^{(c)}$$

$$A_c^{\text{dl}}\left(\mathbf{x}\right)$$

$${\varphi }_c^{\text{dl}}\left(\mathbf{x}\right)$$

$$\frac{g^{dl}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y})}{h^{bm}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y})}={\left[\frac{\sqrt{C}}{\sqrt{2\pi T}}\right]}^C\times {\rho }^{bm}\times Q^{bm}\times \prod _{c=1}^C{e}^{-T{\Phi }_c^{\text{bm}}},$$

$${\rho }^{bm}:={\rho }^{bm}({\mathbf{x}}^{(1)},\cdots ,{\mathbf{x}}^{(C)})=e^{-\frac{C{\sigma }^2}{2T}},{\sigma }^2=C^{-1}\sum _{c=1}^C\Vert {\mathbf{x}}^{(c)}-\bar{\mathbf{x}}{\Vert }^2,$$

$$Q^{bm}={\mathbf{E}}_{\overline{\mathbb{W}}}\left(E^{bm}\right),$$

$$E^{bm}$$

$$p_{T,c}^{\text{dl}}\left(\mathbf{y}|{\mathbf{x}}^{(c)}\right)$$

$$U_2\le E^{bm}$$

$${\sigma }^2\le -\frac{2T\log U_1}{C}.$$

$$h^{\text{bm}}\left({\mathbf{x}}^{(1)},\cdots ,{\mathbf{x}}^{(C)},\mathbf{y}\right)$$

$$h^{\text{ou}}\left({\mathbf{x}}^{(1)},\cdots ,{\mathbf{x}}^{(C)},\mathbf{y}\right)\propto \left[\prod _{c=1}^C{f}_c\left({\mathbf{x}}^{(c)}\right)\right]\text{etr}\left[-\frac{1}{2}{\left[\mathbf{y}-\tilde{\mathbf{x}}\right]}^{\otimes 2}\mathbf{D}\right]$$

$$\mathbf{D}$$

$${\mathbf{V}}_c$$

$$\tilde{\mathbf{x}}$$

$${\mathbf{m}}_c$$

$$\mathrm{d}{\mathbf{X}}_t^{(c)}=\nabla A_c^{\text{ou}}\left({\mathbf{X}}_t^{(c)}\right)\mathrm{d}t+\mathrm{d}{\mathbf{W}}_t^{(c)},{\mathbf{X}}_0^{(c)}\sim f_c^2\left(\mathbf{x}\right)$$

$$A_c^{\text{ou}}\left(\mathbf{x}\right)=-\frac{({\hat{\mathbf{\mu }}}_c-\mathbf{x}{)}^{tr}{\hat{\mathbf{\Lambda }}}_c({\hat{\mathbf{\mu }}}_c-\mathbf{x})}{2}.$$

$${\varphi }_c^{\text{ou}}\left(\mathbf{x}\right)$$

$${\varphi }_c^{\text{dl}}\left(\mathbf{x}\right)-{\varphi }_c^{\text{ou}}\left(\mathbf{x}\right)\ge {\Phi }_c^{\text{ou}}.$$

$${\rho }^{ou}$$

$${\mathbf{M}}_{1,c}$$

$${\mathbf{M}}_{2,c}$$

$$\mathbf{H}$$

$$\frac{g^{dl}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y})}{h^{\text{ou}}\left({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(C)},\mathbf{y}\right)}\propto {\rho }^{ou}({\mathbf{x}}^{(1)},\cdots ,{\mathbf{x}}^{(C)})\times Q^{ou}\times \prod _{c=1}^C{e}^{-T{\Phi }_c^{\text{ou}}}$$

$$Q^{ou}={\mathbf{E}}_{\overline{\mathbb{O}}}\left(E^{ou}\right)$$

$$E^{ou}:=\prod _{c=1}^C\left[\exp \left\{-{\int }_0^T\left({\varphi }_c^{\text{dl}}\left({\mathbf{x}}_t^{(c)}\right)-{\varphi }_c^{\text{ou}}\left({\mathbf{x}}_t^{(c)}\right)-{\Phi }_c^{\text{ou}}\right)\mathrm{d}t\right\}\right].$$

$$U_2\le \exp \left[-\sum _{c=1}^C{\int }_0^T\left({\varphi }_c^{\text{dl}}\left({\mathbf{x}}_t^{(c)}\right)-{\varphi }_c^{\text{ou}}\left({\mathbf{x}}_t^{(c)}\right)-{\Phi }_c^{\text{ou}}\right)\mathrm{d}t\right]$$

$${\tilde{h}}^{ou}=h^{ou}\cdot {\rho }^{ou}=\left[\prod _{c=1}^C{f}_c\left({\mathbf{x}}^{(c)}\right)\right]\text{etr}\left[-\frac{1}{2}{\left[\mathbf{y}-\tilde{\mathbf{x}}\right]}^{\otimes 2}\mathbf{D}\right]\cdot {\rho }^{ou}({\mathbf{x}}^{(1)},\cdots ,{\mathbf{x}}^{(C)})$$

Full text

\authornames

Dai, Pollock, Roberts

\authorone

[University of Essex]Hongsheng Dai \authortwo[University of Warwick]Murray Pollock \authorthree[University of Warwick]Gareth Roberts \addressoneDepartment of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester, CO4 3SQ \emailonehdaia@essex.ac.uk \addresstwoDepartment of Statistics, University of Warwick, Gibbet Hill Road, Coventry, CV4 7ES \emailtwom.pollock@warwick.ac.uk; gareth.o.roberts@warwick.ac.uk

\ams

65C05;65C6062C10;65C30

1 Introduction

A common problem arising in statistical inference is the need to unify distributed analyses and inferences on shared parameters from multiple sources, into a single coherent inference. This unification (or what we term ‘fusion’) problem can arise either explicitly due to the nature of a particular application, or artificially as a consequence of the approach a practitioner takes to tackling an application.

Typically there will exist no closed form analytical approach to unifying distributed inferences, and so we focus on a Monte Carlo approach. Stated generally, in this paper we are interested in sampling (without error) the following $d$-dimensional (fusion) target density,

[TABLE]

where each $f_{c}\!\left(\boldsymbol{x}\right)$ ($c\in\{1,\ldots{},C\}$) is a density (up to a multiplicative constant) representing one of the $C$ distributed inferences we wish to unify. Each $f_{c}\!\left(\boldsymbol{x}\right)$ (which we term a sub-posterior) may in practice itself be represented by a Monte Carlo sample ($\tilde{f}_{c}\!\left(\boldsymbol{x}\right)$), and in this paper we assume we are able to sample (directly and exactly) from each $f_{c}\!\left(\boldsymbol{x}\right)$.

Specific examples of this problem arising naturally in an application include: expert elicitation (Berger (1980), Genest and Zidek (1986)), in which the (distributional) views of multiple experts on a topic (or set of parameters) have to be pooled into a single view before a decision-maker can make an informed decision; and, multi-view learning (Zhao et al. (2017), Li et al. (2015)) and meta-analysis (Fleiss (1993), Smith et al. (1995)), in which an interpretation could be that we are synthesising multiple inferences on a particular parameter set (computed on datasets which may or may not be of the same type), but the underlying raw data is not directly available for the unified inference. Obtaining the raw data may itself be an insoluble problem due to reasons including the nature of the original publication, data confidentiality, or simply time and storage constraints.

This fusion problem also arises artificially in a number of settings, particularly within modern statistical methodologies tackling ‘big data’. The computational cost of algorithms such as Metropolis-Hastings, which is an iterative algorithm requiring full access at every iteration to the full dataset, scale poorly with increasing volumes of data unless a modification is found.

One common modification in light of the challenge of big data is to deploy a ‘divide-and-conquer’ approach (or more accurately termed ‘fork-and-join’ approach (Stamatakis and Aberer (2013))). In this setting the full dataset is artificially split into a large number of smaller data sets, inference is then conducted on each smaller data set in isolation, and the resulting inferences are unified (see for instance, Scott et al. (2016), Agarwal and Duchi (2012), Neiswanger et al. (2014), Wang and Dunson (2013), Minsker et al. (2014), Srivastava et al. (2016) and Li et al. (2017)). The rationale for such an approach is that inference on each small data set can be conducted independently, and in parallel, and so one could exploit large clusters of computing cores to greatly reduce the elapsed time to conduct the full inference. The weakness of these approaches is that the fusion of the separately conducted inferences is inexact. It should be noted that divide-and-conquer methodologies will typically have additional constraints due to hardware concerns – such as minimising or removing any communication between computing cores to reduce the effect of latency. We focus in this paper on the general fusion problem, and so do not fully address the problem in the context of big data (to which we return in subsequent work).

The framework and methodology we outline in this paper (Monte Carlo Fusion) for sampling exactly from (1) can be viewed as a simple rejection sampling scheme on an extended space – we develop and sample from efficient proposal densities for (1), the samples from which we retain according to an appropriate acceptance probability. The mathematical complication in this paper is in computing the intractable acceptance probability – which requires the auxiliary simulation of collections of Brownian (or Ornstein-Uhlenbeck) bridges. Our fusion approach provides a principled way to understand the error in existing unification schemes, using a simple linear combination and correction of Monte Carlo samples (analogous to a traditional meta-analysis approach).

The presentation of this paper broadly follows this pedagogy: In Section 2 we present a density on an extended space which admits (1) as a marginal. The remainder of the paper then develops a rejection sampler for the extended density in Section 2. In Section 3 we develop general theory and methodology for sampling (1) based on a collection of independent Brownian bridges. In Section 4 we present a modification of the theory developed in Section 3 using Ornstein-Uhlenbeck bridges, resulting in sampling efficiencies for the particular (common) setting in which the fusion density is believed to be approximately Gaussian. In Section 5 we consider examples of our methodology applied to both light-tailed and heavy-tailed fusion target densities. Finally, in Section 6 we conclude by discussing the exciting new research directions possible using Monte Carlo Fusion. Much of the technical detail in the paper is suppressed for ease of reading, but can be found in the appendices.

2 An extended fusion density

Consider $f\!\left(\boldsymbol{x}\right)$ and $f_{c}\!\left(\boldsymbol{x}\right)$ as described in (1), where $f\!\left(\boldsymbol{x}\right)$ is integrable and we can sample from the density proportional to $f_{c}\!\left(\boldsymbol{x}\right)$.

The following simple observation will form the foundation of our approach: Suppose that $p^{\text{}}_{c}\!\left(\boldsymbol{y}\,\middle|\,\boldsymbol{x}\right)$ is the transition density (with respect to Lebesgue measure) of a Markov chain on $\mathbf{R}^{d}$ with invariant density proportional to $f_{c}^{2}$.

Proposition 2.1

The $(C+1)d$-dimensional (fusion) density proportional to the integrable function

[TABLE]

admits the marginal density $f$ for ${\boldsymbol{y}}$.

The proof of Proposition 2.1 is elementary and is omitted. The statistical interpretation of Proposition 2.1 is simply if it were possible to sample from the $(C+1)d$-dimensional (fusion) density $g$ in (2), then as a by-product we would obtain a draw from our fusion target density $f$ in (1). How to directly sample from (2) is not clear, even if it were possible to simulate from $p_{c}(\cdot|\boldsymbol{x})$. Our strategy instead will be to use rejection sampling. Two rejection sampling methods (which have differing efficiencies) are provided in Sections 3 and 4.

3 A fusion rejection sampler using Brownian bridges

3.1 The methodology

Consider the proposal density for the extended fusion target (2) proportional to the function

[TABLE]

where $\bar{\boldsymbol{x}}=C^{-1}\sum_{c=1}^{C}\boldsymbol{x}^{(c)}$, and $T$ is an arbitrary positive constant.

Simulation from the proposal $h^{bm}$ can be achieved directly. In particular, $\boldsymbol{x}^{(1)},\ldots\boldsymbol{x}^{(C)}$ are first drawn independently from $f_{1},\ldots f_{C}$ respectively, and then $\boldsymbol{y}$ is simply a Gaussian random variable centred on ${\bar{\boldsymbol{x}}}$.

In Proposition 2.1 we presented a general form of the extended fusion target, in which $p^{\text{}}_{c}\!\left(\boldsymbol{y}\,\middle|\,\boldsymbol{x}\right)$ is the transition density (with respect to Lebesgue measure) of a Markov chain on $\mathbf{R}^{d}$ with invariant density proportional to $f_{c}^{2}$. In this paper we set $p^{\text{}}_{c}\!\left(\boldsymbol{y}\,\middle|\,\boldsymbol{x}\right):=p^{\text{dl}}_{T,c}\!\left(\boldsymbol{y}\,\middle|\,\boldsymbol{x}\right)$, the transition density of a double Langevin diffusion for $f_{c}$ (i.e. the transition density of a Langevin diffusion for $f^{2}_{c}$) from $\boldsymbol{x}$ to $\boldsymbol{y}$ over a pre-defined (user-specified) time $T>0$. To distinguish the resulting extended fusion target from the general case, we further denote the extended fusion target by $g^{dl}$. In particular, for all $1\leq c\leq C$, we consider the $d$-dimensional (double) Langevin (DL) diffusion processes $\mathfrak{X}=\{\boldsymbol{X}_{t}^{(c)},t\in[0,T],c=1,\cdots,C\}$, given by

[TABLE]

where $\boldsymbol{W}_{t}^{(c)}$ is a $d$-dimensional Brownian motion, $\nabla$ is the gradient operator over $\boldsymbol{x}$ and

[TABLE]

$\boldsymbol{X}_{t}^{(c)}$ has invariant distribution $f^{2}_{c}\!\left(\boldsymbol{x}\right)$, for any $t\in[0,T]$ (Hansen, 2003). We also impose the following standard regularity property (where div denotes the divergence operator)

{condition}

Define

[TABLE]

There exists constant $\Phi^{\text{bm}}_{c}>-\infty$ such that for all $\boldsymbol{x}$ and each $c\in\{1,\cdots,C\}$, $\phi^{\text{dl}}_{c}\!\left(\boldsymbol{x}\right)\geq\Phi^{\text{bm}}_{c}$.

Then we have the following proposition which gives a rejection sampling method for $g^{dl}(\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(C)},\boldsymbol{y})$.

Proposition 3.1

Under Condition 3.1 we can write

[TABLE]

where

[TABLE]

and

[TABLE]

with ${\overline{\mathbb{W}}}$ denoting the law of $C$ Brownian bridges $\boldsymbol{x}^{(1)}_{t},\cdots,\boldsymbol{x}^{(C)}_{t}$ with $\boldsymbol{x}^{(c)}_{0}=\boldsymbol{x}^{(c)}$ and $\boldsymbol{x}^{(c)}_{T}=\boldsymbol{y}$ in the time interval $[0,T]$ (noting that these Brownian bridges are independent conditional on the starting and ending points) and

[TABLE]

Proof 3.2

From the Dacunha-Castelle representation (Dacunha-Castelle and Florens-Zmirou, 1986) we have that

[TABLE]

The result then follows from (2) and (3) by rearrangement and recalling that $\sum_{c=1}^{C}\|\boldsymbol{y}-\boldsymbol{x}^{(c)}\|^{2}=\sum_{c=1}^{C}\|\boldsymbol{x}^{(c)}-{\bar{\boldsymbol{x}}}\|^{2}+C||\boldsymbol{y}-{\bar{\boldsymbol{x}}}||^{2}$.

Here $\rho^{bm}$ and $Q^{bm}$ are both necessarily bounded by $1$, and when interpreted methodologically (see next section) correspond to separate acceptance steps within our rejection sampling framework. An event of probability $\rho^{bm}$ can be simulated by direct computation, and an event of probability $Q^{bm}$ can be simulated using the extensive efficient methodology on Poisson samplers (using an auxiliary diffusion bridge path-space rejection sampler as developed in (for instance) Beskos et al. (2006a, b, 2008), Chen and Huang (2013), Dai (2014), Pollock (2013), Pollock et al. (2016a) and Pollock et al. (2016b)). Note that there is a trade-off involved in the (user-specified) choice of $T$. For small $T$, $\rho^{bm}$ will likely be small while $Q^{bm}$ is large, whereas for large $T$ the opposite will be true. A small value of $T$ is usually preferred since the computational cost for the diffusion bridge rejection sampling for $Q^{bm}$ is comparatively expensive.

The algorithm for simulating from $f$ (by means of $g$) therefore proceeds as per Algorithm 1 (which we term Monte Carlo Fusion).

3.2 Practical interpretation of the algorithm

Remark 3.3

Adjustment for the simple average of the sample from sub-densities.

In the above algorithm, for all $c$, $\boldsymbol{x}^{(c)}$ is simulated from $f_{c}\!\left(\boldsymbol{x}\right)$ as in other Monte Carlo Fusion algorithms. The proposed combined value $\boldsymbol{y}$, however, is actually generated from a Gaussian distribution with mean $\bar{\boldsymbol{x}}$ and covariance matrix $C^{-1}T\mathbf{I}_{d}$. Therefore, $\boldsymbol{y}$ can be viewed as a simple average of the values $\{\boldsymbol{x}^{(c)},\ c\in C\}$ added to a Gaussian random error term. This algorithm indicates exactly how the simple average of these independent sub-posterior samples can be adjusted as a draw from the target distribution. This adjustment is in the form of the accept/reject step with acceptance probability $\rho^{bm}\times Q^{bm}$. $\square$

Remark 3.4

Implication on Bayesian group decision theory.

In Step 5 of Algorithm 1, we can also write $\log U_{1}\leq-\frac{C\sigma^{2}}{2T}$ as

[TABLE]

Note that $\sigma^{2}$ is actually the sample variance of the simulated starting points $\boldsymbol{x}^{(c)}$s (strictly speaking, it is the sum of variances for each component of $\boldsymbol{x}^{(c)}$). Thus in this step, the acceptance condition (13) implies that $\boldsymbol{y}$ will have a reasonable probability of being accepted as a sample from $f\!\left(\boldsymbol{x}\right)$ only when the variance of $\boldsymbol{x}^{(c)}$s is small enough. This coincides with our intuition in group decision making. For example, the small variance of $\{\boldsymbol{x}^{(c)},\ c\in C\}$ (satisfying condition (13)) means that the decisions/results from each group are similar, so that we can combine these decisions/results in step 7, using (12), in Algorithm 1. On the other hand, if the variance of $\{\boldsymbol{x}^{(c)},\ c\in C\}$ (not satisfying condition (13)) is large, then $\{\boldsymbol{x}^{(c)},\ c\in C\}$ provide contradictory evidence and should usually be rejected. Note that algorithmically, this initial rejection sampling step is efficient since early rejection then avoids the need to carry out the more complicated (and computationally expensive) step 7, via the condition of (12). $\square$

Remark 3.5

An extreme case.

A very interesting extreme case is to choose $T=0$. Then the event (12) will certainly happen, but the event $U_{1}\leq\exp\left(-\frac{C\sigma^{2}}{2T}\right)$ will occur only if $\boldsymbol{x}^{(1)}=\cdots=\boldsymbol{x}^{(C)}=\boldsymbol{y}$. Therefore in such an extreme case, Algorithm 1 (theoretically) draws a sample $\boldsymbol{y}$ from

[TABLE]

which is the target distribution. In practice, however, we have to choose $T>0$, since the independent sub-posterior samples $\boldsymbol{x}^{(c)}$s have zero probability to be the same. $\square$

4 A fusion rejection sampler using Ornstein-Uhlenbeck bridges

4.1 The methodology

In the previous section, the proposal density $h^{bm}$ uses a simple average of the sub-posterior samples $\boldsymbol{x}^{(c)}$ as the mean of the proposal $\boldsymbol{y}$. Another approach is to consider using a weighted average of the sub-posterior samples $\boldsymbol{x}^{(c)}$ as the mean of the proposal $\boldsymbol{y}$. For example, in a typical meta-analysis, a weighted average from different research outputs is typically used as the unification mean, and an individual output with more certainty (or, smaller variance) should consequently have larger weights (Fleiss, 1993).

Denoting $\hat{\boldsymbol{\mu}}_{c}$ and $\hat{\mathbf{\Lambda}}_{c}$ as the mean and the inverse of covariance matrix estimates for distribution $f_{c}\!\left(\boldsymbol{x}\right)$, respectively. We consider the following more general proposal density

[TABLE]

where

[TABLE]

and

[TABLE]

It is also straightforward to simulate from $h^{\text{ou}}\!\left(\cdot\right)$, since $\boldsymbol{x}^{(c)}$s are independently drawn from $f_{c}\!\left(\boldsymbol{x}\right)$ and then $\boldsymbol{y}$ is generated from a Gaussian distribution with mean $\tilde{\boldsymbol{x}}$ and covariance matrix $\mathbf{D}^{-1}$.

Here the vector $\boldsymbol{m}_{c}$ and the matrix $\mathbf{V}_{c}$ are actually the mean and covariance matrix, respectively, for the following OU process at time $T$ conditional on the starting point $\boldsymbol{x}_{0}^{(c)}=\boldsymbol{x}^{(c)}$ (Masuda, 2004),

[TABLE]

with

[TABLE]

We shall also require the following regularity property. {condition}

Define

[TABLE]

For any $c\in\{1,\cdots,C\}$, there exists $\Phi^{\text{ou}}_{c}>-\infty$ such that, for all $\boldsymbol{x}$,

[TABLE]

We now have the following result.

Proposition 4.1

Define function

[TABLE]

with

[TABLE]

Under Condition 4.1 we can write

[TABLE]

where

[TABLE]

with $\overline{\mathbb{O}}$ denoting the law of $C$ OU bridges $\boldsymbol{x}^{(c)}_{t},c=1,\cdots,C$ in time interval $[0,T]$. These $\boldsymbol{x}^{(c)}_{t}$s are independent conditional on the starting point $\boldsymbol{x}^{(c)}$ and common ending point $\boldsymbol{y}$ and

[TABLE]

Proof 4.2

See Appendix.

Since $\rho^{ou}$ and $Q^{ou}$ in (22) are always no more than $1$, we have the following rejection sampling algorithm, Algorithm 2.

In Algorithm 2, $T$ is a tuning parameter as well. If we choose a small value of $T$, the acceptance probability $\rho^{ou}(\boldsymbol{x}^{(1)},\cdots,\boldsymbol{x}^{(C)})$ will be very low. This is noticed by the fact that trace $(\mathbf{H}\mathbf{D}^{-1})=\infty$ and $\mathbf{M}_{1,c}$ and $\mathbf{M}_{2,c}$ are finite matrices, if $T=0$. Therefore in practice, we should choose a reasonably large value $T$. However, if we choose a large $T$, the acceptance probability (25) will be very small. In practice it is usually preferable to choose a small $T$ since the computational cost for the acceptance / rejection step of (25) is much higher.

4.2 Connection with Consensus Monte Carlo

In practice, people may employ an approximate version of Algorithm 2, since we can ignore the condition (25) in the rejection step 7, if we choose a small value $T$. In other words, from (22) of Proposition 4.1, we have that for small $T$ the target $g^{dl}(\cdot)$ is approximately equal to a function proportional to

[TABLE]

Such an approximation will be very good since $Q^{ou}$, the acceptance probability (25), will be very close to $1$ with a small $T$. In other words, we only simulate $\boldsymbol{y}$ from $\tilde{h}^{ou}(\cdot)$ and accept $\boldsymbol{y}$ as a sample approximately from the target distribution $f\!\left(\cdot\right)$.

If we draw $\boldsymbol{x}_{c},c=1,\cdots,C$ independently from each $f_{c}\!\left(\boldsymbol{x}\right)$, respectively. Another naive approach to combine these draws is to use the following linear combination

[TABLE]

In practice, $\hat{\mathbf{\Lambda}}_{c}^{-1}$ can be obtained via preliminary analysis. Such a $\boldsymbol{y}$ can be viewed as a sample approximately from $f\!\left(\boldsymbol{x}\right)$ as well. This is named as Consensus Monte Carlo (CMC) in (Scott et al., 2016).

The following lemma tells us how the simulation of $\boldsymbol{y}$ from $\tilde{h}^{ou}(\cdot)$ is related to the consensus Monte Carlo sample (27).

Lemma 4.3

If $f_{c}\!\left(\boldsymbol{x}\right)$ is a Gaussian distribution and if we choose $T=\infty$, then simulating $\boldsymbol{y}$ from $\tilde{h}^{ou}(\cdot)$ will be the same as the Consensus Monte Carlo method. Both draw samples exactly from the target distribution.

Proof 4.4

See Appendix B.

This lemma tells us why Consensus Monte Carlo does not provide good results. It is because CMC simulates $\boldsymbol{y}$ from $\tilde{h}^{ou}(\cdot)$ with $T=\infty$, however, $T$ should be chosen as a small value to achieve better approximation or even use the exact Algorithm 2 with the diffusion path rejection sampler.

5 Simulation studies

5.1 Distribution with light tails

We consider the target distribution $f(x)\propto e^{-x^{4}/2}$. We choose $C=4$ and $f_{c}\!\left(x\right)=e^{-x^{4}/2C},c=1,\cdots,C$. It is easy to check that $\phi^{dl}(x)=\frac{1}{2}\left(\frac{4x^{6}}{C^{2}}-\frac{6x^{2}}{C}\right)$ satisfies Condition 3.1. On the other hand, for any chosen values $\hat{\boldsymbol{\mu}}_{c}$, $\hat{\mathbf{\Lambda}}_{c}$, $\phi^{\text{ou}}_{c}\!\left(\boldsymbol{x}\right)$ in (19) satisfies Condition 4.1 as well. Therefore both Algorithm 1 and Algorithm 2 can be applied to simulate from $f(x)$.

We compare the following Monte Carlo methods for the estimation of the density function of $f(x)$:

Simulating MC samples directly from $f(x)$ via a simple rejection sampling with standard Gaussian distribution as the proposal;

• 2.

simulating MC samples based on the exact simulation method, Algorithm 1 with $T=1$;

• 3.

simulating MC samples based on the exact simulation method, Algorithm 2 with $T=1$;

• 4.

simulating MC samples based on the Consensus method in Scott et al. (2016);

The density curve estimation results are summarised in Figure 1. Note that all results are based on $10,000$ realisations. The black solid cure (Simulation [1.], the true fitted density curve), the blue solid curve (Simulation [2.]) and the pink solid curve (Simulation [3.]) are all exact algorithms and they are almost identical. The Consensus methods (Simulation [4.] – the red dashed curve, has very large biases. Note that both CMC algorithm and Algorithm 2 use the same values of $\hat{\boldsymbol{\mu}}_{c}$ and $\hat{\mathbf{\Lambda}}_{c}$ based on preliminary analysis.

The running time of the algorithms are presented in Table 1. It seems that Algorithm 2 uses the most system running time. This is because Algorithm 2 has a smaller acceptance probability (about $0.011$ for the path-space rejection sampling (25)) than that in Algorithm 1 (about $0.139$ for the path-space rejection sampling (12)).

Note that Condition 3.1 is usually satisfied in most applications, however, Condition 4.1 only holds when the target $f(x)$ has lighter tails than Gaussian distributions. We present an example in the following section, where only Condition 3.1 is true.

5.2 Beta distribution

Consider the target distribution as the Beta distribution with density $\pi(u)\propto u^{4}(1-u)$, $u\in[0,1]$, i.e. Beta$(5,2)$. To use the proposed algorithms, the support of the target distribution should be in the whole real axis. Therefore we need to use the variable transformation $x=\log(u/(1-u))$ and consider the target distribution as

[TABLE]

We decompose $\pi(x)$ into $C=5$ components,

[TABLE]

Note that for this simple example Condition 3.1 is satisfied but Condition 4.1 is not satisfied. Therefore we compare the following Monte Carlo methods for the estimation of the density function of Beta$(5,2)$:

Simulating MC samples directly from Beta$(5,2)$ via the simple R command, rbeta;

• 2.

simulating MC samples based on the exact simulation method, Algorithm 1 with $T=3$;

• 3.

simulating MC samples based on the Consensus method in Scott et al. (2016), with variable transformation and with the decomposition in (5.2).

For simulations [3.], $\hat{\boldsymbol{\mu}}_{c},\hat{\mathbf{\Lambda}}_{c}$ (also known as the weight of each consensus sample in CMC algorithm) are chosen, respectively, as the estimated mean and inverse of variance of $f_{c}\!\left(\cdot\right)$, as suggested by Scott et al. (2016).

The density curve estimation results are summarized in Figure 2. Note that all results are based on $10,000$ realisations. The black solid cure (Simulation [1.]) and the blue solid curve (Simulation [2.]) are almost identical, since both of them are based on exact simulation methods. Again, Concensus Monte Carlo gives very biased results.

Remark 4. Tuning parameters and density decomposition

We may choose any value $T$ in the proposed algorithms. However, as we mentioned before, value $T$ is a tunning parameter for the efficiency of Algorithm 1, which indeed shown by our simulation results. The CPU running time for simulating 10,000 realisations based on Algorithm 1 for the Beta-example is summarized in Figure 3. The optimum value is about $T=2$. Note that such optimum value can be found in practice via a small amount of preliminary simulations.

In addition, when we split the target $f$ into $f\propto f_{1}\cdots f_{C}$, we actually chose $f_{1}=\cdots=f_{C}=f^{1/C}$, since such decomposition will always give the smallest variation for $\boldsymbol{x}^{(c)},c=1,\cdots,C$, as suggested in Dai (2017).

6 Conclusion

In this paper we have introduced a novel theoretical framework, and direct methodological implementation (Monte Carlo Fusion), to address the common (but challenging) ‘fusion’ problem – unifying distributed analyses and inferences on shared parameters from multiple sources, into a single coherent inference – by viewing it as a simple rejection sampler on an extended space (Section 2), and developing an appropriate sampling mechanism (Section 3).

Our fusion approach is not only the first paper to answer in a principled manner how to combine samples from multiple sources, but (as shown in Section 4), also provides a principled approach to understand the errors that arise in other existing unification schemes (such as those in big data divide-and-conquer approaches). The errors in existing unification schemes can be considerable, even for simple one dimensional unification targets (such as those we consider in Section 5).

Characterising the error in existing unification schemes is possible by setting the (proposal) sampling mechanisms of those schemes within our framework, and finding a representation for the remaining error (which we could remove by further acceptance or rejection). This opens interesting avenues of research in which existing unification schemes are adapted within our framework into efficient proposal mechanisms for our extended fusion target density (2).

A number of avenues to apply our work directly to interesting applications are possible. In addition to those discussed in Section 1 (namely expert elicitation, multi-view learning and meta-analysis), other application areas include Bayesian group decision theory (see Remark 3.4) and Bayesian sensitivity analysis. In the case of Bayesian sensitivity analysis we need to assess a large number of prior distributions, but existing methods address this by using approximations (Tan et al., 2015).

A key avenue for (on-going) future research is to fully explore how to use the Monte Carlo Fusion framework we introduce to modify, and remove approximation, from existing Monte Carlo methods (particularly within the big data setting) that use the divide-and-conquer (fork-and-join) strategies described in Section 1. In the particular setting of big data the unification of distributed inferences is only part of the problem – additional constraints are imposed due to practical computational and hardware concerns (such as avoiding as far as possible communication between computing cores). As such, key future research will focus on how to implement Monte Carlo Fusion with this particular set of constraints (as direct implementation is not possible).

Another interesting big data direction would be to blend the multi-core approach of divide-and-conquer strategies, with state-of-the-art single-core approaches for big data (for instance, Pollock et al. (2016a)).

Acknowledgements

The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme “Scalable inference; statistical, algorithmic, computational aspects (SIN)” when work on this paper was undertaken. MP and GOR were supported by the EPSRC [grant number EP/K014463/1]. GOR was additionally supported by the EPSRC [grant numbers EP/K034154/1, EP/D002060/1].

Appendix A Proof of Proposition 4.1

To prove the proposition, we first introduce a lemma about a density proportional to the following function

[TABLE]

where $p^{\text{ou}}_{T,c}\!\left(\boldsymbol{y}\,\middle|\,\boldsymbol{x}^{(c)}\right)$ is the transition density from $\boldsymbol{x}^{(c)}$ at time [math] to $\boldsymbol{y}$ at time $T$, for the OU process given by (17).

Lemma A.1

The expression of $\tilde{g}^{ou}$.* The formula in (29) can be rewritten as*

[TABLE]

where $\widetilde{\boldsymbol{x}}$ and $\rho^{ou}$ are given by (4.1) and (4.1), respectively. $\square$

Proof A.2

See the supplementary file.

Now we prove Proposition 4.1.

Proof A.3

If we denote $\overline{\mathbb{DL}}$ as the law of $C$ $d$-dimensional Langevin diffusion bridges given in (4), with starting points $\boldsymbol{x}^{(c)}$ and common ending point $\boldsymbol{y}$, the result of Proposition 3.1 can be written as

[TABLE]

where $\vec{\boldsymbol{x}}=\{\boldsymbol{x}^{(c)}_{t},c=1,\cdots,C,t\in[0,T]\}$ are typical diffusion bridge paths, with starting points $\boldsymbol{x}^{(c)}$ and common ending point $\boldsymbol{y}$.

If we consider a very special case where each $f_{c}\!\left(\boldsymbol{x}\right)$ is a Gaussian density $\exp(A^{\text{ou}}_{c}\left(\boldsymbol{x}\right))$, we immediately have that the target $g^{dl}$ becomes

[TABLE]

In addition, the proposal density $h^{bm}$ in Propostion 3.1 becomes

[TABLE]

and further the rejection sampling ratio in Proposition 3.1 becomes

[TABLE]

[TABLE]

Now we have,

[TABLE]

Lemma A.1 and expressions of $E^{bm}$ and $E_{*}^{ou}$ immediately gives

[TABLE]

*where $E^{ou}$ given in (24). *

Appendix B Proof of Lemma 4.3

We consider the formula of $\rho^{ou}$ in (4.1). If $T=\infty$, we have

[TABLE]

and

[TABLE]

Therefore

[TABLE]

and further $\rho^{ou}(\boldsymbol{x}^{(1)},\cdots,\boldsymbol{x}^{(C)})$ becomes a value not depending on $\boldsymbol{X}_{0}^{(1:C)}$ at all, as $T\rightarrow\infty$. Therefore with $T=\infty$ the density function $\tilde{h}^{ou}(\cdot)$ becomes

[TABLE]

with

[TABLE]

Then we can generate $\boldsymbol{y}$ as, with some standard Gaussian random errors $\boldsymbol{\epsilon}$,

[TABLE]

where $\stackrel{{\scriptstyle distr.}}{{=}}{}$ means ‘equal in distribution’ and $\boldsymbol{\epsilon}_{c},c=1,\cdots,C$ means $C$ independent standard normal vectors.

By noticing that $\hat{\boldsymbol{\mu}}_{c}+\hat{\mathbf{\Lambda}}_{c}^{-1/2}\boldsymbol{\epsilon}_{c}$ has the same distribution as $\boldsymbol{x}^{(c)}$ if $f_{c}\!\left(\boldsymbol{x}\right)$ is a Gaussian distribution, the lemma is proved.