{'fraction_non_alphanumeric': 0.035718982181603, 'fraction_numerical': 0.01285423104740614, 'mean_word_length': 4.535668789808917, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New Estimation Approaches for the Hierarchical Linear Ballistic Accumulator Model

The Linear Ballistic Accumulator (Brown&Heathcote, 2008) model is used as a measurement tool to answer questions about applied psychology. The analyses based on this model depend upon the model selected and its estimated parameters. Modern approaches use hierarchical Bayesian models and Markov chain Monte-Carlo (MCMC) methods to estimate the posterior distribution of the parameters. Although there are several approaches available for model selection, they are all based on the posterior samples produced via MCMC, which means that the model selection inference inherits the properties of the MCMC sampler. To improve on current approaches to LBA inference we propose two methods that are based on recent advances in particle MCMC methodology; they are qualitatively different from existing approaches as well as from each other. The first approach is particle Metropolis-within-Gibbs; the second approach is density tempered sequential Monte Carlo. Both new approaches provide very efficient sampling and can be applied to estimate the marginal likelihood, which provides Bayes factors for model selection. The first approach is usually faster. The second approach provides a direct estimate of the marginal likelihood, uses the first approach in its Markov move step and is very efficient to parallelize on high performance computers. The new methods are illustrated by applying them to simulated and real data, and through pseudo code. The code implementing the methods is freely available.


Introduction
The Linear Ballistic Accumulator (LBA; Brown and Heathcote, 2008) provides a tractable model of decision making which is simpler than some other models of choice response time because it eliminates complexities such as competition between alternatives (Brown and Heathcote, 2005;Ratcliff, 1978;Ratcliff and Rouder, 1998), and passive decay of evidence (Ratcliff and Smith, 2004;Usher and McClelland, 2001). The model's simplicity allows analytic solutions for choices between any number of alternatives. Like other evidence accumulation models, the LBA model is used to address theoretical and applied questions about human cognition, both in the general population and in clinical groups (for reviews, see e.g., Donkin and Brown, 2018;Ratcliff et al., 2016).
When used as a psychometric tool, key inferences are drawn in two ways: either from parameter estimates or by comparing different versions of the LBA model estimated from the same data. These inferences rely on accurate parameter estimates and valid model selection procedures, which can be difficult problems. Most modern applications of the model use hierarchical structures estimated in a Bayesian framework, with posterior distributions over the parameters estimated using Markov chain Monte Carlo (MCMC). In almost all Bayesian applications of the LBA model, the MCMC uses Metropolis steps, with proposals generated by differential evolution (DE-MCMC: Turner et al., 2013). 1 After sampling from the posterior distribution, inferences about model selection are almost always carried out by estimating a marginal likelihood, or some quantity that behaves approximately like the marginal likelihood, from the MCMC samples. Commonly used model selection metrics include the deviance information criterion (DIC: Spiegelhalter et al., 2014) and the Watanabe (or "widely applicable") information criterion (WAIC: Watanabe, 2010). Model selection using Bayes factors or estimated Bayes factors requires estimating the marginal likelihood of each competing model.  estimate the marginal likelihood by generating a large number of samples from the prior. This method requires specialised computing hardware (a general purpose graphical processing unit) to be computationally feasible even in smaller, nonhierarchical applications. However, obtaining reliable and stable estimation of the marginal likelihood of the hierarchical LBA model remains challenging. Gronau et al. (2019) and Evans and Annis (2019) propose estimating the marginal likelihood by bridge sampling and thermodynamic integration respectively, from the MCMC samples. These two methods are promising even for the hierarchical LBA model, but depend on the quality of the MCMC samples. Both Gronau et al. (2019) and Evans and Annis (2019) use the DE-MCMC algorithm to sample from the posterior distribution of the parameters, which often still suffers from the usual problems associated with random walk samplers for the hierarchical LBA model with a large number of parameters, including a high autocorrelation between samples, and slow or uncertain convergence. When the DE-MCMC samples provide an imperfect representation of the posterior, the subsequent model selection methods using bridge sampling or thermodynamic integration will give incorrect estimates. Both bridge sampling and thermodynamic integration would benefit from a more reliable and efficient sampling algorithm to obtain reliable estimates of the marginal likelihood for model comparison. Furthermore, bridge sampling uses an iterative algorithm to estimate the marginal likelihood; Gronau et al. (2019) noted that in rare cases, a very good starting value for the marginal likelihood is crucial for the algorithm to converge. Thermodynamic integration also requires the MCMC algorithm to sample efficiently at different temperatures, which can be difficult to tune. Section 3.3 further discusses the thermodynamic integration method. It is clear that more efficient and robust sampling methods could be beneficial in many ways.
Our paper makes three substantive contributions. First, we allow the individual level parameters to be correlated in the prior by reparameterising them. An essential component to this parameterisation is to first transform the individual level parameters so they can take all values on the real line which makes it straightforward to specify their joint distribution as an unconstrained multivariate normal with full covariance matrix structure. Previous approaches followed Turner et al. (2013) and assumed that the prior joint distribution of these parameters are uncorrelated truncated (positive only) univariate normal distributions. The new assumptions and parameterisation have the twin advantages of more accurately reflect prior knowledge, and increasing sampling efficiency. Furthermore, Section 4 shows that using a hierarchical LBA model with uncorrelated univariate normal distributions for the transformed random effects can result in overconfidence in estimation precision and underestimation of the magnitude of the individual differences. That section also explains why it is difficult or maybe impossible to efficiently apply the existing DE-MCMC sampler to the new parameterisation of the hierarchical LBA model.
The second and third contributions propose two new methods for estimating the improved LBA model that are more efficient than the DE-MCMC sampler used for the estimation of the LBA since Turner et al. (2013). The first method is based on the particle Metropolis within Gibbs (PMwG) approach of Gunawan et al. (2017). It defines a target posterior density on an augmented space that includes the standard model parameters as well as multiple copies of the individual random effects ("particles") and whose marginal density is the joint posterior density of the parameters and random effects. Section 4 shows that for estimating the LBA model, the PMwG sampler is an important alternative to the DE-MCMC sampler because it converges much more reliably and is much more statistically efficient.
The second estimation method is based on a version of sequential Monte Carlo (SMC; Del Moral et al., 2006) that is an alternative to all the MCMC approaches mentioned above, including DE-MCMC and PMwG. Our approach builds on the work by Neal (2001), Del Moral et al. (2006 and Duan and Fulop (2015) by first drawing samples from an easily-generated distribution, such as the multivariate normal prior, and then moving those samples sequentially towards the posterior distribution. We call our algorithm density tempered sequential Monte Carlo (DT-SMC); it uses three main steps to transition from one intermediate density to the next: a reweighting step that moves the particles from one tempered density to the next; a resampling step that eliminates particles with low weights; and the Markov step that applies several iterates of the new PMwG sampler to help ensure that the particles represent the tempered target density adequately.
Section 4 shows that PMwG is much faster than density tempered SMC; it is also easily implemented with very modest computational resources, e.g. personal computers. However, density tempered SMC is easier to parallelize than PMwG, so that it is likely to be faster when powerful computing resources are available. The density tempered SMC explores the parameter space more efficiently when the target posterior distribution is multimodal; such distributions are usually due to multimodal priors with small sample sizes. Both PMwG and density tempered SMC can be used to estimate the marginal likelihood; the marginal likelihood can be obtained from the density tempered SMC algorithm with negligible extra cost, while the output from the PMwG sampler can be used by both thermodynamic integration and bridge sampling. We believe that it is important to present both the PMwG and density tempered SMC algorithms to give users the ability to estimate LBA models using both modest and considerable computational resources; in addition, future improvements in both bridge sampling and thermodynamic integration may make it attractive to estimate the marginal likelihood by first running PMwG. However, it is outside the scope of the article to carefully compare the performance of all current approaches for estimating the marginal likelihood.
The rest of the paper is organised as follows. Section 2 describes the Linear Ballistic Accumulator model; Section 3 presents the Bayesian estimation methodologies; Section 4 discusses the estimation results where the two new methods are applied to simulated and real data; and Section 5 concludes. The paper has several appendices which contain some further implementation details and technical results. An online supplement containing some further empirical and technical results and code applying the two estimation methods to an example dataset is available at osf.io/5b4w3.

The Linear Ballistic Accumulator (LBA) model
To more precisely discuss the algorithms with the updated model specification, we use a slightly different notation for the LBA than previous literature has used. Usually, the LBA model represents a choice between C alternatives (C = 2, 3, ...) using C different evidence accumulators, one for each response choice. Each accumulator begins with an independent amount of starting evidence k c which is sampled independently for each accumulator from a continuous uniform distribution k c ∼ U (0, A). The evidence for accumulator c increases at a drift rate d c which is sampled independently for each accumulator from a normal distribution with mean v c and standard deviation s, so d c ∼ N (v c , s), although other non-normal distributions are possible (Terry et al., 2015). To satisfy the scaling conditions of the model, it is common to set the variance of the sampled drift rates to one, i.e., s = 1; however, see also Donkin et al. (2009). Each accumulator gathers evidence until one accumulator reaches a response threshold b. The LBA model assumes that the observed response time RT is the sum of the decision time, plus some extra time τ for the non-decision process such as stimulus encoding and motor execution. For simplicity, τ is usually assumed to be constant across trials. Thus, the final observed Let T c = (b − k c )/d c + τ be the time for accumulator c to reach the threshold b. Brown and Heathcote (2008) derive the cumulative distribution function of T c as and its density as above, φ and Φ are the density and cumulative distribution functions of the standard normal distribution, respectively, and It follows that the joint density of the response choice RE and response time RT at (RE, RT ) = (c, t), given the values b, A, v, s, τ, is defined as (1) Appendix E gives technical details of the joint density of the response choice RE and response time RT . The supplement at osf.io/5b4w3 shows that the LBA density Eq. (1) is bounded. If a subject makes N independent decisions, with choices RE i , i = 1 . . . N, and corresponding response times RT i , the density of RE := (RE 1 , . . . , RE N ) and RT := (RT 1 , . . . , RT N ) is In applications of the LBA model, it is possible that the parameters (b, A, etc.) vary over different conditions of the experiment, and sometimes also across the different accumulators.
The setup for the hierarchical LBA model is motivated by data collected from the decisions of 19 young subjects and first presented by Forstmann et al. (2008). The participants were asked to decide, repeatedly, whether a cloud of semi-randomly moving dots appeared to move to the left or to the right. Before each decision trial, subjects were instructed about what quality of their decisionmaking they should emphasise. For some trials, they were asked to respond as accurately as possible, for other trials they were asked to respond at their own pace, and for other trials they were asked to respond as quickly as possible. We label these conditions, in order: "accuracy emphasis" (condition 1); "neutral emphasis" (condition 2); and "speed emphasis" (condition 3). The different conditions were randomly mixed from trial to trial, with the subjects cued by a word which appeared on screen before each decision stimulus. Let Z be the number of conditions in the experiment, with Z = 3 here. Each subject made N z = 280 decisions for z = 1, ..., Z (N = 840 trials in total). See Forstmann et al. (2008) for more details on the procedure and the data, including the associated neuroimaging measurements, which are not considered here.
To model the differences between the three conditions in the experiment, we follow Forstmann et al. (2008) and allow different threshold parameters b (1) , b (2) and b (3) for the accuracy, neutral and speed conditions, respectively. RE i,j and RT i,j denote the ith response from the jth subject. Following Forstmann et al., we also collapse data across right-moving and left-moving stimuli, and so we index means of the drift rate distributions as v (1) and v (2) for the accumulator corresponding to incorrect and correct response choice, respectively. We assume that the standard deviation of the drift rate distribution is always s = 1. Together, these assumptions imply that each subject j has the vector of individual-level parameters or random effects With the usual assumptions of independence, the conditional density of all the observations is (2) Each of the individual fixed effects is restricted to be positive (it is theoretically possible, but psychologically implausible, for the means of the drift rate distributions, v j , to be negative). Respecting this, Turner et al. (2013) specified uncorrelated truncated normal distributions for each of these individual level parameters, and this has become standard in hierarchical applications of the LBA since then. Turner et al. (2013) and  also found that the posterior distributions of the individual random effects are highly correlated. Despite this, it has been standard practice to model the random effects as being a priori independent.
To improve both the computational efficiency of the algorithms and inference accuracy, we use instead a hierarchical model based on a multivariate nor-mal distribution of log-transformed random effects, with an explicitly-estimated variance; implying that the group distribution on the non-transformed scale is multivariate log-normal. For each subject j = 1, . . . , S, we define the vector of random effects Throughout, we use short-hand such as α 1,j to refer to the first element of the log-transformed parameter vector for participant j, and α j to refer to the corresponding random effects vector for participant j. The dimension of α j is D α = 7 here. Not all users may want to constrain the mean drift rates to be strictly positive, and so of course those may be left out of the log transformation in some cases.
To account for the dependence between the random effects, the prior distribution of the vector α j is modelled as There are a number of priors in the literature available for the parameters µ and Σ. We take the normal N (0, I Dα ) prior for µ, and the marginally noninformative prior of Huang and Wand (2013) for Σ, where v α , A 1 ,...,A Dα are positive scalars and diag (1/a 1 , ..., 1/a Dα ) is a diagonal matrix with diagonal elements 1/a 1 , ..., 1/a Dα . The notation IW (a, A) means an inverse Wishart distribution with degrees of freedom a and scale matrix A and the notation IG(a, b) means an inverse Gamma distribution with scale parameter a and shape parameter b. We choose this prior for Σ because it leads to psychologically plausible marginal prior distributions for the elements of the covariance matrix. Huang and Wand show that Eq. (5) induces half-t (v α , A d ) distributions for each standard deviation term in Σ and setting v α = 2 leads to marginally uniform distributions for all the correlation terms in Σ. In our application, we set v α = 2 and A d = 1 for all d = 1, . . . , D α . These prior densities cover most possible values in practice, and are relatively non-informative. The specification we have used implies that the distribution for the random effects vector exp(α j )|µ, Σ is a multivariate log-normal distribution with mean and covariance matrix given bỹ

Bayesian Estimation
This section develops efficient Bayesian inference for the hierarchical LBA model described in Section 2. We use the particle MCMC approach of Gunawan et al. (2017) and also develop a density tempered SMC approach which relies on particle MCMC.
Let θ ∈ Θ ⊂ R d θ be the vector of unknown group-level parameters, d θ the dimension of the parameters, and p(θ) be the prior for θ, where R m is m dimensional Euclidean space for a positive integer m. Let y j be the vector of observations for the jth subject, and define y = y 1:S = (y 1 , ..., y S ) as the vector of observations for all S subjects. Let α j ∈ χ ⊂ R Dα be the vector of individual-level parameters (random effects) for subject j, and p(α j |θ) its density under the group-level distribution. Θ and χ are themselves Euclidean spaces in all cases we consider. Now define α = α 1:S = (α 1 , ..., α S ) as the vector of all individual-level parameters, whose dimension is D α .

Particle Markov chain Monte Carlo (PMCMC)
The particle Metropolis with Gibbs (PMwG) sampler of Gunawan et al. (2017) is used for MCMC sampling. This sampler defines a target distribution on an augmented space that includes the model parameters and multiple copies of the individual random effects ("particles").
Let m j α j |θ, y j ; j = 1, ..., S be a family of proposal densities that is used to approximate the conditional densities p α j |θ, y j ; j = 1, ..., S . Appendix D gives the technical assumptions required for these proposal densities.
Let α r j be the rth sample from the proposal density m j α j |θ, y j for subject j. Define α 1:R 1:S := α 1:R 1 , ..., α 1:R S and α 1:R j := α 1 j , ..., α R j . Then the joint density of the particles α 1:R 1:S based on these proposals, and conditional on θ and y, is Let k = (k 1 , ..., k S ), with each k j ∈ {1, ..., R}, α k 1:S = α k1 1 , ..., α k S S be a vector of all selected individual random effects, and α is a collection of all particles excluding the selected individual random effects with α The augmented target density is defined as To understand the role of the index vector k = (k 1 , ..., k S ), we note that the density p θ, α k 1:S y) in Eq. (12) indicates that of all the random effect replicates α 1:R 1:S it is α k 1:S = (α k1 1 , ..., α k S S ) that is generated from the posterior p (θ, α 1:S | y) (after the sampler has converged to the target distribution). Gunawan et al. (2017) show more formally that the marginal density of θ and α 1:S with respect to the joint density p R θ, α 1:R 1:S , k y) is p (θ, α 1:S | y) and give convergence results for the PMwG sampler.

