Likelihood approximation networks (LANs) for fast inference of simulation models in cognitive neuroscience
ABC methods apply when one has access to a parametric stochastic simulator ( also referred to as generative model ) , but , unlike the usual setting for statistical inference , no access to an explicit mathematical formula for the likelihood of observations given the simulator’s parameter setting . While likelihood functions for a stochastic stimulator can in principle be mathematically derived , this can be exceedingly challenging even for some historically famous models such as the Ornstein–Uhlenbeck ( OU ) process ( Lipton and Kaushansky , 2018 ) and may be intractable in many others . Consequently , the statistical community has increasingly developed ABC tools that enable posterior inference of such ‘likelihood-free’ stochastic models while completely bypassing any likelihood derivations ( Cranmer et al . , 2020 ) . Given a parametric stochastic simulator model ℳ and dataset 𝐱 , instead of exact inference based on pℳ⁢ ( θ|𝐱 ) , these methods attempt to draw samples from an approximate posterior p~ℳ⁢ ( θ|𝐱 ) . Consider the following general equation:p~ℳ⁢ ( θ|𝐱 ) ∝∫Kh⁢ ( ||sℳ-s𝐱|| ) ⁢pℳ⁢ ( sℳ|θ ) ⁢𝑑sℳ⁢π⁢ ( θ ) ∝p~ℳ⁢ ( 𝐱|θ ) ⁢π⁢ ( θ ) where sℳ refers to sufficient statistics ( roughly , summary statistics of a dataset that retain sufficient information about the parameters of the generative process ) . Kh ( ||sℳ−sx|| ) refers to a kernel-based distance measure/cost function , which evaluates a probability density function for a given distance between the observed and expected summary statistics ||sℳ-s𝐱|| . The parameter h ( commonly known as bandwidth parameter ) modulates the cost gradient . Higher values of h lead to more graceful decreases in cost ( and therefore a worse approximation of the true posterior ) . By generating simulations , one can use such summary statistics to obtain approximate likelihood functions , denoted as p~ℳ⁢ ( 𝐱|θ ) , where approximation error can be mitigated by generating large numbers of simulations . The caveat is that the amount of simulation runs needed to achieve a desired degree of accuracy in the posterior can render such techniques computationally infeasible . With a focus on amortization , our goal is to leverage some of the insights and developments in the ABC community to develop NN architectures that can learn approximate likelihoods deployable for any inference scenario ( and indeed any inference method , including MCMC , variational inference , or even MLE ) without necessitating repeated training . We next describe our particular approach and return to a more detailed comparison to existing methods in the 'Discussion' section . In this section , we outline our approach to amortization of computational costs of large numbers of simulations required by traditional ABC methods . Amortization approaches incur a one-off ( potentially quite large ) simulation cost to enable cheap , repeated inference for any dataset . Recent research has led to substantial developments in this area ( Cranmer et al . , 2020 ) . The most straightforward approach is to simply simulate large amounts of data and compile a database of how model parameters are related to observed summary statistics of the data ( Mestdagh et al . , 2019 ) . This database can then be used during parameter estimation in empirical datasets using a combination of nearest-neighbor search and local interpolation methods . However , this approach suffers from the curse of dimensionality with respect to storage demands ( a problem that is magnified with increasing model parameters ) . Moreover , its reliance on summary statistics ( Sisson et al . , 2018 ) does not naturally facilitate flexible reuse across inference scenarios ( e . g . , hierarchical models , multiple conditions while fixing some parameters across conditions , etc . ) . To fulfill all desiderata outlined in the introduction , we focus on directly encapsulating the likelihood function over empirical observations of a simulation model so that likelihood evaluation is substantially cheaper than constructing ( empirical ) likelihoods via model simulation online during inference . Such empirical observations could be , for example , trial-wise choices and reaction times ( RTs ) . This strategy then allows for flexible reuse of such approximate likelihood functions ℓ^⁢ ( θ|𝐱 ) in a large variety of inference scenarios applicable to common experimental design paradigms . Specifically , we encapsulate ℓ^⁢ ( θ|𝐱 ) as a feed-forward NN , which allows for parallel evaluation by design . We refer to these networks as likelihood approximation networks ( LANs ) . Figure 1 spells out the setting ( A ) and usage ( B ) of such a method . The LAN architectures used in this paper are simple , small in size , and are made available for download for local usage . While this approach does not allow for instantaneous posterior inference , it does considerably reduce computation time ( by up to three orders of magnitude; see ‘Run time’ section ) when compared to approaches that demand simulations at inference . Notably , this approach also substantially speeds up inference even for models that are not entirely likelihood free but nevertheless require costly numerical methods to obtain likelihoods . Examples are the full-DDM with inter-trial variability in parameters ( for which likelihoods can be obtained via numerical integration , as is commonly done in software packages such as HDDM but with substantial cost ) , but also other numerical methods for generalized DDMs ( Shinn et al . , 2020; Drugowitsch , 2016 ) . At the same time , we maintain the high degree of flexibility with regards to deployment across arbitrary inference scenarios . As such , LANs can be treated as a highly flexible plug-in to existing inference algorithms and remain conceptually simple and lightweight . Before elaborating our LAN approach , we briefly situate it in the context of some related work . One branch of literature that interfaces ABC with deep learning attempts to amortize posterior inference directly in end-to-end NNs ( Radev et al . , 2020b; Radev et al . , 2020a; Papamakarios and Murray , 2016; Papamakarios et al . , 2019a; Gonçalves et al . , 2020 ) . Here , NN architectures are trained with a large number of simulated datasets to produce posterior distributions over parameters , and once trained , such networks can be applied to directly estimate parameters from new datasets without the need for further simulations . However , the goal to directly estimate posterior parameters from data requires the user to first train a NN for the very specific inference scenario in which it is applied empirically . Such approaches are not easily deployable if a user wants to test multiple inference scenarios ( e . g . , parameters may vary as a function of task condition or brain activity , or in hierarchical frameworks ) , and consequently , they do not achieve our second desideratum needed for a user-friendly toolbox , making inference cheap and simple . We return to discuss the relative merits and limitations of these and further alternative approaches in the 'Discussion' section . Formally , we use model simulations to learn a function fΦ⁢ ( x , θ ) , where fΦ ( . ) is the output of a NN with parameter vector Φ ( weights , biases ) . The function fΦ⁢ ( 𝐱 , θ ) is used as an approximation ℓ^⁢ ( θ|𝐱 ) of the likelihood function ℓ⁢ ( θ|𝐱 ) . Once learned , we can use such f as a plug-in to Bayes' rule:p^ ( θ|x ) ∝∏i=1Nℓ^ ( θ|xi ) p ( θ ) To perform posterior inference , we can now use a broad range of existing MC or MCMC algorithms . We note that , fulfilling our second desideratum , an extension to hierarchical inference is as simple as plugging in our NN into a probabilistic model of the formp^⁢ ( θ|𝐱 ) ∝∏j=1J∏i=1Nℓ^⁢ ( θj|xj⁢i ) ⁢p⁢ ( θj|α , β ) ⁢p⁢ ( α , β|γ ) where α , β are generic group-level global parameters , and γ serves as a generic fixed hyperparameter vector . We provide proofs of concepts for two types of LANs . While we use MLPs and convolutional neural networks ( CNNs ) , of conceptual importance is the distinction between the two problem representations they tackle , rather than the network architectures per se . The first problem representation , which we call the pointwise approach , considers the functions fΦ⁢ ( x|θ ) , where θ is the parameter vector of a given stochastic simulator model , and x is a single datapoint ( trial outcome ) . The pointwise approach is a mapping from the input dimension |Θ|+|x| , where | . | refers to the cardinality , to the one-dimensional output . The output is simply the log-likelihood of the single datapoint x given the parameter vector θ . As explained in the next section , for this mapping we chose simple MLPs . The second problem representation , which we will refer to as the histogram approach , instead aims to learn a function fΦ ( . |θ ) , which maps a parameter vector θ to the likelihood over the full ( discretized ) dataspace ( i . e . , the likelihood of the entire RT distributions at once ) . We represent the output space as an outcome histogram with dimensions n×m ( where in our applications n is the number of bins for a discretization of reaction times , and m refers to the number of distinct choices ) . Thus , our mapping has input dimension |Θ| and output dimension |m×n| . Representing the problem this way , we chose CNNs as the network architecture . In each of the above cases , we pursued the architectures that seemed to follow naturally , without any commitment to their optimality . The pointwise approach operates on low-dimensional inputs and outputs . With network evaluation speed being of primary importance to us , we chose a relatively shallow MLP to learn this mapping , given that it was expressive enough . However , when learning a mapping from a parameter vector θ to the likelihood over the full dataspace , as in the histogram approach , we map a low-dimensional data manifold to a much higher dimensional one . Using an MLP for this purpose would imply that the number of NN parameters needed to learn this function would be orders of magnitude larger than using a CNN . Not only would this mean that forward passes through the network would take longer , but also an increased propensity to overfit on an identical data budget . The next two sections will give some detail regarding the networks chosen for the pointwise and histogram LANs . Figure 1 illustrates the general idea of our approach while Figure 2 gives a conceptual overview of the exact training and inference procedure proposed . These details are thoroughly discussed in the last section of the paper . As a first approach , we use simple MLPs to learn the likelihood function of a given stochastic simulator . The network learns the mapping {θ , 𝐱}→log⁡ℓΦ⁢ ( θ|𝐱 ) , where θ represents the parameters of our stochastic simulator , 𝐱 are the simulator outcomes ( in our specific examples below , 𝐱 , refers to the tuple ( r⁢t , c ) of reaction times and choices ) , and log⁡ℓ⁢ ( θ|𝐱 ) is the log-likelihood of the respective outcome . The trained network then serves as the approximate likelihood function log⁡ℓ^Φ⁢ ( θ|𝐱 ) . We emphasize that the network is trained as a function approximator to provide us with a computationally cheap likelihood estimate , not to provide a surrogate simulator . Biological plausibility of the network is therefore not a concern when it comes to architecture selection . Fundamentally this approach attempts to learn the log-likelihood via a nonlinear regression , for which we chose an MLP . Log-likelihood labels for training were derived from empirical likelihoods functions ( details are given in the 'Materials and methods' section ) , which in turn were constructed as kernel density estimates ( KDEs ) . The construction of KDEs roughly followed Turner et al . , 2015 . We chose the Huber loss function ( details are given in the 'Materials and methods' section ) because it is more robust to outliers and thus less susceptible to distortions that can arise in the tails of distributions . Our second approach is based on a CNN architecture . Whereas the MLP learned to output a single scalar likelihood output for each datapoint ( ‘trial’ , given a choice , reaction time , and parameter vector ) , the goal of the CNN was to evaluate , for the given model parameters , the likelihood of an arbitrary number of datapoints via one forward pass through the network . To do so , the output of the CNN was trained to produce a probability distribution over a discretized version of the dataspace , given a stochastic model and parameter vector . The network learns the mapping θ→logℓΦ ( θ| . ) . As a side benefit , in line with the methods proposed by Lueckmann et al . , 2019; Papamakarios et al . , 2019b , the CNN can in fact act as a surrogate simulator; for the purpose of this paper , we do not exploit this possibility however . Since here we attempt to learn distributions , labels were simple binned empirical likelihoods , and as a loss function we chose the Kullback Leibler ( KL ) divergence between the network's output distribution and the label distribution ( details are given in the 'Materials and methods' section ) . Most of the specifics regarding training procedures are discussed in detail in the 'Materials and methods' section; however , we mention some aspects here to aid readability . We used 1 . 5 M and 3 M parameter vectors ( based on 100k simulations each ) to train the MLP and CNN approach , respectively . We chose these numbers consistently across all models and in fact trained on less examples for the MLP only due to RAM limitations we faced on our machines ( which in principle can be circumvented ) . These numbers are purposely high ( but in fact quite achievable with access to a computing cluster , simulations for each model was on the order of hours only ) since we were interested in a workable proof of concept . We did not investigate the systematic minimization of training data; however , some crude experiments indicate that a decrease by an order of magnitude did not seriously affect performance . We emphasize that this is in line with the expressed philosophy of our approach . The point of amortization is to throw a lot of resources at the problem once so that downstream inference is made accessible even on basic setups ( a usual laptop ) . In case simulators are prohibitively expensive even for reasonably sized computer clusters , minimizing training data may gain more relevance . In such scenarios , training the networks will be very cheap compared to simulation time , which implies that retraining with progressively more simulations until one observes asymptotic test performance is a viable strategy . We chose variations of SSMs common in the cognitive neurosciences as our test bed ( Figure 3 ) . The range of models we consider permits great flexibility in allowable data distributions ( choices and response times ) . We believe that initial applications are most promising for such SSMs because ( 1 ) analytical likelihoods are available for the most common variants ( and thus provide an upper-bound benchmark for parameter recovery ) and ( 2 ) there exist many other interesting variants for which no analytic solution exists . We note that there is an intermediate case in which numerical methods can be applied to obtain likelihoods for a broader class of models ( e . g . , Shinn et al . , 2020 ) . These methods are nevertheless computationally expensive and do not necessarily afford rapid posterior inference . Therefore , amortization via LANs is attractive even for these models . Figure 3 further outlines this distinction . We emphasize that our methods are quite general , and any model that generates discrete choices and response times from which simulations can be drawn within a reasonable amount of time can be suitable to the amortization techniques discussed in this paper ( given that the model has parameter vectors of dimension roughly lt15 ) . In fact , LANs are not restricted to models of reaction time and choice to begin with , even though we focus on these as test beds . As a general principle , all models tested below are based on stochastic differential equations ( SDEs ) of the following form:d⁢𝐗t=a⁢ ( t , x ) ⁢d⁢t+b⁢ ( t , x ) ⁢d⁢𝐁t , 𝐗0=wwhere we are concerned with the probabilistic behavior of the particle ( or vector of particles ) 𝐗 . The behavior of this particle is driven by a⁢ ( t , x ) , an underlying drift function , b⁢ ( t , x ) , an underlying noise transformation function , Bt , an incremental noise process , and X0=w , a starting point . Of interest to us are specifically the properties of the first-passage-time-distributions ( FPTD ) for such processes , which are needed to compute the likelihood of a given response time/choice pair {r⁢t , c} . In these models , the exit region of the particle ( i . e . , the specific boundary it reaches ) determines the choice , and the time point of that exit determines the response time . The joint distribution of choices and reaction times is referred to as a FPTD . Given some exit region ℰ , such FPTDs are formally defined asfℰ ( t ) =p ( infτ{ ( Xτ∈ℰ ) }=t ) In words , a first passage time , for a given exit region , is defined as the first time point of entry into the exit region , and the FPTD is the probability , respectively , of such an exit happening at any specified time t ( Ratcliff , 1978; Feller , 1968 ) . Partitioning the exit region into subsets ℰ1 , … . , ℰn ( e . g . , representing n choices ) , we can now define the set of defective distributions{fℰ1⁢ ( t;θ ) , … , fℰN⁢ ( t;θ ) }where θ∈Θ describes the collection of parameters driving the process . For every ℰi , ∫[0 , +∞]fℰi⁢ ( t;θ ) ⁢𝑑t=P⁢ ( 𝐗⁢exits into⁢ℰi ) =P⁢ ( i⁢gets chosen ) {fℰ𝒾 , … , fℰ𝓃} jointly define the FPTD such that∑i=1n∫[0 , +∞]fℰi⁢ ( t;θ ) ⁢𝑑t=1 These functions fℰi , jointly serve as the likelihood function s . t . ℓ⁢ ( θ;{r⁢t , c} ) =fℰc⁢ ( t;θ ) For illustration , we focus the general model formulation above to the standard DDM . Details regarding the other models in our test bed are relegated to the 'Materials and methods' section . To obtain the DDM from the general equation above , we set a⁢ ( t , x ) =v ( a fixed drift across time ) , b⁢ ( t , x ) =1 ( a fixed noise variance across time ) , and Δ⁢𝐁∼𝒩⁢ ( 0 , Δ⁢t ) . The DDM applies to the two alternative decision cases , where decision corresponds to the particle crossings of an upper or lower fixed boundary . Hence , ℰ1={ℝ≥a} and ℰ2={ℝ≤-a} , where a is a parameter of the model . The DDM also includes a normalized starting point w ( capturing potential response biases or priors ) , and finally a nondecision time τ ( capturing the time for perceptual encoding and motor output ) . Hence , the parameter vector for the DDM is then θ= ( v , a , w , τ ) . The SDE is defined asd⁢𝐗τ+t=v⁢d⁢t+d⁢𝐖 , 𝐗τ=w The DDM serves principally as a basic proof of concept for us , in that it is a model for which we can compute the exact likelihoods analytically ( Feller , 1968; Navarro and Fuss , 2009 ) . The other models chosen for our test bed systematically relax some of the fixed assumptions of the basic DDM , as illustrated in Figure 3 . We note that many new models can be constructed from the components tested here . As an example of this modularity , we introduce inhibition/excitation to the race model , which gives us the leaky competing accumulator ( LCA ) ( Usher and McClelland , 2001 ) . We could then further extend this model by introducing parameterized bounds . We could introduce reinforcement learning parameters to a DDM ( Pedersen and Frank , 2020 ) or in combination with any of the other decision models . Again we emphasize that while these diffusion-based models provide a large test bed for our proposed methods , applications are in no way restricted to this class of models .
Across epochs of training , both training and validation loss decrease rapidly and remain low ( Figure 4A ) suggesting that overfitting is not an issue , which is sensible in this context . The low validation loss further shows that the network can interpolate likelihoods to specific parameter values it has not been exposed to ( with the caveat that it has to be exposed to the same range; no claims are made about extrapolation ) . Indeed , a simple interrogation of the learned likelihood manifolds shows that they smoothly vary in an interpretable fashion with respect to changes in generative model parameters ( Figure 4B ) . Moreover , Figure 4C shows that the MLP likelihoods mirror those obtained by KDEs using 100 , 000 simulations , even though the model parameter vectors were drawn randomly and thus not trained per se . We also note that the MLP likelihoods appropriately filter out simulation noise ( random fluctuations in the KDE empirical likelihoods across separate simulation runs of 100⁢K samples each ) . This observation can also be gleaned from Figure 4C , which shows the learned likelihood to sit right at the center of sampled KDEs ( note that for each subplot 100 such KDEs were used ) . As illustrated in the Appendix , these observations hold across all tested models . One perspective on this is to consider the MLP likelihoods as equivalent to KDE likelihoods derived from a much larger number of underlying samples and interpolated . The results for the CNN ( not shown to avoid redundancy ) mirror the MLP results . Finally , while Figure 4 depicts the learned likelihood for the simple DDM for illustration purposes , the same conclusions apply to the learned manifolds for all of the tested models ( as shown in Appendix 1—figures 1–6 ) . Indeed inspection of those manifolds is insightful for facilitating interpretation of the dynamics of the underlying models , how they differ from each other , and the corresponding RT distributions that can be captured . Figure 7 shows that both the MLP and CNN methods consistently yield very good to excellent parameter recovery performance for the LC model , with parameter-wise regression coefficients globally above R2>0 . 9 . As shown in Figure 8 , parameter recovery for the Weibull model is less successful , however , particularly for the Weibull collapsing bound parameters . The drift parameter v , the starting point bias w , and the nondecision time are estimated well; however , the boundary parameters a , α , and β are less well recovered by the posterior mean . Judging by the parameter recovery plot , the MLP seems to perform slightly less well on the boundary parameters when compared to the CNN . To interrogate the source of the poor recovery of α and β parameters , we considered the possibility that the model itself may have issues with identifiability , rather than poor fit . Figure 8 shows that indeed , for two representative datasets in which these parameters are poorly recovered , the model nearly perfectly reproduces the ground truth data in the posterior predictive RT distributions . Moreover , we find that whereas the individual Weibull parameters are poorly recovered , the net boundary B⁢ ( t ) is very well recovered , particularly when evaluated within the range of the observed dataset . This result is reminiscent of the literature on sloppy models ( Gutenkunst et al . , 2007 ) , where sloppiness implies that various parameter configurations can have the same impact on the data . Moreover , two further conclusions can be drawn from this analysis . First , when fitting the Weibull model , researchers should interpret the bound trajectory as a latent parameter rather than the individual α and β parameters per se . Second , the Weibull model may be considered as viable only if the estimated bound trajectory varies sufficiently within the range of the empirical RT distributions . If the bound is instead flat or linearly declining in that range , the simple DDM or LC models may be preferred , and their simpler form would imply that they would be selected by any reasonable model comparison metric . Lastly , given our results the Weibull model could likely benefit from reparameterization if the desire is to recover individual parameters rather than the bound trajectory B⁢ ( t ) . Given the common use of this model in collapsing bound studies ( Hawkins et al . , 2015 ) and that the bound trajectories are nevertheless interpretable , we leave this issue for future work . The Appendix shows parameter recovery studies on a number of other stochastic simulators with non-analytical likelihoods , described in the 'Test bed' section . The appendices show tables of parameter-wise recovery R2 for all models tested . In general , recovery ranges from good to excellent . Given the Weibull results above , we attribute the less good recovery for some of these models to identifiability issues and specific dataset properties rather than to the method per se . We note that our parameter recovery studies here are in general constrained to the simplest inference setting equivalent to a single-subject , single-condition experimental design . Moreover , we use uninformative priors for all parameters of all models . Thus , these results provide a lower bound on parameter recoverability , provided of course that the datasets were generated from the valid parameter ranges on which the networks were trained; see Section 0 . 10 for how recovery can benefit from more complex experimental designs with additional task conditions , which more faithfully represents the typical inference scenario deployed by cognitive neuroscientists . Lastly , some general remarks about the parameter recovery performance . A few factors can negatively impact how well one can recover parameters . First , if the model generally suffers from identifiability issues , the resulting tradeoffs in the parameters can lead the MCMC chain to get stuck on the boundary for one or more parameters . This issue is endemic to all models and unrelated to likelihood-free methods or LANs , and should at best be attacked at the level of reparameterization ( or a different experimental design that can disentangle model parameters ) . Second , if the generative parameters of a dataset are too close to ( or beyond ) the bounds of the trained parameter space , we may also end with a chain that gets stuck on the boundary of the parameter space . We confronted this problem by training on parameter spaces that yield response time distributions that are broader than typically observed experimentally for models of this class , while also excluding obviously defective parameter setups . Defective parameter setups were defined in the context of our applications as parameter vectors that generate data that never allow one or the other choice to occur ( as in p⁢ ( c ) ≪1100 , 000 , data that concentrates more than half of the reaction times within a single 1 ms bin and data that generated mean reaction times beyond 10 s ) . These guidelines were chosen as a mix of basic rationale and domain knowledge regarding usual applications of DDMs to experimental data . As such , the definition of defective data may depend on the model under consideration . A major motivation for this work is the amortization of network training time during inference , affording researchers the ability to test a variety of theoretically interesting models for linking brain and behavior without large computational cost . To quantify this advantage , we provide some results on the posterior sampling runtimes using ( 1 ) the MLP with slice sampling ( Neal , 2003 ) and ( 2 ) CNN with iterated importance sampling . The MLP timings are based on slice sampling ( Neal , 2003 ) , with a minimum of n=2000 samples . The sampler was stopped at some n≥2000 , for which the Geweke statistic ( Geweke , 1992 ) indicated convergence ( the statistic was computed once every 100 samples for n≥2000 ) . Using an alternative sampler , based on differential evolution Markov chain Monte Carlo ( DEMCMC ) and stopped when the Gelman–Rubin R^<1 . 1 ( Gelman and Rubin , 1992 ) , yielded very similar timing results and was omitted in our figures . For the reported importance sampling runs , we used 200K importance samples per iteration , starting with γ values of 64 , which was first reduced to 1 where in iteration i , γi=642i-1 , before a stopping criterion based on relative improvement of the confusion metric was used . Figure 9A shows that all models can be estimated in the order of hundreds of seconds ( minutes ) , comprising a speed improvement of at least two orders of magnitude compared to traditional ABC methods using KDE during inference ( i . e . , the PDA method motivating this work; Turner et al . , 2015 ) . Indeed , this estimate is a lower bound on the speed improvement: we extrapolate only the observed difference between network evaluation and online simulations , ignoring the additional cost of constructing and evaluating the KDE-based likelihood . We decided to use this benchmark because it provides a fairer comparison to more recent PDA approaches in which the KDE evaluations can be sped up considerably ( Holmes , 2015 ) . Notably , due to its potential for parallelization ( especially on GPUs ) , our NN methods can even induce performance speedups relative to analytical likelihood evaluations . Indeed , Figure 9B shows that as the dataset grows runtime is significantly faster than even a highly optimized cython implementation of the Navarro Fuss algorithm ( Navarro and Fuss , 2009 ) for evaluation of the analytic DDM likelihood . This is also noteworthy in light of the full-DDM ( as described in the 'Test bed' section ) , for which it is currently common to compute the likelihood term via quadrature methods , in turn based on repeated evaluations of the Navarro Fuss algorithm . This can easily inflate the evaluation time by 1–2 orders of magnitude . In contrast , evaluation times for the MLP and CNN are only marginally slower ( as a function of the slightly larger network size in response to higher dimensional inputs ) . We confirm ( omitted as separate figure ) from experiments with the HDDM Python toolbox that our methods end up approximately 10-50 times faster for the full-DDM than the current implementation based on numerical integration , maintaining comparable parameter recovery performance . We strongly suspect there to be additional remaining potential for performance optimization . One of the principal benefits of LANs is that they can be directly extended – without further training – to arbitrary hierarchical inference scenarios , including those in which ( 1 ) individual participant parameters are drawn from group distributions , ( 2 ) some parameters are pooled and others separated across task conditions , and ( 3 ) neural measures are estimated as regressors on model parameters ( Figure 10 ) . Hierarchical inference is critical for improving parameter estimation particularly for realistic cognitive neuroscience datasets in which thousands of trials are not available for each participant and/or where one estimates impacts of noisy physiological signals onto model parameters ( Wiecki et al . , 2013; Boehm et al . , 2018; Vandekerckhove et al . , 2011; Ratcliff and Childers , 2015 ) . To provide a proof of concept , we developed an extension to the HDDM Python toolbox ( Wiecki et al . , 2013 ) , widely used for hierarchical inference of the DDM applied to such settings . Lifting the restriction of previous versions of HDDM to only DDM variants with analytical likelihoods , we imported the MLP likelihoods for all two-choice models considered in this paper . Note that GPU-based computation is supported out of the box , which can easily be exploited with minimal overhead using free versions of Google’s Colab notebooks . We generally observed GPUs to improve speed approximately fivefold over CPU-based setups for the inference scenarios we tested . Preliminary access to this interface and corresponding instructions can be found at https://github . com/lnccbrown/lans/tree/master/hddmnn_tutorial ( copy archived at swh:1:rev:e3369b9df138c75d0e490be0c48c53ded3e3a1d6 ) ; Fengler , 2021 . Figure 11 shows example results from hierarchical inference using the LC model , applied to synthetic datasets comprising 5 and 20 subjects ( a superset of participants ) . Recovery of individual parameters was adequate even for five participants , and we also observe the expected improvement of recovery of the group-level parameters μ and σ for 20 participants . Figure 12 shows an example that illustrates how parameter recovery is affected when a dataset contains multiple experimental conditions ( e . g . , different difficulty levels ) . It is common in such scenarios to allow task conditions to affect a single ( or subset ) -model parameter ( in the cases shown: v ) , while other model parameters are pooled across conditions . As expected , for both the full-DDM ( A ) and the Levy model ( B ) , the estimation of global parameters is improved when increasing the number of conditions from 1 to 5 to 10 ( left to right , where the former are subsets of the latter datasets ) . These experiments confirm that one can more confidently estimate parameters that are otherwise difficult to estimate such as the noise α in the Levy model and s⁢v the standard deviation of the drift in the full-DDM . Both of these experiments provide evidence that our MLPs provide approximate likelihoods which behave in accordance with what is expected from proper analytical methods , while also demonstrating their robustness to other samplers ( i . e . , we used HDDM slice samplers without further modification for all models ) . We expect that proper setting of prior distributions ( uniform in our examples ) and further refinements to the slice sampler settings ( to help mode discovery ) can improve these results even further . We include only the MLP method in this section since it is most immediately amenable to the kind of trial-by-trial-level analysis that HDDM is designed for . We plan to investigate the feasibility of including the CNN method into HDDM in future projects .
There are several limitations of the methods presented in this paper , which we hope to address in future work . While allowing for great flexibility , the MLP approach suffers from the drawback that we do not enforce ( or exploit ) the constraint that ℓ^⁢ ( θ|𝐱 ) is a valid probability distribution , and hence the networks have to learn this constraint implicitly and approximately . Enforcing this constraint has the potential to improve estimation of tail probabilities ( a known issue for KDE approaches to ABC more generally; Turner et al . , 2015 ) . The CNN encapsulation exploits the fact that ∫𝒳ℓ^⁢ ( θ|𝐱 ) ⁢𝑑𝐱=1 , however , makes estimation of trial-by-trial effects more resource hungry . We plan to investigate the potential of the CNN for trial-by-trial estimation in future research . A potential solution that combines the strengths of both the CNN and MLP methods is to utilize mixture density networks to encapsulate the likelihood functions . We are currently exploring this avenue . Mixture density networks have been successfully applied in the context of ABC ( Papamakarios and Murray , 2016 ) ; however , training can be unstable without extra care ( Guillaumes , 2017 ) . Similarly , invertible flows ( Rezende and Mohamed , 2015 ) and/or mixture density networks Bishop , 1994 may be used to learn likelihood functions ( Papamakarios et al . , 2019b; Lueckmann et al . , 2019 ) ; however , the philosophy remains focused on distributions of summary statistics for single datasets . While impressive improvements have materialized at the intersection of ABC and deep learning methods ( Papamakarios et al . , 2019a; Greenberg et al . , 2019; Gonçalves et al . , 2020 ) ( showing some success with posterior amortization for models up to 30 parameters , but restricted to a local region of high posterior density in the resulting parameter space ) , generally less attention has been paid to amortization methods that are not only of case-specific efficiency but sufficiently modular to serve a large variety of inference scenarios ( e . g . , Figure 10 ) . This is an important gap that we believe the popularization of the powerful ABC framework in the domain of experimental science hinges upon . A second and short-term avenue for future work is the incorporation of our presented methods into the HDDM Python toolbox ( Wiecki et al . , 2013 ) to extend its capabilities to a larger variety of SSMs . Initial work in this direction has been completed , the alpha version of the extension being available in the form of a tutorial at https://github . com/lnccbrown/lans/tree/master/hddmnn_tutorial . Our current training pipeline can be further optimized on two fronts . First , no attempt was made to minimize the size of the network needed to reliably approximate likelihood functions so as to further improve computational speed . Second , little attempt was made to optimize the amount of training provided to networks . For the models explored here , we found it sufficient to simply train the networks for a very large number of simulated datapoints such that interpolation across the manifold was possible . However , as model complexity increases , it would be useful to obtain a measure of the networks’ uncertainty over likelihood estimates for any given parameter vector . Such uncertainty estimates would be beneficial for multiple reasons . One such benefit would be to provide a handle on the reliability of sampling , given the parameter region . Moreover , such uncertainty estimates could be used to guide the online generation of training data to train the networks in regions with high uncertainty . At the intersection of ABC and NNs , active learning has been explored via uncertainty estimates based on network ensembles ( Lueckmann et al . , 2019 ) . We plan to additionally explore the use of Bayesian NNs , which provide uncertainty over their weights , for this purpose ( Neal , 1995 ) . One more general shortcoming of our methods is the reliance on empirical likelihoods for training , which in turn are based on a fixed number of samples across parameter vectors , just as the PDA method proposed by Turner et al . , 2015 . Recently , this approach has been criticized fundamentally on grounds of producing bias in the generated KDE-based likelihood estimates ( van Opheusden et al . , 2020 ) . A reduction of the approximate likelihood problem to one of inverse binomial sampling was proposed ( van Opheusden et al . , 2020 ) , which will generate unbiased likelihood estimates . To address this concern , we will investigate adaptive strategies for the selection of the simulations count n . We however highlight two points here that aim to put the promise of unbiased likelihoods in perspective . First , our networks add interpolation to the actual estimation of a likelihood . Likelihoods close in parameter space therefore share information that translates into an effectively higher simulation count than the 100k chosen to construct each empirical likelihood used for training . Quantifying this benefit precisely we leave for future research; however , we suspect , as suggested by Figure 4 , that it may be substantial . Second , while we generally acknowledge that bias in the tails remains somewhat of an issue in our approach , resolution is at best partial even in the proposed methods of van Opheusden et al . , 2020 . For the estimation of parameters for which a given datapoint is extremely unlikely ( i . e . , the data is generally unlikely under the model ) , the authors suggest to threshold the simulation count so that their algorithm is guaranteed to stop . This effectively amounts to explicitly allowing for bias again . As another alternative , the authors suggest to introduce a lapse rate in the generative model , which the LAN approach can accommodate as well . However , the introduction of a lapse rate does not deal with tail events directly either , but rather assumes that tail events are unrelated to the process of interest . This in turn will render a lower but fixed number of simulations N feasible for training LANs as well . This is notwithstanding the desirable minimization of simulation times even for high likelihood events , especially when trial-wise simulations are in fact necessary ( which tends to be in cases where amortization with LANs is a priori not a good computational strategy to begin with ) . Hence , although the inverse binomial sampling approach is elegant conceptually , excessive computation remains an issue when we need accurate estimates of the probability of actual tail events . Generally , however , we maintain it is desirable and important for future work to make use of the otherwise great potential of adaptive sampling to minimize total computational cost . Furthermore , we relegate to future research proper exploitation of the fact that LANs are by design differentiable in the parameters . We are currently working on an integration of LANs with Tensorflow probability ( Abadi et al . , 2016 ) , utilizing autograd to switch our MCMC method to the gradient-based NUTS sampler ( Hoffman and Gelman , 2014 ) . The main benefits of this sampler are robust mixing behavior , tolerance for high levels of correlations in the parameter space , while at the same time maintaining the ability to sample from high-dimensional posteriors . High level of correlations in posteriors is traditionally an Achilles' heel of the otherwise robust coordinate-wise slice samplers . DEMCMC and iterated Importance samplers are somewhat more robust in this regards; however , both may not scale efficiently to high-dimensional problems . Robustness concerns aside , initial numerical results additionally show some promising further speedups . Another important branch for future work lies in the utilization of LANs for model comparison . Initial results are promising in that we obtained satisfactory model recovery using the deviance information criterion ( DIC ) used for model selection in the standard HDDM package . However , this issue demands much more attention to evaluate other model selection metrics and extensive further numerical experiments , which we relegate to future work . Lastly , in contrast to the importance sampler driving the posterior inference for the CNN , we believe that some of the performance deficiencies of the MLP are the result of our MCs not having converged to the target distribution . A common problem seems to be that the sampler hits the bounds of the constrained parameter space and does not recover from that . As shown in Figures 7 and 8 , even ostensibly bad parameter recoveries follow a conceptual coherence and lead to good posterior predictive performance . We therefore may be underreporting the performance of the MLP and plan to test the method on an even more comprehensive suite of MCMC samplers , moreover including thus far neglected potential for reparameterization .
The Race model departs from previous model formulations in that it has a particle for each of N choice options instead of a single particle representing the evidence for one option over another . The function fEi⁢ ( t , θ ) now represents the probability of particle i to be the first of all particle to cross the bound a at time t . We consider race models for which the drift and starting point can vary for each particle separately . Treating the boundary as a constant a leaves us with a parameter vector θ= ( v1 , . . . , vn , a , w1 , . . . , wn , ndt ) . The SDE is defined for each particle separately ( or in vector form ) asd⁢𝐗τ+ti=vi⁢d⁢t+d⁢𝐖 , 𝐗0i=…=𝐗τi=wi These models represent the most straightforward extension to a multichoice scenario . In this section , we clarify a few strengths and weaknesses of the two presented methods and their respective use cases . First , representing the likelihood function datapoint-wise as an MLP output , or globally via the CNN output histogram , affects the potential for parallelization . As exploited by the choice of sampler , the CNN is very amenable to parallelization across parameters since inputs are parameter tuples only . Since the output is represented as a global likelihood histogram , the dataset likelihood is computed as the summation of the elementwise multiplied of bin-log-likelihoods , with a correspondingly binned dataset ( counts over bins ) . This has the highly desirable property of making evaluation cost ( time ) independent of dataset size . While the MLP in principle allows parallel processing of inputs , the datapoint-wise representation of input values ( {θ , x} ) makes the potential for cross-parameter parallelization dependent on dataset sizes . While a single evaluation of the CNN is more costly , cross-parameter batch processing can make it preferable to the MLP . Second , the CNN has an advantage during training , where the representation of the output as a softmax layer , and corresponding training via minimization of the KL divergence , provides a more robust training signal to ensure probability distributions compared to the purely local one in which the MLP learns a scalar likelihood output as a simple regression problem . Third , and conversely , the MLP formulation is more natural for trial-wise parameter estimates since the histogram representations may be redundant in case datapoints are in fact evaluated one by one ( given datapoint-wise parameters induced by trial-by-trial effects on parameter vectors ) . Give equivalent success in learning likelihoods , we see potential for speedup when using the pointwise approach in this case . In principle , both approaches however allow one to estimate the impact of trial-wise regressors on model parameters during inference , without further training . It is , for example , common in the cognitive neuroscience literature to allow the cross-trial time-course of EEG , fMRI , or spike signals to be modeled as a trial-by-trial regressor on model parameters of , for example , DDMs ( Wiecki et al . , 2013; Frank et al . , 2015; Cavanagh et al . , 2011; Herz et al . , 2016; Pedersen and Frank , 2020 ) . Another relevant example is the incorporation of latent learning dynamics . If a subject’s choice behavior is driven by reinforcement learning across stimuli , we can translate this into trial-by-trial effects on the parameter vectors of a generative process model ( Pedersen and Frank , 2020 ) . These applications are implicitly enabled at no extra cost with the MLP method , while the trial-by-trial split multiplies the necessary computations for the CNN by the number N of datapoints when compared to scenarios that only need dataset-wise parameters . We stress again , however , that in general both the CNN and the MLP can directly be used for hierarchical inference scenarios . The preceding discussion pertains to further potential for optimization and relative strengths , not categorical potential for application to a given scenario . With respect to the latter , both methods are essentially equal .
In cognitive neuroscience , computational modeling can formally adjudicate between theories and affords quantitative fits to behavioral/brain data . Pragmatically , however , the space of plausible generative models considered is dramatically limited by the set of models with known likelihood functions . For many models , the lack of a closed-form likelihood typically impedes Bayesian inference methods . As a result , standard models are evaluated for convenience , even when other models might be superior . Likelihood-free methods exist but are limited by their computational cost or their restriction to particular inference scenarios . Here , we propose neural networks that learn approximate likelihoods for arbitrary generative models , allowing fast posterior sampling with only a one-off cost for model simulations that is amortized for future inference . We show that these methods can accurately recover posterior parameter distributions for a variety of neurocognitive process models . We provide code allowing users to deploy these methods for arbitrary hierarchical model instantiations without further training .
Cognitive neuroscience studies the links between the physical brain and cognition . Computational models that attempt to describe the relationships between the brain and specific behaviours quantitatively is becoming increasingly popular in this field . This approach may help determine the causes of certain behaviours and make predictions about what behaviours will be triggered by specific changes in the brain . Many of the computational models used in cognitive neuroscience are built based on experimental data . A good model can predict the results of new experiments given a specific set of conditions with few parameters . Candidate models are often called ‘generative’: models that can simulate data . However , typically , cognitive neuroscience studies require going the other way around: they need to infer models and their parameters from experimental data . Ideally , it should also be possible to properly assess the remaining uncertainty over the parameters after having access to the experimental data . To facilitate this , the Bayesian approach to statistical analysis has become popular in cognitive neuroscience . Common software tools for Bayesian estimation require a ‘likelihood function’ , which measures how well a model fits experimental data for given values of the unknown parameters . A major obstacle is that for all but the most common models , obtaining precise likelihood functions is computationally costly . In practice , this requirement limits researchers to evaluating and comparing a small subset of neurocognitive models for which a likelihood function is known . As a result , it is convenience , rather than theoretical interest , that guides this process . In order to provide one solution for this problem , Fengler et al . developed a method that allows users to perform Bayesian inference on a larger number of models without high simulation costs . This method uses likelihood approximation networks ( LANs ) , a computational tool that can estimate likelihood functions for a broad class of models of decision making , allowing researchers to estimate parameters and to measure how well models fit the data . Additionally , Fengler et al . provide both the code needed to build networks using their approach and a tutorial for users . The new method , along with the user-friendly tool , may help to discover more realistic brain dynamics underlying cognition and behaviour as well as alterations in brain function .