---
abstract: |
    Sequential Monte Carlo has become a standard tool for Bayesian Inference of complex models. This approach can be computationally demanding, especially when initialized from the prior distribution. On the other hand, deterministic approximations of the posterior distribution are often available with no theoretical guaranties. We propose a bridge sampling scheme starting from such a deterministic approximation of the posterior distribution and targeting the true one. The resulting Shortened Bridge Sampler (SBS) relies on a sequence of distributions that is determined in an adaptive way.

    We illustrate the robustness and the efficiency of the methodology on a large simulation study. When applied to network datasets, SBS inference leads to different statistical conclusions from the one supplied by the standard variational Bayes approximation.
author:
- Sophie Donnet
- Stéphane Robin
bibliography:
- 'biblio.bib'
date: 'Received: date / Accepted: date'
subtitle: Using deterministic approximations to accelerate SMC for posterior sampling
title: 'Shortened Bridge Sampler: '
---

[example.eps]{} gsave newpath 20 20 moveto 20 220 lineto 220 220 lineto 220 20 lineto closepath 2 setlinewidth gsave .4 setgray fill grestore stroke grestore

Introduction \[sec:Intro\]
==========================

In Bayesian statistics, except in a restricted number of conjugate models, the posterior distribution does not have a close form and requires the use of approximation methods. The 1990’s have witnessed the developments of powerful stochastic methods and simulation-based algorithms able to perform Bayesian statistical inference on complex statistical models.

Among them, Monte Carlo Markov Chains (MCMC) –whose principle is to generate a Markov Chain such that its ergodic distribution is the posterior distribution (see for instance [@Robert:2005; @RoC09] for an introduction)– have been successfully applied to many problems, as attested by the countless publications in that domain. However, MCMCs suffer from several drawbacks. First of all, it is difficult to assess whether the chain has reached its ergodic distribution or not. Secondly, if the distribution of interest is highly multi-modal, MCMC algorithms can be trapped in local modes. More generally, when the space of parameters to be explored is of high dimension, MCMC algorithms will have difficulties in reaching their equilibrium distribution within a reasonable computational time.

Recently, population based Monte Carlo methods have proved their efficiency and robustness in front of high dimensional and multimodal spaces. In a few words, population based Monte Carlo algorithms generate a large sample of parameters with a tractable distribution and update the importance sampling weights at each iteration, in order to finally match the distribution of interest. Among population based Monte Carlo, Sequential Monte Carlo (SMC) is a method combining parameters sampling and resampling. More precisely, a sequence of distributions of interest is designed, such that the first one is simple (i.e. easy to sample from) and the last one is the posterior distribution. This sequence of distributions defines the iterations of the algorithm. Then, at the first iteration, a sample of parameters is simulated with the first distribution. In the following iterations, the parameters are stochastically moved, weighted and resampled to follow the current distribution. The true posterior distribution is reached at the last iteration. The sequence of distributions can be dynamically designed. Primarily developed in the context of filtering problems (see [@Doucet2001]), they have been extended to the general problem of posterior sampling by [@DelMoral2006]. In comparison with MCMC methods, SMC does not require any burn-in period or convergence diagnostic. In addition, whereas computing marginal likelihood (for model comparison) has always been a challenging issue when using MCMC, SMC supplies an unbiased estimator of this quantity as a by-product of the algorithm. For all these reasons, SMC has proved its superiority over MCMC for complex models.

In the recent years, particular fields (such as genomics or network analysis to name but a few) brought news statistical problems involving an increasing amount of data or statistical models with a large number of parameters. In such cases, not only MCMC but also population Monte Carlo algorithms have reached their limitations, requiring unreasonable computational time to explore the posterior distribution. To deal with such difficulties, deterministic approximations of the posterior distribution through optimization mathematical tools – such as variational approximation ([@WaJ08; @Blei2016]), Expectation-Propagation ([@Minka:2001]) or Integrated nested Laplace approximation ([@Rue09]) for instance– have been proposed. These methods have the great advantage to be computationally light and can handle large data. However, their theoretical properties and accuracy is still under study. In particular, we do know that variational approximations can supply underestimated posterior variances (see for instance [@CoM07] for a large illustration of this phenomena on the Probit model).

On the one hand, SMC supplies a sample from the exact posterior distribution but can require unacceptable computational time. On the other hand, deterministic approximations (optimal in a sense to be determined) of the posterior distribution are fast methods but non-exact. One may therefore be tempted to take advantage of the two approaches in a combined strategy. The idea of combining variational Bayes inference with SMC is actually not new. [@RAJ15] split the data into block and compute the posterior distribution of $\theta$ given each block. They use a variational argument to propose the product of this partial posterior as a proxy for the true posterior. Focusing on Gaussian mixtures, [@McGPTAK16] consider online-inference and propose an sequential sampling scheme where, for each new batch of data, the variational approximation is iteratively updated and used as a proposal. [@Naesseth2017] use a SMC approach to get an improved, but still biased, variational approximation. Our approach is different from all these. Our main idea is to design a bridge sampling from the approximated posterior distribution to the true posterior distribution, the transfer from the approximate to the exact distribution being performed with an SMC algorithm ([@DelMoral2006]). The sampling method we propose can be considered from two points of view : either SMC is seen as a tool to correct the approximate distribution, or the approximate posterior distribution is seen as a mean to drastically accelerate the SMC procedure.

Adopting the last perspective plunges the problematic into in a larger topic. Indeed, in general, for any challenging statistical model at stake, there exists a frequentist solution, suppling a point estimation of the parameter of interest. In the Bayesian practice, this point estimation is standardly set as an initial value in the MCMC algorithm, thus hoping a decrease of the computational time. In a SMC strategy starting with a sample from the prior distribution, such an initial point value is meaningless. We claim that the posterior distribution can be reached in a reduced computational time if the bridge sampling scheme starts from an approximate posterior distribution based on that point estimate. We therefore refer to the proposed sampling method as the Shortened Bridge Sampler ([SBS]{}).

The paper is organized as follows. Section \[sec:Algo\] is dedicated to the description of the methodology. We first remind the principle of importance sampling, then introduce the sampling path and expose the algorithm in Subsection \[subsec:algo\]. Its robustness and efficiency are illustrated on several simulated experiments in Section \[sec:Simul\]. In Subsection \[sec:log\], the logit regression serves as a toy example to illustrate the computational time reduction and to test the robustness of the method with respect to the quality of the deterministic approximation of the posterior distribution. The Latent Class Analysis model (Subsection \[subsec:LCA\]) is exploited to illustrate the relevance of our methodology on a mixture model, in particular, we propose a new strategy to tackle the label switching issue. On the Stochastic Block Model (SBM) with covariates (Subsection \[subsec:SBMreg\]), we compare our strategy with the Variational Bayesian one in terms of model selection and model averaging. Finally, real datasets of social networks with covariates are presented in Section \[sec:Illust\] : we stress the new insights brought by a “correction” of the Variational posterior approximation by the SMC strategy in term of both model averaging and significance of the parameters. Perspectives are discussed in Section \[sec:Discuss\].