The Conditional Monte Carlo (MC) Algorithm
The conditional MC algorithm outlined in Algorithm 1 is an important component of the PMwG sampler and updates R − 1 particles simultaneously, while keeping the particle α k 1:S from the posterior p (θ, α 1:S | y) fixed in p R α 1:R 1:S |θ, y . Hence, the density of all the particles that are generated by the MC algorithm conditional on α k 1:S , θ, y is, ψ α 1:R 1:S y, θ)/ S j=1 m j α kj j |θ, y j which appears in the augmented target density in Eq. (12).
The following simple example illustrates the notation and the target distribution in Eq. (12). Suppose that there are S = 2 individuals and R = 3 particles for each individual. Let α 1 = α 1 1 , α 2 1 , α 3 1 and α 2 = α 1 2 , α 2 2 , α 3 2 be vectors of particles for subjects 1 and 2, respectively. If k 1 = 2 and k 2 = 3, then we define α k 1:2 = α 2 1 , α 3 2 as the vector of selected individual random effects for subjects 1 and 2, α is the collection of all particles excluding the selected random effects for subject 1, and α (−k2) 2 = α 1 2 , α 2 2 is the collection of all particles excluding the selected random effects for subject 2. Then the particles α k 1:2 are from the posterior p (θ, α 1:2 | y) and the rest of the particles α (−k) 1:2 from the proposal distributions m j α j |θ, y j for j = 1, 2. The index k indicates the α generated from the posterior p(θ, α 1:S |y) and is an important element of the PMwG sampler described by Algorithm 2.
, and nor- Note that step 2 is easily parallelized across the r = 1, ..., R particles and across the j = 1, ..., S subjects, or even both. There is no dependence between these parallel computations within each step.
Algorithm 2 describes the PMwG sampling scheme for the hierarchical LBA model defined in Section 2. The sampler starts at an initial set of parameters θ = (µ, Σ) and random effects α 1:S . We now explain one iteration of the PMwG algorithm. Steps (2a)-(2c) of the algorithm sample the group-level parameters of the LBA model using Gibbs steps conditional on the selected particles α k 1:S from previous iteration.
The PMwG sampler is applied in three stages to improve its effectiveness: burnin, adaptation, and sampling stages. The burnin stage allows the Markov chain to move from its initial position, which is randomly drawn from the prior, to the typical set of the posterior, i.e. the region of greatest posterior concentration (Betancourt, 2018). The adaptation stage draws samples from a reasonable approximation to the posterior distribution and uses those samples to construct improved proposal distributions for the sampling stage. These adapted proposal distributions are further adapted in the sampling stage and allow for very efficient sampling. Appendix A discusses the practical implementation of the PMwG sampler.
We take p 0 (θ, α) := p (θ) p (α|θ), because in the current LBA model, it is both easy to generate from and evaluate the prior densities p (θ) and p (α|θ). With this choice, At the initial temperature, the particle cloud θ at iteration p − 1 is an estimate of ξ ap−1 (θ, α|y). The transition from the particle cloud estimate of ξ ap−1 (θ, α|y) to the particle cloud estimate of ξ ap (θ, α|y) is implemented by first reweighting to obtain the updated weights W (p) (15) We now follow Del Moral et al. (2012) and select the next value of a p to target a pre-defined effective sample size, ESS T 2 . We do so by evaluating the ESS over a grid of points a 1:G,p of potential a p values and select as a p the value of a j,p whose ESS is the closest to ESS T . After reweighting, the effective sample size (ESS) is close to ESS T . To eliminate particles with low weight and replicate particles with larger weights, θ To improve the approximation of the particle cloud to ξ ap , we carry out L Markov move steps for each particle, using a Markov kernel K ξa p that has ξ ap as its invariant density. This Markov move step increases particle diversity, and in particular makes identical particles (which are produced during re-sampling) different from each other.
The Markov kernel K ξa p is constructed based on the PMwG algorithm. The augmented tempered target density at SMC step p is defined as ξ ap θ, α 1:R 1:S , k|y := where Eq. (11) gives ψ α 1:R 1:S |θ, y ; we note that it is possible to use different proposal densities m (p) j (·) for different SMC steps p. Using the same derivation as in Gunawan et al. (2017), we can show that the marginal density of θ and α 1:S w.r.t. ξ ap θ, α 1:R 1:S , k|y is ξ ap (θ, α 1:S |y). The augmented tempered density involves the term ψ α 1:R 1:S |θ, y / S j=1 m j α kj j |θ, y j , which is the density under ξ ap of all particles that are generated by the MC algorithm conditional on α k 1:S , θ, y . The conditional MC algorithm is similar to the one given in Algorithm 1, except that the density tempered SMC version adopts the tempered conditional density p (y|θ, α) ap instead of p (y|θ, α). Therefore, the Markov move step is based on the PMwG sampling scheme in Algorithm 2, except that instead of p R , we have augmented tempered target densities ξ ap .
Algorithm 3 describes the density tempered SMC algorithm. Steps (1), (2a)-(2d) are standard and apply to any model with slight modification.
Step (2e) 2 ESS measures variability in the weights, and is defined as and varies between 1 and M . A low value of ESS indicates that the weights are concentrated on only a few particles.
performs M parallel PMwG algorithm (Algorithm 2) L times for each temperature, except that in Step 2 of that algorithm we sample α (−k) 1:S ∼ ξ ap · |k, α k 1:S , θ, y using the conditional MC procedure in Algorithm 1 with the likelihood p (y|θ, α) replaced by p (y|θ, α) ap . This is an attractive and important feature as it allows the algorithm to fully use the computational power of modern graphical processing units (GPUs) with thousands of parallel cores. Appendix B discusses the tuning parameters and the proposal densities in the density tempered SMC algorithm.
Algorithm 3 The density tempered SMC algorithm 1. Set p = 0 and generate θ

