A tutorial on bridge sampling.

The marginal likelihood plays an important role in many areas of Bayesian statistics such as parameter estimation, model comparison, and model averaging. In most applications, however, the marginal likelihood is not analytically tractable and must be approximated using numerical methods. Here we provide a tutorial on bridge sampling (Bennett, 1976; Meng & Wong, 1996), a reliable and relatively straightforward sampling method that allows researchers to obtain the marginal likelihood for models of varying complexity. First, we introduce bridge sampling and three related sampling methods using the beta-binomial model as a running example. We then apply bridge sampling to estimate the marginal likelihood for the Expectancy Valence (EV) model-a popular model for reinforcement learning. Our results indicate that bridge sampling provides accurate estimates for both a single participant and a hierarchical version of the EV model. We conclude that bridge sampling is an attractive method for mathematical psychologists who typically aim to approximate the marginal likelihood for a limited set of possibly high-dimensional models.

Bayesian statistics has become increasingly popular in mathematical psychology Andrews and Baguley (2013), Bayarri et al. (2016), Poirier (2006), Vanpaemel (2016), Verhagen et al. (2015), Wetzels et al. (2016). The Bayesian approach is conceptually simple, theoretically coherent, and easily applied to relatively complex problems. These problems include, forinstance, hierarchical modeling Matzke et al. (2015), Matzke and Wagenmakers (2009), Rouder and Lu (2005), Rouder et al. (2005), Rouder et al. (2007) or the comparison of non-nested models Lee (2008), Pitt et al. (2002), Shiffrin et al. (2008). Three major applications of Bayesian statistics concern parameter estimation, model comparison, and Bayesian model averaging. In all three areas, the marginal likelihood – that is, the probability of the observed data given the model of interest – plays a central role (see also Gelman & Meng, 1998).

First, in parameter estimation, we consider a single model and aim to quantify the uncertainty for a parameter of interest after having observed the data . This is realized by means of a posterior distribution that can be obtained using Bayes theorem:

Here, the marginal likelihood of the data ensures that the posterior distribution is a proper probability density function (PDF) in the sense that it integrates to 1. This illustrates why in parameter estimation the marginal likelihood is referred to as a normalizing constant.

Second, in model comparison, we consider () competing models, and are interested in the relative plausibility of a particular model () given the prior model probability and the evidence from the data (see three special issues on this topic in the Journal of Mathematical Psychology : Mulder & Wagenmakers, 2016; Myung, Forster, & Browne, 2000; Wagenmakers & Waldorp, 2006). This relative plausibility is quantified by the so-called posterior model probability of model given the data  (Berger & Molina, 2005) : where the denominator is the sum of the marginal likelihood times the prior model probability of all models. In model comparison, the marginal likelihood for a specific model is also referred to as the model evidence (Didelot, Everitt, Johansen, & Lawson, 2011), the integrated likelihood (Kass & Raftery, 1995), the predictive likelihood of the model (Gamerman & Lopes, 2006, chapter 7), the predictive probability of the data (Kass & Raftery, 1995), or the prior predictive density (Ntzoufras, 2009). Note that conceptually the marginal likelihood of Eq. (2) is the same as the marginal likelihood of Eq. (1). However, for the latter equation we droped the model index because in parameter estimation we consider only one model.

If only two models and are considered, Eq. (2) can be used to quantify the relative posterior model plausibility of model compared to model . This relative plausibility is given by the ratio of the posterior probabilities of both models, and is referred to as the posterior model odds:

Eq. (3) illustrates that the posterior model odds are the product of two factors: The first factor is the ratio of the prior probabilities of both models—the prior model odds. The second factor is the ratio of the marginal likelihoods of both models—the so-called Bayes factor Etz and Wagenmakers (in press), Jeffreys (1961), Ly et al. (2016a), Ly et al. (2016b), Robert (2016). The Bayes factor plays an important role in model comparison and is referred to as the “standard Bayesian solution to the hypothesis testing and model selection problems” (Lewis & Raftery, 1997, p. 648) and “the primary tool used in Bayesian inference for hypothesis testing and model selection” (Berger, 2006, p. 378).

Third, the marginal likelihood plays an important role in Bayesian model averaging (BMA; Hoeting, Madigan, Raftery, & Volinsky, 1999) where aspects of parameter estimation and model comparison are combined. As in model comparison, BMA considers several models; however, it does not aim to identify a single best model. Instead it fully acknowledges model uncertainty. Model averaged parameter inference can be obtained by combining, across all models, the posterior distribution of the parameter of interest weighted by each model’s posterior model probability, and as such depends on the marginal likelihood of the models. This procedure assumes that the parameter of interest has identical interpretation across the different models. Model averaged predictions can be obtained in a similar manner.

A problem that arises in all three areas – parameter estimation, model comparison, and BMA – is that an analytical expression of the marginal likelihood can be obtained only for certain restricted examples. This is a pressing problem in Bayesian modeling, and in particular in mathematical psychology where models can be non-linear and equipped with a large number of parameters, especially when the models are implemented in a hierarchical framework. Such a framework incorporates both commonalities and differences between participants of one group by assuming that the model parameters of each participant are drawn from a group-level distribution (for advantages of the Bayesian hierarchical framework see Ahn, Krawitz, Kim, Busemeyer, & Brown, 2011; Navarro, Griffiths, Steyvers, & Lee, 2006; Rouder & Lu, 2005; Rouder et al., 2005; Rouder, Lu, Morey, Sun, & Speckman, 2008; Scheibehenne & Pachur, 2015; Shiffrin et al., 2008; Wetzels, Vandekerckhove, Tuerlinckx, & Wagenmakers, 2010). For instance, consider a four-parameter Bayesian hierarchical model with four group-level distributions each characterized by two parameters and a group size of 30 participants; this then results in 30  4 individual-level parameters and 2  4 group-level parameters for a total of 128 parameters. In sum, even simple models quickly become complex once hierarchical aspects are introduced and this frustrates the derivation of the marginal likelihood.

To overcome this problem, several Monte Carlo sampling methods have been proposed to approximate the marginal likelihood. In this tutorial we focus on four such methods: the bridge sampling estimator (Bennett, 1976, Chapter 5 of Chen, Shao, & Ibrahim, 2012; Meng & Wong, 1996), and its three commonly used special cases—the naive Monte Carlo estimator, the importance sampling estimator, and the generalized harmonic mean estimator (for alternative methods see Gamerman & Lopes, 2006, Chapter 7; and for alternative approximation methods relevant to model comparison and BMA see Carlin & Chib, 1995; Green, 1995).1 As we will illustrate throughout this tutorial, the bridge sampler is accurate, efficient, and relatively straightforward to implement (e.g., DiCiccio, Kass, Raftery, & Wasserman, 1997; Frühwirth-Schnatter, 2004; Meng & Wong, 1996).

