Abstract

We describe a Bayesian filtering scheme for nonlinear state-space models in continuous time. This scheme is called Generalised Filtering and furnishes posterior (conditional) densities on hidden states and unknown parameters generating observed data. Crucially, the scheme operates online, assimilating data to optimize the conditional density on time-varying states and time-invariant parameters. In contrast to Kalman and Particle smoothing, Generalised Filtering does not require a backwards pass. In contrast to variational schemes, it does not assume conditional independence between the states and parameters. Generalised Filtering optimises the conditional density with respect to a free-energy bound on the model's log-evidence. This optimisation uses the generalised motion of hidden states and parameters, under the prior assumption that the motion of the parameters is small. We describe the scheme, present comparative evaluations with a fixed-form variational version, and conclude with an illustrative application to a nonlinear state-space model of brain imaging time-series.

1. Introduction

This paper is about the inversion of dynamic causal models based on nonlinear state-space models in continuous time. These models are formulated in terms of ordinary or stochastic differential equations and are ubiquitous in the biological and physical sciences. The problem we address is how to make inferences about the hidden states and unknown parameters generating data, given only observed responses and prior beliefs about the form of the underlying generative model, and its parameters. The parameters here include quantities that parameterise the model’s equations of motion and control the amplitude (variance or inverse precision) of random fluctuations. If we consider the parameters and precisions as separable quantities, model inversion represents a triple estimation problem. There are relatively few schemes in the literature that can deal with problems of this sort. Classical filtering and smoothing schemes such as those based on Kalman and Particle filtering (e.g., [1, 2]) deal only with estimation of hidden states. Recently, we introduced several schemes that solve the triple estimation problem, using variational or ensemble learning [35]. Variational schemes of this sort simplify the problem by assuming conditional independence among sets of unknown quantities, usually states, parameters and precisions. This is a natural partition, in terms of time-varying hidden states and time-invariant parameters and precisions. The implicit mean-field approximation leads to schemes that optimise the posterior or conditional density on time-varying hidden states, while accumulating the sufficient statistics necessary to optimise the conditional density of parameters and precisions, after all the data have been observed.

In this paper, we dispense with the mean-field approximation and treat all unknown quantities as conditionally dependent variables, under the prior constraint that the changes in parameters and precisions are very small. This constraint is implemented by representing all unknown variables in generalised coordinates of motion, which allows one to optimise the moments of the joint posterior as data arrive. The resulting scheme enables an efficient assimilation of data and the possibility of online and real-time deconvolution. We refer to this Bayesian filtering in generalised coordinates as Generalised Filtering (GF). Furthermore, by assuming a fixed form for the conditional density (the Laplace assumption) one can reduce the triple estimation problem to integrating or solving a set of relatively simple ordinary differential equations. In this paper, we focus on GF under the Laplace assumption.

We have previously described Variational filtering in [3] and smoothing in [5] for probabilistic inference on the hidden states of generative models based upon stochastic differential equations (see also [68], for recent and related advances). Variational filtering and its fixed form variant (Dynamic Expectation Maximization) outperform extended Kalman filtering and Particle filtering (provided the true posterior density is unimodal) in terms of the accuracy of estimating hidden states [3]. This improvement rests upon using generalised coordinates of motion. In other words, instead of just trying to estimate the conditional density of hidden states, one optimises the conditional density on their generalised motion to arbitrarily high order. The implicit modelling of generalised states has several fundamental advantages. First, it can accommodate temporal correlations in random fluctuations on the hidden states (e.g., fractal time and spectra in biological systems [9]). In other words, random terms in the model’s stochastic differential equations can have analytic autocovariance functions, whose smoothness can be estimated. This allows one to eschew standard Weiner assumptions, which is important in many realistic settings, particularly in the analysis of biological time-series. Second, generalised states afford a very simple filtering scheme that is formally distinct from Kalman and Particle filtering. In brief, Variational filtering uses a gradient descent on a free-energy bound on the model’s (negative) log-evidence. This means that filtering can be reduced to the solution of differential equations that necessarily entail continuity constraints on the conditional estimates. Clearly, the free-energy is a functional of time, which means that the gradient descent has to “hit a moving target.” Generalised coordinates finesse this problem by placing the gradient descent in a frame of reference that moves with the conditional expectation or mean (see [4] for details). This is heuristically related to the separation of temporal scales in centre manifold theory [10], where the motion of free-energy minima (the centre manifold) enslaves the dynamics of the gradient descent.

Variational filtering of this sort is fundamentally different in its mathematical construction from conventional schemes like Kalman filtering because of its dynamical formulation. It can be implemented without any assumptions on the form of the conditional density by using an ensemble of “particles” that are subject to unit (standard) Wiener perturbations. The ensuing ensemble density tracks the conditional mean of the hidden states and its dispersion encodes conditional uncertainty. Variational filtering can be further simplified by assuming the ensemble density (conditional density) is Gaussian, using the Laplace assumption. Crucially, under this approximation, the conditional covariance (second moment of the conditional density) becomes an analytic function of the conditional mean. In other words, only the mean per se has to be optimized. This allows one to replace an ensemble of particles, whose motion is governed by stochastic differential equations, with a single ordinary differential equation describing the motion of the conditional mean. The solution of this ordinary differential equation corresponds to the D-step in Dynamic Expectation Maximisation (DEM). DEM comprises additional E (expectation) and M (maximization) steps that optimise the conditional density on parameters and precisions after the D (deconvolution) step has been solved. Iteration of these steps proceeds until convergence, in a way that is formally related to conventional variational schemes (cf. [5, 11]).