Estimating the Marginal Likelihood
This section shows show how to estimate the marginal likelihood with negligible post-processing cost from the intermediate outputs of density tempered SMC using both the standard method (Del Moral et al., 2006) as well as thermodynamic integration (TI).

Standard density tempered SMC estimation of the marginal likelihood
The marginal likelihood since Z a0 = 1 and p (y) = Z a P . The particle cloud θ

Thermodynamic Integration Estimators from density tempered SMC Output
Gelman and Meng (1998) and Friel and Pettitt (2008) show how to compute the marginal likelihood of the data given the model using ideas from thermodynamic integration or path sampling, which relies on sampling from the posterior at different temperatures (called "power posteriors" or "tempered posteriors").
The thermodynamic identity (Friel and Pettitt, 2008;Friel et al., 2014) is The log of the marginal likelihood is the integral of E ξa p (log {p (y|θ, α)}) over the tempering sequence a p , where a p moves from 0 to 1. The derivation of the thermodynamic identity is in Friel and Pettitt (2008) and Friel et al. (2014). For each value of the tempering sequence a p , a sample from ξ ap (θ, α|y) can be used to estimate E ξa p (log {p (y|θ, α)}). Appendix C describes first and second order quadrature approximations to the integral in Eq. (17), labelled in the results below as T I 1 and T I 2 , respectively. We present two approximations as T I 2 corrects for bias that can be introduced in the method of estimating T I 1 .
There are some advantages in estimating the marginal likelihood by thermodynamic integration using density tempered SMC output: (i) density tempered SMC produces the M triples θ 1:M for each value of the tempering sequence a p , p = 0, ..., P ; this means that estimating the marginal likelihood via thermodynamic integration does not incur any extra computational costmarginal likelihood can be calculated by TI directly from the density tempered SMC outputs; (ii) assessing the convergence of density tempered SMC is much less of an issue than that of MCMC methods; and density tempered SMC is easily parallelizable for each a p in the tempering sequence; (iii) the number of tempering steps and the tempering sequence are chosen adaptively to target a pre-defined effective sample size of the SMC samples.