From the approximate posterior distribution to the true posterior distribution {#sec:Algo}
==============================================================================

Let us first introduce some notations. ${\boldsymbol{Y}}$ denotes the observations, $\ell({\boldsymbol{Y}}|{\theta})$ is the likelihood function with ${\theta}\in \Theta$ the unknown parameters and $\pi({\theta})$ is the prior distribution on ${\theta}$. The Bayesian inference is based on the posterior distribution: $$\label{eq:Bayes}
 p({\theta}| {\boldsymbol{Y}}) = \frac{\ell({\boldsymbol{Y}}| {\theta}) \pi({\theta})}{p({\boldsymbol{Y}})}.$$ where $p({\boldsymbol{Y}})$ is the marginal likelihood defined as: $$\label{eq:marg}
 p({\boldsymbol{Y}}) = \int \ell({\boldsymbol{Y}}| {\theta}) \pi({\theta}) d{\theta}$$ and is required in the Bayesian model selection procedure.

In what follows, ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ is an approximate posterior distribution on ${\theta}$. *We assume that ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ can be easily intensively simulated and that the density function of ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ has an explicit expression* The aim of this paper is to propose a way to use such an approximate posterior to actually sample from the true posterior.

Note that, in general, complex statistical models are written as hierarchical models and involve latent variables ${\boldsymbol{Z}}$ (see Sections \[subsec:LCA\] and \[subsec:SBMreg\]). In such cases, the distributions of interest are the joint distribution $p({\boldsymbol{Z}},{\theta}| {\boldsymbol{Y}})$ or the marginal one $p({\theta}| {\boldsymbol{Y}})$. For the sake of simplicity, we chose to present the method without latent variables but obviously, all the following results and algorithms can be extended to this situation, replacing ${\theta}$ by $({\boldsymbol{Z}},{\theta})$. A substantial part of the simulations presented in Section \[sec:Simul\] is devoted to such models.

A first approach: Importance Sampling {#sec:is}
-------------------------------------

A first naive approach to use ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ consists in resorting to a simple importance sampling (IS) strategy, that is to say sampling $({\theta}^m)_{m = 1\dots M}$ from ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ and weighting the sample by $W^m \propto {\ell({\boldsymbol{Y}}| {\theta}^m) \pi({\theta}^m)}/{{{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta}^m)}$. However, this strategy is obviously naive for several reasons. First of all, there is no guarantee that the support of the approximate distribution includes the support of the true distribution, the contrary is even observed in practice for the variational approximation for instance [see @CoM07; @WaT05]. As a consequence, the posterior sample obtained through such an importance sampling strategy would be restricted to the support of ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ which can be strictly included in the support of $p(\cdot | {\boldsymbol{Y}})$. Secondly, if ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ and ${p(\cdot | {\boldsymbol{Y}})}$ are very different, the sample will be degenerated, meaning that very few particles will have a non-negligible weight. This results in a small Effective Sample Size ($ESS$), which the algorithm we propose aims at keeping along iterations. In such situations, there is no hope to efficiently sample using ’one-step’ IS but the principle can be used iteratively to progressively shift from the initial proposal to the true posterior distribution.

A path sampling between the approximate and the true posterior distributions
----------------------------------------------------------------------------

The main idea of this paper is to take advantage of the deterministic approximation ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ of the posterior distribution to accelerate SMC procedure, or, inversely to transform sequentially a sample from the deterministic approximated posterior distribution into a sample from the true posterior distribution.

Sequential Monte Carlo samplers generate samples from a sequence of intermediate distributions $({p_h})_{h = 0\dots H}$ where the intermediate distributions $({p_h})_{h = 0\dots H}$ are smooth transitions from a simple distribution $p_0$ to the distribution of interest $p_H = p(\cdot | {\boldsymbol{Y}})$. A classical choice for $({p_h})_{h = 0\dots H}$ [@Neal2001] is to consider: $$\begin{aligned}
\label{eq:pih2}
 {p_h}({\theta}) & \propto & \pi({\theta}) \ell({\boldsymbol{Y}}| {\theta}) ^{\rho_h}\end{aligned}$$ where $\rho_0 = 0$, $\rho_H = 1$, thus slowly moving from the prior distribution to the posterior by progressively integrating the data ${\boldsymbol{Y}}$ through the likelihood function. In this paper, we propose an alternative scheme moving smoothly from the approximate posterior distribution ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ to the true ${p(\cdot | {\boldsymbol{Y}})}$. The path is thus defined by: $$\begin{aligned}
\label{eq:pih}
 {p_h}({\theta}) & \propto & {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}(p( {\theta}| {\boldsymbol{Y}})) ^{\rho_h}\nonumber\\
 & \propto & {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}( \pi({\theta}) \ell({\boldsymbol{Y}}| {\theta})) ^{\rho_h}.\end{aligned}$$ where, $\rho_0 = 0$, $\rho_H = 1$. In a few words, we start from the easy-to-sample distribution $ {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})$ and progressively replace it with the true posterior distribution, this strategy being known as annealed importance sampling procedure [@Neal2001]. We claim that this scheme significantly reduces the computational time and is robust with respect to $ {{\widetilde{p}}_{{\boldsymbol{Y}}}}$.

Note that if ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ is chosen to be the prior distribution $\pi(\cdot)$, then schemes and are identical.

An alternative strategy would consist in using ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ has an importance sampler in the first iteration of the standard annealing scheme defined in . However, in such a strategy, the approximate distribution ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ is under-exploited, since at the first iteration, the particles are reweighted with $W_0^m \propto {\pi({\theta}^m)}/{{{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta}^m |{\boldsymbol{Y}})}$, thus going back, in practice, to a (possibly truncated) version of the prior distribution. This phenomena will be illustrated in Section \[sec:log\].

To sample from the sequence of distributions $(p_h)_{h = 1, \dots, H}$, we adopt the sequential sampler proposed by [@DelMoral2006] where the annealing coefficients $\rho_h$ will be adjusted dynamically. We describe the algorithm in the following subsection.

Shortened Bridge Sampling Algorithm {#subsec:algo}
-----------------------------------

We now need to design an algorithm sequentially sampling from $ {p_h}({\theta}) \propto {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}( \pi({\theta}) \ell({\boldsymbol{Y}}| {\theta})) ^{\rho_h}$. The last years have witnessed a proliferation of scientific papers dealing with the problem of SMC methods and their applications [see @Doucet2001 for an overview]. In our work, we resort to the algorithm proposed by @DelMoral2006.

Let us introduce the following notations: $$\begin{aligned}
\label{eq:gammah_Zh}
 \gamma_h(\theta) & = & {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-\rho_h} \left[\ell({\boldsymbol{Y}}| {\theta}) \pi({\theta})\right]^{\rho_h},\end{aligned}$$ where $ Z_h = \int \gamma_h(\theta) {\text{ d}}\theta \nonumber $, so that ${p_h}(\theta) = \gamma_h(\theta) / Z_h$ is a probability density. The main idea of [@DelMoral2006] is to plunge the problem of sampling a sequence of distributions defined on a single set $\Theta$ into the standard SMC filtering framework. To that purpose, the sequence $({p_h})_{h = 0\dots H}$ is replaced by a sequence of extended distributions: $$\label{eq:pibarh}
{\overline{p}}_h({{\theta}_{0:h}}) = \frac{{\overline{\gamma}}_h({{\theta}_{0:h}})}{Z_h}$$ with $$\label{eq:pibarh2}
{\overline{\gamma}}_h( {{\theta}_{0:h}}) = \gamma_h({{\theta}_h}) \prod_{k = 1}^{h} L_{k}\left({{\theta}_{k-1}}| {{\theta}_{k}}\right)$$ where ${{\theta}_{0:h}}= ({{\theta}_{0}},\dots, {{\theta}_h}) \in \Theta \times \dots\times \Theta = \Theta^{h+1}$ and $(L_k)_{k = 0,\dots H-1}$ is a sequence of backward kernels satisfying: $$\label{eq:L}
\int L_k\left({{\theta}_{k-1}}| {{\theta}_{k}}\right)d{\theta}_{k-1} = 1, \quad \forall k = 0\dots H-1.$$ Due to Property , the marginal version of ${\overline{p}}_h$ (i.e. when integrating out ${{\theta}_{0}}$ , $\dots$, ${\theta}_{h-1}$) is the distribution of interest ${p_h}$. Once defined the sequence $({\overline{p}}_h)_{h = 0\dots H}$, one may use the original SMC algorithm designed by @Doucet2001 for filtering. At iteration $h$, the SMC sampler involves three steps:

-   *Moving the particles* from ${\theta}_{h-1}$ to ${\theta}_{h}$ using a transition kernel $ K_h({{\theta}_h}| {\theta}_{h-1})$. As a consequence, let $ \eta_{h-1}({\theta}_{0:h-1})$ denote the sampling kernel for $ {\theta}_{0:h-1}$ until iteration $h-1$, $\eta_h$’s expression is: $$\label{eq:eta}
     \eta_h({{\theta}_{0:h}}) = \eta_{h-1}({\theta}_{0:h-1}) K_h({{\theta}_h}| {\theta}_{h-1})$$

-   *Reweighing the particles* in order to correct the discrepancy between the sampling distribution $\eta_h$ and the distribution of interest at iteration $h$, ${\overline{p}}_h$.

-   *Selecting the particles* in order to reduce the variability of the importance sampling weights and avoid degeneracy. In practice the particles will be resampled when the $ESS$ decreases below a pre-specified rate.

#### About the importance weights.

At iteration $h$, the importance sampling weights for $({{\theta}_{0:h}}^m)_{m = 1\dots M}$ are : $\forall m = 1\dots M$, $$\label{eq:w1}
w_h^m = {w_h}({{\theta}_{0:h}}^m) = \frac{{\overline{\gamma}}_h({{\theta}_{0:h}}^m)}{\eta_h({{\theta}_{0:h}}^m)}$$ in their unnormalized version. $(W^m_h)_{m = 1\dots M }$ denotes the normalized weights, i.e. $$\label{eq:wW}
{W_h}^m = \frac{{w_h^m}}{\sum_{m' = 1}^M w^{m'}_h}, \quad \forall m = 1\dots M$$ Equations (\[eq:pibarh\]-\[eq:pibarh2\]-\[eq:eta\]-\[eq:w1\]) imply a recurrence formula for the weight of any particle ${{\theta}_{0:h}}$: $$\label{eq:w}
{w_h}({{\theta}_{0:h}}) = w_{h-1}({\theta}_{0:h-1}) {\widetilde{w}}_{h-1:h}({\theta}_{h-1},{{\theta}_h})$$ where the incremental weight $ {\widetilde{w}}_{h-1:h}({\theta}_{h-1},{{\theta}_h})$ is equal to: $$\label{eq:wtilde}
 {\widetilde{w}}_{h-1:h}({\theta}_{h-1},{{\theta}_h}) = \frac{L_h({\theta}_{h-1} | {\theta}_h)}{ K_h({{\theta}_h}| {\theta}_{h-1})}\frac{\gamma_h({{\theta}_h})}{\gamma_{h-1}({\theta}_{h-1})}$$

#### About the transition kernel $K_h$.

As, at this step, the target distribution is ${p_h}$, it seems natural to choose $K_h({\theta}_h | {\theta}_{h-1})$ as a Monte Carlo Markov Chain (MCMC) kernel with ${p_h}({\theta}) \propto {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}( \ell({\boldsymbol{Y}}| {\theta}) \pi({\theta}) ) ^{\rho_h}$ as stationary distribution. Following [@DelMoral2006], we choose the backward kernel: $$\label{eq:L2}
L_h({\theta}_{h-1} | {\theta}_h) = \frac{ K_h({\theta}_h | {\theta}_{h-1}) p_{h}({\theta}_{h-1}) }{{p_h}({\theta}_h)}$$ which satisfies Property and enables us to rewrite the weight increment $ {\widetilde{w}}_{h-1:h}({\theta}_{h-1},{{\theta}_h})$ appearing in and defined in as $$\label{eq:wtilde2}
 {\widetilde{w}}_{h-1:h}({\theta}_{h-1},{\theta}_h) = \frac{\gamma_{h}({\theta}_{h-1}) }{\gamma_{h-1}({\theta}_{h-1})} = \left[\alpha(\theta_{h-1})\right]^{\rho_h-\rho_{h-1}}$$ where $$\label{eq:alpha}
\alpha(\theta) = \frac{\ell({\boldsymbol{Y}}| {\theta}) \pi({\theta})}{{{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta}| {\boldsymbol{Y}})}.$$ In what follows, we denote $\alpha_h = \alpha({\theta}_h)$.

Using this particular backward kernel has two major consequences. First it is not required having an explicit expression for the transition kernel $K_h({\theta}_h | {\theta}_{h-1})$, which is quite welcome for MCMC kernels. Secondly, examining equations and , one may notice that the weight for a particle ${{\theta}_{0:h}}$ does not depend on $\theta_h$ but only on ${\theta}_{0:h-1}$. As a consequence, the weights of the particles ${{\theta}_{0:h}}$ can be computed before they are simulated and for any new ${p_h}$

#### Adaptive design of $(\rho_h)_{h = 0\dots H}$.

As a consequence of this last remark, we are able to design an adaptive strategy for $({\rho_h})_{h = 0, \dots H}$ [as in @Schafer2013; @Jasra2011]. Indeed, being able to compute the weights of the up-coming particles for any new $\rho_h$, we can increase $\rho_h$ until the quality of the sample (measured through an indicator computed from the weights) decreases for the next distribution. In practice, following [@Zhou2016], we use the conditional Effective Sampling Size ($cESS$) to measure the quality of $p_ {h-1}$ as an importance sampler when estimating an expectation against ${p_h}$. It is defined as: $$\begin{aligned}
\label{eq:cESS}
 cESS & = & \left[\sum_{m = 1}^M M W_{h-1}^m \left(\frac{ {\widetilde{w}}_{h-1:h}^m}{\sum_{m = 1}^M MW_{h-1}^m {\widetilde{w}}_{h-1:h}^m} \right)^2 \right]^{-1} \\
 & = & \frac{M \left(\sum_{m = 1}^M W_{h-1}^m {\widetilde{w}}_{h-1:h}^m\right)^2}{\sum_{m = 1}^M W_{h-1}^m ( {\widetilde{w}}_{h-1:h}^m)^2}, \end{aligned}$$ becoming $$\begin{aligned}
\label{eq:cESS2}
&&cESS\left(\rho_h;\rho_{h-1}, (W_{h-1}^m,\alpha_{h-1}^m )_{m  \leq  M}\right) =cESS_{h-1}(\rho)\nonumber\\
&&
= \frac{M \left(\sum_{m = 1}^M W_{h-1}^m (\alpha_{h-1}^m)^{\rho_h -\rho_{h-1}}\right)^2}{\sum_{m = 1}^M W_{h-1}^m (\alpha_{h-1}^m)^{2(\rho_h -\rho_{h-1})}}. \end{aligned}$$ where $\alpha_{h-1}^m = \alpha({\theta}_{h-1}^m)$ has been defined in equation . If $\rho_h = \rho_{h-1}$ , $cESS$ si maximal (equal to $M$, the number of particles). As $\rho_h$ increases, the discrepancy between $p_ {h-1}$ and ${p_h}$ increases and so the quality of $p_ {h-1}$ as an importance sampling distribution when estimating an expectation against ${p_h}$ decreases and so does $cESS$. As a consequence, our strategy to find the next $\rho_h$ is to set: $$\rho_ h = 1 \wedge \sup_{\rho}\left\{\rho > \rho_{h-1}, cESS_{h-1}(\rho) \geq \tau_1 M \right\}$$

#### Selection of the particles.

In order to prevent a degeneration of the particle approximation, we use a standard resampling of the particles. In other words, if the variance of weights $(W_{h}^m)_{m = 1\dots M}$ is too high (or in other words, if the $ESS$ is too small), we resample the particles using a multinomial distribution, thus discarding the particles with low weights and duplicating the particles with high weights.

#### Sampling algorithm.

Finally, we propose the following Shorten Bridge Sampling [SBS]{}algorithm adapted to the sequence (\[eq:pih\]).

------------------------------------------------------------------------

**[SBS]{}algorithm**

------------------------------------------------------------------------

1.  Set $(\tau_1,\tau_2) \in [0,1]^2$, $\rho_0 = 0$.

2.  *At iteration $0$* , sample $(\theta^m_{0})_{m = 1\dots M}$ from the approximate distribution ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$. $ \forall m = 1\dots M$, set: $$w_{0}^m = 1, \quad W_{0}^m = \frac{1}{M}, \quad \alpha_{0}^m = \frac{\ell({\boldsymbol{Y}}| {\theta}^m_{0}) \pi({\theta}^m_{0})}{{{\widetilde{p}}_{{\boldsymbol{Y}}}}(\theta^m_{0})}$$

3.  *At iteration $h$*: starting from $({\theta}_{h-1}^m, W_{h-1}^m,\alpha_{h-1}^m)_{m = 1\dots M}$

    1.  Find $\rho_h$ such that: $$\rho_ h = 1 \wedge \sup_{\rho}\left\{\rho > \rho_{h-1},cESS_{h-1}(\rho) \geq \tau_1 M \right\},$$

    2.  $\forall m = 1\dots M $, compute ${w_h^m}= w^m_{h-1}\left(\alpha_{h-1}^m\right)^{\rho_h - \rho_{h-1}}$ and ${W_h^m}= {{w_h^m}}\left/{\sum_{m' = 1}^M w_h^{m'}}\right.$

    3.  Compute $$ESS_{h} = \frac{\left(\sum_{m = 1}^M{W_h^m}\right)^2 }{\sum_{m = 1}^M ({W_h^m})^2} \in [1, M]$$ If $ESS_{h} < \tau_2\, M$, resample the particles $$\begin{array}{ccl}
         ({\theta}_{h-1}^m)' & \sim_{i.i.d} & \sum_{m = 1}^M {W_h^m}\delta_{\{ {\theta}^m_{h-1}\}}
         ({{\theta}_h^m})\\
         {{\theta}_{h-1}^m}& \leftarrow & ({\theta}_{h-1}^m)'\\
         w_{h}^m &  \leftarrow  & 1 \\
         W_{h}^m & \leftarrow  & 1/M
         \end{array}
         \quad \forall m = 1\dots M$$

    4.  $\forall m = 1\dots M$, : propagate the particle ${{\theta}_h^m}\sim K_h( \cdot | {\theta}_{h-1}^m) $ where $K_h$ is a MCMC kernel with ${p_h}({\theta}) \propto {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}( \ell({\boldsymbol{Y}}| {\theta}) \pi({\theta}) ) ^{\rho_h}$ as an invariant distribution and compute: $$\alpha^m_{h} = \frac{\ell({\boldsymbol{Y}}| {\theta}^m_{h}) \pi({\theta}^m_{h})}{{{\widetilde{p}}_{{\boldsymbol{Y}}}}(\theta^m_{h}| {\boldsymbol{Y}})}$$

4.  If $\rho_h = 1$, stop. If $\rho_h < 1$ return to $1$.

------------------------------------------------------------------------

Let $\phi$ be a function defined in $\Theta$. The study of the statistical properties of $\sum_{m=1}^ M W^H_m \phi(\theta^{H}_m)$ as an estimator of $\mathbf{E}[\phi(\theta)|{\boldsymbol{Y}}]$ is a difficult task due to the sampling and resampling steps in the algorithm. However, many results can be found in the literature [see @Doucet2009 and references there in]. First of all, $\sum_{m=1}^ M W^H_m \phi(\theta^{H}_m)$ is known to be strongly convergent. Moreover, following [@DelMoral2006], a Central Limit Theorem can be obtained. Besides, in addition to these asymptotic properties, it is possible to control the mean-square error of the estimator for a given number of particles $M$, provided additional assumptions on $\phi$ . Results of convergence were also provided by [@delmoral2012_conv] for adaptive sequential Monte Carlo algorithms.

Estimation of the marginal likelihood 
--------------------------------------

With respect to MCMC strategies, Annealing Importance Sampling and SMC have the great advantage to supply good estimators of the marginal likelihood. Indeed, as proved by [@DelMoral2006], a non-biaised estimator of the marginal likelihood derives as a by-product of SMC. Moreover, the path sampling identity also provides an estimate of the marginal likelihood, as detailed hereafter.

Let us recall that $Z_ h = \int_{\Theta} \gamma_h({\theta}) d\theta$. Following [@DelMoral2006] and using the notations introduced below, the ratio of the marginal likelihoods ${Z_{h}} / {Z_{h-1}}$ is estimated by: $$\widehat{\frac{Z_{h}}{Z_{h-1}}} = \sum_{m = 1}^M {W_h^m}{\widetilde{w}}_{h-1:h}^m.$$ and $$\label{eq:pZhat}
\widehat{\frac{Z_{H}}{Z_{0}}} = \prod_{h = 1}^H \widehat{\frac{Z_{h}}{Z_{h-1}}} = \prod_{h = 1}^H \sum_{m = 1}^M {W_h^m}{\widetilde{w}}_{h-1:h}^m$$ is an unbiased estimator of $Z_H/Z_0$. Now, $Z_H =  p({\boldsymbol{Y}})$ and $Z_0 = 1 $.

Another estimate is given by the path sampling identity. Indeed, under non-restrictive regularity assumptions, the following equality holds: $$\label{eq:Int}
\log p(Y) - \log Z_0 = \int_{0}^{1} \mathbb{E}_{p_{\rho}}\left[\frac{d \log \gamma_{\rho}(\cdot)}{d \rho} \right]d\rho$$ where $\gamma_\rho(\theta) = {{\widetilde{p}}_{{\boldsymbol{Y}}}}(\theta)^{1-\rho}(\ell({\boldsymbol{Y}}| \theta) \pi(\theta))^{\rho}$, $p_\rho()$ is the associated probability density distribution and, in our geometric path sampling: $$\frac{d \log \gamma_{\rho}(\cdot)}{d \rho} = \log \frac{\ell({\boldsymbol{Y}}| {\theta}) \pi({\theta})}{{{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})} = \log \alpha(\theta)$$ An elementary trapezoidal scheme and Monte Carlo approximations of the expectations involved in lead to the following approximation of the marginal likelihood: $$\label{eq:pZhat2}
\widehat{\widehat{\log p({\boldsymbol{Y}})}} = \sum_{h = 1}^H \frac{\rho_h - \rho_{h-1}}{2} (U^M_h + U^M_{h-1})$$ where $U^M_h = \widehat{ \mathbb{E}}_{p_{\rho_h}}\left[ \log \alpha(\theta)\right] = \sum_{m = 1}^M W_{h}^m \log \alpha^m_{h}
 $.

Note that as suggested in [@Zhou2016], we noticed on simulation studies that the two estimators behave similarly in our examples. A precise comparison of the two estimators is out of the scope of this paper.

Simulation study \[sec:Simul\]
==============================

We now present a large simulation study. The goal of the study is to assess the fact that the proposed [SBS]{}algorithm –that combines an optimization-based approximation of the posterior distribution with a SMC sampler– drastically decreases the computational time with respect to a classical annealing-scheme or, equivalently, that the approximated posterior distribution can be corrected into the true posterior distribution at a low computational cost.

Logistic regression {#sec:log}
-------------------

This first model is used as a toy example to illustrate the efficiency and the robustness of our methodology. Let $(Y_1, \dots, Y_n)$ be a set a $n$ independent observations with values in $\{0,1\}$. Any individual observation $i$ is described by a vector ${\boldsymbol{x}}_i \in \mathbb{R}^p$ of $p$ covariates and we consider the logistic regression model: $$P(Y_i = 1) = 1-P(Y_i = 0) = \frac{e^{{\boldsymbol{x}}_i^t \theta}}{1+e^{{\boldsymbol{x}}_i^t \theta} }$$

We generate a simulated dataset with a randomly chosen regression matrix $X$ and the following parameters: $$n = 200, \quad p = 4, \quad \theta = (0.5, -0.6, 0, -1).$$

Setting a Gaussian prior distribution on $\theta \in \mathbb{R}^p$, $\theta \sim \mathcal{N}(0,\sigma^2\boldsymbol{I}_p),$ with $\sigma^2 = 100$, it is natural to propose a Gaussian approximation of the posterior distribution ${{\widetilde{p}}_{{\boldsymbol{Y}}}}: = \mathcal{N}(\widehat{\mu},\widehat{\Sigma})$. We first computed the Gaussian variational Bayes estimator and we obtained: $$\widehat{\mu}^{VB} = (0.398, -0.643, -0.280, -0.847)$$ and $$ \^[VB]{} = 10\^[-3]{} (

[rrrr]{} 23.3 & -1.2 & 0.1 & 1.1\
-1.2 & 20.6 & -0.1 & 2.0\
0.1 & -0.1 & 24.0 & 1.4\
1.1 & 2.0 & 1.4 & 25.5

) . $$ We also considered an approximate Gaussian posterior distribution based on the maximum likelihood estimator (ML) $\widehat{\theta}^{ML}$. The ML and its asymptotic variance $ \widehat{\Sigma}^{ML}$ are obtained with the -function . To test the robustness of our bridge sampling with respect to ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$, we also consider an artificially increased (respectively decreased) variance. In addition, we consider a distribution centered on an aberrant value with a small variance, thus leading to five different ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$: $$\begin{aligned}
\widetilde{p}_{{\boldsymbol{Y}},1} & = & \mathcal{N}(\widehat{\mu}^{VB}, \widehat{\Sigma}^{VB})\\
\widetilde{p}_{{\boldsymbol{Y}},2} & = & \mathcal{N}(\widehat{\theta}^{ML}, \widehat{\Sigma}^{ML})\\
 \widetilde{p}_{{\boldsymbol{Y}},3} & = & \mathcal{N}(\widehat{\mu}^{VB}, {\text{diag}}(\widehat{\Sigma}^{VB} )/5)\\
 \widetilde{p}_{{\boldsymbol{Y}},4} & = & \mathcal{N}(\widehat{\mu}^{VB}, {\text{diag}}(\widehat{\Sigma}^{VB} )\times10)\\
 \widetilde{p}_{{\boldsymbol{Y}},5} & = & \mathcal{N}(\widehat{\mu}^{VB} + 0.5, {\text{diag}}(\widehat{\Sigma}^{VB} )/5)
 \end{aligned}$$ $\widetilde{p}_{{\boldsymbol{Y}},1}$, $\widetilde{p}_{{\boldsymbol{Y}},2}$, $\widetilde{p}_{{\boldsymbol{Y}},3}$ $\widetilde{p}_{{\boldsymbol{Y}},4}$ and $\widetilde{p}_{{\boldsymbol{Y}},5}$ are plotted in Figure \[fig:logreg:approx\] (for $\theta_2$).

For each $(\widetilde{p}_{{\boldsymbol{Y}},k})_{k = 1\dots 5}$, we sample the posterior distribution using three methods.

-   [CBS]{}refers to a Classical Bridge Sampling $\pi(\theta)\ell({\boldsymbol{Y}}| \theta)^ {\rho_h}$, sequentially sampled with a SMC algorithm. The strategy serves as a reference to be compared to the other ones.

-   [CBS$+$IS]{}refers to the same annealing scheme as [CBS]{}but the first sample $(\theta_0^m)_{m=1,\dots,M}$ is generated with $\widetilde{p}_{{\boldsymbol{Y}}}$ and the adequate weights are computed.

-   Finally, we use [SBS]{}, described in section \[subsec:algo\] corresponding to the annealing scheme: $ {{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}( \pi({\theta}) \ell({\boldsymbol{Y}}| {\theta})) ^{\rho_h}$.

#### Tunings.

In each case, the SMC is performed using $M = 10000$ particles, $\tau_1 = 0.9$ and $\tau_2 = 0.8$. The kernel $K_h$ is a composed of $B = 5$ iteration of a standard Metropolis-Hastings kernel, proposing a new parameter as: $ \theta^c \sim \frac{1}{3} \sum_{i = 1}^3 \mathcal{N}(\theta^{\ell-1}, \rho_i \times \widehat{\Sigma}^{ML})$ where $\widehat{\Sigma}^{ML}$ is the asymptotic variance of the maximum likelihood estimator and $(\rho_1,\rho_2,\rho_3) = (1,0.1,10)$.

. \[fig: log reg post\]

#### Results.

The results are plotted in Figures \[fig: log reg post\] and \[fig:logreg:approx\] (right). For the five $(\widetilde{p}_{{\boldsymbol{Y}},k})_{k = 1 \dots 5}$, the posterior distributions given by our algorithm [SBS]{}(Figure \[fig: log reg post\], right) are confounded and stick to the one obtained from the reference algorithm [CBS]{}, thus illustrating the practical robustness of our new bridge sampler even. More precisely, even if the approximated posterior distribution has an underestimated variance (which is known to be the case for the Variational Bayes estimator in this case), our methodology will supply a sample from the true posterior distribution. Moreover, the algorithm is also robust when ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ is piked around on an absurd value (see results for $\widetilde{p}_{{\boldsymbol{Y}},5})$).

On the contrary, using algorithm [CBS$+$IS]{}can be a bad idea. Indeed, if the approximation $\widetilde{p}_{{\boldsymbol{Y}}}$ has a wider support than the true posterior, the algorithm will perform well and provide a sample from the right posterior distribution: as can be seen in Figure \[fig: log reg post\], left frame, the black and the blue curves are indistinguishable. However, if $\widetilde{p}_{{\boldsymbol{Y}}}$ has an underestimated variance ($\widetilde{p}_{{\boldsymbol{Y}},1}$, $\widetilde{p}_{{\boldsymbol{Y}},2}$, $\widetilde{p}_{{\boldsymbol{Y}},4}$ and $\widetilde{p}_{{\boldsymbol{Y}},5}$ ), the standard sampling strategy will lead to a false posterior distribution (see the orange, red, purple and green curves in the left frame of Figure \[fig: log reg post\]). Note that, there is no algorithmic indicator detecting such a bad behavior.

To compare the computational times of the various strategies, we can have a look at the number of iterations required for the sequences $(\rho_h)_{h \geq 0}$ to reach $1$. These sequences $(\rho_h)_{h \geq 0}$ are plotted in Figure \[fig:logreg:approx\] (right). We only plot the curves for the combinations “algorithm/approximated posterior $\widetilde{p}_{{\boldsymbol{Y}}}$” leading to the true posterior distribution. As expected, [CBS]{}is the most time consuming (see black dotted curve). As an indicative basis, on this example, the $30$ iterations of [CBS]{}requires roughly $5$ minutes on an using six cores. With $\widetilde{p}_{{\boldsymbol{Y}},3}$ (increased variance), [SBS]{}and [CBS$+$IS]{}finish in quite comparable computational time with slightly better results for our methodology. [SBS]{}requires the same computational time with $\widetilde{p}_{{\boldsymbol{Y}},4}$ (small variance). [SBS]{}clearly outperforms [CBS]{}. Finally, considering an extreme case where the approximation distribution $\widetilde{p}_{{\boldsymbol{Y}},5}$ is concentrated around an aberrant value (which is unlikely to be the case in practice), the [SBS]{}and [CBS]{}have comparable computation times.

As a conclusion, [CBS$+$IS]{}can not be used is general cases since it can supply a wrong approximation of the posterior distribution. This is due to the fact that, at the first iteration of [CBS$+$IS]{}, the sample must be of the prior distribution. Using IS, i.e. simulating with $\widetilde{p}_{{\boldsymbol{Y}},i}$ and assessing weight can lead to a sample from the truncated prior distribution $\pi({\theta})\mathbbm{1}_{\widetilde{p}_{{\boldsymbol{Y}}}({\theta})>0}$. This simulation from the wrong distribution at the first step of the sequential importance sampler is not corrected in the following iterations. Besides, we are not able to detect such a phenomena.

On the contrary, our new bridge sampler behaves well, whatever ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$. An under-evaluated variance in ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ is not a limit to the use of our scheme. The gain in computational time depends obviously on the distance between ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ and $p(\cdot | {\boldsymbol{Y}})$ but is expected to be drastic when ${{\widetilde{p}}_{{\boldsymbol{Y}}}}$ comes from deterministic approximations of the posterior (Variational Bayes, Expectation Propagation, etc).

Latent Class Analysis model {#subsec:LCA}
---------------------------

In this section we consider the latent class analysis (LCA) model. On this model, we focus on the label switching issue and show that the strategy we propose can tackle this difficulty.

### Model and prior distribution

LCA is a mixture model for multivariate binary observations such as the correct or incorrect answers submitted during an exam [@Barthol11], the symptoms presented by persons with major depressive disorder [@Garrett00] or a disability index recorded by long-term survey [@erosheva2007]. Let $({\boldsymbol{Y}}_i)_{i = 1,\dots, n} = (Y_{i1},\dots, Y_{iq})_{i = 1,\dots, n}$ be $n$ i.i.d observations where $\forall (i,j)$, $Y_{ij} \in \{0,1\}$. $i$ and $j$ are respectively the individual and the response indices. Each ${\boldsymbol{Y}}_i$ is assumed to arise from the following mixture model: $$\mathbb{P}({\boldsymbol{Y}}_i = (y_{i1},\dots, y_{iq})) = \sum_{k = 1}^g \pi_k \prod_{j = 1}^q \gamma_{kj}^{y_{ij}} (1-\gamma_{kj})^{1-y_{ij}}$$ where $\pi_k$ represents the proportion of the $k$-th component ($\sum_{k = 1}^g \pi_k = 1$) and $\gamma_{kj}$ is the success probability for the $j$-th response in the $k$-th group. The model is equivalently written as $$\begin{array}{cclr}
Y_{ij} | Z_i & \sim & \mathcal{B}(\gamma_{Z_i \, j}), & \quad \quad \forall i = 1,\dots n, j = 1,\dots, q\\
\mathbb{P}(Z_{i} = k) & = & \pi_k, & \quad \quad \forall i = 1,\dots n, k = 1,\dots, g
\end{array}$$ where ${\boldsymbol{Z}}= (Z_1,\dots, Z_n)$ is a latent random vector.

We set the following standard exchangeable prior distributions on ${\boldsymbol{\pi}}= (\pi_1,\dots, \pi_g)$ and ${\boldsymbol{\gamma}}= (\gamma_{kj})_{k = 1,\dots,q, j = 1,\dots,q}$: $$\label{eq:LCAprior}
\begin{array}{cccl}
(\pi_1,\dots, \pi_g) & \sim & & \mathcal{D}ir(d, \dots d)\\
(\gamma_{kj})_{k = 1,\dots,g, j = 1,\dots,q} & \sim & _{i.i.d} & \mathcal{B}eta(a,b) 
\end{array}$$

### Posterior distribution and label switching

As for any mixture model, the posterior distribution should reproduce the likelihood invariance under permutation of the mixture indices. In other words, the posterior distribution is multi-modal, each mode corresponding to a permutation of the index of the mixture components. In such cases, it is well documented that MCMC algorithms often fail into exploring the various modes of the posterior distribution. Note that the label switching issue arises not only when it comes out to sample the posterior distribution but also for evidence approximation in model selection [see for instance the introduction of @Lee2016 and references inside]. In this section, we illustrate the fact that a simple solution –based on the Variational Bayes posterior approximation and our sampling algorithm [SBS]{}– can be proposed to handle the label switching problem.

On this model, it is easy to derive a mean-field variational approximation of the posterior distribution, resulting into a posterior approximation of the form: $$\begin{aligned}
\label{eq:LCApostVB}
 &&{{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB}}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}) =f_{ \mathcal{D}ir(\tilde \delta_1, \dots \tilde \delta_g)}({\boldsymbol{\pi}})\nonumber \\
 && \times \prod_{k = 1,j=1}^{g,q}   f_{ \mathcal{B}eta(\tilde \alpha_{kj},\tilde \beta_{kj})} (\gamma_{kj})   \prod_{i = 1}^n \prod_{k = 1}^ q (\tau_{ik})^{\mathbbm{1}_{Z_{i} = k}} \end{aligned}$$ Details can be found in [@White13] and the algorithm is implemented in the corresponding R-package . However, contrary to the true posterior distribution $(p(\cdot | {\boldsymbol{Y}})$, ${{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB}}$ is not exchangeable. Moreover, this posterior approximation is known to be excessively concentrated around one mode. As a consequence, we can presume (and we illustrate it on the following numerical experiments) that our sampling algorithm [SBS]{}starting from $ {{\widetilde{p}}_{{\boldsymbol{Y}}}}$ and using a standard Gibbs transition kernel will not be able to propagate particles on the other modes of the posterior distribution.

When talking about a *“standard Gibbs transition kernel”*, we refer to the most naive Gibbs algorithm, sequentially sampling $[{\boldsymbol{Z}}| {\boldsymbol{Y}}, {\boldsymbol{\gamma}}, {\boldsymbol{\pi}}]$, $[{\boldsymbol{\gamma}}| {\boldsymbol{Y}}, {\boldsymbol{Z}}{\boldsymbol{\pi}}]$ and $[{\boldsymbol{\pi}}| {\boldsymbol{Z}}]$. Using the expressions for $ {{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB}}$ and $p_{\rho_n}({\boldsymbol{Z}},{\boldsymbol{\gamma}},{\boldsymbol{\pi}})$, these three distributions are conjugate. We stick to this MCMC kernel, and do not introduce any modification to *force or prevent* the label switching phenomena at the propagation step.

$ {{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB}}$ being unsatisfying from the label switching perspective, we introduce its so called *symmetrized version*, forcing the invariance by permutation. More precisely, let $\mathcal{S}_g$ be the set of all the permutations of $\{1,\dots,g\}$, we define $$\label{eq:LCApostVBsym0}
{{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB.Sym}}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}) = \sum_{\sigma \in \mathcal{S}_q} {{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB.Sym}}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}, \sigma)$$ where $$\begin{aligned}
\label{eq:LCApostVBsym}
 &&{{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB.Sym}}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}, \sigma) 
  = \frac{1}{g !} f_{ \mathcal{D}ir(\tilde \delta_{\sigma(1)}, \dots \tilde \delta_{\sigma(g)})}({\boldsymbol{\pi}})
 \nonumber\\
 & & \times \prod_{k = 1}^g \prod_{j = 1}^ q f_{ \mathcal{B}eta(\tilde \alpha_{\sigma(k)j},\tilde \beta_{\sigma(k)j})} (\gamma_{kj}) 
 \prod_{i = 1}^n \prod_{k = 1}^ q (\tau_{i\sigma(k)})^{\mathbbm{1}_{Z_{i} = k}}.\nonumber\\ \end{aligned}$$