In this work, we retain the Laplace approximation to the conditional density but dispense with the mean-field approximation; in other words, we do not assume conditional independence between the states, parameters, and precisions. We implement this by absorbing parameters and precisions into the hidden states. This means that we can formulate a set of ordinary differential equations that describe the motion of time-dependent conditional means and implicitly the conditional precisions (inverse covariances) of all unknown variables. This furnishes (marginal) conditional densities on the parameters and precisions that are functionals of time. The associated conditional density of the average parameters and precisions over time can then be accessed using Bayesian parameter averaging. Treating time-invariant parameters (and precisions) as states rests on modelling their motion. Crucially, we impose prior knowledge that this motion is zero, leading to a gradient descent on free-energy, which is very smooth (cf. the use of filtering as a “second-order” technique for learning parameters [12]). In brief, the resulting scheme assimilates evidence in the generalised motion of data to provide a time-varying density on all the model’s unknown variables, where the marginal density on the parameters and precisions converges slowly to a steady-state solution. The scheme can be iterated until the time or path-integral of free-energy (free-action) converges and may represent a useful and pragmatic (online) alternative to variational schemes.

This paper comprises four sections. In the first, we describe the technical details of Generalised Filtering from the first principles. This section starts with the objective (to maximize the path-integral of a free-energy bound on a model’s log-evidence). It ends with set of ordinary differential equations, whose solution provides the conditional moments of a Gaussian approximation to the conditional density we seek. The second section reviews a generic model that embodies both dynamic and structural (hierarchical) constraints. We then look at Generalised Filtering from the first section, under this model. The third section presents comparative evaluations of GF using a simple linear convolution model, which is a special case of the model in Section 2. These evaluations are restricted to a comparison with DEM because DEM is formally the closest scheme to GF and has been compared with Extended Kalman filtering and Particle filtering previously [4]. In the final section, we apply Generalised Filtering to a neurobiological problem, namely, inferring the parameters and hidden physiological states generating a time-series of functional magnetic resonance imaging (fMRI) data from the human brain. We use this to illustrate the effect of the mean-field assumption implicit in DEM and establish the face validity of Generalised Filtering in terms of known neurobiology. We conclude with a brief discussion.

2. Generalised Filtering

In this section, we present the conceptual background and technical details behind Generalised Filtering, which (in principle) can be applied to any nonlinear state-space or dynamic causal model formulated with stochastic differential equations. Given the simplicity of the ensuing scheme, we also take the opportunity to consider state-dependant changes in the precision of random fluctuations. This represents a generalisation of our previous work on dynamic causal models and will be exploited in a neurobiological context, as a metaphor for attention (Feldman et al.; in preparation). However, we retain a focus on cascades of state-space models, which we have referred to previously as hierarchical dynamic models [13].

2.1. Filtering from Basic Principles

Given a model and generalised sensor data comprising real values, their velocity, acceleration, jerk, and so forth, we want to evaluate the log-evidence integrated over the time that data are observed (cf. [14, 15]): Generally, this path-integral cannot be evaluated directly; however, one can induce an upper bound that can be evaluated with a recognition density on the causes (i.e., states and parameters) of the data. We will see later that these causes comprise time-varying states and slowly varying parameters . This bound is based on free-energy, which is constructed by simply adding a nonnegative term to the (negative) log-evidence [11, 16, 17]. The resulting integral is called free-action because it is a path-integral of free-energy (see also [18]) By Gibb’s inequality the Kullback-Leibler divergence is greater than zero, with equality when is the true conditional density. In this case, (negative) free-action becomes accumulated log-evidence . Minimising this bound, by optimising the recognition density at each time point, makes it an approximate conditional density on the causes and makes free-action a bound approximation to the accumulated log-evidence. The approximate conditional density can then be used for inference on the states or parameters of a particular model, and the accumulated log-evidence can be used to compare different models (e.g., [19]).