Estimating Marginal Likelihood using the PMwG output
Appendix G discusses the estimation of marginal likelihood using the PMwG output.

Application to Simulated Data
This section applies the PMwG and density tempered SMC methods to fit the hierarchical LBA model specified in Section 2; the data simulated from an LBA model mimicking the conditions inspired by the experiment of Forstmann et al. (2008). The three generated datasets have S = 19, 50, and 100 subjects, and N = 1000 trials (N z ≈ 333 trials in each condition) for each subject for each value of S; the generated data have nonzero correlations between the random effects which reflects both the improved parametric specification of the LBA model as well as plausible psychological assumptions about individual differences. To generate the simulated data, we used group-level parameters which matched those estimated for the real data reported by Forstmann et al. (2008). This defined a multivariate normal distribution for the log-transformed random effects. The random effects for the simulated data sets were sampled randomly from this multivariate normal.
A total of 10, 000 draws were obtained using the PMwG sampler for the subsequent analysis of the posterior distribution. See Appendix A for a discussion on the tuning parameters, the proposal densities, and the implementation of the PMwG sampler. We assess the performance of the sampler using the "integrated autocorrelated time" (IACT), which measures the inefficiency of the sampling scheme in terms of the multiple of its iterates that are required to obtain the same variance as an independent sampling scheme, e.g. if IACT = 10, then we need 10 times as many iterates as an independent scheme, so a larger value indicates poorer performance. The IACT of a scalar parameter ω is defined as (Chib and Greenberg, 1996) IACT where ρ ω (t) is the lag-t autocorrelation of the iterates of ω for the underlying stationary Markov chain. Table 1 reports the estimated IACT values for the group level parameters from the PMwG sampler, with the estimates obtained using the CODA R package of Plummer et al. (2006) 3 . The results suggest that the PMwG sampler is efficient because all its IACT values are small. Similar conclusions can be drawn from the IACT values for the random effects, as shown in the online supplement at osf.io/5b4w3. Table 1: Estimated Inefficiency Factors (IACT) of the LBA parameters using the PMwG method with S = 19, 50, 100 subjects and N = 1, 000 in the simulation study. The order of the random effect parameters in the covariance matrix Σ is We used 10 independent runs with M = 250 samples each for the density tempered SMC method to generate 2, 500 samples of the LBA individual random effects and parameters. The independent samplers mean that 10 independent estimates of the marginal likelihood are also obtained, which allows a rough estimate of the sampling variability in the marginal likelihood. Sampling error in the marginal likelihood is important in inference, but often overlooked. See Appendix B for a discussion of the tuning parameters, the proposal densities, and the implementation of the density tempered SMC sampler. Table 2 shows the wall-clock computation time to run both the PMwG and the density tempered SMC methods which are based on a Matlab implementation running on 28 CPU-cores. The running time for PMwG includes the time taken for all three stages. The table shows that PMwG is much faster than density tempered SMC. In general, the PMwG method can be used with minimal computational resources, e.g. a personal computer; however, density tempered SMC is easier to parallelize than PMwG, so it is likely to be faster than PMwG if there is access to many more CPU cores and the model or data are large. The optimal number of CPU cores required for density tempered SMC is equal to the number of SMC samples M , which means the properties of the sampler can be easily tuned to provide maximum parallel efficiency on a large range of hardware.  Figure 1 plots the posterior distributions estimated by PMwG (in blue) and density tempered SMC (in red) for the population mean parameters (µ), with the vertical lines showing the true parameter values. Figure 2 shows similar plots for the variances (the diagonal elements of Σ) of the random effects. In both figures, the top, middle and bottom panels correspond to S = 19, S = 50 S = 100 subjects, respectively.  The two figures show that the PMwG and density tempered SMC estimators give very similar results-suggesting that the two algorithms are correctly implemented. The posterior densities narrow as the number of subjects S increases, and the posterior modes of the parameters are consistent with the true data generating values. Each panel in the two figures also includes a third posterior distribution (in yellow), which illustrates the effect of making the standard hierarchical LBA assumption that the random effects are independent, using the same priors for the group level parameters as in the correlated case discussed above, and using PMwG to sample.
The results for this simplified, uncorrelated, LBA model differ systematically to the LBA model proposed here. The posterior mean estimates are more certain, that is, the distributions are more peaked (Figure 1), and the posterior variance estimates are smaller (Figure 2). This suggests that using a model with independent distributions for the random effects, as is standard in the literature, can -if there really is correlation in the data -lead to unwarranted overconfidence in estimation precision, and under-estimation of the magnitude of individual differences.
We estimated the marginal likelihood for each simulated experiment (S = 19, S = 50, and S = 100 subjects) using the outputs of the density tempered SMC sampler and the three estimators: the standard density tempered SMC method, and the thermodynamic integration estimators T I 1 and T I 2 discussed in Section 3.3. Table 3 reports the logs of the marginal likelihood estimates (with the standard error in brackets); the standard errors are obtained using ten replicates for each estimator. The three methods agree very closely, and the standard errors for each method are very small, suggesting that the log of the marginal likelihoods are estimated accurately.  Figure 3 shows the kernel density estimates of the marginal posterior densities for the three parameters µ 1 , µ 2 and µ 3 which govern the decision threshold in the accuracy, neutral, and speed emphasis conditions (respectively) of the simulated experiment. There is considerable overlap between the two marginal distributions representing the different threshold parameters for the accuracy and neutral conditions. If these data were from a real (not simulated) experiment, this might be interpreted as evidence that the participants in the experiment failed to distinguish between those two conditions; that they did (or could) not adopt different decision-making styles when asked to. This is exactly the kind of model selection question facing researchers using the LBA model. Section 4.2 demonstrates how to use the marginal likelihood estimates obtained using density tempered SMC to further investigate whether participants adopted different decision-making thresholds in the three conditions. In that section, we compare the unrestricted (three parameter) model against restricted models having shared threshold parameters. An alternative approach to this question could enforce order constraints on the estimated random effects. This is easily accomplished by parameter transformations -define the speed-emphasis threshold as usual, and then estimate parameters for the differences between that threshold and the neutral condition threshold, and between the neutral and accuracy condition thresholds. These incremental parameters can be constrained to be positive by estimating their logarithms, as with the other parameters, and inferential tests made by comparing the estimated increments with zero.

