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.

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 y. This is realized by means of a posterior distribution that can be obtained using Bayes theorem: Here, the marginal likelihood of the data p(y) 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 m (m ∈ N) competing models, and are interested in the relative plausibility of a particular model M i (i ∈ {1, 2, . . . , m}) given the prior model probability and the evidence from the data y (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 p(M i | y) of model M i given the data y (Berger & Molina, 2005): where the denominator is the sum of the marginal likelihood times the prior model probability of all m 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) Equation 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 & Wagenmakers, in press;Jeffreys, 1961;Ly, Verhagen, & Wagenmakers, 2016a, 2016bRobert, 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, 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 y given a specific model of interest M, and is defined as the integral of the likelihood over the prior: with θ a vector containing the model parameters. Equation 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.
Equation 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 k = 2 out of n = 10 true/false questions correctly. Assume that the number of correct answers follows a binomial distribution, that is, k ∼ Binomial(n, θ) with θ ∈ (0, 1), where θ represents the latent probability for Prior and posterior distribution for the rate parameter θ from the beta-binomial model.
The Beta(1, 1) prior on the rate parameter θ is represented by the dotted line; the Beta(3, 9) posterior distribution is represented by the solid line and was obtained after having observed 2 correct responses out of 10 trials.
answering any one question correctly. The probability mass function (PMF) of the binomial distribution is given by: where k, n ∈ Z ≥0 , and k ≤ n. 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, θ ∼ Uniform(0, 1) which can equivalently be written in terms of a beta distribution θ ∼ Beta(1, 1). This prior distribution is represented by the dotted line in Figure 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 Beta(α, β) distribution (α, β > 0) is given by: where B(α, β) is the beta function that is defined as: Γ(α+β) , and Γ(n) = (n − 1)! for n ∈ N. 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 k and n, 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 k = 2 and n = 10 of our example, we obtain: p(k = 2 | n = 10) = 1/11 ≈ 0.0909. 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 (Equation 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 (Equation 6) and Bayes theorem, we obtain: which we recognize as the PDF of the Beta(k + 1, n − k + 1) distribution. Thus, if we assume a uniform prior on θ and observe k = 2 correct responses out of n = 10 trials, we obtain a Beta(3, 9) distribution as posterior distribution. This distribution is represented by the solid line in Figure 1. In general, if k | n, θ ∼ Binomial(n, θ) and θ ∼ Beta(1, 1), then θ | n, k ∼ Beta(k + 1, n − k + 1).

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 & Handscomb, 1964;Raftery & Banfield, 1991 This yields the naive Monte Carlo estimatorp 1 (y): Running example. To obtain the naive Monte Carlo estimate of the marginal likelihood in our running example, we need N samples from the Beta(1, 1) prior distribution for θ. For illustrative purposes, we limit the number of samples to 12 whereas in practice one should take N to be very large. We obtain the following samples:

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 & 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 Equation 4 is large. This is realized by using samples from a so-called importance density g IS (θ) 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, Equation 4 is used as starting point which is then extended by the importance density g IS (θ): This yields the importance sampling estimatorp 2 (y): samples from the importance density 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, Matzke, & Wagenmakers, 2015). The latter criterion ensures that values in the tails of the distribution cannot misleadingly dominate the estimate (Neal, 2001). 3