In order to keep the conjugacy properties in conditional distributions we use our algorithm [SBS]{}to sequentially sample: $$\begin{aligned}
\label{eq:LCApostVBsym2}
 p_{\rho_h}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}, \sigma ) \propto [{{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB.Sym}}({\boldsymbol{Z}},{\boldsymbol{\gamma}}, {\boldsymbol{\pi}}, \sigma)]^ {1-\rho_h} && \nonumber  \\
  \times [ \ell({\boldsymbol{Y}}| {\boldsymbol{Z}}, {\boldsymbol{\gamma}},{\boldsymbol{\pi}})p({\boldsymbol{Z}}| {\boldsymbol{\pi}}) \pi({\boldsymbol{\pi}},{\boldsymbol{\gamma}})]^ {\rho_h}&&
 \end{aligned}$$ Once $\sigma$ has been integrated out, our sequential sampling scheme starts at $({{\widetilde{p}}_{{\boldsymbol{Y}}}}^{\mathit{VB.Sym}})$ given in equation and terminates at the true posterior distribution.

### Simulation design and test {#subsubsec:criteria}

This example aims at proving that our strategy supplies a sample from the true posterior distribution, even in the difficult framework of mixture models. To assess the validity of our method, we propose using the following testing procedure.

#### Testing procedure.