Crucially, the free-energy can be evaluated easily because it is a function of and a generative model entailed by The free-energy has been expressed here in terms of , the negentropy of , and an energy expected under . In physics, is called Gibb’s energy and is a log-probability that reports the joint surprise about data and their causes. If we assume that is Gaussian (the Laplace assumption), then we can express free-energy in terms of its sufficient statistics, the mean and covariance of the recognition density where . Here and throughout, subscripts denote derivatives. We can now minimise free-action with respect to the conditional precisions (inverse covariances) by solving to give Equation (2.5) shows that the precision is an analytic function of the mean, which means all we have worry about is the (approximate) conditional mean. One can see this clearly by eliminating the conditional covariance to express the free-energy as a function of (and only of) the conditional means The conditional means, which minimise free-energy, are the solutions to the following ordinary differential equations. For the generalised states these equations are where is a derivative matrix operator with identity matrices above the leading diagonal, such that . Here and throughout, we assume that all gradients are evaluated at the mean; here . The stationary solution of (2.7), in a frame of reference that moves with the generalised motion of the mean, minimises free-energy and action. This can be seen by noting that the variation of free-action with respect to the solution is zero: This ensures that when Gibb’s energy is minimized, the mean of the motion is the motion of the mean, that is, . For slowly varying parameters this motion disappears and we can use the following scheme: Here, the solution minimises free-energy, under constraint that the motion of the mean is small; . This can be seen by noting Equations (2.7) and (2.9) prescribe recognition (filtering) dynamics for expected states and parameters respectively. The dynamics for states can be regarded as a gradient descent in a frame of reference that moves with the expected motion (cf. a moving target). Conversely, the dynamics for the parameters can be thought of as a gradient descent that resists transient fluctuations in free-energy with a damping term . This instantiates our prior belief that fluctuations in the parameters are small, where can be regarded as a prior precision. Figure 1 shows the implicit kernel (solid line) that is applied to fluctuating free-energy gradients to produce changes in the conditional mean. It can be seen that the height of the kernel is . We find that using ensures stable convergence in most situations. This renders the integral of the kernel over the time-series about one eighth (see Appendix A for a discussion of this assimilation scheme in relation to classical methods).

The solutions of (2.7) and (2.9) produce the time-dependent moments of an approximate Gaussian conditional density on the unknown states and parameters. These can be used for inference on states (or parameters) directly or for model comparison using the (bound on) accumulated log-evidence (see (2.3) and (2.6)). Because the conditional means of the parameters change slowly over time (unlike their conditional precisions), we can summarise the conditional estimates with the conditional density on the average over time using Bayesian parameter averaging Here, and are the conditional moments of the average while and are the mean and precision of the marginal conditional density on the parameters at each point in time. This summary will be overconfident because it neglects conditional dependencies due to temporally correlated random terms (and the use of generalised data). However, it usefully connects the conditional density with posterior distributions estimated in conventional schemes that treat the parameters as static.

In Generalised Filtering, changes in the conditional uncertainty about the parameters are modelled explicitly as part of a time-varying conditional density on states and parameters. In contrast, variational schemes optimise a conditional density on static parameters, that is, that is not a functional of time (although fluctuations in the conditional dependence between states and parameters are retained as mean-field terms that couple their respective marginals). This difference is related to a formal difference between the free-energy bound on log-evidence in variational schemes and free-action. Variational schemes return the free-energy of the accumulated data (e.g., see [3, equation ()]). Conversely, in Generalised Filtering, the free-energy is accumulated over time. This means that the variational free-energy bounds the log-evidence of accumulated data, while free-action bounds the accumulated log-evidence of data These are not the same because evidence is a nonlinear function of the data (see, Appendix B for an example). This distinction does not mean that one form of evidence accumulation is better than another; although it raises interesting issues for model comparison that will be pursued elsewhere (Li et al; in preparation). For the moment, we will approximate the free-action associated with the conditional densities from variational schemes by assuming that , where and (Appendix B). This allows us to compare fixed-form schemes with (DEM) and without (GF) the mean field assumption, in terms of free-action.

2.2. Summary

In summary, we have derived recognition or filtering dynamics for expected states and parameters (in generalised coordinates of motion), which cause data. The solutions to these equations minimise free-action (at least locally) and therefore minimise a bound on the accumulated evidence for a model of how we think the data were caused. This minimisation furnishes the conditional density on the unknown variables in terms of conditional means and precisions. The precise form of the filtering depends on the energy associated with a particular generative model. In what follows, we examine the forms entailed by hierarchical dynamic models.

3. Hierarchical Dynamic Models

In this section, we review the form of models that will be used in subsequent sections to evaluate Generalised Filtering. Consider the state-space model Using and unit noise , this model can also be written as The nonlinear functions represent a mapping from hidden states to observed data and the equations of motion of the hidden states. These equations are parameterised by . The states are variously referred to as sources or causes, while the hidden states meditate the influence of causes on data and endow the system with memory. We assume that the random fluctuations are analytic, such that the covariance of is well defined. Unlike our previous treatment of these models, we allow for state-dependent changes in the amplitude of random fluctuations. These effects are meditated by the vector and matrix functions and respectively, which are parameterised by first and second-order parameters .

Under local linearity assumptions, the generalised motion of the data and hidden states can be expressed compactly as where the generalised predictions are (with ) Gaussian assumptions about the random fluctuations prescribe a generative model in terms of a likelihood and empirical priors on the motion of hidden states These probability densities are encoded by their covariances or precisions with precision parameters that control the amplitude and smoothness of the random fluctuations. Generally, the covariances factorise into a covariance proper and a matrix of correlations among generalised fluctuations that encode their autocorrelation [4, 20]. In this paper, we will deal with simple covariance functions of the form . This renders the precision parameters log-precisions.

