Global Consensus Monte Carlo

The goal of the present paper is to approximate (1). We take an approach that has also been developed in contemporaneous work by Vono et al., 2019a , although their objectives were somewhat different. Alongside the variable of interest $Z$, we introduce a collection of $b$ instrumental variables each also defined on $\mathsf{E}$, denoted by $X_{1:b}$. On the extended state space $\mathsf{E}\times\mathsf{E}^{b}$, we define the probability density function $\tilde{\pi}_{\lambda}$ by

where for each $j\in\{1,\dots,b\}$, $\{K_{j}^{(\lambda)}:\lambda\in\mathbb{R}_{+}\}$ is a family of Markov transition densities on $\mathsf{E}$. Defining

the density of the $Z$-marginal of $\tilde{\pi}_{\lambda}$ may be written as

Here, we may view each $f_{j}^{(\lambda)}$ as a smoothed form of $f_{j}$, with $\pi_{\lambda}$ being the corresponding smoothed form of the target density (1).

The role of $\lambda$ is to control the fidelity of $f_{j}^{(\lambda)}$ to $f_{j}$, and so we assume the following in the sequel.

For example, Assumption 1 implies that $\pi_{\lambda}$ converges in total variation to $\pi$ by Scheffé’s lemma (Scheffé,, 1947), and therefore $\pi_{\lambda}(\varphi)\to\pi(\varphi)$ for bounded $\varphi:\mathsf{E}\to\mathbb{R}$. Assumption 1 is satisfied for essentially any kernel that may be used for kernel density estimation; here, $\lambda$ takes a similar role to that of the smoothing bandwidth.

On a first reading one may assume that the $K_{j}^{(\lambda)}$ are chosen to be independent of $j$; for example, with $\mathsf{E}=\mathbb{R}$ one could take $K_{j}^{(\lambda)}(z,x)=\mathcal{N}(x;z,\lambda)$. We describe some considerations in choosing these transition kernels in Section LABEL:subsec:choosingTransitionDensities of the supplement, and describe settings in which choosing these to differ with $j$ may be beneficial.

The instrumental hierarchical model is presented diagrammatically in Figure 1. The variables $X_{1:b}$ may be seen as ‘proxies’ for $Z$ associated with each of the data subsets, which are conditionally independent given $Z$ and $\lambda$. Loosely speaking, $\lambda$ represents the extent to which we allow the local variables $X_{1:b}$ to differ from the global variable $Z$. In terms of computation, it is the separation of $Z$ from the subsets of the data $\mathbf{y}_{1:b}$, given $X_{1:b}$ introduced by the instrumental model, that can be exploited by distributed algorithms.

This approach to constructing an artificial joint target density is easily extended to accommodate random effects models, in which the original statistical model itself contains local variables associated with each data subset. These variables may be retained in the resulting instrumental model, alongside the local proxies $X_{1:b}$ for $Z$. A full description of the resulting model is presented in Rendell, (2020).

The framework we describe is motivated by concepts in distributed optimisation, a connection that is also explored in the contemporaneous work of Vono et al., 2019a . The global consensus optimisation problem (see, e.g., Boyd et al.,, 2011, Section 7) is that of minimising a sum of functions on a common domain, under the constraint that their arguments are all equal to some global common value. If for each $j\in\{1,\ldots,b\}$ one uses the Gaussian kernel density $K_{j}^{(\lambda)}(z,x)=\mathcal{N}(x;z,\lambda)$, then taking the negative logarithm of (2) gives

where $C$ is a normalising constant. Maximising $\pi(z)$ is equivalent to minimising this function under the constraint that $z=x_{j}$ for $j\in\{1,\dots,b\}$, which may be achieved using the alternating direction method of multipliers (Bertsekas and Tsitsiklis,, 1989). Specifically, (5) corresponds to the use of $1/\lambda$ as the penalty parameter in this procedure.

There are some similarities between this framework and Approximate Bayesian Computation (ABC; see Marin et al.,, 2012, for a review of such methods). In both cases one introduces a kernel that can be viewed as acting to smooth the likelihood. In the case of (2) the role of $\lambda$ is to control the scale of smoothing that occurs in the parameter space; the tolerance parameter used in ABC, in contrast, controls the extent of a comparable form of smoothing in the observation (or summary statistic) space.