We introduce the following property from which we derive a negative criterion, in the sense that if the obtained distribution is the true posterior distribution, it must satisfy the criterion.

\[lem:criterion\] Let $\Phi: \Theta \mapsto \mathbb{R}$ be such that $\exists \; \Psi$ verifying $ H: {\theta}\mapsto \left(\Phi({\theta}), \Psi({\theta})\right)$ is injective and continuously differentiable.Assume that $$({\theta}^\star, {\boldsymbol{Y}}^{\star}) \sim \pi({\theta}^\star)\ell({\boldsymbol{Y}}^\star | {\theta}^\star).$$ Given ${\boldsymbol{Y}}^\star$, let $p_{{\boldsymbol{Y}}^\star}(\Phi(\theta))$ be any probability distribution on $\Phi(\Theta)$. Let $U({\theta}^\star, {\boldsymbol{Y}}^\star,\Phi,p_{{\boldsymbol{Y}}^\star})$ be the following statistic: $$\label{eq:crit-theo}
 U({\theta}^\star, {\boldsymbol{Y}}^\star,\Phi,p_{{\boldsymbol{Y}}^\star}) = \mathbb{E}_{p_{{\boldsymbol{Y}}^\star}} \left[ \mathbbm{1}_{\Phi({\theta}) < \Phi({\theta}^{\star})}\right].$$ Then, if $