Running example.
To obtain the importance sampling estimate of the marginal likelihood in our running example, we first need to choose an importance density g IS (θ). 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 that regions with low likelihood values are eliminated, thereby increasing the accuracy and efficiency of the estimator. 3 To illustrate the need for an importance density with fatter tails than the posterior distribution, imagine you sample from the tail region of an importance density with thinner tails. In this case, the numerator in Equation 8 would be substantially larger than the denominator resulting in a very large ratio. Since this specific ratio is only one component of the sum displayed in Equation 8, this component would dominate the importance sampling estimate. Hence, thinner tails of the importance density run the risk of producing unstable estimates across repeated computations. In fact, the estimator may have infinite variance (e.g., Ionides, 2008;Owen & Zhou, 2000). density is quantified by a mixture weight γ that ranges between 0 and 1. The larger γ, the higher the influence of the uniform density resulting in a less peaked distribution with thick tails. If γ = 1, the beta mixture density simplifies to the uniform distribution on [0, 1]; 4 and if γ = 0, 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 Beta(3, 9) 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 Beta(3, 9) posterior distribution and obtain: 5 Note that here we use θ * i to refer to the i th sample from the posterior distribution to distinguish it from the previously usedθ i -the i 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θ * = 0.232 and a variance of s 2 θ * = 0.014. Third, knowing that, if X ∼ Beta(α, β), then E(X) = α/(α + β) and V (X) = αβ/ (α + β) 2 (α + β + 1) , we obtain the following method of moment estimates for α and Using a mixture weight on the uniform component of γ = 0.30-a choice that was made to ensure that, visually, the tails of the importance density are clearly thicker than the tails of  Figure 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 {θ 1 ,θ 2 , . . . ,θ 12 } randomly drawn from our beta mixture importance density. Available at https://tinyurl.com/yc7ho7hr under CC license https://creativecommons.org/licenses/by/2.0/.
As is evident from the figure, the beta mixture importance density resembles the posterior distribution, but has fatter tails.
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 N = 12 samples for θ from our beta mixture importance density results in: These samples are represented by the gray dots in Figure 3.
The final step is to compute the average adjusted likelihood for the 12 samples using Equation 8. This yields the importance sampling estimate of the marginal likelihood as: = 0.0827.

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 Equation 4 is large by using an importance density g IS (θ) (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: Note that the generalized harmonic mean estimator is a more stable version of the harmonic mean estimator (Newton & Raftery, 1994). A problem of the harmonic mean estimator is that it is dominated by the samples that have small likelihood values.
which is used to define the generalized harmonic mean estimatorp 3 (y) (Gelfand & Dey, 1994) as follows:p samples from the posterior distribution (9) 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 Equation 9 is the reciprocal of the ratio in Equation 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 Equation 9. Thus, in the case of the generalized harmonic mean estimator, a suitable importance density should (1) have thinner tails than the posterior distribution (Newton & Raftery, 1994;DiCiccio et al., 1997), 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.

Running example.
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 N = 12 samples from the posterior distribution. Reusing the samples from the last section, we obtain: The result of this transformation is that the samples range across the entire real line instead of the (0, 1) interval only. We obtain: ξ Density Figure 4 . Illustration of the generalized harmonic mean estimator for the beta-binomial model. The solid line represents the probit-transformed Beta(3, 9) posterior distribution that was obtained after having observed 2 correct responses out of 10 trials, and the dashed line represents the importance density N ξ; µ = −0.793, σ = 0.423/1.5 . The gray dots represent the 12 probit-transformed samples {ξ * 1 , ξ * 2 , . . . , ξ * 12 } randomly drawn from the Beta(3, 9) 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 probittransformed posterior samples ξ * j . Using the method of moments, we obtain as estimateŝ µ = −0.793 andσ = 0.423. 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μ = −0.793 and standard deviationσ = 0.423/1.5. 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 Figure   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 θ * j or the probit-transformed samples ξ * j . 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, Equation 9 becomes: 8 probit-transformed samples from the posterior distribution (10) 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 M 1 and M 2 (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 8 A detailed explanation is provided in the appendix. Note that using the original posterior samples θ * j would involve transforming the importance density (e.g., the normal density on ξ) to the (0, 1) interval.
for example Overstall & Forster, 2010). This version is based on the following identity: where g(θ) is the so-called proposal distribution and h(θ) the so-called bridge function.
Multiplying both sides of Equation 11 by the marginal likelihood p(y) results in: .
The marginal likelihood can now be approximated using: samples from the posterior distribution (12) Equation 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 Equation 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 s 1 = N 1 N 2 +N 1 , s 2 = N 2 N 2 +N 1 , and C a constant; its particular value is not required because h(θ) is part of both the numerator and the denominator of Equation 12, and therefore the constant C 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 (Equation 16).
Equation 13 shows that the optimal bridge function depends on the marginal likelihood p(y) 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 (Equation 13) in Equation 12 (Meng & Wong, 1996). The formula to approximate the marginal likelihood on iteration t + 1 is then specified as follows: wherep(y) (t) denotes the estimate of the marginal likelihood on iteration t of the iterative scheme. Note that Equation 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., Equation 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., Equation 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., Equation 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 1 s 2p (y) (t) . 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., Equation 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., Equation 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., Equation 14), the ratio in the numerator will go to 1/s 1 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 Equation 14 would be zero.
Equation 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 N 2 samples from the proposal distribution g(θ) and N 1 samples from the posterior distribution p(θ|y). Second, for all N 2 samples from the proposal distribution, we have to evaluate l 2,i . 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 l 1,j for all N 1 samples from the posterior distribution. This is analogous to evaluating l 2,i . Fourth, we have to determine the constants s 1 and s 2 that only depend on N 1 and N 2 . Fifth, we need an initial guess of the marginal likelihoodp 4 (y). 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 (Equation 15).
A second challenge is how to use the N 1 samples from the posterior distribution. One option is to use all N 1 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 2N 1 samples from the posterior distribution.
The first N 1 of the total of 2N 1 samples are used for fitting the proposal distribution and the remaining N 1 samples are used in the iterative scheme (i.e., Equation 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 Running example. To obtain the bridge sampling estimate of the marginal likelihood in the beta-binomial example, we follow the eight steps illustrated in Figure 5: 1. We draw 2N 1 = 24 samples from the Beta(3, 9) posterior distribution for θ.
We obtain the following sample of 24 values: .22, 0.16, 0.09, 0.35, 0.06, 0.27, 0.26, 0.41, 0.20, 0.43, 0.21, 0.12, 9 In case the posterior samples are obtained via MCMC sampling using multiple chains, we use the first half of the iterations per chain for fitting the proposal distribution and the second half of the iterations per chain for the iterative scheme.
10 For an explanation of where the name "bridge" comes from see https://osf.io/9jzm3/. Note that the first 12 samples equal the ones used in the last section, whereas the last 12 samples were obtained from drawing again 12 values from the Beta(3, 9) 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 g(θ). 11 3. We transform the first batch of N 1 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 probittransforming the posterior samples, that is, ξ * j = Φ −1 (θ * j ) with j ∈ {1, 2, . . . , 12}. We obtain: 6. We calculate l 2,i for all N 2 samples from the proposal distribution.
This step involves assessing the value of the unnormalized posterior and the proposal distribution for all N 2 samples from the proposal distribution. As in the running example for the generalized harmonic mean estimator, we obtain the unnormalized 11 There exist several candidates for the proposal distribution. Alternative proposal distributions are, for example, the importance density that we used for the importance sampling estimator or for the generalized harmonic mean estimator, or the analytically derived Beta(3, 9) posterior distribution.
posterior as: p k = 2 | n = 10, Φ ξ i φ ξ i , where φ ξ i comes from using the change-of-variable method (see running example for the generalized harmonic mean estimator and the appendix 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). 7. We transform the second batch of N 1 posterior samples.
As in step 2, we use the probit transformation and obtain: This is analogous to step 6. 9. We run the iterative scheme (Equation 15) until our predefined tolerance criterion is reached.
In addition, we need to calculate l 2,i (i ∈ {1, 2, . . . , 12}) for all samples from the proposal distribution, and l 1,j (j ∈ {1, 2, . . . , 12}) for the second batch of the probit-transformed samples from the posterior distribution. Here we show how to calculate l 2,1 and l 1,1 using the first sample from the proposal distribution and the first sample of the second batch of the posterior samples, respectively: A better initial guess can be obtained from, for example, the importance sampling estimator or the generalized harmonic mean estimator explained in the previous sections. In our experience, however, usually the exact choice of the initial value does not seem to influence the convergence of the bridge sampler much. Forp 4 (k = 2 | n = 10) (t+1) , we then get: Usingp(y) (0) = 0, we obtain as updated estimate of the marginal likelihoodp 4 (k = 2 | n = 10) (1) = 0.0908. 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 aŝ p 4 (k = 2 | n = 10) (5) = 0.0902.

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 g(θ) are independent, whereas the MCMC samples from the posterior distribution p(θ|y) may be autocorrelated.
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 g(θ), we use the N 2 samples forθ i from the proposal distribution. To evaluate the variance and expected value with respect to the posterior distribution, we use the second batch of N 1 samples θ * j 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 Equation 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 f 1 (θ) and the denominator of both f 1 (θ) and f 2 (θ), we use the bridge sampling estimate as 13 We estimate the spectral density at frequency zero by fitting an autoregressive model using the spectrum0.ar() function as implemented in the coda R package (Plummer, Best, Cowles, & Vines, 2006). normalizing constant.
Note that, under the assumption that the bridge sampling estimatorp 4 (y) is an unbiased estimator of the marginal likelihood p(y), the square root of the expected relative mean-squared error (Equation 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 & Lientz, 1970;Wagenmakers, Lodewyckx, Kuriyal, & Grasman, 2010;Wetzels, Grasman, & Wagenmakers, 2010).
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, Damasio, Tranel, & Damasio, 1997). The four decks differ in their Table 1 Summary of the Payoff Scheme of the Traditional IGT as Developed by Bechara et al. (1994) 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.