The instrumental model described forms the basis of our proposed global consensus framework; ‘global consensus Monte Carlo’ is correspondingly the application of Monte Carlo methods to form an approximation of $\pi_{\lambda}$. We focus here on the construction of a Metropolis-within-Gibbs Markov kernel that leaves $\tilde{\pi}_{\lambda}$ invariant. If $\lambda$ is chosen to be sufficiently small, then the $Z$-marginal $\pi_{\lambda}$ provides an approximation of $\pi$. Therefore given a chain with values denoted $(Z^{i},X_{1:b}^{i})$ for $i\in\{1,\ldots,N\}$, an empirical approximation of $\pi$ is given by

where $\delta_{z}$ denotes the Dirac measure at $z$.

The Metropolis-within-Gibbs kernel we consider utilises the full conditional densities

for $j\in\{1,\ldots,b\}$, and

where (7) follows from the mutual conditional independence of $X_{1:b}$ given $Z$. Here we observe that $K_{j}^{(\lambda)}(z,x_{j})$ simultaneously provides a pseudo-prior for $X_{j}$ and a pseudo-likelihood for $Z$.

We define $M_{1}^{(\lambda)}$ to be a $\tilde{\pi}_{\lambda}$-invariant Markov kernel that fixes $z$; we may write

where for each $j$, $P_{j,z}^{(\lambda)}(x_{j},\cdot)$ is a Markov kernel leaving (7) invariant. We similarly define $M_{2}^{(\lambda)}$ to be a $\tilde{\pi}_{\lambda}$-invariant Markov kernel that fixes $x_{1:b}$,

where $P_{x_{1:b}}^{(\lambda)}(z,\cdot)$ is a Markov kernel leaving (8) invariant.

Using these Markov kernels we construct an MCMC kernel that leaves $\tilde{\pi}_{\lambda}$ invariant; we present the resulting sampling procedure as Algorithm 1. In the special case in which one may sample exactly from the conditional distributions (7)–(8), this algorithm takes the form of a Gibbs sampler.

The interest from a distributed perspective is that the full conditional density (7) of each $X_{j}$, for given values $x_{j}$ and $z$, depends only on the $j$th block of data (through the partial likelihood $f_{j}$) and may be computed on the $j$th machine. Within Algorithm 1, the sampling of each $X_{j}^{i}$ from $P_{j,Z^{i-1}}^{(\lambda)}(X_{j}^{i-1},\cdot)$ may therefore occur on the $j$th machine; these $X_{1:b}^{i}$ may then be communicated to a central machine that draws $Z^{i}$.

Our approach has particular benefits when sampling exactly from (7) is not possible, in which case Algorithm 1 takes a Metropolis-within-Gibbs form. One may choose the Markov kernels $P_{j,z}^{(\lambda)}$ to comprise multiple iterations of an MCMC kernel leaving (7) invariant; indeed multiple computations of each $f_{j}$ (and therefore multiple accept/reject steps) may be conducted on each of the $b$ nodes, without requiring communication between machines. This stands in contrast to more straightforward MCMC approaches directly targeting $\pi$, in which such communication is required for each evaluation of (1), and therefore for every accept/reject step. Similar to such a ‘direct’ MCMC approach, each iteration of Algorithm 1 requires communication to and from each of the $b$ machines on which the data are stored; but in cases where the communication latency is high, the resulting sampler will spend a greater proportion of time exploring the state space. This may in turn result in faster mixing (e.g. with respect to wall-clock time). We further discuss this comparison, and the role of communication latency, in Section LABEL:subsec:repeatedMCMCkernels of the supplement.

The setting in which Algorithm 1 describes a Gibbs sampler, with all variables drawn exactly from their full conditional distributions, is particularly amenable to analysis. We provide such a study in Section LABEL:sec:Gaussians of the supplement. This analysis may also be informative about the more general Metropolis-within-Gibbs setting, when $P_{j,z}^{(\lambda)}$ comprises enough MCMC iterations to exhibit good mixing.

For $\varphi:\mathsf{E}\to\mathbb{R}$ we may estimate $\pi(\varphi)$ using (6) as $\pi_{\lambda}^{N}(\varphi)$. The regularisation parameter $\lambda$ here takes the role of a tuning parameter; we can view its effect on the mean squared error of such estimates using the bias–variance decomposition,