p_{{\boldsymbol{Y}}^\star}(\Phi(\theta)) = p(\Phi({\theta})|{\boldsymbol{Y}}^\star)$, $  \forall {\theta}\in \Theta$ then $U({\theta}^\star, {\boldsymbol{Y}}^\star,\Phi,q_{{\boldsymbol{Y}}^\star}) \sim \mathcal{U}_{[0,1]}$.

Note that when $U({\theta}^\star, {\boldsymbol{Y}}^\star,\Phi, p_{{\boldsymbol{Y}}^\star})$ has no explicit expression and when we have access to a sample from $q_{{\boldsymbol{Y}}^\star}$, we can replace $U({\theta}^{\star}, {\boldsymbol{Y}}^{\star },\Phi, p_{{\boldsymbol{Y}}^{\star}})$ by its non-biased and convergent estimator: $$\label{eq:crit-emp}
U_M({\theta}^{\star},{\boldsymbol{Y}}^{\star},\Phi, ({\theta}_m)_{m = 1\dots M}) = \frac{1}{M}\sum_{m = 1}^M \mathbbm{1}_{\Phi(\theta_m) <\Phi({\theta}^{\star}) }$$ where ${\theta}_m \sim_{i.i.d.} p_{{\boldsymbol{Y}}^\star}$ and moreover if $p_{{\boldsymbol{Y}}^\star}(\Phi(\theta)) = p(\Phi({\theta})|{\boldsymbol{Y}}^\star)$ for all ${\theta}\in \Theta$, then $$U_M({\theta}^\star,{\boldsymbol{Y}}^\star,\Phi, ({\theta}_m)_{m = 1\dots M}) \sim \mathcal{U}_{\left\{0,\frac{1}{M},
\frac{2}{M},\dots, 1\right\}}.$$