The goal of this tutorial is to bring the bridge sampling estimator to the attention of mathematical psychologists. We aim to explain this estimator and facilitate its application by suggesting a step-by-step implementation scheme. To this end, we first show how bridge sampling and the three special cases can be used to approximate the marginal likelihood in a simple beta-binomial model. We begin with the naive Monte Carlo estimator and progressively work our way up – via the importance sampling estimator and the generalized harmonic mean estimator – to the most general case considered: the bridge sampling estimator. This order was chosen such that key concepts are introduced gradually and estimators are of increasing complexity and sophistication. The first three estimators are included in this tutorial with the sole purpose of facilitating the reader’s understanding of bridge sampling. In the second part of this tutorial, we outline how the bridge sampling estimator can be used to derive the marginal likelihood for the Expectancy Valence (EV; Busemeyer & Stout, 2002) model—a popular, yet relatively complex reinforcement-learning model for the Iowa gambling task (Bechara, Damasio, Damasio, & Anderson, 1994). We apply bridge sampling to both an individual-level and a hierarchical implementation of the EV model.

Throughout the article, we use the software package R to implement the bridge sampling estimator for the various models. The interested reader is invited to reproduce our results by downloading the code and all relevant materials from our Open Science Framework folder at osf.io/f9cq4.

Four sampling methods to approximate the marginal likelihood

In this section we outline four standard methods to approximate the marginal likelihood. For more detailed explanations and derivations, we recommend Ntzoufras (2009, Chapter 11) and Gamerman and Lopes (2006, Chapter 7); a comparative review of the different sampling methods is presented in DiCiccio et al. (1997). The marginal likelihood is the probability of the observed data given a specific model of interest , and is defined as the integral of the likelihood over the prior: with a vector containing the model parameters. Eq. (4) illustrates that the marginal likelihood can be interpreted as a weighted average of the likelihood of the data given a specific value for where the weight is the a priori plausibility of that specific value. Eq. (4) can therefore be written as an expected value: where the expectation is taken with respect to the prior distribution. This idea is central to the four sampling methods that we discuss in this tutorial.

Introduction of the running example: the beta-binomial model

Our running example focuses on estimating the marginal likelihood for a binomial model assuming a uniform prior on the rate parameter (i.e., the beta-binomial model). Consider a single participant who answered out of true/false questions correctly. Assume that the number of correct answers follows a binomial distribution, that is, with , where represents the latent probability for answering any one question correctly. The probability mass function (PMF) of the binomial distribution is given by: where , and . The PMF of the binomial distribution serves as the likelihood function in our running example.

In the Bayesian framework, we also have to specify the prior distribution of the model parameters; the prior distribution expresses our knowledge about the parameters before the data have been observed. In our running example, we assume that all values of are equally likely a priori. This prior belief is captured by a uniform distribution across the range of , that is, which can equivalently be written in terms of a beta distribution . This prior distribution is represented by the dotted line in Fig. 1. It is evident that the density of the prior distribution equals 1 for all values of . One advantage of expressing the prior distribution by a beta distribution is that its two parameters (i.e., in its general form the shape parameters and ) can be thought of as counts of “prior successes” and “prior failures”, respectively. In its general form, the PDF of a distribution () is given by: where is the beta function that is defined as: , and for .

Fig. 1

Prior and posterior distribution for the rate parameter from the beta-binomial model. The prior on the rate parameter is represented by the dotted line; the posterior distribution is represented by the solid line and was obtained after having observed 2 correct responses out of 10 trials. Available at https://tinyurl.com/yc8bw98v under CC license https://creativecommons.org/licenses/by/2.0/.

Analytical derivation of the marginal likelihood

As we will see in this section, the beta-binomial model constitutes one of the rare examples where the marginal likelihood is analytic. Assuming a general and , we obtain the marginal likelihood as: where we suppress the “model” in the conditioning part of the probability statements because we focus on a single model in this running example. Using and of our example, we obtain: . This value will be estimated in the remainder of the running example using the naive Monte Carlo estimator, the importance sampling estimator, the generalized harmonic mean estimator, and finally the bridge sampling estimator.

As we will see below, the importance sampling, generalized harmonic mean estimator, and bridge sampling estimator require samples from the posterior distribution. These samples can be obtained using computer software such as WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000), JAGS (Plummer, 2003), or Stan (Stan Development Team, 2016), even when the marginal likelihood that functions here as a normalizing constant is not known (Eq. (1)). However, in our running example MCMC samples are not required because we can derive an analytical expression of the posterior distribution for after having observed the data. Using the analytic expression of the marginal likelihood (Eq. (6)) and Bayes theorem, we obtain: which we recognize as the PDF of the distribution. Thus, if we assume a uniform prior on and observe correct responses out of trials, we obtain a distribution as posterior distribution. This distribution is represented by the solid line in Fig. 1. In general, if and , then .

Method 1: the naive Monte Carlo estimator of the marginal likelihood

The simplest method to approximate the marginal likelihood is provided by the naive Monte Carlo estimator Hammersley and Handscomb (1964), Raftery and Banfield (1991). This method uses the standard definition of the marginal likelihood (Eq. (4)), and relies on the central idea that the marginal likelihood can be written as an expected value with respect to the prior distribution, that is, . This expected value of the likelihood of the data with respect to the prior can be approximated by evaluating the likelihood in samples from the prior distribution for and averaging the resulting values. This yields the naive Monte Carlo estimator :

Running example

To obtain the naive Monte Carlo estimate of the marginal likelihood in our running example, we need samples from the prior distribution for . For illustrative purposes, we limit the number of samples to 12 whereas in practice one should take to be very large. We obtain the following samples: where we use the tilde symbol to emphasize that we refer to a sampled value. All sampled values are represented by the gray dots in Fig. 2.

Fig. 2

Illustration of the naive Monte Carlo estimator for the beta-binomial example. The dotted line represents the prior distribution and the solid line represents the posterior distribution that was obtained after having observed 2 correct responses out of 10 trials. The gray dots represent the 12 samples randomly drawn from the prior distribution. Available at https://tinyurl.com/y8uf6t8f under CC license https://creativecommons.org/licenses/by/2.0/.

Following Eq. (7), the next step is to calculate the likelihood (Eq. (5)) for each , and then to average all obtained likelihood values. This yields the naive Monte Carlo estimate of the marginal likelihood:

Method 2: the importance sampling estimator of the marginal likelihood

The naive Monte Carlo estimator introduced in the last section performs well if the prior and posterior distribution have a similar shape and strong overlap. However, the estimator is unstable if the posterior distribution is peaked relative to the prior (e.g., Gamerman & Lopes, 2006; Ntzoufras, 2009). In such a situation, most of the sampled values for result in likelihood values close to zero and contribute only minimally to the estimate. This means that those few samples that result in high likelihood values dominate estimates of the marginal likelihood. Consequently, the variance of the estimator is increased Newton and Raftery (1994), Pajor (2016).2

The importance sampling estimator, on the other hand, overcomes this shortcoming by boosting sampled values in regions of the parameter space where the integrand of Eq. (4) is large. This is realized by using samples from a so-called importance density instead of the prior distribution. The advantage of sampling from an importance density is that values for that result in high likelihood values are sampled most frequently, whereas values for with low likelihood values are sampled only rarely.

To derive the importance sampling estimator, Eq. (4) is used as starting point which is then extended by the importance density : This yields the importance sampling estimator :