which holds exactly when $\operatorname{\mathbb{E}}[\pi_{\lambda}^{N}(\varphi)]=\pi_{\lambda}(\varphi)$. In many practical cases this decomposition will provide a very accurate approximation for large $N$, as the squared bias of $\pi_{\lambda}^{N}(\varphi)$ is typically asymptotically negligible in comparison to its variance.

The decomposition (11) separates the contributions to the error from the bias introduced by the instrumental model and the variance associated with the MCMC approximation. If $\lambda$ is too large, the squared bias term in (11) can dominate while if $\lambda$ is too small, the Markov chain may exhibit poor mixing due to strong conditional dependencies between $X_{1:b}$ and $Z$, and so the variance term in (11) can dominate.

It follows that $\lambda$ should ideally be chosen in order to balance these two considerations. We provide a theoretical analysis of the role of $\lambda$ in Section LABEL:sec:Gaussians of the supplement, by considering a simple Gaussian setting. In particular we show that for a fixed number of data, one should scale $\lambda$ with the number of samples $N$ as $\mathcal{O}(N^{-1/3})$ in order to minimise the mean squared error. We also consider the consistency of the approximate posterior $\pi_{\lambda}$ as the number of data $n$ tends to infinity. Specifically, suppose these are split equally into $b$ blocks; considering $\lambda$ and $b$ as functions of $n$, we find that $\pi_{\lambda}$ exhibits posterior consistency if $\lambda/b$ decreases to $0$ as $n\to\infty$, and that credible intervals have asymptotically correct coverage if $\lambda/b$ decreases at a rate strictly faster than $n^{-1}$.

As an alternative to selecting a single value of $\lambda$, we propose in Section 3 a Sequential Monte Carlo sampler employing Markov kernels formed via Algorithm 1. In this manner a decreasing sequence of $\lambda$ values may be considered, which may result in lower-variance estimates for small $\lambda$ values; we also describe a possible bias correction technique.

As previously mentioned, a Gibbs sampler construction essentially corresponding to Algorithm 1 has independently been proposed by Vono et al., 2019a . Their main objective is to improve algorithmic performance when computation of the full posterior density is intensive by constructing full conditional distributions that are more computationally tractable, for which they primarily consider a setting in which $b=1$. In contrast, we focus on the exploitation of this framework in distributed settings (i.e. with $b>1$), in the manner described in Section 2.2.

The objectives of our algorithm are similar, but not identical, to those of the previously-introduced ‘embarrassingly parallel’ approaches proposed by many authors. These reduce the costs of communication latency to a near-minimum by simulating a separate MCMC chain on each machine; typically, the chain on the $j$th machine is invariant with respect to a ‘subposterior’ distribution with density proportional to $\mu(z)^{1/b}f_{j}(z)$. Communication is necessary only for the final aggregation step, in which an approximation of the full posterior is obtained using the samples from all $b$ chains.

A well-studied approach within this framework is Consensus Monte Carlo (Scott et al.,, 2016), in which the post-processing step amounts to forming a ‘consensus chain’ by weighted averaging of the separate chains. In the case that each subposterior density is Gaussian this approach can be used to produce samples asymptotically distributed according to $\pi$, by weighting each chain using the precision matrices of the subposterior distributions. Motivated by Bayesian asymptotics, the authors suggest using this approach more generally. In cases where the subposterior distributions exhibit near-Gaussianity this performs well, with the final ‘consensus chain’ providing a good approximation of posterior expectations. However, there are no theoretical guarantees associated with this approach in settings in which the subposterior densities are poorly approximated by Gaussians. In such cases consensus Monte Carlo sometimes performs poorly in forming an approximation of the posterior $\pi$ (as in examples of Wang et al.,, 2015; Srivastava et al.,, 2015; Dai et al.,, 2019), and so the resulting estimates of integrals $\pi(\varphi)$ exhibit high bias.

Various authors (e.g. Minsker et al.,, 2014; Rabinovich et al.,, 2015; Srivastava et al.,, 2015; Wang et al.,, 2015) have therefore proposed alternative techniques for utilising the values from each of these chains in order to approximate posterior expectations, each of which presents benefits and drawbacks. For example, Neiswanger et al., (2014) propose a strategy based on kernel density estimation; based on this approach, Scott, (2017) suggests a strategy based on finite mixture models, though notes that both methods may be impractical in high-dimensional settings.