The Expectancy Valence Model
In this section, we describe the EV model (see also Steingroever, Wetzels, & Wagenmakers, 2013a;Steingroever et al., 2014Steingroever et al., , 2016. Originally proposed by Busemeyer and Stout (2002), the EV model is arguably the most popular model for the IGT (for references see Steingroever, Wetzels, & Wagenmakers, 2013a, and for alternative IGT models see Ahn, Busemeyer, Wagenmakers, & Stout, 2008;Dai, Kerestes, Upton, Busemeyer, & Stout, 2015;Steingroever et al., 2014;Worthy, Pang, & Byrne, 2013;Worthy & Maddox, 2014). 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 ∈ {1, 2, 3, 4}, 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 k (t): The weight that participants assign to losses relative to rewards is the attention weight parameter w. A small value of w, that is, w < .5, 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, w > .5, 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 & Stout, 2002).
The EV model further assumes that decision makers use the utility of deck k on trial t, v k (t), to update only the expected utility of deck k, Ev k (t); 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 v k (t) is higher than expected, the expected utility of deck k is adjusted upward. If the experienced utility v k (t) 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, Ev k (0) = 0 (k ∈ {1, 2, 3, 4}).
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 t + 1. 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 θ(t), which also depends on the final model parameter, response consistency c ∈ [−5, 5]: If c is positive, successive choices become less random and more determined by the expected deck utilities; if c 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 [−2, 2] instead of the proposed range [−5, 5] (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 w ∈ [0, 1], which quantifies the weight of losses over rewards; (2) the updating parameter a ∈ [0, 1], which determines the memory for past expectancies; and (3) the response consistency parameter c ∈ [−2, 2], 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 T = 100 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 remainder of this article, we also use this reparameterization.
For each participant, we choose 2N 1 to match the number of samples obtained from Steingroever et al. (2016) which was at least 5, 000; 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 5, 000 additional samples. Steingroever et al. (2016) confirmed convergence of the HMC chains by reporting that allR 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.
3. We transform the first batch of N 1 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

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

5.
We draw N 2 samples from the proposal distribution.
We use the R software to randomly draw N 2 samples from the proposal distribution obtained in step 4. We obtain (ω s,i ,α s,i ,γ s,i ) with i ∈ {1, 2, . . . , N 2 }.
6. We calculate l 2,i for all N 2 samples from the proposal distribution.
This step involves assessing the value of the unnormalized posterior and the proposal distribution for all N 2 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 (Φ(ω s,i ), Φ(α s,i ), Φ(γ s,i )).
Before formalizing the likelihood of the observed choices of participant s, we define the following variables: We define Ch s (t) as a vector containing the sequence of choices made by participant s up to and including trial t, and X s (t) 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 s: Here T is the total number of trials, P r[S k (t + 1)] is the probability of choosing deck k on trial t + 1, and δ k (t + 1) is a dummy variable which is 1 if deck k is chosen on trial t + 1 and 0 otherwise. This is analogous to step 2.
8. We calculate l 1,j for the second batch of N 1 probit-transformed samples from the posterior distribution.
This is analogous to step 6. 9. We run the iterative scheme (Equation 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 Equation 17. The largest coefficient of variation is 2.07% 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). Figure 6 shows the log marginal likelihoods for the 30 participants of Busemeyer and Stout (2002) obtained with bridge sampling (x-axis) and importance sampling reported by Steingroever et al. (2016;y-axis). The two sets of estimates correspond almost perfectly.
15 Note that this measure relates to the marginal likelihoods, not to the log marginal likelihoods. The main diagonal indicates perfect correspondence between the two methods.
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.

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 w, a, and c 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.

Schematic execution of the bridge sampler.
To compute the marginal likelihood, we again follow the steps outlined in Figure 5, with a few minor modifications.
1. For each parameter, that is, all individual-level and group-level parameters, we draw 2N 1 = 60, 000 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 s, s ∈ {1, 2, . . . , 30}, each probit-transformed individual-level parameter (i.e., ω s = Φ −1 (w s ), α s = Φ −1 (a s ), γ s = Φ −1 (c s )) is drawn from a group-level normal distribution characterized by a grouplevel 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 0 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 fourth value yielding 30, 000 samples per chain. This setting resulted in allR 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.
3. We transform the first batch of N 1 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, The group-level mean parameters do not have to be transformed because they already range across the entire real line.

4.
We fit the proposal distribution to the first batch of the N 1 probit-transformed posterior samples.
We use method of moment estimates for the mean vector and the covariance matrix obtained from the first batch of N 1 probit-transformed posterior samples to specify our multivariate normal proposal distribution. 16 We used a model file that is an adapted version of the model file used by Ahn et al. (2011).

We draw N 2 samples from the proposal distribution.
We use the R software to randomly draw N 2 samples from the proposal distribution obtained in step 4. We obtain (ω s,i ,α s,i ,γ s,i ) and (μ ω,i ,τ ω,i ,μ α,i ,τ α,i ,μ γ,i ,τ γ,i 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 individuallevel 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). 7. We follow steps 7 -9, as outlined for the bridge sampler of the individual-level implementation of the EV model.
Assessing the accuracy of our implementation. 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 & Lientz, 1970;Dickey, 1971; 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 M f 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.
According to the Savage-Dickey density ratio test, the Bayes factor for the full model Since we choose θ 0 to be the intersection point of the prior and posterior distribution,  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 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 16.44%.

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 & 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 & Lientz, 1970;Dickey, 1971) 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 17 In general, a large number of approaches for model selection exist which are based on MCMC posterior sampling and some of them are not based on approximating the models' marginal likelihoods (e.g., Ando, 2007;Spiegelhalter, Best, Carlin, & van der Linde, 2002). A comparison of these methods is beyond the scope of this tutorial.
unstable in case the test value falls in the tail of the distribution. If only an individuallevel 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 & 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.