This property enables the elaboration of our checking procedure in $4$ steps.

------------------------------------------------------------------------

**Checking procedure for posterior approximation**

------------------------------------------------------------------------

For a given approximation method ${\mathcal{M}}$ of the posterior distribution. Case (a) corresponds to deterministic approximations and case (b) to stochastic approximations.

1.  Generate $S$ parameters and datasets $({\theta}^{\star s}, {\boldsymbol{Y}}^{\star s})_{s = 1,\dots,S}$ according to the Bayesian model $\pi({\theta}^{\star s})\ell({\boldsymbol{Y}}^{\star s} | {\theta}^{\star s})$.

2.  For each $s = 1\dots S$, from dataset ${\boldsymbol{Y}}^{\star s}$,

    1.  derive the deterministic approximation of the posterior $p_{{\boldsymbol{Y}}^{\star s}}$ using ${\mathcal{M}}$,

    2.  get a sample $({\theta}^{s}_m)_{m = 1\dots M}$ from ${\mathcal{M}}$.

3.  Choose a real-valued function of the parameter $\Phi({\theta})$ and compute

    1.  $ U({\theta}^{\star s}, {\boldsymbol{Y}}^{\star s},\Phi,p_{{\boldsymbol{Y}}^{\star s}})$ using ,

    2.  $U_M({\theta}^{\star s},{\boldsymbol{Y}}^{\star s},\Phi, ({\theta}^{s}_m)_{m = 1\dots M})$ using .