An aggregation procedure proposed by Wang and Dunson, (2013) bears some relation to our proposed framework, being based on the application of Weierstrass transforms to each subposterior density. The resulting smoothed densities are analogous to (3), which represents a smoothed form of the partial likelihood $f_{j}$. As well as proposing an aggregation method based on rejection sampling, the authors propose a technique for ‘refining’ an initial posterior approximation, which may be expressed in terms of a Gibbs kernel on an extended state space. Comparing with our framework, this is analogous to applying Algorithm 1 for one iteration with $N$ different initial values.

A potential issue common to these approaches is the treatment of the prior density $\mu$. Each subposterior density is typically assigned an equal share of the prior information in the form of a fractionated prior density $\mu(z)^{1/b}$, but it is not clear when this approach is satisfactory. For example, suppose $\mu$ belongs to an exponential family; any property that is not invariant to multiplying the canonical parameters by a constant will not be preserved in the fractionated prior. As such if $\mu(z)^{1/b}$ is proportional to a valid probability density function of $z$, then the corresponding distribution may be qualitatively very different to the full prior. Although Scott et al., (2016) note that fractionated priors perform poorly on some examples (for which tailored solutions are provided), no other general way of assigning prior information to each block naturally presents itself. In contrast our approach avoids this problem entirely, with $\mu$ providing prior information for $Z$ at the ‘global’ level.

Finally, we believe this work is complementary to other distributed algorithms lying outside of the ‘embarrassingly parallel’ framework (including Xu et al.,, 2014; Jordan et al.,, 2019), and to approaches that aim to reduce the amount of computation associated with each likelihood calculation on a single node, e.g. by using only a subsample or batch of the data (as in Korattikara et al.,, 2014; Bardenet et al.,, 2014; Huggins et al.,, 2016; Maclaurin and Adams,, 2014).

We present an approach to use many of the particle approximations produced by Algorithm 2. A natural idea is to regress the values of $\pi_{\lambda}^{N}(\varphi)$ on $\lambda$, extrapolating to $\lambda=0$ to obtain an estimate of $\pi(\varphi)$. A similar idea has been used for bias correction in the context of Approximate Bayesian Computation, albeit not in an SMC setting, regressing on the discrepancy between the observed data and simulated pseudo-observations (Beaumont et al.,, 2002; Blum and François,, 2010).

Under very mild assumptions on the transition densities $K_{j}^{(\lambda)}$, $\pi_{\lambda}(\varphi)$ is smooth as a function of $\lambda$. Considering a first-order Taylor expansion of this function, a simple approach is to model the dependence of $\pi_{\lambda}(\varphi)$ on $\lambda$ as linear, for $\lambda$ sufficiently close to $0$. Having determined a subset of the values of $\lambda$ used for which a linear approximation is appropriate (we propose a heuristic approach in Section LABEL:subsec:inclusionProcedure of the supplement) one can use linear least squares to carry out the regression. To account for the SMC estimates $\pi_{\lambda_{p}}^{N}(\varphi)$ having different variances, we propose the use of weighted least squares, with the ‘observations’ $\pi_{\lambda_{p}}^{N}(\varphi)$ assigned weights inversely proportional to their estimated variances; we describe methods for computing such variance estimates in Section 3.2. A bias-corrected estimate of $\pi(\varphi)$ is then obtained by extrapolating the resulting fit to $\lambda=0$, which corresponds to taking the estimated intercept term.

To make this explicit, first consider the case in which $\varphi:\mathsf{E}\to\mathbb{R}$, so that the estimates $\pi_{\lambda}^{N}(\varphi)$ are univariate. For each value $\lambda_{p}$ denote the corresponding SMC estimate by $\eta_{p}\coloneqq\pi_{\lambda_{p}}^{N}(\varphi)$, and let $v_{p}$ denote some proxy for the variance of this estimate. Then for some set of indices $S\coloneqq\{p^{*},\dots,n\}$ chosen such that the relationship between $\eta_{p}$ and $\lambda_{p}$ is approximately linear for $p\in S$, a bias-corrected estimate for $\pi(\varphi)$ may be computed as

where $\tilde{\lambda}_{S}$ and $\tilde{\eta}_{S}$ denote weighted means given by

The formal justification of this estimate assumes that the observations are uncorrelated, which does not hold here. We demonstrate in Section 4 examples on which this simple approach is nevertheless effective, but in principle one could use generalised least squares combined with some approximation of the full covariance matrix of the SMC estimates.