A suitable importance density should (1) be easy to evaluate; (2) have the same domain as the posterior distribution; (3) closely resemble the posterior distribution; and (4) have fatter tails than the posterior distribution Neal (2001), Vandekerckhove et al. (2015). The latter criterion ensures that values in the tails of the distribution cannot misleadingly dominate the estimate (Neal, 2001).3

To obtain the importance sampling estimate of the marginal likelihood in our running example, we first need to choose an importance density . An importance density that fulfills the four above mentioned desiderata is a mixture between a beta density that provides the best fit to the posterior distribution and a uniform density across the range of  (Vandekerckhove et al., 2015). The relative impact of the uniform density is quantified by a mixture weight that ranges between and . The larger , the higher the influence of the uniform density resulting in a less peaked distribution with thick tails. If , the beta mixture density simplifies to the uniform distribution on ;4 and if , the beta mixture density simplifies to the beta density that provides the best fit to the posterior distribution.

In our specific example, we already know that the density is the beta density that provides the best fit to the posterior distribution because this is the analytic expression of the posterior distribution. However, to demonstrate the general case, we show how we can find the beta distribution with the best fit to the posterior distribution using the method of moments. This particular method works as follows. First, we draw samples from our posterior distribution and obtain:5 Note that here we use to refer to the th sample from the posterior distribution to distinguish it from the previously used —the th sample from a distribution other than the posterior distribution, such as a prior distribution or an importance density. Second, we compute the mean and variance of these posterior samples. We obtain a mean of and a variance of .

Third, knowing that, if , then and , we obtain the following method of moment estimates for and : Using a mixture weight on the uniform component of  – a choice that was made to ensure that, visually, the tails of the importance density are clearly thicker than the tails of the posterior distribution – we obtain the following importance density: . This importance density is represented by the dashed line in Fig. 3. The figure also shows the posterior distribution (solid line). As is evident from the figure, the beta mixture importance density resembles the posterior distribution, but has fatter tails.

Fig. 3

Illustration of the importance sampling estimator for the beta-binomial model. The dashed line represents our beta mixture importance density and the solid gray line represents the posterior distribution that was obtained after having observed 2 correct responses out of 10 trials. The gray dots represent the 12 samples randomly drawn from our beta mixture importance density. Available at https://tinyurl.com/yc7ho7hr under CC license https://creativecommons.org/licenses/by/2.0/.

In general, it is advised to choose the mixture weight on the uniform component small enough to make the estimator efficient, yet large enough to produce fat tails to stabilize the estimator. A suitable mixture weight can be realized by gradually minimizing the mixture weight and investigating whether stability is still guaranteed (i.e., robustness analysis).

Drawing samples for from our beta mixture importance density results in: These samples are represented by the gray dots in Fig. 3.

The final step is to compute the average adjusted likelihood for the 12 samples using Eq. (8). This yields the importance sampling estimate of the marginal likelihood as:

Method 3: the generalized harmonic mean estimator of the marginal likelihood

Just as the importance sampling estimator, the generalized harmonic mean estimator focuses on regions of the parameter space where the integrand of Eq. (4) is large by using an importance density  (Gelfand & Dey, 1994).6 However, in contrast to the importance sampling estimator, the generalized harmonic mean estimator requires an importance density with thinner tails for an analogous reason as in importance sampling.

To derive the generalized harmonic mean estimator, also known as reciprocal importance sampling estimator (Frühwirth-Schnatter, 2004), we use the following identity: Rewriting results in: which is used to define the generalized harmonic mean estimator  (Gelfand & Dey, 1994) as follows:

Note that the generalized harmonic mean estimator – in contrast to the importance sampling estimator – evaluates samples from the posterior distribution. In addition, note that the ratio in Eq. (9) is the reciprocal of the ratio in Eq. (8); this explains why the importance density for the generalized harmonic mean estimator should have thinner tails than the posterior distribution in order to avoid inflation of the ratios that are part of the summation displayed in Eq. (9). Thus, in the case of the generalized harmonic mean estimator, a suitable importance density should (1) have thinner tails than the posterior distribution DiCiccio et al. (1997), Newton and Raftery (1994), and as in importance sampling, it should (2) be easy to evaluate; (3) have the same domain as the posterior distribution; and (4) closely resemble the posterior distribution.

To obtain the generalized harmonic mean estimate of the marginal likelihood in our running example, we need to choose a suitable importance density. In our running example, an importance density that fulfills the four above mentioned desiderata can be obtained by following four steps: First, we draw samples from the posterior distribution. Reusing the samples from the last section, we obtain:

Second, we probit-transform all posterior samples (i.e., ).7 The result of this transformation is that the samples range across the entire real line instead of the interval only. We obtain:

These probit-transformed samples are represented by the gray dots in Fig. 4.

Fig. 4

Illustration of the generalized harmonic mean estimator for the beta-binomial model. The solid line represents the probit-transformed posterior distribution that was obtained after having observed 2 correct responses out of 10 trials, and the dashed line represents the importance density . The gray dots represent the 12 probit-transformed samples randomly drawn from the posterior distribution. Available at https://tinyurl.com/yazgk8kj under CC license https://creativecommons.org/licenses/by/2.0/.

Third, we search for the normal distribution that provides the best fit to the probit-transformed posterior samples . Using the method of moments, we obtain as estimates and . Note that the choice of a normal importance density justifies step 2; the probit transformation (or an equivalent transformation) was required to match the range of the posterior distribution to the one of the normal distribution.

Finally, as importance density we choose a normal distribution with mean and standard deviation . This additional division by 1.5 is to ensure thinner tails of the importance density than of the probit-transformed posterior distribution (for a discussion of alternative importance densities see DiCiccio et al., 1997). We decided to divide by 1.5 for illustrative purposes only. Our importance density is displayed in Fig. 4 (dashed line) together with the probit-transformed posterior distribution (solid line).

The generalized harmonic mean estimate can now be obtained using either the original posterior samples or the probit-transformed samples . Here we use the latter ones (see also Overstall & Forster, 2010). Incorporating our specific importance density and a correction for having used the probit-transformation, Eq. (9) becomes:8

For our beta-binomial model, we now obtain the generalized harmonic mean estimate of the marginal likelihood as:

Method 4: the bridge sampling estimator of the marginal likelihood

As became evident in the last two sections, both the importance sampling estimator and the generalized harmonic mean estimator impose strong constraints on the tail behavior of the importance density relative to the posterior distribution to guarantee a stable estimator. Such requirements can make it difficult to find a suitable importance density, especially when a high dimensional posterior is considered. The bridge sampler, on the other hand, alleviates such requirements (e.g., Frühwirth-Schnatter, 2004).

Originally, bridge sampling was developed to directly estimate the Bayes factor, that is, the ratio of the marginal likelihoods of two models and (e.g., Jeffreys, 1961; Kass & Raftery, 1995). However, in this tutorial, we use a version of bridge sampling that allows us to approximate the marginal likelihood of a single model (for an earlier application see for example Overstall & Forster, 2010). This version is based on the following identity: where is the so-called proposal distribution and the so-called bridge function. Multiplying both sides of Eq. (11) by the marginal likelihood results in:

The marginal likelihood can now be approximated using:

Eq. (12) illustrates that we need samples from both the proposal distribution and the posterior distribution to obtain the bridge sampling estimate for the marginal likelihood. However, before we can apply Eq. (12) to our running example, we have to discuss how we can obtain a suitable proposal distribution and bridge function. Conceptually, the proposal distribution is similar to an importance density, should resemble the posterior distribution, and should have sufficient overlap with the posterior distribution. According to Overstall and Forster (2010), a convenient proposal distribution is often a normal distribution with its first two moments chosen to match those of the posterior distribution. In our experience, this choice for the proposal distribution works well for a wide range of scenarios. However, this proposal distribution might produce unstable estimates in case of high-dimensional posterior distributions that clearly do not follow a multivariate normal distribution. In such a situation, it might be advisable to consider more sophisticated versions of bridge sampling (e.g., Frühwirth-Schnatter, 2004; Meng & Schilling, 2002; Wang & Meng, 2016).

Choosing the optimal bridge function

In this tutorial we use the bridge function defined as (Meng & Wong, 1996): where , , and a constant; its particular value is not required because is part of both the numerator and the denominator of Eq. (12), and therefore the constant cancels. This particular bridge function is referred to as the “optimal bridge function” because Meng and Wong (1996, p. 837) proved that it minimizes the relative mean-squared error (Eq. (16)).

Eq. (13) shows that the optimal bridge function depends on the marginal likelihood which is the very entity we want to approximate. We can resolve this issue by applying an iterative scheme that updates an initial guess of the marginal likelihood until the estimate of the marginal likelihood has converged according to a predefined tolerance level. To do so, we insert the expression for the optimal bridge function (Eq. (13)) in Eq. (12) (Meng & Wong, 1996). The formula to approximate the marginal likelihood on iteration is then specified as follows:

where denotes the estimate of the marginal likelihood on iteration of the iterative scheme. Note that Eq. (14) illustrates why bridge sampling is robust to the tail behavior of the proposal distribution relative to the posterior distribution; the difference to the importance sampling and generalized harmonic mean estimator is that, in the case of the bridge sampling estimator, samples from the tail region cannot inflate individual summation terms and thus dominate the estimate. To illustrate this, we consider what happens to the bridge sampling estimator, the importance sampling estimator, and the generalized harmonic mean estimator in case (1) the proposal/importance distribution has fatter tails than the posterior distribution, and (2) the proposal/importance distribution has thinner tails than the posterior distribution (see also Frühwirth-Schnatter, 2004). Specifically, we look at a single term in the respective sums and consider the limit of that term as we move further and further out in the tails. This is insightful since a single term can have a lasting effect on the estimator (e.g., in case a single term in a sum is very large or even infinite).

In case (1) (i.e., the proposal/importance distribution has fatter tails than the posterior), the ratio in the importance sampling estimator (i.e., Eq. (8)) goes to zero as we move further out in the tails. Since samples in the tails may only be obtained occasionally and a zero term in the sum does not inflate the estimate this is not a reason for concern. In contrast, when we consider the ratio in the generalized harmonic mean estimator (i.e., Eq. (9)), we see that the ratio goes to infinity as we move further out in the tails. Even if this occurs only very rarely, this is an issue since the resulting value will dominate the estimate. Consequently, the resulting estimator may have a large variance since samples from the tail regions may be obtained only occasionally across repeated applications. For the bridge sampling estimator (i.e., Eq. (14)), we need to consider the ratio in the numerator and denominator. The ratio in the numerator will go to zero and the ratio in the denominator will go to . Hence, both of these ratios are bounded and will not inflate the two sums, hence also not the resulting estimate.

In case (2) (i.e., the proposal/importance distribution has thinner tails than the posterior), the ratio in the importance sampling estimator (i.e., Eq. (8)) goes to infinity as we move further out in the tails, inflating the estimate. In contrast, when we consider the ratio in the generalized harmonic mean estimator (i.e., Eq. (9)), we see that the ratio goes to zero. As explained above, this is not a reason for concern. These considerations explain why in importance sampling, the importance distribution should have fatter tails than the posterior whereas for the generalized harmonic mean estimator, it should have thinner tails. For the bridge sampling estimator (i.e., Eq. (14)), the ratio in the numerator will go to and the ratio in the denominator will go to zero. Again, both of these ratios are bounded making the bridge sampling estimator more robust to the tail behavior than the other two estimators. This of course assumes that not all terms in the denominator (for case (2)) and the numerator (for case (1)) will be zero, that is, the proposal and the posterior distribution have sufficient overlap. In the extreme scenario of no overlap the bridge sampling estimate is not defined because both sums of Eq. (14) would be zero.

Extending the numerator of the right side of Eq. (14) with , and the denominator with , and subsequently defining and , we obtain the formula for the iterative scheme of the bridge sampling estimator at iteration  (Meng & Wong, 1996, p. 837).

Eq. (15) suggests that, in order to obtain the bridge sampling estimate of the marginal likelihood, a number of requirements need to be fulfilled. First, we need samples from the proposal distribution and samples from the posterior distribution . Second, for all samples from the proposal distribution, we have to evaluate . This involves obtaining the value of the unnormalized posterior (i.e., the product of the likelihood times the prior) and of the proposal distribution for all samples. Third, we evaluate for all samples from the posterior distribution. This is analogous to evaluating . Fourth, we have to determine the constants and that only depend on and . Fifth, we need an initial guess of the marginal likelihood . Since some of these five requirements can be obtained easier than others, we will point out possible challenges.

A first challenge is that using a suitable proposal distribution may involve transforming the posterior samples. Consequently, we have to determine how the transformation affects the definition of the bridge sampling estimator for the marginal likelihood(Eq. (15)).

A second challenge is how to use the samples from the posterior distribution. One option is to use all samples for both fitting the proposal distribution and for computing the bridge sampling estimate. However, Overstall and Forster (2010) showed that such a procedure may result in an underestimation of the marginal likelihood. To obtain more reliable estimates they propose to divide the posterior samples in two parts; the first part is used to obtain the best-fitting proposal distribution, and the second part is used to compute the bridge sampling estimate. Throughout this tutorial, we use two equally large parts. In the remainder we therefore state that we draw samples from the posterior distribution. The first of the total of samples are used for fitting the proposal distribution and the remaining samples are used in the iterative scheme (i.e., Eq. (15)).9

To summarize, the discussion of the requirements and challenges encountered in bridge sampling illustrated that the bridge sampling estimator imposes less strict requirements on the proposal distribution than the importance sampling and generalized harmonic mean estimator and allows for an almost automatic application due to the default choice of the bridge function.10

To obtain the bridge sampling estimate of the marginal likelihood in the beta-binomial example, we follow the eight steps illustrated in Fig. 5 :

Fig. 5

Schematic illustration of the steps involved in obtaining the bridge sampling estimate of the marginal likelihood. Available at https://tinyurl.com/y7b2kze7 under CC license https://creativecommons.org/licenses/by/2.0/.

We draw samples from the posterior distribution for .

We obtain the following sample of 24 values:

Note that the first samples equal the ones used in the last section, whereas the last 12 samples were obtained from drawing again 12 values from the posterior distribution for .

We choose a proposal distribution.

Here we opt for an approach that can be easily generalized to models with multiple parameters and select a normal distribution as the proposal distribution .11

We transform the first batch of posterior samples.