4.  Compare the empirical distribution

    1.  of $U({\theta}^{\star s}, {\boldsymbol{Y}}^{\star s},\Phi,p_{{\boldsymbol{Y}}^{\star s}}))_{s=1\dots S}$ to the uniform distribution on $[0, 1]$,

    2.  of $\left(U_M({\theta}^{\star s},{\boldsymbol{Y}}^{\star s},\Phi, ({\theta}^{s}_m)_{m = 1\dots M}\right)_{s=1\dots S}$ to the uniform distribution on $\{1, \dots, \frac{1}{M}\}$.

------------------------------------------------------------------------

The comparison can be performed through graphical tools – for instance the empirical Cumulative Distribution Function (cdf) – or through a statistical test such as the discrete goodness of fit test, discrete version of the Kolmogorov-Smirnov (KS) test. If we reject the hypothesis $$H_0 = \left\{U_M({\theta}^{\star},{\boldsymbol{Y}}^{\star  },\Phi, ({\theta}_m)_{m = 1\dots M}) \sim \mathcal{U}_{\left\{0,\frac{1}{M}, \frac{2}{M},\dots, 1\right\}} \right\},$$ we conclude that the obtained sample is not distributed from the true posterior distribution.

In our case, several methods ${\mathcal{M}}$ will be considered: the variational Bayes approximation (VB), the [SBS]{}starting from VB and [SBS]{}starting from a symmetrized version of VB.

#### Simulation design.

We set $n = 100$ individuals, $p = 10$ observations by individual and $g = 2$ groups in the mixture. Referring to equation , the hyper-parameters are set to $\delta = a = b = 2$. $S = 500$ datasets ${\boldsymbol{Y}}^{\star s}$ and are simulated. For each simulated dataset, we run the VB algorithm implemented in the BayesLCA package [@White13]. The sampling algorithms [SBS]{}starting from ${\widetilde{p}}^\mathit{VB}_{{\boldsymbol{Y}}^{\star s} }$ and ${\widetilde{p}}_{{\boldsymbol{Y}}^{\star s} }^\mathit{VB.Sym}$ respectively are implemented with $M = 5000$, $\tau_1 = 0.9$ and $ \tau_2 = 0.9$. The $K_h$ kernel is standard Gibbs algorithm of length $B = 5$.

### Results

On Figure \[Fig:post\_pi1\] (left panel), we plot the estimated posterior density (for a arbitrarily chosen dataset) obtained from the Variational approximation, its symmetrized version, the [SBS]{}starting from with VB and the [SBS]{}stating from with VB symmetrized. As expected, the variational approximation underestimates the posterior variance. [SBS]{}succeeds in inflating this posterior variance. This phenomenon can be observed on $ {\Phi_1(\theta) =}|\pi_1-\pi_2|$ (left) and $ \Phi_2(\theta) = \pi_1$. About $\pi_1$, we observe that using only the classical variational Bayes estimator –highly concentrated on the MAP– we are enable to explore the other modes. However, if forcing the exploration of the different modes, we are able to observe the full posterior distribution, charging all the modes.

On Figure \[Fig:post\_pi1\] (right panel), we plot the empirical cdf’s (ecdf) for $U_M$. The results for $\Phi_1$ ($\Phi_2$ respectively) are on the left (right respectively). The ecdf for $\left(U({\theta}^{\star s}, {\boldsymbol{Y}}^{\star s},\Phi_l, {\widetilde{p}}^{\mathit{VB}}_{{\boldsymbol{Y}}^{\star s}})\right)_{s = 1\dots S}$ is in blue and the ecdf for $\left(U({\theta}^{\star s}, {\boldsymbol{Y}}^{\star s},\Phi_l, {\widetilde{p}}^{\mathit{VB Sym}}_{{\boldsymbol{Y}}^{\star s}})\right)_{s = 1\dots S}$ is in purple. The ecdf obtained from $({\theta}^{s}_m)_{m = 1\dots M}$ and $({\theta}^{s, \mathit{Sym}}_m)_{m = 1\dots M}$ sampled by algorithm [SBS]{}starting from ${\widetilde{p}}^\mathit{VB}_{{\boldsymbol{Y}}^{\star s} }$ ( respectively ${\widetilde{p}}_{{\boldsymbol{Y}}^{\star s} }^\mathit{VB.Sym}$) are plotted in red (green respectively).

The phenomena observed before on a single dataset are confirmed here on the $500$ datasets. On $|\pi_1-\pi_2|$, which is insensitive to label switching, the non-symmetrized and symmetrized versions of the algorithms give equivalent results (curves red/green, and blue/purple can not be distinguished), which is not the case on $\pi_1$ (which is actually non-identifiable because of label switching). On $|\pi_1-\pi_2|$, ${\widetilde{p}}^\mathit{VB}_{{\boldsymbol{Y}}^{\star s} }$ and ${\widetilde{p}}_{{\boldsymbol{Y}}^{\star s} }^\mathit{VB.Sym}$ are different from the true posterior distribution but our two algorithms starting respectively from ${\widetilde{p}}^\mathit{VB}_{{\boldsymbol{Y}}^{\star s} }$ and ${\widetilde{p}}_{{\boldsymbol{Y}}^{\star s} }^\mathit{VB.Sym}$ supply a sample from the true posterior distribution.

On $\pi_1$, we observe that starting [SBS]{}from ${\widetilde{p}}^\mathit{VB}_{{\boldsymbol{Y}}^{\star s} }$ clearly leads to the wrong posterior. The equality of the true posterior distribution with the one obtained via [SBS]{}starting from VB.Sym is not rejected.

[lcc]{} & $|\pi_1-\pi_2|$ & $\pi_1$\
\
VB & $4.497e^{-06}$ & $<2.2e^{-16}$\
Symmetrized VB & $1.563e^{-05}$ & $1.431e^{-09}$\
[SBS]{}with VB & $0.596$ & $<2.2e^{-16}$\
[SBS]{}with Symmetrized VB& $0.567$ & $0.903$\

Stochastic block models with covariates {#subsec:SBMreg}
---------------------------------------

#### Model.

As a last example, we consider the combination of stochastic block-model [@NoS01 SBM] and logistic regression (shortened as ’[SBM-reg]{}’ in the sequel) considered in [@LRO15]. This model aims at deciphering some residual structure in an observed network once accounted for the effect of some edge covariates. The model is as follows. Consider a set of $n$ nodes; for each pair ($1 \leq i < j \leq n$) of nodes, we observe a $p$-dimensional covariates vector ${\boldsymbol{x}}_{ij}$. Likewise in SBM, we further assume that each node belongs to one among $g$ groups and we denote $Z_i$ the (unobserved) group where node $i$ is affected; ${\boldsymbol{\pi}}= (\pi_k)_k$ denotes the vector of group proportions. The model states that the edges of the observed binary undirected network ${\boldsymbol{Y}}= (Y_{ij})$ are drawn independently conditionally on the set of latent variables ${\boldsymbol{Z}}= (Z_i)$ as Bernoulli variables: $$(Y_{ij} | Z_i, Z_j, {\boldsymbol{\alpha}}, {\boldsymbol{\beta}}) \sim {\mathcal{B}}(p_{ij}), 
\quad 
{\text{logit}}(p_{ij}) = {\boldsymbol{x}}_{ij}^\intercal {\boldsymbol{\beta}}+ \alpha_{Z_i, Z_j}$$ where ${\boldsymbol{\alpha}}= (\alpha_{kl})$ stands for the matrix of between-group effects (analogous to the between-group connection probabilities from SBM, in logit scale) and ${\boldsymbol{\beta}}= (\beta_\ell)_{\ell =1,\dots, p}$ for the vector of regression coefficients. As for the priors, $\pi$ has a Dirichlet distribution, both $\alpha$ and $\beta$ are Gaussian. When considering model selection or averaging, the number of groups $g$ is supposed to be uniformly distributed among $\{1, \dots, g_{\max}\}$.

#### Bayesian model averaging (BMA).

BMA [@HMR99] is a general principle, which consists in combining the results obtained with several models, rather than to choose the ’best’ one. Among other interests, it allows to account for model uncertainty. We apply this principle to a regression parameter $\beta_\ell$. While model selection consists in choosing $g$ as $\widehat{g} = \arg\max_g p(g | {\boldsymbol{Y}})$ and considering the posterior $p(\beta_\ell | {\boldsymbol{Y}}, g = \widehat{g})$, BMA directly considers the unconditional posterior $$p(\beta_\ell | {\boldsymbol{Y}}) = \sum_g p(g | {\boldsymbol{Y}}) p(\beta_\ell | {\boldsymbol{Y}}, g).$$ In terms of moments, it results in ${\mathbb E}(\beta_\ell | {\boldsymbol{Y}}) = \sum_g p(g | {\boldsymbol{Y}}) {\mathbb E}(\beta_\ell | {\boldsymbol{Y}}, g)$ and ${\mathbb V}(\beta_\ell | {\boldsymbol{Y}}) = {\mathbb V}_{\text{within}}(\beta_\ell | {\boldsymbol{Y}}) + {\mathbb V}_{\text{between}}(\beta_\ell | {\boldsymbol{Y}})$ where ${\mathbb V}_{\text{within}}$ measures the mean variance of the parameter conditionally on $g$ and ${\mathbb V}_{\text{between}}$ is the variance of the parameter due to model uncertainty: $$\begin{aligned}
{\mathbb V}_{\text{within}}(\beta_\ell | {\boldsymbol{Y}}) & = & \sum_g p(g | {\boldsymbol{Y}}) {\mathbb V}(\beta_\ell | {\boldsymbol{Y}}, g), \\
{\mathbb V}_{\text{between}}(\beta_\ell | {\boldsymbol{Y}}) & = & \sum_g p(g | {\boldsymbol{Y}}) \left({\mathbb E}(\beta_\ell | {\boldsymbol{Y}}, g) - {\mathbb E}(\beta_\ell | {\boldsymbol{Y}})\right)^2.\end{aligned}$$

#### VB approximation.

As illustrated in Subsection \[subsec:SBMreg\], the VB approximate posterior is quite accurate for logistic regression and [@GDR12] also proved its empirical accuracy for SBM. A first goal of this simulation study is to check if this accuracy still holds when the two models are combined into the [SBM-reg]{}model. To this aim, we focus on the posterior distribution of the regression parameters. Secondly, we want to check the accuracy of the VB posterior distribution of the number of groups, that can be used either to assess goodness-of-fit or for model averaging ([@LRO15])

#### Simulation design.

We simulate networks with $n \in \{20, 50\}$ nodes according to an [SBM-reg]{}model with ${g_*}\in \{1, 2\}$ groups and $p = 3$ covariates. To apply Property \[lem:criterion\], the parameters are sampled from the prior distribution. $S = 100$ replicates are simulated for each configuration and for each of them , the [SBM-reg]{}models with $g \in \{1, \dots, g_ {\max} = 5\}$ were fitted with the VB algorithm described in [@LRO15]. The [SBS]{}algorithm is then run on each dataset.

#### Results for parameter estimation.

We first consider the posterior distribution of $\beta$ when the number of groups $g$ is known. On Figure \[Fig:SimulSMBreg-beta\], we plot on the left the boxplots for the posterior means $(\widehat{\mathbb{E}}^{VB}(\beta_\ell | {\boldsymbol{Y}}^s))_{s=1\dots 100}$ and $(\widehat{\mathbb{E}}^{{SBS\xspace}}(\beta_\ell | {\boldsymbol{Y}}^s))_{s=1\dots 100}$. The boxplot (over the 100 simulated datasets) of the posterior standard deviations $(\widehat{\sigma}^{VB}(\beta_\ell | {\boldsymbol{Y}}^s))_{s=1\dots 100}$ and $(\widehat{\sigma}^{{SBS\xspace}}(\beta_\ell | {\boldsymbol{Y}}^s))_{s=1\dots 100}$ are on the top-right. We clearly observe that the posterior means provided by VB and [SBS]{}are both accurate and similar, but the VB’s posterior standard deviations (sd) are smaller than [SBS]{}’s posterior standard deviations.

To further assess the quality of the posterior distribution provided by VB and [SBS]{}, we checked Property \[lem:criterion\] for the regression coefficients, which are not subject to label switching. We observe that the ecdf of the [SBS]{}sample is close to uniform, whereas there is a departure for VB. The p-values resulting from KS test of Property \[lem:criterion\] (see Table \[Tab:pvalSBMreg\], upper table) lead to the same conclusions. All these observations concur to show that, although the VB approximate posterior distribution is accurate for logistic regression and SBM separately, it is biased for the [SBM-reg]{}model, and that the proposed [SBS]{}is a way to correct it. As a consequence of this phenomenon, the empirical level of VB’s credibility intervals is equal to 84.75%, which is below the nominal level 95%, whereas [SBS]{}’s credibility intervals almost reach the targeted level (93.75%).

  --------- --------- --------- --------- ---------
                                          
             $g = 1$   $g = 2$   $g = 1$   $g = 2$
     VB       0.027     0.004     0.002     0.077
   [SBS]{}    0.785     0.121     0.839     0.238
    $\;$      $\;$      $\;$      $\;$    
    $\;$      $\;$      $\;$      $\;$    
                                          
                                          
             $g = 1$   $g = 2$   $g = 1$   $g = 2$
     VB       0.017     0.003     0.002     0.079
   [SBS]{}    0.740     0.277     0.778     0.312
  --------- --------- --------- --------- ---------

  : Simulation results for the [SBM-reg]{}model: p-values for the KS test of Property \[lem:criterion\]. \[Tab:pvalSBMreg\]

   ![Simulation results for the [SBM-reg]{}model: VB (white) and [SBS]{}(red) posterior of the regression coefficients ${\boldsymbol{\beta}}= (\beta_\ell)$. Top: posterior mean (left), posterior standard deviation (right); $x$-axis label: ${g_*}. n$ (e.g. ’1.20’ means ${g_*}= 1$, $n= 20$). Bottom: graphical check of Property \[lem:criterion\] for VB (dashed blue) and for [SBS]{}(solid red). Left: $g$ = ${g_*}$, right: with model averaging. []{data-label="Fig:SimulSMBreg-beta"}](Fig7.eps "fig:"){width=".4\columnwidth" height=".4\textwidth"}   ![Simulation results for the [SBM-reg]{}model: VB (white) and [SBS]{}(red) posterior of the regression coefficients ${\boldsymbol{\beta}}= (\beta_\ell)$. Top: posterior mean (left), posterior standard deviation (right); $x$-axis label: ${g_*}. n$ (e.g. ’1.20’ means ${g_*}= 1$, $n= 20$). Bottom: graphical check of Property \[lem:criterion\] for VB (dashed blue) and for [SBS]{}(solid red). Left: $g$ = ${g_*}$, right: with model averaging. []{data-label="Fig:SimulSMBreg-beta"}](Fig8.eps "fig:"){width=".4\columnwidth" height=".4\textwidth"}
  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Simulation results for the [SBM-reg]{}model: VB (white) and [SBS]{}(red) posterior of the regression coefficients ${\boldsymbol{\beta}}= (\beta_\ell)$. Top: posterior mean (left), posterior standard deviation (right); $x$-axis label: ${g_*}. n$ (e.g. ’1.20’ means ${g_*}= 1$, $n= 20$). Bottom: graphical check of Property \[lem:criterion\] for VB (dashed blue) and for [SBS]{}(solid red). Left: $g$ = ${g_*}$, right: with model averaging. []{data-label="Fig:SimulSMBreg-beta"}](Fig9.eps "fig:"){width=".4\columnwidth" height=".4\textwidth"}   ![Simulation results for the [SBM-reg]{}model: VB (white) and [SBS]{}(red) posterior of the regression coefficients ${\boldsymbol{\beta}}= (\beta_\ell)$. Top: posterior mean (left), posterior standard deviation (right); $x$-axis label: ${g_*}. n$ (e.g. ’1.20’ means ${g_*}= 1$, $n= 20$). Bottom: graphical check of Property \[lem:criterion\] for VB (dashed blue) and for [SBS]{}(solid red). Left: $g$ = ${g_*}$, right: with model averaging. []{data-label="Fig:SimulSMBreg-beta"}](Fig10.eps "fig:"){width=".4\columnwidth" height=".4\textwidth"}

#### Results for model selection.

We now consider the posterior distribution of the number of groups $p(g|{\boldsymbol{Y}})$ and its use for model selection. Figure \[Fig:SimulSMBreg-modsel\] provides a comparison of the posterior provided by VB and [SBS]{}. We observe that the VB approximation always results in a more concentrated distribution than [SBS]{}. This behavior can be compared to the under-estimation of the posterior variance of the parameters that we already discussed. To compare the results in terms of model selection we computed the frequency at which the right model is selected (i.e. when $\widehat{g} = {g_*}$) and the mean posterior probability of the ${g_*}$ (see Table \[Tab:SimulSMBreg-Pgstar\]). We observe that VB performs better than [SBS]{}for both criteria. This parallels [@Min05], who shows that the minimization of the Kullback-Leibler (KL) divergence leads to an accurate estimate of the mode, which is convenient for model selection.

   ![Simulation results for the [SBM-reg]{}model: box-plots for the posterior probability $p(g|{\boldsymbol{Y}})$ as a function of $g$. Top $n=20$, bottom: $n=50$. Left: ${g_*}= 1$, right: ${g_*}= 2$. \[Fig:SimulSMBreg-modsel\]](Fig11.eps "fig:"){width=".4\columnwidth" height=".3\textwidth"}   ![Simulation results for the [SBM-reg]{}model: box-plots for the posterior probability $p(g|{\boldsymbol{Y}})$ as a function of $g$. Top $n=20$, bottom: $n=50$. Left: ${g_*}= 1$, right: ${g_*}= 2$. \[Fig:SimulSMBreg-modsel\]](Fig12.eps "fig:"){width=".4\columnwidth" height=".3\textwidth"}
  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Simulation results for the [SBM-reg]{}model: box-plots for the posterior probability $p(g|{\boldsymbol{Y}})$ as a function of $g$. Top $n=20$, bottom: $n=50$. Left: ${g_*}= 1$, right: ${g_*}= 2$. \[Fig:SimulSMBreg-modsel\]](Fig13.eps "fig:"){width=".4\columnwidth" height=".3\textwidth"}   ![Simulation results for the [SBM-reg]{}model: box-plots for the posterior probability $p(g|{\boldsymbol{Y}})$ as a function of $g$. Top $n=20$, bottom: $n=50$. Left: ${g_*}= 1$, right: ${g_*}= 2$. \[Fig:SimulSMBreg-modsel\]](Fig14.eps "fig:"){width=".4\columnwidth" height=".3\textwidth"}

  --------- ----------- ----------- ----------- -----------
                                                
             ${g_*}=1$   ${g_*}=2$   ${g_*}=1$   ${g_*}=2$
     VB         100         10          100         42
   [SBS]{}      46          23          60          36
    $\;$       $\;$        $\;$        $\;$     
    $\;$       $\;$        $\;$        $\;$     
                                                
                                                
             ${g_*}=1$   ${g_*}=2$   ${g_*}=1$   ${g_*}=2$
     VB        0.947       0.138       0.982       0.410
   [SBS]{}     0.435       0.257       0.562       0.387
  --------- ----------- ----------- ----------- -----------

  : Simulation results for the [SBM-reg]{}model: model selection \[Tab:SimulSMBreg-Pgstar\]