Comparing the PMwG sampler to existing estimation methods
The hierarchical LBA model is most often estimated using DE-MCMC as developed by Turner et al. (2013). However, such applications have been restricted to an LBA model which assumes independent truncated normal prior distributions for the random effects, whereas the LBA model defined in Section 2 allows for a (correlated) multivariate normal prior distribution on the log of the random effects. The considerably simpler model estimated via DE-MCMC does not require the estimation of covariance elements. An uncorrelated model with D α random effects per person has 2D α group-level parameters (a mean and a standard deviation for each random effect). In contrast, the current LBA model has D α +(D α +1)D α /2 parameters (a mean and standard deviation for each random effect, plus all pairwise covariances). Appendix F shows that the PMwG sampler agrees with the DE-MCMC sampler, when estimating the simplified LBA model without correlations.
Estimating the LBA model with non-zero correlations in the prior may be impossible with the current DE-MCMC sampler and will require extensive modifications at a minimum. The DE-MCMC sampler generates the group level mean and standard deviation for each of the random effects blocked in pairs. This blocking breaks down when the off-diagonal covariances in the prior are non-zero, as they are highly correlated with the (diagonal) variance parameters, and so some other blocking scheme will need to be devised. It is not clear what that scheme should be. Further, DE-MCMC generates proposals by taking linear combinations of existing samples. When the correlations are non-zero, this means generating covariance matrix proposals, but the differential evolution approach does not ensure that the proposals will be valid covariance matrices (positive definite). It is also standard to set the tuning parameter in the DE-MCMC algorithm inversely proportional to the number of parameters. This is likely to make the method very inefficient given the much larger number of unknown parameters when a full covariance matrix is used in the prior (scales with D 2 α ). It is possible that these difficulties for DE-MCMC may be overcome with its further development. However, such an extension of the DE-MCMC approach is beyond the scope of this paper.
In contrast, PMwG scales well as the number of parameters increases because it uses Gibbs sampling to generate the group-level parameters. In demonstrating this scaling property, Wall et al. (2019) extend the PMwG approach to estimate a high-dimensional hierarchical LBA model with S = 110 subjects, N = 1, 350 trials, and D α = 30 individual random effects parameters-giving a covariance matrix with 435 free parameters. Gunawan et al. (2019) extend the hierarchical LBA model to allow the individual level parameters of each subject to change over blocks of the trials, thus extending the PMwG approach to estimate timevarying LBA models. Both these extensions show that the methods proposed in this article allow exploration of important psychological questions that are at present neglected due to computational intractability.