Since we use a normal proposal distribution, we have to transform the posterior samples from the rate scale to the real line so that the range of the posterior distribution matches the range of the proposal distribution. This can be achieved by probit-transforming the posterior samples, that is, with . We obtain:

We fit the proposal distribution to the first batch of probit-transformed posterior samples.

We use the method of moment estimates and from the first batch of probit-transformed posterior samples to obtain our proposal distribution .

We draw samples from the proposal distribution.

We obtain:

We calculate for all samples from the proposal distribution.

This step involves assessing the value of the unnormalized posterior and the proposal distribution for all samples from the proposal distribution. As in the running example for the generalized harmonic mean estimator, we obtain the unnormalized posterior as: , where comes from using the change-of-variable method (see running example for the generalized harmonic mean estimator and the Appendix A, Appendix B, Appendix C, Appendix D, Appendix E for details). Thus, as in the case of the generalized harmonic mean estimator, the uniform prior on translates to a standard normal prior on . The values of the proposal distribution can easily be obtained (for example using the R software).

We transform the second batch of posterior samples.

As in step 2, we use the probit transformation and obtain:

We calculate for the second batch of probit-transformed samples from the posterior distribution.

This is analogous to step 6.

We run the iterative scheme (Eq. (15)) until our predefined tolerance criterion is reached.

As tolerance criterion we choose . This requires an initial guess for the marginal likelihood which we set to 0.12

The simplicity of the beta-binomial model allows us to calculate the bridge sampling estimate by hand. To determine according to Eq. (15), we need to calculate the constants and . Since , we obtain: . In addition, we need to calculate () for all samples from the proposal distribution, and () for the second batch of the probit-transformed samples from the posterior distribution. Here we show how to calculate and using the first sample from the proposal distribution and the first sample of the second batch of the posterior samples, respectively:

For , we then get:

Using , we obtain as updated estimate of the marginal likelihood . This iterative procedure has to be repeated until our predefined tolerance criterion is reached. For our running example, this criterion is reached after five iterations. We now obtain the bridge sampling estimate of the marginal likelihood as .

Interim summary

So far we used the beta-binomial model to illustrate the computation of four different estimators of the marginal likelihood. These four estimators were discussed in order of increasing sophistication, such that the first three estimators provided the proper context for understanding the fourth, most general estimator—the bridge sampler. This estimator is the focus in the remainder of this tutorial. The goal of the next sections is to demonstrate that bridge sampling is particularly suitable to estimate the marginal likelihood of popular models in mathematical psychology. Importantly, bridge sampling may be used to obtain accurate estimates of the marginal likelihood of hierarchical models (for a detailed comparison of bridge sampling versus its special cases see Frühwirth-Schnatter, 2004; Sinharay & Stern, 2005).

Assessing the accuracy of the bridge sampling estimate

In this section we show how to quantify the accuracy of the bridge sampling estimate. A straightforward approach would be to apply the bridge sampling procedure multiple times and investigate the variability of the marginal likelihood estimate. In practice, however, this solution is often impractical due to the substantial computational burden of obtaining the posterior samples and evaluating the relevant quantities in the bridge sampling procedure.

Frühwirth-Schnatter (2004) proposed an alternative approach that approximates the estimator’s expected relative mean-squared error: The derivation of this approximate relative mean-squared error by Frühwirth-Schnatter takes into account that the samples from the proposal distribution are independent, whereas the MCMC samples from the posterior distribution may be autocorrelated. The approximate relative mean-squared error is given by: where , , denotes the variance of with respect to the proposal distribution (the variance is defined analogously), and corresponds to the normalized spectral density of the autocorrelated process at the frequency 0.

In practice, we approximate the unknown variances and expected values by the corresponding sample variances and means. Hence, for evaluating the variance and expected value with respect to , we use the samples for from the proposal distribution. To evaluate the variance and expected value with respect to the posterior distribution, we use the second batch of samples from the posterior distribution which we also use in the iterative scheme for computing the marginal likelihood. Because the posterior samples are obtained via an MCMC procedure and are hence autocorrelated, the second term in Eq. (17) is adjusted by the normalized spectral density (for details see Frühwirth-Schnatter, 2004).13 To evaluate the normalized posterior density which appears in the numerator of and the denominator of both and , we use the bridge sampling estimate as normalizing constant.

Note that, under the assumption that the bridge sampling estimator is an unbiased estimator of the marginal likelihood , the square root of the expected relative mean-squared error (Eq. (16)) can be interpreted as the coefficient of variation (i.e., the ratio of the standard deviation and the mean; Brown, 1998). In the remainder of this article, we report the coefficient of variation to quantify the accuracy of the bridge sampling estimate.

Case study: bridge sampling for reinforcement learning models

In this section, we illustrate the computation of the marginal likelihood using bridge sampling in the context of a published data set (Busemeyer & Stout, 2002) featuring the Expectancy Valence (EV) model—a popular reinforcement learning (RL) model for the Iowa gambling task (IGT; Bechara et al., 1994). We first introduce the task and the model, and then use bridge sampling to estimate the marginal likelihood of the EV model implemented in both an individual-level and a hierarchical Bayesian framework. For the individual-level framework, we compare estimates obtained from bridge sampling to importance sampling estimates published in Steingroever, Wetzels, and Wagenmakers (2016). For the hierarchical framework, we compare our results to estimates from the Savage–Dickey density ratio test Dickey (1971), Dickey and Lientz (1970), Wagenmakers et al. (2010), Wetzels, Grasman et al. (2010).

The Iowa gambling task

In this section we describe the IGT (see also Steingroever, Pachur, Šmíra, & Lee, submitted for publication; Steingroever, Wetzels, Horstmann, Neumann, & Wagenmakers, 2013; Steingroever, Wetzels, & Wagenmakers, 2013a; Steingroever, Wetzels, & Wagenmakers, 2013b; Steingroever, Wetzels, & Wagenmakers, 2014; Steingroever et al., 2016). Originally, Bechara et al. (1994) developed the IGT to distinguish decision-making strategies of patients with lesions to the ventromedial prefrontal cortex from the ones of healthy controls (see also Bechara, Damasio, Tranel, & Anderson, 1998; Bechara, Damasio, Damasio, & Lee, 1999; Bechara, Tranel, & Damasio, 2000). During the last decades, the scope of application of the IGT has increased tremendously covering clinical populations with, for example, pathological gambling (Cavedini, Riboldi, Keller, D’Annucci, & Bellodi, 2002), obsessive–compulsive disorder (Cavedini, Riboldi, D’Annucci, Belotti, Cisima, & Bellodi, 2002), psychopathic tendencies (Blair, Colledge, & Mitchell, 2001), and schizophrenia Bark et al. (2005), Martino et al. (2007).

The IGT is a card game that requires participants to choose, over several rounds, cards from four different decks in order to maximize their long-term net outcome Bechara et al. (1994), Bechara et al. (1997). The four decks differ in their payoffs, and two of them result in negative long-term outcomes (i.e., the bad decks), whereas the remaining two decks result in positive long-term outcomes (i.e., the good decks). After each choice, participants receive feedback on the rewards and losses (if any) associated with that card, as well as their running tally of net outcomes over all trials so far. Unbeknownst to the participants, the task (typically) contains 100 trials.