Given this generative model, we can now write down the energy as a function of the conditional expectations, in terms of a log-likelihood and log-priors on the motion of hidden states and the parameters (ignoring constants): This energy has a simple quadratic form, where the auxiliary variables are prediction errors for data, the motion of hidden states and parameters respectively. The predictions of the states are and the predictions of the parameters are their prior expectations. Equation (3.6) assumes flat priors on the states and that priors on the parameters are Gaussian , where This assumption allows us to express the learning scheme in (2.9) succinctly as , where . Equation (3.6) provides the form of the free-energy and its gradients required for filtering, where (with a slight abuse of notation), Note that the constant in the free-energy just includes because we have used Gaussian priors. The derivatives are provided in Appendix C. These have simple forms that comprise only quadratic terms and trace operators.

3.1. Hierarchical Forms

We next consider hierarchical forms of this model. These are just special cases of Equation (3.1), in which we make certain conditional independences explicit. Although they may look more complicated, they are simpler than the general form above. They are useful because they provide an empirical Bayesian perspective on inference [21, 22]. Hierarchical dynamic models have the following form Again, are continuous nonlinear functions and is a prior mean on the causes at the highest level. The random terms are conditionally independent and enter each level of the hierarchy. They play the role of observation noise at the first level and induce random fluctuations in the states at higher levels. The causes link levels, whereas the hidden states link dynamics over time. In hierarchical form, the output of one level acts as an input to the next. This input can enter nonlinearly to produce quite complicated generalised convolutions with “deep” (i.e., hierarchical) structure [22]. The energy for hierarchical models is (ignoring constants)

This is exactly the same as (3.6) but now includes extra terms that mediate empirical (structural) priors on the causes . These are induced by the conditional independence assumptions in hierarchical models. Note that the data enter the prediction errors at the lowest level such that  .

3.2. Summary

In summary, hierarchical dynamic models are nearly as complicated as one could imagine; they comprise causal and hidden states, whose dynamics can be coupled with arbitrary (analytic) nonlinear functions. Furthermore, the states can be subject to random fluctuations with state-dependent changes in amplitude and arbitrary (analytic) autocorrelation functions. A key aspect is their hierarchical form that induces empirical priors on the causes that link successive levels and complement the dynamic priors afforded by the model’s equations of motion (see [13] for more details). These models provide a form for the free-energy and its gradients that are needed for filtering, according to (3.8) (see Appendix C for details). We now evaluate this scheme using simulated and empirical data.

4. Comparative Evaluations

In this section, we generate synthetic data using a simple linear convolution model used previously to cross-validate Kalman filtering, Particle filtering, Variational filtering and DEM [3, 4]. Here we restrict the comparative evaluations to DEM because it is among the few schemes that can solve the triple estimation problem in the context of hierarchical models, and provides a useful benchmark in relation to other schemes. Our hope was to show that the conditional expectations of states, parameters and precisions were consistent between GF and DEM but that the GF provided more realistic conditional precisions that are not confounded by the mean-field assumption in implicit in DEM. To compare DEM and GF, we used the a model based on (3.9) that is specified with the functions

Here, the parameters comprise the state and other matrices. In this standard state-space model, noisy input perturbs hidden states, which decay exponentially to produce an output that is a linear mixture of hidden states. Our example used a single input, two hidden states and four outputs. This model was used to generate data for the examples below. This entails the integration of stochastic differential equations in generalised coordinates, which is relatively straightforward (see [4, Appendix  2]). We generated data over 32 time bins, using innovations sampled from Gaussian densities with the log-precisions of eight and six for observation noise and state noise respectively. We used orders of generalised motion in all simulations and all random fluctuations were smoothed with a Gaussian kernel whose standard deviation was one quarter of a time bin.

When generating data, we used a deterministic Gaussian bump function centred on . However, when inverting the model, this cause was unknown and was subject to mildly informative shrinkage priors with zero mean and unit precision; . These were implemented by fixing the log-precision . This model and an example of the data it generates are shown in Figure 2. The format of Figure 2 will be used in subsequent figures and shows the data and hidden states in the top panels and the causes below. The upper left panel shows the simulated data in terms of the output due to hidden states (coloured lines) and observation noise (red lines). The (noisy) dynamics of the hidden states are shown in the upper right panels, which are the response to the cause or input on the lower left. The model is shown as a Bayesian dependency graph on the lower right.

These data were then subject to GF and DEM to recover the conditional densities on the hidden states, unknown cause, parameters and log-precisions. For both schemes, we used uninformative priors on four parameters; , and set the remainder to their true value, with infinitely precise priors. The four parameters included two from the equations of motion and two from the observer function . The log-precision priors were , where . In effect, model inversion or filtering entails estimating the hidden states and unknown inputs perturbing these states while, at the same time, identifying the parameters underlying the influence of the hidden states on their motion (and the system’s response) and the precision of both observation and state noise. In addition, we estimated the smoothness of the random fluctuations (see Appendix D) but placed very tight priors on the smoothness parameters (whose prior expectations were the true values) to ensure a fair comparison with DEM. Note that we are trying to infer the inputs to the system, which makes this a particularly difficult problem. This inference is precluded in conventional filtering, which usually treats inputs as noise (or as known).