Although it does not seem to hamper model selection, the biased estimation of the posterior $p(g|{\boldsymbol{Y}})$ may have undesired consequences when used for model averaging. To illustrate this point, we simply computed the empirical coverage of credibility intervals for each $\beta_\ell$ after model averaging. The mean coverage across simulation condition and covariate index $\ell$ for VB (85.8%) is still below the nominal level, whereas this of [SBS]{}(93.25%) is close to 95%. Figure \[Fig:SimulSMBreg-beta\] (bottom right) also shows that the distribution of the ecdf after model averaging is almost confounded with the uniform for [SBS]{}, whereas it still displays a significant bias for VB. The p-values for the KS test of Property \[lem:criterion\] (see Table \[Tab:pvalSBMreg\], bottom) lead to the same conclusion.

Illustrations on network datasets {#sec:Illust}
=================================

#### Network analysis with [SBM-reg]{}.

To illustrate the use of the proposed sampling algorithm, studied a series of examples analyzed by [@LRO15] with an [SBM-reg]{}model. The main two purposes of such an analysis is ($i$) to estimate the effect ${\boldsymbol{\beta}}$ of the covariates and ($ii$) to assess the goodness-of-fit of the model based on the covariates. Task ($ii$) is achieved by computing the posterior probability for the SBM part of the model to involve only $g=1$ class, that is $p(g = 1 | {\boldsymbol{Y}})$. A low value of this probability is an indication for a residual structure in the network.

We refer to [@LRO15] for the presentation of the data. We considered the datasets (networks) refereed to as Florentine (business), Florentine (marriage), Karate, Tree and Blog. Their respective sizes range from few tens to few hundreds nodes and their densities from 1% to 50%. Note that the numerical results presented here for the VB inference slightly differ, as we kept all nodes from each graph whereas [@LRO15] removed all isolated nodes.

For each of these datasets, we fitted an [SBM-reg]{}model with $g = 1 \dots g_{\max}$ groups with a VBEM algorithm to obtain the Gaussian VB approximate distribution for ${{\widetilde{p}}_{{\boldsymbol{Y}}}}^{VB}({\boldsymbol{\beta}})$. We also run the proposed [SBS]{}with $M = 1000$ particles and obtained a weighted sample from $\widehat{p}_{{\boldsymbol{Y}}}^{{SBS\xspace}}({\boldsymbol{\beta}})$. To compare the posterior distributions of ${\boldsymbol{\beta}}$, we adopted the Bayesian model averaging principle described in Section \[subsec:SBMreg\]. For each dataset $s$ and each covariate $\ell$, we first computed the ratio between the VB and [SBS]{}posterior standard deviation (sd) $\sqrt{{\mathbb V}^{VB}(\beta_j) / {\mathbb V}^{{SBS\xspace}}(\beta_j)}$. For each configuration, we also computed the ratios $\widetilde{{\boldsymbol{\beta}}}^{VB}_j / \sqrt{\widetilde{\Sigma}^{VB}_{jj}}$ and $\widetilde{{\boldsymbol{\beta}}}^{{SBS\xspace}}_j / \sqrt{\widetilde{\Sigma}^{{SBS\xspace}}_{jj}}$, which are typically used to evaluate the effect of the covariates. Figure \[Fig:IllustPostBeta\], left, shows that the VB approximation tends to under estimate the posterior variance of the parameter. As expected, Figure \[Fig:IllustPostBeta\], right, shows that it yields in over-estimating the significance of the effect of the covariates.

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Posterior moments of the regression coefficients for the datasets from [@LRO15]. Left: ratio between the VB and [SBS]{}sd as a function of the network size (color = network). Right: ratio between the posterior mean and the posterior sd ($x$: VB, $y$: [SBS]{}, signed-log scale, color: network, horizontal and vertical lines: 2.5% and 97.5% $\mathcal{N}(0, 1)$ quantiles) \[Fig:IllustPostBeta\]](Fig15.eps "fig:"){width=".4\columnwidth"}   ![Posterior moments of the regression coefficients for the datasets from [@LRO15]. Left: ratio between the VB and [SBS]{}sd as a function of the network size (color = network). Right: ratio between the posterior mean and the posterior sd ($x$: VB, $y$: [SBS]{}, signed-log scale, color: network, horizontal and vertical lines: 2.5% and 97.5% $\mathcal{N}(0, 1)$ quantiles) \[Fig:IllustPostBeta\]](Fig16.eps "fig:"){width=".4\columnwidth"}
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

To further illustrate these differences, we studied the posterior distribution of the three regression coefficients used in the ’Tree’ example. The regression coefficient are associated with the genetic, geographic and taxonomic distance between the tree species. The results given in Table \[Tab:IllustPostBeta\] indicate again that the posterior sd provided by [SBS]{}are all larger than these resulting from VB. We observe that, because of the high concentration of the posterior distribution of $g$, VB strongly under-estimates the variance due to model uncertainty. A result is that the [SBS]{}posterior standard deviations are almost twice larger than the VB posterior standard deviations. As a consequence, the significance ratio is reduced, and the influence of both the genetic and the geographic distance ($\beta_1$ and $\beta_2$) turn out to be more questionable according to [SBS]{}than to VB.

  -------------------- ----------- ----------- -----------
                        $\beta_1$   $\beta_2$   $\beta_3$
  post. mean            4.62e-05      0.23        -0.9
  post. within var.     2.24e-10     0.0432      0.00175
  post. between var.    5.55e-17    1.18e-06    2.42e-07
  posterior sd           1.5e-05      0.208      0.0418
  ratio                   3.09        1.11        -21.5
  $\;$                    $\;$        $\;$        $\;$
  $\;$                    $\;$        $\;$        $\;$
                                               
                        $\beta_1$   $\beta_2$   $\beta_3$
  post. mean            4.13e-05      0.355      -0.906
  post. within var.     1.09e-09      0.219      0.00889
  post. between var.    3.99e-12     0.0019      0.00281
  posterior sd          3.31e-05      0.47        0.108
  ratio                   1.25        0.755       -8.38
  -------------------- ----------- ----------- -----------

  : Posterior moments of the regression coefficients. \[Tab:IllustPostBeta\]

For the goodness-of-fit study, we compared the values of $\widetilde{p}_{{\boldsymbol{Y}}}^{VB}(1)$ and $\widehat{p}_{{\boldsymbol{Y}}}^{{SBS\xspace}}(1)$ (Table \[Tab:IllustPostPrM1\]). Except in the most uncertain case (Karate), the posterior probabilities are similar and lead to the same conclusion about the existence of a residual structure in the network.

              $\widetilde{p}_{{\boldsymbol{Y}}}^{VB}(1)$   $\widehat{p}_{{\boldsymbol{Y}}}^{{SBS\xspace}}(1)$
  ---------- -------------------------------------------- ----------------------------------------------------
  Marriage                $9.54 \; 10^{-1}$                                      $1.00$
  Business                $7.04 \; 10^{-1}$                                      $1.00$
  Karate                  $2.56 \; 10^{-1}$                                $7.07 \; 10^{-3}$
  Tree                   $4.83 \; 10^{-153}$                              $1.06 \; 10^{-161}$
  Blog                   $8.63 \; 10^{-174}$                              $4.04 \; 10^{-290}$

  : Posterior probability for the [SBM-reg]{}model with only one class. \[Tab:IllustPostPrM1\]

Discussion and perspectives {#sec:Discuss}
===========================

In this paper, we presented a simple strategy to combine the strength of deterministic approximations of the posterior distribution with sequential Monte Carlo samplers. We illustrated the efficiency of our approach and its robustness with respect to the deterministic approximation on a large simulation study. Its application on network datasets stresses the fact that the well-known underestimation of the posterior variance by the variational approximation can be easily corrected, sometimes leading to different statistical conclusions. Besides, if dependencies between parameters have been neglected in the deterministic posterior approximation, they will be recovered by the sequential sampling.

Our approach is not restricted to the case where a standard deterministic posterior approximation can be derived (such as Variational Bayes, Laplace or Expectation Propagation estimate). Any point estimate can be used to design a rough posterior (using a Gaussian or a log-Gaussian seems to be the simplest solution) and serves as an accelerator of the sampling sequence. This strategy is different from an empirical Bayes strategy, the point estimate being only used to explore more efficiently the posterior distribution and not to elicit a prior distribution. The method is not as sensible as standard Importance Sampling to an eventual under-evaluation of the approximate posterior variance : even with a too narrow approximation of the posterior distribution, the algorithm is able to get back to the true posterior variance.

SMC directly supplies a final population of particles arising from the true posterior distribution, as opposed to MCMC strategies, whose convergence is difficult to assess. The proposed SBS algorithm is adaptive in the sense that the sequence ${{\widetilde{p}}_{{\boldsymbol{Y}}}}({\theta})^ {1-{\rho_h}}(p( {\theta}| {\boldsymbol{Y}})) ^{\rho_h}$ is determined on the fly in an automatic way. Furthermore, the algorithm path (summarized by the sequence $\rho_h$) is an indicator of the quality of the deterministic posterior distribution used to initiate the bridge sampling.

A natural extension of the present work is its adaptation to Approximate Bayesian Computation (ABC) context for models with no explicit likelihood, following [@DelMoral2012]. The difficulty will arise from the specification of the distributions sequence.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors are thankful to Nicolas Chopin for fruitful discussions.