A crucial aspect of the IGT is whether and to what extent participants eventually learn to prefer the good decks because only choosing from the good decks maximizes their long-term net outcome. The good decks are typically labeled as decks C and D, whereas the bad decks are labeled as decks A and B. Table 1 presents a summary of the traditional payoff scheme as developed by Bechara et al. (1994). This table illustrates that decks A and B yield high constant rewards, but even higher unpredictable losses: hence, the long-term net outcome is negative. Decks C and D, on the other hand, yield low constant rewards, but even lower unpredictable losses: hence, the long-term net outcome is positive. In addition to the different payoff magnitudes, the decks also differ in the frequency of losses: decks A and C yield frequent losses, while decks B and D yield infrequent losses.

Table 1

Summary of the Payoff Scheme of the Traditional IGT as Developed by Bechara et al. (1994).

	Deck A	Deck B	Deck C	Deck D
	Bad deckwith infrequent losses	Bad deckwith frequent losses	Good deckwith infrequent losses	Good deckwith frequent losses
Reward/trial	100	100	50	50
Number of losses/10 cards	5	1	5	1
Loss/10 cards	−1250	−1250	−250	−250
Net outcome/10 cards	−250	−250	250	250

The expectancy valence model

In this section, we describe the EV model (see also Steingroever et al., 2013a; Steingroever et al. (2014), Steingroever et al. (2016), Steingroever et al. (submitted for publication)). Originally proposed by Busemeyer and Stout (2002), the EV model is arguably the most popular model for the IGT (for references see Steingroever et al., 2013a, and for alternative IGT models see Ahn, Busemeyer, Wagenmakers, & Stout, 2008; Dai, Kerestes, Upton, Busemeyer, & Stout, 2015; Steingroever et al., 2014; Worthy & Maddox, 2014; Worthy, Pang, & Byrne, 2013). The model formalizes participants’ performance on the IGT through the interaction of three model parameters that represent distinct psychological processes. The first model assumption is that after choosing a card from deck k, k , on trial t, participants compute a weighted mean of the experienced reward W(t) and loss L(t) to obtain the utility of deck k on trial t, v : The weight that participants assign to losses relative to rewards is the attention weight parameter w. A small value of w, that is, , is characteristic for decision makers who put more weight on the immediate rewards and can thus be described as reward-seeking, whereas a large value of w, that is, , is characteristic for decision makers who put more weight on the immediate losses and can thus be described as loss-averse Ahn et al. (2008), Busemeyer and Stout (2002).

The EV model further assumes that decision makers use the utility of deck k on trial t, , to update only the expected utility of deck k, ; the expected utilities of the unchosen decks are left unchanged. This updating process is described by the Delta learning rule, also known as the Rescorla–Wagner rule (Rescorla & Wagner, 1972) : If the experienced utility is higher than expected, the expected utility of deck k is adjusted upward. If the experienced utility is lower than expected, the expected utility of deck k is adjusted downward. This updating process is influenced by the second model parameter—the updating parameter a. This parameter quantifies the memory for rewards and losses. A value of a close to zero indicates slow forgetting and weak recency effects, whereas a value of a close to one indicates rapid forgetting and strong recency effects. For all models, we initialized the expectancies of all decks to zero, (k ). This setting reflects neutral prior knowledge about the payoffs of the decks.

In the next step, the model assumes that the expected utilities of each deck guide participants’ choices on the next trial . This assumption is formalized by the softmax choice rule, also known as the ratio-of-strength choice rule (Luce, 1959) : The EV model uses this rule to compute the probability of choosing each deck on each trial. This rule contains a sensitivity parameter that indexes the extent to which trial-by-trial choices match the expected deck utilities. Values of close to zero indicate random choice behavior (i.e., strong exploration), whereas large values of indicate choice behavior that is strongly determined by the expected utilities (i.e., strong exploitation).

The EV model uses a trial-dependent sensitivity parameter , which also depends on the final model parameter, response consistency : If is positive, successive choices become less random and more determined by the expected deck utilities; if is negative, successive choices become more random and less determined by the expected deck utilities, a pattern that is clearly non-optimal. We restricted the consistency parameter of the EV model to the range instead of the proposed range  (Busemeyer & Stout, 2002). This modification improved the estimation of the EV model and prevented the choice rule from producing numbers that exceed machine precision (see also Steingroever et al., 2014).

In sum, the EV model has three parameters: (1) the attention weight parameter , which quantifies the weight of losses over rewards; (2) the updating parameter , which determines the memory for past expectancies; and (3) the response consistency parameter , which determines the balance between exploitation and exploration.

Data

We applied bridge sampling to a data set published by Busemeyer and Stout (2002). The data set consists of 30 healthy participants each contributing IGT card selections (see Busemeyer and Stout for more details on the data sets).14

Application of bridge sampling to an individual-level implementation of the EV model

In this section we describe how we use bridge sampling to estimate the marginal likelihood of an individual-level implementation of the EV model. This implementation estimates model parameters for each participant separately. Accordingly, we also obtain a marginal likelihood of the EV model for every participant.

Schematic execution of the bridge sampler

To obtain the bridge sampling estimate of the marginal likelihood for each participant, we follow the steps outlined in Fig. 5.

For each participant , , we proceed as follows:

For each parameter, we draw samples from the posterior distribution.

Since Steingroever et al. (2016) already fit an individual-level implementation of the EV model separately to the data of each participant in Busemeyer and Stout (2002), we reuse their posterior samples (see Steingroever et al., 2016, for details on the prior distributions and model implementation). Note that they parameterized the model not in terms of , but in terms of , and in the remainder of this article, we also use this reparameterization.

For each participant, we choose to match the number of samples obtained from Steingroever et al. (2016) which was at least 5000; however, whenever this number of samples was insufficient to ensure convergence of the Hamiltonian Monte Carlo (HMC) chains, Steingroever et al. (2016) repeated the fitting routine with 5000 additional samples. Steingroever et al. (2016) confirmed convergence of the HMC chains by reporting that all statistics were below 1.05.

We choose a proposal distribution.

We generalize our approach from the running example and use a multivariate normal distribution as a proposal distribution.

We transform the first batch of posterior samples.

Since we use a multivariate normal distribution as a proposal distribution, we have to transform all posterior samples to the real line using the probit transformation, that is, , , , .

We fit the proposal distribution to the first batch of probit-transformed posterior samples.

We use method of moment estimates for the mean vector and the covariance matrix obtained from the first batch of probit-transformed posterior samples to specify our multivariate normal proposal distribution.

We draw samples from the proposal distribution.

We use the R software to randomly draw samples from the proposal distribution obtained in step 4. We obtain with .

We calculate for all samples from the proposal distribution.

This step involves assessing the value of the unnormalized posterior and the proposal distribution for all samples from the proposal distribution. Before we can assess the value of the unnormalized posterior (i.e., the product of the likelihood and the prior), we have to derive how our transformation in step 3 affects the unnormalized posterior.