Real Data Application
We applied PMwG and density tempered SMC to the behavioural data first presented by Forstmann et al. (2008), and introduced in Section 2. PMwG and density tempered SMC were run with the same settings as in Appendix A and Appendix B, respectively. The wall-clock computation times to run PMwG and density tempered SMC were around 30 and 138 minutes, respectively, using a Matlab implementation of the algorithm and 28 CPU-cores. Table 4 shows the estimated group-level parameters from both methods, along with the estimated posterior standard deviations and the IACT inefficiency factors for PMwG. All the IACT values are again small, indicating that the chains mixed well and that the performance of the sampler did not deteriorate markedly when moving from simulated data to real data. The posterior mean estimates and the associated posterior standard deviations from PMwG and density tempered SMC are very close to each other, providing an accuracy check on both samplers. Figure 4 shows the kernel density estimates of marginal posterior densities for the group mean parameters, using samples from PMwG (blue) and density tempered SMC (red); the figure also shows the estimated posterior distributions for the reduced model, with zero prior correlation between the individual level parameters estimated using PMwG (yellow). Figure 5 shows the corresponding results for the between-subject variance parameters instead of the means, i.e., the diagonal elements of Σ. Both figures demonstrate that the density tempered SMC estimates are very close to the PMwG estimates, for all the parameters for the LBA model with a full prior covariance matrix for the random effects. However, similarly to the simulation study, the posterior estimates from the simplified LBA model (with zero covariance) are quite different to the posterior estimates from the LBA model with a full covariance matrix. The differences are particularly pronounced for the variance parameters in Figure 5. In all cases, Table 4: Posterior means (with posterior standard deviations in brackets) and IACT of the LBA group-level parameters of the full model estimated using PMwG and density tempered SMC (DT-SMC) for the data from Forstmann et al. (2008). The order of the random effect parameters in the covariance matrix the zero-covariance reduced model estimates are much smaller for the betweensubjects variances, and often also have narrower posterior distributions. This matches the results from the simulation study above, and is consistent with the hypothesis that the real data includes non-zero between-subject correlations, and failing to take this into account in the prior leads to unwarranted overconfidence in the posteriors. The model setup in Section 2 allows us to estimate the correlation matrix Γ between individual level parameters by using the standard transformation of the estimated covariance matrix. Table 5 shows that the threshold parameters for all three conditions are highly correlated: Γ b (1) , b (2) = .  parameters -presumably because of trade-offs between explaining the same RT as either composed of more or less decision time vs. non-decision time. The mean drift rate for the accumulator corresponding to the correct response is not highly correlated with other individual level parameters. The magnitudes of these correlations emphasize the importance of explicitly modelling the covariance matrix, rather than forcing it to have zero correlations. Table 6 summarizes the posterior mean estimates of the parameters on the original (not logarithmic) scale obtained using Eq. (6). These values suggest that there may only be slight differences in the estimates of the threshold parameters between the three conditions: accuracy, neutral, and speed (μ 1 ,μ 2 , andμ 3 ). To investigate this, we estimated a restricted model with two threshold parameters by combining the accuracy and neutral conditions, as well as a more restricted model with a single, shared, threshold parameter for all three conditions. We used density tempered SMC, specified as above, to estimate the marginal likelihood for each model. Table 7 reports the estimated log marginal likelihoods (with standard errors in brackets) for the three models. The differences between the log marginal likelihoods are much larger than the standard errors, and also large relative to the scales usually used to judge statistical re- liability, e.g. the corresponding Bayes factors for the model comparisons are all much larger than 10 6 . The results favour the unrestricted model for these data, which also supports the analyses by Forstmann et al. (2008). An important caveat for this analysis is that the marginal likelihoods -for any modeldepend on the prior. For example, another approach to the problem would set up the random effects using an intercept plus effect coding, and specify a prior distribution for the effect sizes. Further work will be required to investigate the sensitivity of model selection outcomes to different choices of prior distribution.