Figure 3 shows the conditional estimates of the states during the first iteration (pass through the time series) of the GF scheme. Here the predicted response and prediction error are shown on the upper left while the conditional expectations of the hidden states and cause are shown on the upper right and lower left respectively. For the hidden states (upper right) the conditional means are depicted by blue lines and the 90% conditional confidence regions by grey areas. These are sometimes referred to as “tubes”. Here, the confidence tubes are based upon the marginal conditional density of the states. The key thing to observe here is that the conditional confidence increases with time. This reflects smooth increases in the conditional log-precisions (lower right panel) as they assimilate gradients and are drawn from their initial value (prior expectation or four) toward the true values (of eight and six). Note further how the parameter estimates show similar drifts; however, there is a marked movement towards the true values (from the prior expectation of zero) when their role is disclosed by the perturbation at around fifteen time bins.

The dynamics of the conditional states are prescribed by (2.7) and the slower assimilation dynamics of the conditional parameters and log-precisions by (2.9). By passing repeatedly through the time-series, the parameter and log-precision estimates converge to their stationary solution, at which point free-action stops increasing. As they become better estimates, the estimates of the states become more veridical and confident. By about sixteen iterations we get the estimates shown in Figure 4 (each iteration takes a second or so on a modern PC). Figure 4 uses the same format as Figure 3 but replaces the time-dependent evolution of the parameters and log-precisions with the conditional estimates of their Bayesian average (lower right; see (2.11)). Here, we see that the confidence tubes on the hidden states have shrunk to the extent we can be very confident about the conditional means. In this figure, the confidence tube around the cause is shown and suggests there is much less confidence here; despite the fact that the estimates are remarkably close to the true values (shown as broken grey lines) during the onset and offset of the true cause.

The conditional estimates of the parameters show good agreement with the true values, with the exception of the first parameter of the equation of motion . Crucially, this parameter has the highest conditional uncertainty, shown in terms of 90% confidence intervals (red bars). Happily the true values lie within or near these intervals, with a slight overconfidence for the parameters of the observer function (second pair). These estimates of conditional uncertainty should now be compared with the equivalent results from the DEM scheme.

Figure 5 shows the results of DEM using exactly the same format as Figure 4. The first thing to note is the remarkable similarity between the conditional expectations, both for the states and parameters, which are virtually indistinguishable. The most poignant difference between the two schemes lies in the confidence intervals: Although the results of the DEM scheme look much “tighter” in terms of the confidence tubes on states, they fail to contain the true values at all times. This implicit overconfidence problem is even more acute for the parameter estimates that have very tight posterior precisions, most notably on the first parameter that was identified with low confidence by GF. This example illustrates how the mean-field assumption implicit in DEM can manifest in terms of overconfidence; and how relaxing this assumption with GF partly resolves this problem. We are not claiming that this is a ubiquitous behaviour but are just demonstrating that such differences are easy to find.

The free-action bound on accumulated log-evidence is shown in Figure 6 as a function of (twenty four) iterations of the GF and DEM schemes. The first thing to note is that GF furnishes a much tighter bound than DEM. This is not surprising because DEM optimises the free-energy of the accumulated data under mean-field assumptions, not the accumulated free-energy. However, it does speak to a profound difference in the bounds on evidence that might have important implications for Bayesian model comparison (this will be pursued in the work of Li et al.—in preparation). GF extremises free-action until the last few iterations, when there is a paradoxical (if small) reversal. We suppose that this is due to the first-order approximations employed during integration of the scheme (see Appendix C). This effect can be mitigated by using smaller step-sizes during solution of the recognition dynamics and by increasing the precision of the random fluctuation in parameters: in practice, we increase with each iteration, in proportion to the length of the time-series). The negative free-action increases for the first few iterations of DEM and then decreases to well below its starting point. A quantitative examination of the contributions to free-action suggests that this decrease is largely due to the high conditional log-precisions estimated by DEM (upper right panel).

It can be seen that DEM overestimates the precision of both observation and state noise while GF overestimates observation noise but underestimates state noise. Both schemes are overconfident about their estimate, in that the true values lie outside the 90% confidence intervals (red bars). These confidence intervals are based on accumulating the conditional precisions at each time step. For DEM, this accumulation is an integral part of optimisation whereas for GF it rests on the Bayesian parameter averaging of time-dependent precisions. These are shown on the lower right in terms of the corresponding confidence regions (grey areas). This panel shows a mild contraction of the confidence tube for the precision of state noise, when the hidden states are changing the most (shortly after the cause arrives). This is sensible because state noise is on the motion of hidden sates. A similar but more pronounced effect is seen in the equivalent confidence tubes for the parameters (lower left). Here, all the parameters estimates enjoy a transient decrease in conditional uncertainty during the perturbation because there is more information in the data about their putative role at these times.

4.1. Summary