First, we derive how our transformation affects the likelihood. To evaluate the likelihood, we need to transform the probit-transformed samples back to the original parameter scales. That is, we evaluate the likelihood for . Before formalizing the likelihood of the observed choices of participant , we define the following variables: We define as a vector containing the sequence of choices made by participant up to and including trial , and as a vector containing the corresponding sequence of net outcomes. We now obtain the following expression for the likelihood of the observed choices of participant :

Here is the total number of trials, is the probability of choosing deck on trial , and is a dummy variable which is 1 if deck is chosen on trial and 0 otherwise.

Second, we have to derive how our transformation affects the priors on each EV model parameter to yield priors on the probit-transformed model parameters. Since Steingroever et al. (2016) used independent uniform priors on we obtain standard normal priors on the probit-transformed model parameters (see beta-binomial example and Appendix D for an explanation).

We transform the second batch of posterior samples.

This is analogous to step 2.

We calculate for the second batch of probit-transformed samples from the posterior distribution.

This is analogous to step 6.

We run the iterative scheme (Eq. (15)) until our predefined tolerance criterion is reached.

We use a tolerance criterion and initialization analogous to the running example. Once convergence is reached, we receive an estimate of the marginal likelihood for each participant, and derive the coefficient of variation for each participant using Eq. (17). The largest coefficient of variation is suggesting that the bridge sampler has low variance.15

Assessing the accuracy of our implementation

To assess the accuracy of our implementation, we compared the marginal likelihood estimates obtained with our bridge sampler to the estimates obtained with importance sampling (Steingroever et al., 2016). Fig. 6 shows the log marginal likelihoods for the 30 participants of Busemeyer and Stout (2002) obtained with bridge sampling (-axis) and importance sampling reported by Steingroever et al. (2016; -axis). The two sets of estimates correspond almost perfectly. These results indicate a successful implementation of the bridge sampler. Thus, this section emphasizes that both the importance sampler and bridge sampler can be used to estimate the marginal likelihood for the data of individual participants. However, when we want to estimate the marginal likelihood of a Bayesian hierarchical model, it may be difficult to find a suitable importance density. The bridge sampler, on the other hand, can be applied more easily and more efficiently.

Fig. 6

Comparison of the log marginal likelihoods obtained with bridge sampling (-axis) and importance sampling reported by Steingroever et al. (2016) (-axis). The main diagonal indicates perfect correspondence between the two methods. Available at https://tinyurl.com/yac3o8qs under CC license https://creativecommons.org/licenses/by/2.0/.

Application of bridge sampling to a hierarchical implementation of the EV model

In this section we illustrate how bridge sampling can be used to estimate the marginal likelihood of a hierarchical EV model. This hierarchical implementation assumes that the parameters , , and from each participant are drawn from three separate group-level distributions. This model specification hence incorporates both the differences and the similarities between participants. We illustrate this application using again the Busemeyer and Stout (2002) data set, and assume that these participants constitute one group.

To compute the marginal likelihood, we again follow the steps outlined in Fig. 5, with a few minor modifications.

For each parameter, that is, all individual-level and group-level parameters, we draw samples from the posterior distribution.

To obtain the posterior samples, we fit a hierarchical Bayesian implementation of the EV model to the Busemeyer and Stout (2002) data set using the software JAGS (Plummer, 2003).16 We assume that, for each participant , , each probit-transformed individual-level parameter (i.e., , , ) is drawn from a group-level normal distribution characterized by a group-level mean and standard deviation parameter. For all group-level mean parameters we assume a standard normal distribution, and for all group-level standard deviation parameters a uniform distribution ranging from to 1.5. For a detailed explanation of the hierarchical implementation of the EV model, see Wetzels, Vandekerckhove et al. (2010).

To reach convergence and reduce autocorrelation, we collect two MCMC chains, each with 120,000 samples from the posterior distributions after having excluded the first 30,000 samples as burn-in. Out of these 120,000 samples per chain, we retained every th value yielding 30,000 samples per chain. This setting resulted in all statistics below 1.05 suggesting that all chains have successfully converged from their starting values to their stationary distributions.

We choose a proposal distribution.

We use a multivariate normal distribution as a proposal distribution.

We transform the first batch of posterior samples.

As before, we ensure that the range of the posterior distribution matches the range of the proposal distribution by using the probit transformation, that is, , , , , , and , . The group-level mean parameters do not have to be transformed because they already range across the entire real line.

We fit the proposal distribution to the first batch of the probit-transformed posterior samples.

We use method of moment estimates for the mean vector and the covariance matrix obtained from the first batch of probit-transformed posterior samples to specify our multivariate normal proposal distribution.

We draw samples from the proposal distribution.

We use the R software to randomly draw samples from the proposal distribution obtained in step 4. We obtain and with and .

We calculate for all samples from the proposal distribution.

This step involves assessing the value of the unnormalized posterior and the proposal distribution for all samples from the proposal distribution. The unnormalized posterior is defined as:

, where refers to all choices of subject , to the net outcomes that subject experienced on trials , to the th sample from the proposal distribution for the individual-level parameters of subject , and to the th sample from the proposal distribution for all group-level parameters (e.g., ).

The likelihood function for a given participant is the same as in the individual case. However, for each participant we now have to add besides the prior on the individual-level parameters also the prior on the group-level parameters. The product of the likelihood and the priors gives the unnormalized posterior density (see Appendix E for more details).

We follow steps 7–9, as outlined for the bridge sampler of the individual-level implementation of the EV model.

To investigate the accuracy of our implementation, we compare Bayes factors obtained with bridge sampling to Bayes factors obtained from the Savage–Dickey density ratio test (Dickey (1971), Dickey and Lientz (1970); for a tutorial, see Wagenmakers et al., 2010). The Savage–Dickey density ratio is a simple method for computing Bayes factors for nested models. We artificially create three nested models by taking the full EV model in which all parameters are free to vary, and then restricting one of the three group-level mean parameters, that is, , , or , to a predefined value. For these values we choose the intersection point of the prior and posterior distribution of each group-level mean parameter. To obtain these intersection points, we fit the full EV model and then use a nonparametric logspline density estimator (Stone, Hansen, Kooperberg, Truong, et al., 1997). The obtained values are presented in Table 2. Since we compare the full model to each restricted model, we obtain three Bayes factors.

Table 2

Bayes factors comparing the full EV model to the restricted EV models, log marginal likelihoods, and coefficient of variation (with respect to the marginal likelihood) expressed as a percentage.

Model	Bayes Factor	Log Marginal Likelihood	CV[%]
Full model	–	−3800.434	10.13
Restricted at μ_ω=−0.334	1.202	−3800.618	16.44
Restricted at μ_α=−0.604	1.052	−3800.484	9.71
Restricted at μ_γ=0.92	1.068	−3800.500	12.03

According to the Savage–Dickey density ratio test, the Bayes factor for the full model versus a specific restricted model can be obtained by dividing the height of the prior density at the predefined parameter value by the height of the posterior at the same location:

Since we choose to be the intersection point of the prior and posterior distribution, equals 1. This Savage–Dickey Bayes factor of 1 indicates that the marginal likelihood under the full model equals the marginal likelihood under the restricted model. Fig. 7 illustrates the Savage–Dickey Bayes factor comparing the full model to the model assuming fixed to .

Fig. 7