Appendix A
The Bridge Sampling Estimator as a General Case of Methods 1 -3 In this section we show that the naive Monte Carlo, the importance sampling, and the generalized harmonic mean estimators are special cases of the bridge sampling estimator under specific choices of the bridge function h(θ) and the proposal distribution g(θ). 18 An overview is provided in Table A1.
and g(θ) = p(θ) and g(θ) = g IS (θ) Generalized harmonic mean and g(θ) = g IS (θ) Note. p(θ) is the prior distribution, gIS(θ) is the importance density, p(θ|y) is the posterior distribution, g(θ) is the proposal distribution, h(θ) is the bridge function, and C is a constant. The last column shows the bridge function needed to obtain the special cases.
To prove that the bridge sampling estimator reduces to the naive Monte Carlo estimator, consider bridge sampling, choose the prior distribution as the proposal distribution (i.e., g(θ) = p(θ)), and specify the bridge function as h(θ) = 1/g(θ). Inserting these specifications into Equation 12 yields: which is equivalent to the naive Monte Carlo estimator shown in Equation 7.
To prove that the bridge sampling estimator reduces to the importance sampling estimator, consider bridge sampling, choose the importance density as the proposal distribution (i.e., g(θ) = g IS (θ)), and specify the bridge function as h(θ) = 1/g(θ) . Inserting these specifications into Equation 12 yields: which is equivalent to the importance sampling estimator shown in Equation 8.
To prove that the bridge sampling estimator reduces to the generalized harmonic mean estimator, consider bridge sampling, choose the importance density as the proposal distribution (i.e., g(θ) = g IS (θ)), and specify the bridge function as h(θ) = 1/(p(y | θ) p(θ)).