In this section, we have tried to illustrate some of the basic features of Generalised Filtering and provide some comparative evaluations using an established and formally similar variational scheme (DEM). In this example, the estimates of the conditional means were very similar. The main difference emerged in the estimation of posterior confidence intervals and the behaviour of the free-action bound on accumulated log-evidence. These differences are largely attributable to the mean-field approximation inherent in DEM and related variational schemes. In the next section, we turn to a more complicated (and nonlinear) model to show that GF can recover causal structure from data, which DEM fails to disclose.

5. An Empirical Illustration

In this section, we turn to a model that is more representative of real-world applications and involves a larger number of states, whose motion is coupled in a nonlinear fashion. This model and the data used for its inversion have been presented previously in a comparative evaluation of variational filtering and DEM. Here, we use it to illustrate that the GF scheme operates with nonlinear models and to provide a face validation in this context. This validation rests upon analysing data from a part of the brain known to be functionally selective for visual motion processing [23]. We hoped to show that GF could establish a significant response to visual motion using evoked brain imaging responses. This example was chosen because inference about brain states from noninvasive neurophysiologic observations is an important issue in cognitive neuroscience and functional imaging (e.g., [2426]), and GF may play a useful role in stochastic dynamic causal modelling.

5.1. The Biophysical Model

We used a hemodynamic model of brain responses to explain evoked neuronal activity that has been described extensively in previous communications (e.g., [27, 28]). In brief, neuronal activity causes an increase in a vasodilatory signal that is subject to autoregulatory feedback. Blood flow responds in proportion to this signal and causes change in blood volume and deoxyhaemoglobin content. The observed signal is a nonlinear function of volume and deoxyhaemoglobin. These dynamics are modelled by the differential equations

In this model, changes in vasodilatory signal are elicited by neuronal input . Relative oxygen extraction is a function of flow, where is resting oxygen extraction fraction and outflow is a function of volume with Grubb's exponent . A description of the parameters of this model and their assumed values (prior expectations) are provided in Table 1. Blood flow, volume, and deoxyhaemoglobin concentration are all nonnegative quantities. One can implement this formal constraint with the transformation . Under this transformation the differential equations above can be written as

This allows us to formulate the model in terms of hidden states with unbounded support (i.e., the conditional means can be positive or negative) to give the following functions (see (3.9)):

This model represents a multiple-input, single-output model with four hidden states. The parameters of interest here were that couple three distinct neuronal responses to the vasodilatory signal . These evoked responses correspond to neuronal activity elicited experimentally, as described next:

5.2. Data and Preprocessing