Prior and posterior distribution of the group-level mean in the Busemeyer and Stout (2002) data set. The figure shows the posterior distribution (solid line) and the prior distribution (dotted line). The gray dot indicates the intersection of the prior and the posterior distributions, for which the Savage–Dickey Bayes factor equals . Available at https://tinyurl.com/y7cyxclq under CC license https://creativecommons.org/licenses/by/2.0/.

The computation of the three bridge sampling Bayes factors, on the other hand, works as follows: First, we follow the steps outlined above to obtain the bridge sampling estimate of the full EV model. Second, we obtain the bridge sampling estimate of the marginal likelihood for the three restricted models. This requires adapting the steps outlined above to each of the three restricted models. Lastly, we use the first equality in Eq. (19) to obtain the three Bayes factors.

The Bayes factors derived from bridge sampling are reported in Table 2. It is evident that Bayes factors derived from bridge sampling closely approximate the Savage–Dickey Bayes factors of 1. These results suggest a successful implementation of the bridge sampler. This is also reflected by the close match between the log marginal likelihoods of the four models presented in the third column of Table 2.

Finally, we confirm that the bridge sampler has low variance; the coefficient of variation with respect to the marginal likelihood of the full model and the three restricted models ranges between 9.71 and .

Discussion

In this tutorial, we explained how bridge sampling can be used to estimate the marginal likelihood of popular models in mathematical psychology. As a running example, we used the beta-binomial model to illustrate step-by-step the bridge sampling estimator. To facilitate the understanding of the bridge sampler, we first discussed three of its special cases—the naive Monte Carlo estimator, the importance sampling estimator, and the generalized harmonic mean estimator. Consequently, we introduced key concepts that became gradually more complicated and sophisticated. In the second part of this tutorial, we showed how bridge sampling can be used to estimate the marginal likelihood of both an individual-level and a hierarchical implementation of the Expectancy Valence (EV; Busemeyer & Stout, 2002) model—a popular reinforcement-learning model for the Iowa gambling task (IGT; Bechara et al., 1994). The running example and the application of bridge sampling to the EV model demonstrated the positive aspects of the bridge sampling estimator, that is, its accuracy, reliability, practicality, and ease-of-implementation DiCiccio et al. (1997), Frühwirth-Schnatter (2004), Meng and Wong (1996).

The bridge sampling estimator is superior to the naive Monte Carlo estimator, the importance sampling estimator, and the generalized harmonic mean estimator for several reasons. First, Meng and Wong (1996) showed that, among the four estimators discussed in this article, the bridge sampler presented in this article minimizes the mean-squared error because it uses the optimal bridge function. Second, in bridge sampling, choosing a suitable proposal distribution is much easier than choosing a suitable importance density for the importance sampling estimator or the generalized harmonic mean estimator because bridge sampling is more robust to the tail behavior of the proposal distribution relative to the posterior distribution. This advantage facilitates the application of the bridge sampler to higher dimensional and complex models. This characteristic of the bridge sampler combined with the popularity of higher dimensional and complex models in mathematical psychology suggests that bridge sampling can advance model comparison exercises in many areas of mathematical psychology (e.g., reinforcement-learning models, response time models, multinomial processing tree models, etc.). Third, bridge sampling is relatively straightforward to implement. In particular, our step-by-step procedure can be easily applied to other models with only minor changes of the code (i.e., the unnormalized posterior and potentially the proposal function have to be adapted). In our opinion, this is one of the most attractive features of bridge sampling: It is an accurate yet very generic method. Exploiting this generic characteristic, we have implemented the bridge sampling procedure in the bridgesampling R package (Gronau, Singmann, & Wagenmakers, 2017) in order to maximize its accessibility.

Despite the numerous advantages of the bridge sampler, the take-home message of this tutorial is not that the bridge sampler should be used blindly. There exist a large variety of methods to approximate the marginal likelihood that differ in their efficiency.17 The most appropriate method optimizes the trade-off between accuracy and implementation effort. This trade-off depends on a number of aspects such as the complexity of the model, the number of models under consideration, the statistical experience of the researcher, and the time available. This suggests that the choice of the method should be reconsidered each time a marginal likelihood needs to be obtained. Obviously, when the marginal likelihood can be determined analytically, bridge sampling is not needed at all. If the goal is to compare (at least) two nested models, the Savage–Dickey density ratio test Dickey (1971), Dickey and Lientz (1970) might be a better alternative. Note, however, that this requires an approximation of the marginal posterior density of one or more parameters which can be unstable in case the test value falls in the tail of the distribution. If only an individual-level implementation of a model is used, importance sampling may be easier to implement and may require less computational effort. This presupposes that one can find a proposal distribution with fatter tails than the posterior which may not always be trivial (even in an individual-level case). If the goal is to obtain the marginal likelihood of a large number of relatively simple models, the product space or reversible jump method (RJMCMC) might be more appropriate Carlin and Chib (1995), Green (1995), Lodewyckx et al. (2011). In contrast to bridge sampling, implementations of these methods tend to be problem-specific rather than generic (but see Lunn, Best, & Whittaker, 2009). If a researcher with a limited programming background and/or little time resources wants to conduct a model comparison exercise, rough approximations of the Bayes factor, such as the Bayesian information criterion, might be more suitable (Schwarz, 1978). It should be kept in mind, however, that this approximation assumes a certain prior structure that may not respect the knowledge or information that researchers have at their disposal. On the other hand, a researcher with an extensive background in programming and mathematical statistics might consider using path sampling—a generalization of bridge sampling (Gelman & Meng, 1998).

To conclude, in this tutorial we showed that bridge sampling offers a reliable and easy-to-implement approach to estimating a model’s marginal likelihood. Bridge sampling can be profitably applied to a wide range of problems in mathematical psychology involving parameter estimation, model comparison, and Bayesian model averaging.

Table A.1

Summary of the bridge sampling estimators for the marginal likelihood, and its special cases: the naive Monte Carlo, importance sampling, and generalized harmonic mean estimator.

Method	Estimator	Samples	Bridge Function h(θ)
Bridge sampling	1N_2∑i=1N_2p(y∣θ~_i)p(θ~_i)h(θ~_i)1N_1∑j=1N_1h(θj∗)g(θj∗)	θ~_i∼g(θ)	h(θ)=CN_1N_2+N_1p(y∣θ)p(θ)+N_2N_2+N_1p(y)g(θ)
		θj∗∼p(θ∣y)
Naive Monte Carlo	1N∑i=1Np(y∣θ~_i)	θ~_i∼p(θ)	h(θ)=1g(θ) and g(θ)=p(θ)
Importance sampling	1N∑i=1Np(y∣θ~_i)p(θ~_i)g_IS(θ~_i)	θ~_i∼g_IS(θ)	h(θ)=1g(θ) and g(θ)=g_IS(θ)
Generalized harmonic mean	1N∑i=1Ng_IS(θi∗)p(y∣θi∗)p(θi∗)^−1	θi∗∼p(θ∣y)	h(θ)=1p(y∣θ)p(θ) and g(θ)=g_IS(θ)

Note. is the prior distribution, is the importance density, is the posterior distribution, is the proposal distribution, is the bridge function, and is a constant. The last column shows the bridge function needed to obtain the special cases.