Appendix B
Bridge Sampling Implementation: Avoiding Numerical Issues In order to avoid numerical issues, we can rewrite Equation 15 in the following way: l 2,i s 1 l 2,i +s 2p4 (y) (t) l * is a constant which we can choose in a way that keeps the terms in the sums manageable.
Hence, we can run the iterative scheme with respect tor which is more convenient because it keeps the terms in the sums manageable and multiply the result by exp(l * ) to obtain the estimate of the marginal likelihood or, equivalently, we can take the logarithm of the result and add l * to obtain an estimate of the logarithm of the marginal likelihood.
Appendix C

Correcting for the Probit Transformation
In this section we describe how the probit transformation affects our expression of the generalized harmonic mean estimator (Equation 9) to yield Equation 10. Recall that we derived the generalized harmonic mean estimator using the following equality: For practical reasons, in the running example, we used a normal distribution on ξ as importance density. This ξ was defined as the probit transform of θ (i.e, ξ = Φ −1 (θ)). In particular, the normal importance density was given by 1 σ φ ξ−μ σ . Note that this importance density is a function of ξ, whereas the general importance density g IS in Equation C1 is specified in terms of θ. Therefore, to include our specific importance density into Equation C1, we need to write it in terms of θ. This yields 1 σ φ Φ −1 (θ)−μ σ 1 φ(Φ −1 (θ)) , where the latter factor comes from applying the change-of-variable method. Replacing g IS (θ) in Equation C1 by this expression, results in: Rewriting results in: which can be approximated as: samples from the posterior distribution probit-transformed samples from the posterior distribution (C3) which shows that the generalized harmonic estimate can be obtained using the samples from the posterior distribution, or the probit-transformed ones. In the online-provided code, we use the latter approach (see also Overstall & Forster, 2010). Note the in our running example, ∀ξ * j : p Φ ξ * j = 1.
Appendix D

Details on the Application of Bridge Sampling to the Individual-Level EV Model
In this section, we provide more details on how we obtained the unnormalized posterior distribution for a specific participant s, s ∈ {1, 2, . . . , 30}. Since we focus on one specific participant, we drop the subscript s in the remainder of this section. As explained in Appendix B, we run the iterative scheme with respect tor to avoid numerical issues.