Data were acquired from a normal subject at 2-Tesla using a Magnetom VISION (Siemens, Erlangen) whole body MRI system, during a visual attention study. Contiguous multislice images were obtained with a gradient echo-planar sequence (TE = 40 ms; TR = 3.22 seconds; matrix size = , voxel size  mm). Four consecutive hundred-scan sessions were acquired, comprising a sequence of 10-scan blocks under five conditions. The first was a dummy condition to allow for magnetic saturation effects. In the second, Fixation, subjects viewed a fixation point at the centre of a screen. In an Attention condition, subjects viewed 250 dots moving radially from the centre at 4.7 degrees per second and were asked to detect changes in radial velocity. In no attention, the subjects were asked simply to view the moving dots. In another condition, subjects viewed stationary dots. The order of the conditions alternated between Fixation and visual stimulation. In all conditions subjects fixated the centre of the screen. No overt response was required in any condition and there were no actual speed changes. The data were analysed using a conventional SPM analysis (http://www.fil.ion.ucl.ac.uk/spm). The activity from extrastriate cortex (motion-sensitive area V5 [23]) was summarised using the principal local eigenvariate of a region centred on the maximum of a contrast testing for the effect of visual motion. The first 256 samples from this regional response were used for Generalised Filtering and DEM.

The three potential causes of neuronal activity were encoded as box-car functions corresponding to the presence of a visual stimulus, motion in the visual field, and attention. These stimulus functions constitute the priors on the three causes in the model. The associated parameters, encode the degree to which these experimental effects induce hemodynamic responses. Given we selected a motion-selective part of the brain; one would anticipate that the conditional probability that exceeds zero would be large. The regional data were subject to GF using the model in (5.3) and informative shrinkage priors on all but the neuronal coupling parameters (see Table 1 for the prior means). We used mildly informative priors for the coupling parameters and similarly for the log-precisions . Otherwise, the analysis was identical to the analysis of the simulated data of the previous section (including unit prior precisions on the causes).

The ensuing conditional means and 90% confidence regions for the causal and hidden states are shown in Figure 7. The dynamics of inferred activity, flow and other biophysical states are physiologically plausible. For example, activity-dependent changes in flow are around , producing about a 3% change in fMRI signal. The conditional estimates of the neuronal coupling parameters, are shown on the lower right. As anticipated, we can be almost certain that neuronal activity encoding visual motion induces a response. Interestingly, both visual stimulation per se and motion appear to elicit responses in this area. This was somewhat unexpected because this area is meant to be selective for motion processing. The interpretation of these results is confounded by profound negative correlations between the visual and motion coupling parameters as evidenced by the conditional covariance among the parameters in Figure 8 (these correlations mean that one can infer confidently the sum of the two parameters but not their differences). This figure shows the conditional means of the eight (Bayesian average) parameters (upper left) and their conditional covariance (upper right). The first five (biophysical) parameters have been estimated with a high conditional confidence while the neuronal parameters are less precise and are interdependent. These conditional dependencies arise from the experiential design, which is highly inefficient for disambiguating between visually evoked responses and motion-specific responses. This is because there were only three blocks of static stimuli, compared to twelve blocks of visual motion. We confirmed this by simulating data using the conditional estimates of the parameters and precision from the current analysis but enforcing motion-specific responses by making , and . The conditional estimates were virtually identical to those from the empirical analysis in Figure 8 (results not shown). This highlights the importance of interpreting not just the marginal conditional density of a single parameter but the conditional density over all parameters estimated.

Figure 8 also shows the time-dependent changes in conditional confidence during the experimental session for the parameters (lower left) and log-precisions (lower right). We have focused on the conditional precision of the visual coupling parameter (grey area) to highlight the increase in precision during the two periods of stationary visual stimulation. These are the only times that the data inform the strength of this parameter in a way that is not confounded by the simultaneous presence of motion in the visual stimulus. Similar transient increases in the conditional precision of the log-precision estimates are seen at the onset and offset of visual stimulation on the lower right.

A detailed summary of the hemodynamics is shown in Figure 9. This figure focuses on the first 120 time bins and plots the hemodynamic states in terms of the conditional mean of . Each time bin corresponds to 3.22 seconds. The hidden states are shown in the upper panel, overlaid on periods (grey) of visual motion. It can be seen that a mixture of neuronal activity (the conditional estimate of shown in the lower panel), induces a transient burst of vasodilatory signal (blue), which is suppressed rapidly by the resulting increase in flow (green). The increase in flow dilates the venous capillary bed to increase volume (red) and dilute deoxyhaemoglobin (cyan). The concentration of deoxyhaemoglobin determines the measured activity. Interestingly, the underlying neuronal activity appears to show adaptation during visual stimulation and a marked offset transient in nearly all the epochs shown. Note that the conditional densities of the hemodynamic states are non-Gaussian (i.e., lognormal) despite the Laplace assumption entailed by the filtering scheme This is an example of how nonlinear models, under Gaussian assumptions, can be used to model non-Gaussian densities.

5.3. Comparison with DEM

Finally, we analysed the same data using DEM. The results are shown in Figure 10 using the same format as Figure 8. It is immediately obvious that DEM has failed to detect visual motion-dependent responses in this brain area. The coupling parameter estimates are small, both quantitatively and statistically (their confidence intervals include zero). This failure is also evident in the (negative) free-action bound on accumulated log-evidence in Figure 11. This shows that GF provides a much tighter bound relative to DEM. Again, we do not mean to suggest that this is a generic behaviour of variational schemes, just that one can find examples where the mean-field assumption can have a profound effect on inference. Having said this, it took some time to find the length of the time-series and priors on the log-precisions that revealed this difference. In most of our analyses, the conditional estimates from GF and DEM were very similar. The behaviour of the GF free-action over iterations is interesting and speaks to the important point that online evidence accumulation leads to an efficient inversion scheme for long time-series. This is because the conditional moments of the parameters and precisions are near optimal at the end of long sequences (as in conventional assimilation schemes). This is manifest by the convergence of free-action within a handful of iterations (four in the example here, where each iteration took about ten seconds). This computationally efficient form of data assimilation was one of the pragmatic motivations for developing Generalised Filtering.

5.4. Summary

As noted in [4], it is perhaps remarkable that so much conditional information about the underlying neuronal and hemodynamics can be extracted from a single time-series, given only the functional form of its generation. This speaks to the power of generative modelling, in which constraints on the form of the model allow one to access hidden quantities. To date, dynamic causal models of neuronal systems, measured using fMRI or electroencephalography (EEG), have used known, deterministic causes and have ignored state-noise (see, [25, 26] for important exceptions). One of the motivations for Generalised Filtering was to develop a computational efficient inference scheme that can deal with correlated state-noise, of the sort seen in biological time-series.

6. Conclusion

In this paper, we have introduced Generalised Filtering, an online Bayesian scheme for inverting generative models cast as stochastic differential equations in generalised coordinates of motion. This scheme is based upon a path-integral optimisation of free-energy, where free-energy bounds the log-evidence for a model. Under a Laplace approximation to the true posterior density on the model’s unknown variables, one can formulate deconvolution or model inversion as a set of ordinary differential equations, whose solution provides their conditional mean (which implicitly prescribes their conditional precision). Crucially, this density covers not only time-varying hidden states but also parameters, and precisions that change slowly. We have seen that its performance is consistent with equivalent fixed-form variational schemes (Dynamic Expectation Maximisation) that entail the extra assumption that the states, parameters and precisions are conditionally independent.

Although not emphasised in this paper, the basic approach on which Generalised Filtering is based was developed with neurobiological implementation in mind. In other words, we have tried to construct a scheme that could be implemented by the brain in a neurobiologically plausible fashion. This was one of the primary motivations for a dynamical optimisation of the parameter and precision estimates. In future communications, we will focus on the neurobiological interpretation of Generalised Filtering and how it might relate to the optimisation of synaptic activity, efficacy, and gain during perceptual inference in the brain. Our particular focus here will be on state-dependent changes in precision as a model of visual attention (Feldman et al; in preparation). In this context, the recognition dynamics entailed by optimisation can be regarded as simulations of neuronal responses to sensory inputs.

In a more practical setting, this sort of filtering may find a useful role, not only in data analysis but also in online applications, such as speech recognition or active noise cancellation. Indeed, we have already used DEM to infer the hidden states of chaotic systems (hierarchically coupled Lorentz attractors) that were used to simulate bird songs [29]. Applications of Generalised Filtering to similar time-series from chaotic systems may finesse model optimisation in a variety of systems. This is clearly speculative but highlights the potential importance of assimilating data to make inference dynamically, even when the unknown quantities do not change with time.

Appendices

A. Parameter Optimisation, Newton’s Method, and Stability

There is a close connection between the updates implied by (2.9) and Newton’s method for optimisation. Consider the update under a local linearisation [30], assuming that As time proceeds, the change in generalised mean becomes The first line means the motion cancels itself and becomes zero, while the change in the conditional mean becomes a classical Newton update.

An intuition about the need for a high prior precision on the fluctuations of the model parameters can be motivated by a linear stability analysis of the associated recognition dynamics (see (2.9)), with Jacobian in (A.1) For the ensuing dynamics to converge exponentially to the conditional mean (when ) we require . In other words, the prior precision on the motion of parameters should be greater than the (twice the root) conditional precision of the parameters per se. Otherwise, the conditional estimate will exhibit (damped) oscillations.

B. Free-Energy and Action under Mean-Field Approximations

In this appendix, we compare the free-energy and action under the mean-field approximation assumed in variational schemes, where The mean-field approximation finesses evaluation of the free-energy of the accumulated data by making the (implausible) assumption that for discrete samples at . This means that, under the Laplace assumption, Note that (under flat priors on the parameters) the key difference between free-energy and action lies in the contribution from conditional uncertainty (precision) about the parameters. The free-energy contains the log-determinant of the sum of precisions (cf. the precision of the Bayesian parameter average in (2.11)), while the free-action contains the sum of the log-determinant of precisions.

C. Integrating Recognition Dynamics (Filtering)

Filtering involves integrating the ordinary differential equations (2.7) and (2.9) to optimise the conditional means. We can simplify the numerics for hierarchical dynamic models by first collapsing the hierarchy and then collapsing over causal and hidden states This gives a simple form for the Gibbs energy that comprises a log-likelihood and log-prior with the following integration (Generalised Filtering) scheme:

This system can be solved (integrated) using a local linearisation [30] with updates over time steps , where the filter’s Jacobian. Note that we have omitted terms that mediate changes in the motion of state estimates due to changes in parameter estimates. This is because changes in parameter estimates are negligible at the time scale of changes in state estimates. The requisite gradients (evaluated at the conditional expectation) are, with a slight abuse of notion when dealing with derivatives w.r.t. vectors

The corresponding curvatures are (neglecting second-order terms involving states and parameters and second-order derivatives of the conditional entropy)

Finally, the conditional precision and its derivatives are given by the curvature of Gibb’s energy

Note that we have simplified the numerics here by neglecting conditional dependencies between the precisions and the states or parameters. This is easy to motivate because one is not interested in the conditional precision of the precisions but in the (conditional expectation of the) precisions per se; cf. the mean-field approximation implicit in variational approximations. We have also ignored terms due to state-dependent noise, which are not called on in this paper. Finally, we find that the integration of (C.3) is much more stable (in the first few iterations) if we make . Again, this rests on a mean-field like approximation, where we ignore the effects of rapid fluctuations in the states on the entropy of the parameters (but not vice versa).

These equations may look complicated but can be evaluated automatically using numerical derivatives. All the simulations in this paper used just one routine—spm_LAP.m. All the user has to supply are equations defining the generative model. Demonstrations of this scheme are available as part of the SPM software (http://www.fil.ion.ion.ucl.ac.uk.com/spm; DEM_demo.m) and reproduce the examples in the main text.

D. Optimising Smoothness Hyperparameters

When including a hyperparameterisation of the smoothness of the random fluctuations, encoded by the precision matrix on generalised motion , one needs the following derivatives: where is the order of generalised motion. We have not included these derivatives above, because we used assumed values for smoothness in this paper. However, our software implementation of this scheme estimates smoothness by default.

Software Note

The schemes described in this paper are implemented in Matlab code and are available freely http://www.fil.ion.ucl.ac.uk.com/spm. A DEM toolbox provides several demonstrations from a graphical user interface. These demonstrations reproduce the figures of this paper (see spm_LAP.m and ancillary routines).

Acknowledgments

This work was funded by the Wellcome Trust and supported by a grant from the China Scholarship Council (CSC). The author would like to thank Marcia Bennett for help preparing this paper.