In the more general case where $\varphi:\mathsf{E}\to\mathbb{R}^{d}$ for $d>1$, we propose simply evaluating (13) for each component of this quantity separately, which corresponds to fitting an independent weighted least squares regression to each component. This facilitates the use of the variance estimators described in the following section, though in principle one could use multivariate weighted least squares or other approaches.

We propose the weighted form of least squares here since, as the values of $\lambda$ used in the SMC procedure approach zero, the estimators generated may increase in variance: partly due to poorer mixing of the MCMC kernels as previously described, but also due to the gradual degeneracy of the particle set. In order to estimate the variances of estimates generated using SMC, several recent approaches have been proposed that allow this estimation to be carried out online by considering the genealogy of the particles. Using any such procedure, one may estimate the variance of $\pi_{\lambda}^{N}(\varphi)$ for each $\lambda$ value considered by Algorithm 2, with these values used for each $v_{p}$ in (13).

Within our examples, we use the estimator proposed by Lee and Whiteley, (2018), which for fixed $N$ coincides with an earlier proposal of Chan and Lai, (2013) (up to a multiplicative constant). Specifically, after each step of the SMC sampler we compute an estimate of the asymptotic variance of each estimate $\pi_{\lambda_{p}}^{N}(\varphi)$; that is, the limit of $N\operatorname{var}[\pi_{\lambda_{p}}^{N}(\varphi)]$ as $N\to\infty$. While this is not equivalent to computing the true variance of each estimate, for fixed large $N$ the relative sizes of these asymptotic variance estimates should provide a useful indicator of the relative variances of each estimate $\pi_{\lambda}^{N}(\varphi)$. In Section LABEL:subsec:gaussianSMCexample of the supplement we show empirically that inversely weighting the SMC estimates according to these estimated variances can result in more stable bias-corrected estimates as the particle set degenerates. We also explain in Section LABEL:subsec:stoppingRule how these estimated variances can be used within a rule to determine when to terminate the algorithm.

The asymptotic variance estimator described is consistent in $N$. However, if in practice resampling at the $p$th time step causes the particle set to degenerate to having a single common ancestor, then the estimator evaluates to zero, and so it is impossible to use this value as the inverse weight $v_{p}$ in (13). Such an outcome may be interpreted as a warning that too few particles have been used for the resulting SMC estimates to be reliable, and that a greater number should be used when re-running the procedure.

To compare the posterior approximations formed by our global consensus framework with those formed by some embarrassingly parallel approaches discussed in Section 2.4, we conduct a simulation study based on a simple model. Let $\mathcal{LN}(x;\mu,\sigma^{2})$ denote the density of a log-normal distribution with parameters $(\mu,\sigma^{2})$; that is,

One may consider a model with prior density $\mu(z)=\mathcal{LN}(z;\mu_{0},\sigma_{0}^{2})$ and likelihood contributions $f_{j}(z)=\mathcal{LN}(\log(\mu_{j});\log(z),\sigma_{j}^{2})$ for $j\in\{1,\dots,b\}$. This may be seen as a reparametrisation of the Gaussian model in Section LABEL:sec:Gaussians of the supplement, in which each likelihood contribution is that of a data subset with a Gaussian likelihood. This convenient setting allows for the target distribution $\pi$ to be expressed analytically. For the implementation of the global consensus algorithm, we choose Markov transition kernels given by $K_{j}^{(\lambda)}(z,x)=\mathcal{LN}(x;\log(z),\lambda)$ for each $j\in\{1,\dots,b\}$, which satisfy Assumption 1; this allows for exact sampling from all the full conditional distributions.