Conclusions
Based on recent advances in particle MCMC, the article develops two new estimation approaches for the Linear Ballistic Accumulator model of Brown and Heathcote (2008); Particle Metropolis within Gibbs and density tempered SMC. We show that PMwG and density tempered SMC perform well for both simulated and real data. The new methods are alternatives to the existing approach that is based on MCMC with proposals generated by differential evolution (Turner et al., 2013) and provide important advantages. Density tempered SMC is extremely well-suited to parallelisation on high-performance computers, which is likely to be an advantage in future work with large-scale models and data. Although both PMwG and DE-MCMC are also both appropriate for parallelisation in high-performance computing environments, DE-MCMC requires more frequent dependence between the multiple chains, which limits its efficiency when parallelised. We also found that the PMwG sampler is much faster than density tempered SMC and can be easily implemented with minimal computational resources, e.g., personal computers. Furthermore, density tempered SMC provides an estimate of the marginal likelihood at little extra cost, and hence it can be readily used for model selection via Bayes factors.
Another important contribution of our work is to explicitly model the full covariance structure of the prior for the random effects. Like all plausible cognitive models, there are substantial correlations between the individual level parameters of the LBA model: subjects with a large decision threshold also tend to have a large starting point distribution, etc. In previous applications of the model, these prior correlations were set to zero, with the group-level distributions treated a priori as independent. Despite making this a priori assumption of independence, the resulting posterior samples always exhibited strong posterior correlations between individual level parameters. Explicitly allowing non-zero correlations in the prior, as we have done, provides better estimates of the parameters and their variances, and improves computational efficiency.
The computational flexibility of the new methods allows exploration of important psychological questions which have hitherto been neglected, due to statistical intractability. For example, it is well known that there can be substantial sequential effects in decision-making data: both response choices and response times tend to be positively autocorrelated. All applications of the LBA modeland indeed, almost all decision-making models -have ignored these sequential effects, treating the data as i.i.d. and attributing the effects of any sequential dependence to error terms. Both our approaches allow tractable extensions that explicitly take into account of within subject dependence and other interesting sequential effects such as parameter evolution due to fatigue or learning. We are investigating these models in ongoing work.
To aid researchers in adopting the proposed methods, we provide scripts that implement both the PMwG and density tempered SMC methods as applied to the real data from Forstmann et al. (2008); see osf.io/5b4w3 for more details.
choosing efficient proposal densities. To simplify this choice, we develop efficient proposal densities in three stages: burnin, initial adaptation, and sampling. In the burnin and the initial adaptation stages, the proposal density for subject j is a mixture over the prior group-level distribution for the random effects, p (α j |θ), and a normal distribution centred on the previous sample for the random effect N α j ; α is the previous iterate α kj j for the individual jth random effect. In practice, we should use a larger number of particles in the burnin and initial adaptation stages than in the sampling stages.
This proposal can be made more flexible by using Σ instead of Σ, where 0 < < 1 is a scale factor. In this paper, we set = 1. It is necessary to reduce the scale factor when we have larger number of random effects. Wall et al. (2019) estimate a high-dimensional hierarchical LBA model with S = 110 subjects, N = 1, 350 trials, and D α = 30 individual random effects parameters and set = 0.1. We find empirically in our examples that this strategy works well in practice.
In the sampling stage, we use the posterior MCMC draws (α 1:S , θ) from the initial adaptation stage to adaptively build more efficient proposal densities m j α j |θ, y j , for each subject j = 1, ..., S. This usually allows the use of a much smaller number of particles. We first transform the posterior draws of the parameters Σ so that they all lie on the real line. The covariance matrix Σ is reparameterised in terms of its Cholesky factorisation Σ = LL T , where L is a lower triangular matrix. We also apply a log transformation for the diagonal elements of L, while the subdiagonal elements of L are unrestricted. For each subject, we fit a normal distribution to the vectors formed by joining the posterior draws of α j with (µ, L) and obtain the conditional distribution g (α j |µ, L) ∼ N α j ; µ j,prop , Σ j,prop for j = 1, ..., S. The efficient proposal density for subject j is then the two component mixture In this paper, the number of particles in the PMwG method is set to R burnin = R adapt = R sampling = 100; as discussed above, it is sometimes helpful to use more particles in the burnin and initial adaptation stages. We also set w mix = 0.9.
We used 500 burnin iterations, 500 adaptation iterations, and 10000 iterations for sampling. Following Hesterberg (1995), including the prior density p (α j |θ) in Eq. (A.1) and Eq. (A.2) ensures that the importance weights are bounded because it is straightforward to show that the density p y j |θ, α j is bounded. This ensures that the sampler is ergodic; see the online supplement at osf.io/5b4w3 for further details.

Appendix B. Tuning Parameters and Proposal Densities for Density tempered SMC
Density tempered SMC has three tuning parameters: the number of particles R, the number of Markov move steps L, and the number of SMC samples M . The bigger the number of SMC samples M , the better the approximation to the posterior density p (θ, α|y). By using the result of Del Moral et al. (2006), density tempered SMC provides consistent inference for the posterior density p (θ, α|y) as the number of annealed samples M goes to infinity, for any given number of particles R. The L Markov moves in step (2e) in Algorithm 3 help to diversify the collection of parameters and random effects after the resampling in step (2d) so that they better approximate the tempered target density. In addition, we can obtain more accurate marginal likelihood estimates with larger number of SMC samples M , the number of particles R, and the number of Markov move steps L.
The prior density p (α j |θ) is an efficient proposal for the random effects for each subject when the tempering value a p is small as it dominates the tempered density which is quite flat because p y j |θ, α j ap ≈ 1. In our application we use the prior as a proposal density when a p < 0.1. It is often also adequate to use a smaller number of particles R and smaller number of Markov moves L when a p is small. When a p is larger than 0.1, we first fit a normal distribution -in the same manner as for the adaptive proposal densities in the PMwG algorithm -to the current transformed particle cloud (µ, L) for j = 1, ..., S. This gives the conditional distribution g (α j |µ, L) ∼ N α j ; µ j,prop , Σ j,prop at each stage of the SMC process. We then use the two component mixture given in Eq. (A.2) as a proposal density and set the mixture weight to w mix = 0.9. The number of MC samples and Markov moves were set to R = 100 and L = 10 respectively. We set ESS T = 0.8M ; i.e., we target an effective sample size of 80% of the maximum SMC sample size. Unlike the PMwG algorithm, density tempered SMC does not require an initial adaptation stage to construct the efficient proposal density. Instead, the proposal densities of the random effects for each subject are obtained from the current particle cloud at each stage of the SMC process.