As a toy example to illustrate the effects of non-Gaussian partial likelihoods we consider a case in which $f_{j}(z)=\mathcal{LN}(\log(\mu_{j});\log(z),1)$ for each $j$, and $\mu(z)=\mathcal{LN}(z;0,25)$. Here we took $b=32$, and selected the location parameters $\mu_{j}$ as i.i.d. samples from a standard normal distribution. We ran Global Consensus Monte Carlo (GCMC) using the Gibbs sampler form of Algorithm 1, for values of $\lambda$ between $10^{-5}$ and $10$. For comparison we also drew samples from each subposterior distribution as defined in Section 2.4, combining the samples using various approaches. These are the Consensus Monte Carlo (CMC) averaging of Scott et al., (2016); the nonparametric density product estimation (NDPE) approach of Neiswanger et al., (2014); and the Weierstrass rejection sampling (WRS) combination technique of Wang and Dunson, (2013), using their R implementation (https://github.com/wwrechard/weierstrass). In each case we ran the algorithm $25$ times, drawing $N=10^{5}$ samples.

To demonstrate the role of $\lambda$ in the bias–variance decomposition (11), Table 1 presents the means and standard deviations of estimates of $\pi(\varphi)$, for various test functions $\varphi$. In estimating the first moment of $\pi$, GCMC generates a low-bias estimator when $\lambda$ is chosen to be sufficiently small; however, as expected, the variance of such estimators increases when very small values of $\lambda$ are chosen. While the other methods produce estimators of reasonably low variance, these exhibit somewhat higher bias. For CMC the bias is especially pronounced when estimating higher moments of the posterior distribution, as exemplified by the estimates of the fifth moment. Note however that high biases are also introduced when using Global Consensus Monte Carlo with large values of $\lambda$ (as seen here with $\lambda=10$), for which $\pi_{\lambda}$ is a poor approximation of $\pi$.

Also of note are estimates of $\int\log(z)\pi(z)\mathop{}\!\text{d}z$, corresponding to the mean of the Gaussian model of which this a reparametrisation. While Global Consensus Monte Carlo performs well across a range of $\lambda$ values, the other methods perform less favourably; Consensus Monte Carlo produces an estimate that is incorrect by an order of magnitude. While this could be solved by a simple reparametrisation of the problem in this case, in more general settings no such straightforward solution may exist.

In Section LABEL:subsec:additionalLognormalExample of the supplement we present second example based on a log-normal model, demonstrating the robustness of these methods to permutation and repartitioning of the data.

Binary logistic regression models are commonly used in settings related to marketing. In web design for example, A/B testing may be used to determine which content choices lead to maximised user interaction, such as the user clicking on a product for sale.

We assume that we have a data set of size $n$ formed of responses $\eta_{\ell}\in\{-1,1\}$, and vectors $\xi_{\ell}\in\{0,1\}^{d}$ of binary covariates, where $\ell\in\{1,\dots,n\}$. The likelihood contribution of each block of data then takes the form $f_{j}(z)=\prod_{\ell}S(\eta_{\ell}z^{\mathsf{T}}\xi_{\ell})$, $z\in\mathbb{R}^{d}$, where the product is taken over those indices $\ell$ included in the $j$th block of data, and $S$ denotes the logistic function, $S(x)=(1+\exp(x))^{-1}$.

For the prior $\mu$, we use a product of independent zero-mean Gaussians, with standard deviation $20$ for the parameter corresponding to the intercept term, and $5$ for all other parameters. For the Markov transition densities in GCMC, we use multivariate Gaussian densities: $K_{j}^{(\lambda)}(z,x)=\mathcal{N}(x;z,\lambda I)$ for each $j\in\{1,\dots\ b\}$.

We investigated several such simulated data sets and the efficacy of various approaches in approximating the true posterior $\pi$. To illustrate the bias–variance trade-off described in Section 2.3, in the presentation of these results we focus on the estimation of the posterior first moment; denoting the identity function on $\mathbb{R}^{d}$ by $\operatorname{Id}$, we may write this as $\pi(\operatorname{Id})$. While our global consensus approach was consistently successful in forming estimators with low mean squared error in each component, in low-dimensional settings the application of Consensus Monte Carlo often resulted in marginal improvements. However, in many higher-dimensional settings, the estimators resulting from CMC exhibited relatively large biases.

We present here an example in which the $d$ predictors correspond to $p$ binary input variables, their pairwise products, and an intercept term, so that $d=1+p+\binom{p}{2}$. In settings where the interaction effects corresponding to these pairwise products are of interest, the dimensionality $d$ of the space can be very large compared to $p$. We used a simulated data set with $p=20$ input variables, resulting in a parameter space of dimension $d=211$. The data comprise $n=80{,}000$ observations, split into $b=8$ equally-sized blocks. Each observation of the $20$ binary variables was generated from a Bernoulli distribution with parameter $0.1$, and for each vector of covariates, the response was generated from the correct model, for a fixed underlying parameter vector $z^{*}$.