Appendix C. Details on the Thermodynamic Integration
The first order quadrature approximation (T I 1 ) to the integral in Eq. (17) of the main text is based on the trapezoidal rule and is Discretising the tempering sequence a p and using the trapezoidal rule introduces some bias to the estimate in Eq. (C.1). Friel et al. (2014) propose using the corrected trapezium rule method of Atkinson and Han (2004) to reduce the bias of T I 1 . The corrected trapezium rule approximates the integral of a function f between points a and b as where c ∈ [a, b] and f (·) and f (·) are the first and second derivative of f . The first term in Eq. (C.2) is the usual trapezium rule approximation used in T I 1 .
The second derivative f (c) in the second term of Eq. (C.2) can be approximated as Using the results in Eq. (C.4), the second order quadrature approximation (T I 2 ) to the integral in Eq. (17) of the main text is where both the expectation E ξa p (log {p (y|θ, α)}) and variance V ξa p (log {p (y|θ, α)}) can be estimated using the density tempered SMC output at the tempering value a p . See Appendix C for more details. Estimating the marginal likelihood using thermodynamic integration requires careful consideration of three tuning issues: (i) the algorithm that samples from the tempered posterior ξ ap (θ, α|y), for p = 0, ..., P ; (ii) the number of tempering steps (P ); and (iii) the tempering sequence a p for p = 0, ..., P . Friel and Pettitt (2008) use standard MCMC algorithms, such as Gibbs and Metropolis-Hastings, because they deal with tractable likelihoods; Evans and Annis (2019) use DE-MCMC to sample from the tempered target posterior. It is important to make sure that the MCMC sampler used converges for each value of the tempering sequence a p , p = 0, ..., P to ensure the accuracy of the marginal likelihood estimates. However, it is difficult in general to assess whether the chains mix adequately and converge to the invariant tempered target density. Furthermore, standard random walk Metropolis-Hastings and DE-MCMC algorithms usually suffer from high autocorrelations between samples, and slow or uncertain convergence for models with a large number of parameters. The MCMC convergence problems can sometimes be solved by increasing the number of MCMC samples, but this leads to increased computational workload and the MCMC algorithms suffer from limited parallelizability. Friel and Pettitt (2008), Xie et al. (2010), and Evans and Annis (2019) use the tempering sequence a p = ((p − 1)/(P − 1)) 1/0.3 , for p = 0, ..., P . These tempering sequences place more computational effort on temperatures near 0, where the tempered posterior changes rapidly. They also fix the tempering steps before running the algorithm by experimenting with different numbers of tempering steps and then choose the tempering step P that gives the lowest standard error of the log of marginal likelihood estimates. In contrast, our approach through density tempered SMC gives a principled, but perhaps not optimal, sequence of tempering steps.

Appendix D. Assumptions for the Proposal Densities
We define the support of the posterior and the proposal densities as S θ j := α j ∈ χ α : p α j |θ, y j > 0 and Q θ j := α j ∈ χ α : m j α j |θ, y j > 0 .
and assume that S θ j ⊆ Q θ j for any θ ∈ R d θ and j = 1, ..., S. This ensures that the m j α j |θ, y j can be used as proposal densities to approximate p α j |θ, y j .
Appendix E. Further details on the joint density of response times and response choice We derive in detail the joint density in Eq. (1) of the response choice RE and response time RT . Let µ 1 (·) be the counting measure on the set {2, ..., C}, i.e. µ 1 (c) = 1 for c = 2, . . . , C, and µ 2 (·) be the usual Lebesgue measure on the line. The joint probability for RE = c and RT ∈ (t, t + dt) is P RE = c, RT ∈ (t, t + dt) = P RE = c P RT ∈ (t, t + dt)|RE = c = P (T k > t, k = c)P T c ∈ (t, t + dt) = k =c (1 − F k (t))f c (t)µ 1 (c)µ 2 (dt).
This implies that Eq. (1) is the joint density of RE and RT with respect to the product measure µ 1 ⊗ µ 2 .

Appendix F. Comparing PMwG and DE-MCMC for the simplified LBA model
This section compares the performance of PMwG and DE-MCMC applied to a simplified hierarchical LBA model that assumes independent normal distributions for the log-transformed individual random effects using the same simulated dataset with S = 100 and N = 1000 trials given in Section 4.1. This is a standard assumption for DE-MCMC applications of LBA, and simplifies model estimation considerably because only the diagonal elements of Σ are estimated. For each random effect component in α j , we define the independent univariate normal distributions α dj ∼ N µ d , σ 2 d , d = 1, ..., D α . We use the same priors for the group level parameters defined in Section 2.
For DE-MCMC we used 15 chains and ran each for 5, 000 iterations, discarding the first 2, 500 iterations from each chain, and thinning by keeping only every 10th draws; we obtain a total of 3, 750 samples. To match with DE-MCMC sampler, we generated 4, 000 draws using the PMwG sampler for the analysis of the posterior distribution. The posterior distributions from the two samplers agreed closely, which is to be expected with the large (100 participants, each with 1, 000 trials) and clean (synthetic data generated without noise) sample. Figure F.6 shows the trace plots of the iterates of two group-level LBA parameters estimated using PMwG and DE-MCMC, respectively. It is clear that even with substantial thinning, the DE-MCMC samples do not mix as well as the (un-thinned) PMwG samples. Similar plots were obtained for the other group level parameters. Since the DE-MCMC is based on multiple interacting chains, the DE-MCMC samples is plotted by taking the mean across iterates from each chain. The impression of poorer mixing for the DE-MCMC sampler was confirmed by calculating IACT values for all parameters. Those calculations showed higher inefficiency factors for the 10x thinned DE-MCMC samples than for the un-thinned PMwG samples (median IACT for parameters was 2.10, compared with 1.19 for the PMwG samples). We also calculated IACT on the chain formed by taking the mean across iterates from each chain. This likely provides an underestimate of inefficiency experienced in practice, where the function calculated from samples is often not summed over chains.
Appendix G. Estimating Marginal Likelihood using the PMwG output Gronau et al. (2019) and Evans and Annis (2019) use the posterior samples obtained from the DE-MCMC sampler to estimate the marginal likelihood by bridge sampling and thermodynamic integration, respectively. Section 4 and Appendix F show that the PMwG sampler is more reliable and efficient than the DE-MCMC sampler. Therefore, it is instructive to use the PMwG output to estimate the marginal likelihood by bridge sampling and thermodynamic integration.
For thermodynamic integration, we first run the PMwG sampler for each value of the tempering sequence a p . The posterior samples at each tempering sequence a p is the estimate of E pa p (θ,α|y) (log {p (y|θ, α)}) and these samples can be used to obtain the marginal likelihood using thermodynamic integration methods described in Appendix C.
We can similarly first obtain the posterior samples using the PMwG; and then use these posterior samples to obtain the marginal likelihood estimate by the bridge sampling method as described in Gronau et al. (2019).