Home | About | Journals | Submit | Contact Us | Français |

**|**Front Hum Neurosci**|**v.5; 2011**|**PMC3096853

Formats

Article sections

Authors

Related links

Front Hum Neurosci. 2011; 5: 39.

Published online 2011 May 2. doi: 10.3389/fnhum.2011.00039

PMCID: PMC3096853

Christoph Mathys,^{1,}^{2,}^{*} Jean Daunizeau,^{1,}^{3} Karl J. Friston,^{3} and Klaas E. Stephan^{1,}^{3}

Edited by: Sven Bestmann, University College London, UK

Reviewed by: Quentin Huys, University College London, UK; Robert C. Wilson, Princeton University, USA

*Correspondence: Christoph Mathys, Laboratory for Social and Neural Systems Research, Department of Economics, University of Zurich, Bluemlisalpstrasse 10, 8006 Zurich, Switzerland. e-mail: hc.zhte@syhtamhc

Received 2011 January 1; Accepted 2011 March 29.

Copyright © 2011 Mathys, Daunizeau, Friston and Stephan.

This is an open-access article subject to a non-exclusive license between the authors and Frontiers Media SA, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and other Frontiers conditions are complied with.

This article has been cited by other articles in PMC.

Computational learning models are critical for understanding mechanisms of adaptive behavior. However, the two major current frameworks, reinforcement learning (RL) and Bayesian learning, both have certain limitations. For example, many Bayesian models are agnostic of inter-individual variability and involve complicated integrals, making online learning difficult. Here, we introduce a generic hierarchical Bayesian framework for individual learning under multiple forms of uncertainty (e.g., environmental volatility and perceptual uncertainty). The model assumes Gaussian random walks of states at all but the first level, with the step size determined by the next highest level. The coupling between levels is controlled by parameters that shape the influence of uncertainty on learning in a subject-specific fashion. Using variational Bayes under a mean-field approximation and a novel approximation to the posterior energy function, we derive trial-by-trial update equations which (i) are analytical and extremely efficient, enabling real-time learning, (ii) have a natural interpretation in terms of RL, and (iii) contain parameters representing processes which play a key role in current theories of learning, e.g., precision-weighting of prediction error. These parameters allow for the expression of individual differences in learning and may relate to specific neuromodulatory mechanisms in the brain. Our model is very general: it can deal with both discrete and continuous states and equally accounts for deterministic and probabilistic relations between environmental events and perceptual states (i.e., situations with and without perceptual uncertainty). These properties are illustrated by simulations and analyses of empirical time series. Overall, our framework provides a novel foundation for understanding normal and pathological learning that contextualizes RL within a generic Bayesian scheme and thus connects it to principles of optimality from probability theory.

Learning can be understood as the process of updating an agent's beliefs about the world by integrating new and old information. This enables the agent to exploit past experience and improve predictions about the future; e.g., the consequences of chosen actions. Understanding how biological agents, such as humans or animals, learn requires a specification of both the computational principles and their neurophysiological implementation in the brain. This can be approached in a bottom-up fashion, building a neuronal circuit from neurons and synapses and studying what forms of learning are supported by the ensuing neuronal architecture. Alternatively, one can choose a top-down approach, using generic computational principles to construct generative models of learning and use these to infer on underlying mechanisms (e.g., Daunizeau et al., 2010a). The latter approach is the one that we pursue in this paper.

The laws of inductive inference, prescribing an optimal way to learn from new information, have long been known (Laplace, 1774, 1812). They have a unique mathematical form, i.e., it has been proven that there is no alternative formulation of inductive reasoning that does not violate either consistency or common sense (Cox, 1946). Inductive reasoning is also known as Bayesian learning because the requisite updating of conditional probabilities is described by Bayes’ theorem. Since a Bayesian learner processes information optimally, it should have an evolutionary advantage over other types of agents, and one might therefore expect the human brain to have evolved such that it implements an ideal Bayesian learner (Geisler and Diehl, 2002). Indeed, there is substantial evidence from studies on various domains of learning and perception that human behavior is better described by Bayesian models than by other theories (e.g., Kording and Wolpert, 2004; Bresciani et al., 2006; Yuille and Kersten, 2006; Behrens et al., 2007; Xu and Tenenbaum, 2007; Orbán et al., 2008; den Ouden et al., 2010). However, there remain at least three serious difficulties with the hypothesis that humans act as ideal Bayesian learners. The first problem is that in all but the simplest cases, the application of Bayes’ theorem involves complicated integrals that require burdensome and time-consuming numerical calculations. This makes online learning a challenging task for Bayesian models, and any evolutionary advantage conferred by optimal learning might be outweighed by these computational costs. A second and associated problem is how ideal Bayesian learning, with its requirement to evaluate high-dimensional integrals, would be implemented neuronally (cf. Yang and Shadlen, 2007; Beck et al., 2008; Deneve, 2008). The third difficulty is that Bayesian learning constitutes a normative framework that prescribes how information *should* be dealt with. In reality, though, even when endowed with equal prior knowledge, not all agents process new information alike. Instead, even under carefully controlled conditions, animals and humans display considerable inter-individual variability in learning (e.g., Gluck et al., 2002; Daunizeau et al., 2010b). Despite previous attempts of Bayesian models to deal with individual variability (e.g., Steyvers et al., 2009; Nassar et al., 2010), the failure of orthodox Bayesian learning theory to account for these individual differences remains a key problem for understanding (mal)adaptive behavior of humans. Formal and mechanistic characterizations of this inter-subject variability are needed to comprehend fundamental aspects of brain function and disease. For example, individual differences in learning may result from inter-individual variability in basic physiological mechanisms, such as the neuromodulatory regulation of synaptic plasticity (Thiel et al., 1998), and such differences may explain the heterogeneous nature of psychiatric diseases (Stephan et al., 2009).

These difficulties have been avoided by descriptive approaches to learning, which are not grounded in probability theory, notably some forms of reinforcement learning (RL), where agents learn the “value” of different stimuli and actions (Sutton and Barto, 1998; Dayan and Niv, 2008). While RL is a wide field encompassing a variety of schemes, perhaps the most prototypical and widely used model is that by Rescorla and Wagner (1972). In this description, predictions of value are updated in relation to the current prediction error, weighted by a learning rate (which may differ across individuals and contexts). The great advantage of this scheme is its conceptual simplicity and computational efficiency. It has been applied empirically to various learning tasks (Sutton and Barto, 1998) and has played a major role in attempts to explain electrophysiological and functional magnetic resonance imaging (fMRI) measures of brain activity during reward learning (e.g., Schultz et al., 1997; Montague et al., 2004; O'Doherty et al., 2004; Daw et al., 2006). Furthermore, its non-normative descriptive nature allows for modeling aberrant modes of learning, such as in schizophrenia or depression (Smith et al., 2006; Murray et al., 2007; Frank, 2008; Dayan and Huys, 2009). Similarly, it has found widespread use in modeling the effects of neuromodulatory transmitters, such as dopamine, on learning (e.g., Yu and Dayan, 2005; Pessiglione et al., 2006; Doya, 2008).

Despite these advantages, RL also suffers from major limitations. On the theoretical side, it is a heuristic approach that does not follow from the principles of probability theory. In practical terms, it often performs badly in real-world situations where environmental states and the outcomes of actions are not known to the agent, but must also be inferred or learned. These practical limitations have led some authors to argue that Bayesian principles and “structure learning” are essential in improving RL approaches (Gershman and Niv, 2010). In this article, we introduce a novel model of Bayesian learning that overcomes the three limitations of ideal Bayesian learning discussed above (i.e., computational complexity, questionable biological implementation, and failure to account for individual differences) and that connects Bayesian learning to RL schemes.

Any Bayesian learning scheme relies upon the definition of a so-called “generative model,” i.e., a set of probabilistic assumptions about how sensory signals are generated. The generative model we propose is inspired by the seminal work of Behrens et al. (2007) and comprises a hierarchy of states that evolve in time as Gaussian random walks, with each walk's step size determined by the next highest level of the hierarchy. This model can be inverted (fitted) by an agent using a mean-field approximation and a novel approximation to the conditional probabilities over unknown quantities that replaces the conventional Laplace approximation. This enables us to derive closed-form update equations for the posterior expectations of all hidden states governing contingencies in the environment. This results in extremely efficient computations that allow for real-time learning. The form of these update equations is similar to those of Rescorla–Wagner learning, providing a Bayesian analogon to RL theory. Finally, by introducing parameters that determine the nature of the coupling between the levels of the hierarchical model, the optimality of an update is made conditional upon parameter values that may vary from agent to agent. These parameters encode prior beliefs about higher-order structure in the environment and enable the model to account for inter-individual (and inter-temporal intra-individual) differences in learning. In other words, the model is capable of describing behavior that is subjectively optimal (in relation to the agent's prior beliefs) but objectively maladaptive. Importantly, the model parameters that determine the nature of learning may relate to specific physiological processes, such as the neuromodulation of synaptic plasticity. For example, it has been hypothesized that dopamine, which regulates plasticity of glutamatergic synapses (Gu, 2002), may encode the precision of prediction errors (Friston, 2009). In our model, this precision-weighting of prediction errors is determined by the model's parameters (cf. Figure Figure4).4). Ultimately, our approach may therefore be useful for model-based inference on subject-specific computational and physiological mechanisms of learning, with potential clinical applications for diagnostic classifications of psychiatric spectrum disorders (Stephan et al., 2009).

To prevent any false expectations, we would like to point out that this paper is of a purely theoretical nature. Its purpose is to introduce the theoretical foundations and derivation of our model in detail and convey an understanding of the phenomena it can capture. Owing to length restrictions and to maintain clarity and focus, several important aspects cannot be addressed in this paper. For example, it is beyond the scope of this initial theory paper to investigate the numerical exactness of our variational inversion scheme, present applications of model selection, or present evidence that our model provides for more powerful inference on learning mechanisms from empirical data than other approaches. These topics will be addressed in forthcoming papers by our group.

This paper is structured as follows: First, we derive the general structure of the model, level by level, and consider both its exact and variational inversion. We then present analytical update equations for each level of the model that derive from our variational approach and a quadratic approximation to the variational energies. Following a structural interpretation of the model's update equations in terms of RL, we present simulations in which we demonstrate the model's behavior under different parameter values. Finally, we illustrate the generality of our model by demonstrating that it can equally deal with (i) discrete and continuous environmental states, and (ii) deterministic and probabilistic relations between environmental and perceptual states (i.e., situations with and without perceptual uncertainty).

An overview of our generative model is given by Figures Figures11 and and2.2. To model learning in general terms, we imagine an agent who receives a sequence of sensory inputs *u*^{(1)}, *u*^{(2)},…,*u*^{(n)}. Given a generative model of how the agent's environment generates these inputs, probability theory tells us how the agent can make optimal use of the inputs *u*^{(1)},…,*u*^{(k−1)} and any further “prior” information to predict the next input *u*^{(k)}. The generative model we introduce here is an extension of the model proposed by Daunizeau et al. (2010b) and also draws inspiration from the work by Behrens et al. (2007). Our model is very general: it can deal with states and inputs that are discrete or continuous and uni- or multivariate, and it equally accounts for deterministic and probabilistic relationships between environmental events and perceptual states (i.e., situations with and without perceptual uncertainty). However, to keep the mathematical derivation as simple as possible, we initially deal with an environment where the sensory input *u*^{(k)}{0, 1} on trial *k* is of a binary form; note that this can be readily extended to much more complex environments and input structures. In fact, at a later stage we will also deal with stochastic mappings (i.e., perceptual uncertainty; cf. Eq. 45) and continuous (real-valued) inputs and states (cf. Eq. 48).

For simplicity, imagine a situation where the agent is only interested in a single (binary) state of its environment; e.g., whether it is light or dark. In our model, the environmental state *x*_{1} at time *k*, denoted by ${x}_{1}^{(k)}\{0,1\},$, causes input *u*^{(k)}. Here, ${x}_{1}^{(k)}$ could represent the on/off state of a light switch and *u* the sensation of light or darkness. (For simplicity, we shall often omit the time index *k*). In the following, we assume the following form for the likelihood model:

$$p(u|{x}_{1})={(u)}^{{x}_{1}}{(1-u)}^{1-{x}_{1}}$$

(1)

In other words: *u*=*x*_{1} for both *x*_{1}=1 and *x*_{1}=0 (and *vice versa*). This means that knowing state *x*_{1} allows for an accurate prediction of input *u*: when the switch is on, it is light; when it is off, it is dark. The deterministic nature of this relation does not affect the generality of our argument and will later be replaced by a stochastic mapping when dealing with perceptual uncertainty below.

Since *x*_{1} is binary, its probability distribution can be described by a single real number, the state *x*_{2} at the next level of the hierarchy. We then map *x*_{2} to the probability of *x*_{1} such that *x*_{2}=0 means that *x*_{1}=0 and *x*_{1}=1 are equally probable. For *x*_{2}→∞ the probability for *x*_{1}=1 and *x*_{1}=0 should approach 1 and 0, respectively. Conversely, for *x*_{2}→ the probabilities for *x*_{1}=1 and *x*_{1}=0 should approach 0 and 1, respectively. This can be achieved with the following empirical (conditional) prior density:

$$p({x}_{1}|{x}_{2})=s{({x}_{2})}^{{x}_{1}}{(1-s({x}_{2}))}^{1-{x}_{1}}=\text{Bernoulli(x1;s(x2))}$$

(2)

where *s*(·) is a sigmoid (softmax) function:

$$s(x)\stackrel{\text{def}}{=}\frac{1}{1+\text{exp}(-x)}$$

(3)

Put simply, *x*_{2} is an unbounded real parameter of the probability that *x*_{1}=1. In our light/dark example, one might interpret *x*_{2} as the tendency of the light to be on.

For the sake of generality, we make no assumptions about the probability of *x*_{2} except that it may change with time as a Gaussian random walk. This means that the value of *x*_{2} at time *k* will be normally distributed around its value at the previous time point, ${x}_{2}^{(k-1)}$:

$$p({x}_{2}^{(k)}|{x}_{2}^{(k-1)},{x}_{3}^{(k)})=N({x}_{2}^{(k)};{x}_{2}^{(k-1)},\text{exp(\kappa x3(k)+\omega ))}$$

(4)

Importantly, the dispersion of the random walk (i.e., the variance exp(κ*x*_{3}+ω) of the conditional probability) is determined by the parameters κ and ω (which may differ across agents) as well as by the state *x*_{3}. Here, this state determines the log-volatility of the environment (cf. Behrens et al., 2007, 2008). In other words, the tendency *x*_{2} of the light switch to be on performs a Gaussian random walk with volatility exp(κ*x*_{3}+ω). Introducing ω allows for a volatility that scales independently of the state *x*_{3}. Everything applying to *x*_{2} now equally applies to *x*_{3}, such that we could add as many levels as we please. Here, we stop at the fourth level, and set the volatility of *x*_{3} to *, a constant parameter (which may again differ across agents):*

$$p({x}_{3}^{(k)}|{x}_{3}^{(k-1)},\mathrm{)=N({x}_{3}^{(k)};{x}_{3}^{(k-1)},)}$$

(5)

Given full priors on the parameters, i.e., *p*(κ, ω, ), we can now write the full generative model

$$\begin{array}{l}p({u}^{(k)},{x}_{1}^{(k)},{x}_{2}^{(k)},{x}_{3}^{(k)},{x}_{2}^{(k-1)},{x}_{3}^{(k-1)},\kappa ,\omega ,)\hfill & =p({u}^{(k)}|{x}_{1}^{(k)})p({x}_{1}^{(k)}|{x}_{2}^{(k)})p({x}_{2}^{(k)}|{x}_{2}^{(k-1)},{x}_{3}^{(k)},\kappa ,\omega )\hfill \\ \mathrm{\mathrm{p({x}_{3}^{(k)}|{x}_{3}^{(k-1)},\mathrm{)p({x}_{2}^{(k-1)},{x}_{3}^{(k-1)})p(\kappa ,\omega ,)}}}\hfill \end{array}$$

(6)

Given priors on the initial state $p({x}_{2}^{(0)},{x}_{3}^{(0)})$ the generative model is defined for all times *k* by recursion to *k*=1. Inverting this model corresponds to optimizing the posterior densities over the unknown (hidden) states *x*={*x*_{1}, *x*_{2}, *x*_{3}} and parameters χ={κ, ω, }. This corresponds to perceptual inference and learning, respectively. In the next section, we consider the nature of this inversion or optimization.

It is instructive to consider the factorization of the generative density

$$\begin{array}{l}p({u}^{(k)},{x}_{1}^{(k)},{x}_{2}^{(k)},{x}_{3}^{(k)},{x}_{2}^{(k-1)},{x}_{3}^{(k-1)},\chi )\hfill \\ =p({u}^{(k)},{x}_{1}^{(k)},{x}_{2}^{(k)},{x}_{3}^{(k)},\chi |{x}_{2}^{(k-1)},{x}_{3}^{(k-1)})p({x}_{2}^{(k-1)},{x}_{3}^{(k-1)})\hfill \end{array}$$

(7)

In this form, the Markovian structure of the model becomes apparent: the joint probability of the input and the states at time *k* depends only on the states at the immediately preceding time *k*−1. It is the probability distribution of these states that contains the information conveyed by previous inputs ${u}^{(1\dots k-1)}\stackrel{\text{def}}{=}({u}^{(1)},\dots ,{u}^{(k-1)})$; i.e.:

$$p({x}_{2}^{(k-1)},{x}_{3}^{(k-1)})=p({x}_{2}^{(k-1)},{x}_{3}^{(k-1)}|{u}^{(\mathrm{1...}k-1)})$$

(8)

By integrating out ${x}_{2}^{(k-1)}$ and ${x}_{3}^{(k-1)}$ we obtain the following compact form of the generative model at time *k*:

$$\begin{array}{l}\underset{-\infty}{\overset{\infty}{\int}}\underset{-\infty}{\overset{\infty}{\int}}p({u}^{k},\text{x1(k),x2(k),x3(k),\chi |x2(k\u22121),x3(k\u22121))p(x2(k\u22121),x3(k\u22121)|u(1\u2026k\u22121))dx2k\u22121dx2(k\u22121)=p(u(k),x(k),\chi |u(1\u2026k\u22121))}\hfill \end{array}$$

(9)

Once *u*^{(k)} is observed, we can plug it into this expression and obtain

$$p({x}^{(k)},\chi |{u}^{(\mathrm{1...}k)})$$

(10a)

This is the quantity of interest to us, because it describes the posterior probability of the time-dependent states, *x*^{(k)} (and time-independent parameters, χ), in the agent's environment. This is what the agent infers (and learns), given the history of previous inputs. Computing this probability is called model inversion: unlike the likelihood $p({u}^{(k)}|{x}^{(k)},\text{\chi ,u(1\u2026k\u22121))}$ the posterior does not predict data (*u*^{(k)}) from hidden states and parameters but predicts states and parameters from data.

In the framework we introduce here, we model the individual variability between agents by putting delta-function priors on the parameters:

$$\begin{array}{l}p({x}^{(k)},\chi |{u}^{(\mathrm{1...}k)})=p({x}^{(k)}|\chi ,{u}^{(\mathrm{1...}k)})p(\chi |{u}^{(\mathrm{1...}k)})\hfill \\ p(\chi |{u}^{(\mathrm{1...}k)})=\delta (\chi -{\chi}_{a})\hfill \end{array}$$

(10b)

where χ* _{a}* are the fixed parameter values that characterize a particular agent at a particular time (e.g., during the experimental session). This corresponds to the common distinction between states (as variables that change quickly) and parameters (as variables that change slowly or not at all). In other words, we assume that the timescale at which parameter estimates change is much larger than the one on which state estimates change, and also larger than the periods during which we observe agents in a single experiment. This is not a strong assumption, given that the model has multiple hierarchical levels of states that give the agent the required flexibility to adapt its beliefs to a changing environment. In effect, this approach gives us a family of Bayesian learners whose (slowly changing) individuality is captured by their priors on the parameters. For other examples where subjects’ beliefs about the nature of their environment were modeled as priors on parameters, see Daunizeau et al. (2010b) and Steyvers et al. (2009).

In principle, model inversion can proceed in an online fashion: By (numerical) marginalization, we can obtain the (marginal) posteriors $p({x}_{2}^{(k)}|{u}^{(1\dots k)})$ and $p({x}_{3}^{(k)}|{u}^{(1\dots k)})$; this is the approach adopted by Behrens et al. (2007), allowing one to compute *p*(*u*^{(k+1)}, *x*^{(k+1)}, χ|*u*^{(1…k)}) according to Eq. 6, and subsequently *p*(*u*^{(k+1)}, χ|*u*^{(1…k+1)}) once *u*^{(k+1)} becomes known, and so on. Unfortunately, this (exact) inversion involves many complicated (non-analytical) integrals for every new input, rendering exact Bayesian inversion unsuitable for real-time learning in a biological setting. If the brain uses a Bayesian scheme, it is likely that it relies on some sufficiently accurate, but fast, approximation to Eq. 9: this is approximate Bayesian inference. As described in the next section, a generic and efficient approach is to employ a mean-field approximation within a variational scheme. This furnishes an efficient solution with biological plausibility and interpretability.

Variational Bayesian (VB) inversion determines the posterior distributions $p({x}^{(k)},\text{\chi |u1\u2026k)}$ by maximizing the log-model evidence. The log-evidence corresponds to the negative surprise about the data, given a model, and is approximated by a lower bound, the negative free energy. Detailed treatments of the general principles of the VB procedure can be found in numerous papers (e.g., Beal, 2003; Friston and Stephan, 2007). The approximations inherent in VB enable a computationally efficient inversion scheme with closed-form single-step probability updates from trial to trial. In particular, VB can incorporate the so-called *mean-field approximation* which turns the joint posterior distribution into the product of approximate marginal posterior distributions:

$$\begin{array}{l}p({x}^{(k)},\text{\chi |u(1\u2026k))=p(x(k)|\chi ,u(1\u2026k))p(\chi |u1\u2026k)}p({x}^{(k)}|\chi ,\text{u(1\u2026k))\u2248i=1nq^(xi(k))}q({x}_{i}^{(k)})\approx \widehat{q}({x}_{i}^{(k)})\hfill \hfill \hfill \end{array}$$

(11)

Based on this assumption, the variational maximization of the negative free energy is implemented in a series of variational updates for each level *i* of the model separately. The second line in Eq. 11 represents the mean-field assumption (factorization of the posterior), while the third line reflects the fact that we assume a fixed form *q*(·) for the approximate marginals $\widehat{q}()$. We make minimal assumptions about the form of the approximate posteriors by following the maximum entropy principle: given knowledge of, or assumptions about, constraints on a distribution, the least arbitrary choice of distribution is the one that maximizes entropy (Jaynes, 1957). To keep the description of the posteriors simple and biologically plausible, we take them to be characterized only by their first two moments; i.e., by their mean and variance. At the first level, we have a binary state *x*_{1} with a mean *μ*_{1}=*p*(*x*_{1}=1). Under this constraint, the maximum entropy distribution is the Bernoulli distribution with parameter *μ*_{1} (where the variance *μ*_{1}(1 − *μ*_{1}) is a function of the mean):

$$\begin{array}{l}p({x}_{1}^{(k)}|\chi ,{u}^{(\mathrm{1...}k)})\approx q({x}_{1}^{(k)})\hfill \\ \text{=Bernoulli(x1(k);\mu 1(k))=(\mu 1(k))x1(k)(1\u2212\mu 1(k))1\u2212x1(k)}\hfill \end{array}$$

(12)

At the second and third level, the maximum entropy distribution of the unbounded real variables *x*_{2} and *x*_{3}, given their means and variances, is Gaussian. Note that the choice of a Gaussian distribution for the approximate posteriors is not due simply to computational expediency (or the law of large numbers) but follows from the fact that, given the assumption that the posterior is encoded by its first two moments, the maximum entropy principle prescribes a Gaussian distribution. Labeling the means *μ*_{2}, *μ*_{3} and the variances *σ*_{2}, *σ*_{3}, we obtain

$$p({x}_{2}^{(k)}|\chi ,{u}^{(\mathrm{1...}k)})\approx q({x}_{2}^{(k)})=\mathcal{N}({x}_{2}^{(k)};{\mu}_{2}^{(k)},{\sigma}_{2}^{(k)})$$

(13)

$$p({x}_{3}^{(k)}|\chi ,{u}^{(\mathrm{1...}k)})\approx q({x}_{3}^{(k)})=\mathcal{N}({x}_{3}^{(k)};{\mu}_{3}^{(k)},{\sigma}_{3}^{(k)})$$

(14)

Now that we have approximate posteriors *q* that we treat as known for all but the *i*th level, the next step is to determine the variational posterior $\widehat{q}({x}_{i})$ for this level *i*. Variational calculus shows that given *q*(*x _{j}*) (

$$\widehat{q}({x}_{i}^{(k)})=\frac{1}{{Z}_{i}}\text{exp}(I({x}_{i}^{(k)}))$$

(15)

${\mathcal{Z}}_{i}$ is a normalization constant that ensures that the integral (or sum, in the discrete case) of $\widehat{q}$ over *x* equals unity. Under our generative model, the variational energy is

$$\begin{array}{c}I({x}_{i}^{(k)})=\underset{{X}_{\backslash i}}{\int}q({x}_{\backslash i}^{(k)})\text{lnp(u(k),xi(k),xi(k),\chi |u(1...k\u22121))dxi(k)}& q({x}_{\backslash i})={\displaystyle \underset{q({x}_{j})}{j\ne i}}\end{array}$$

(16)

where *x _{\i}* denotes all

To compute $\widehat{q}({x}_{i}^{(k)})$, we need $q({x}_{\backslash i}^{(k)})$ and therefore the sufficient statistics λ_{\i}={*μ*_{\i},σ_{\i}} for the posteriors at all but the *i*th level. One could try to extract them from $\widehat{q}({x}_{\backslash i}^{(k)})$, but that would constitute a circular problem. We avoid this by exploiting the hierarchical form of the model: for the first level, we use the sufficient statistics of the higher levels from the previous time point *k*−1, since information about input *u*^{(k)} cannot yet have reached those levels. From there we proceed upward through the hierarchy of levels, always using the updated parameters ${\lambda}_{\backslash i}^{(k)}$ for levels lower than the current level and pre-update values ${\lambda}_{\backslash i}^{(k-1)}$ for higher levels. Extending the approach suggested by Daunizeau et al. (2010b), and using power series approximations where necessary, the variational energy integrals can then be calculated for all *x _{i}*, giving

$$I({x}_{1}^{(k)})=\text{ln}\left(\frac{{u}^{(k)}s({\mu}_{2}^{(k-1)})}{(1-{u}^{(k)})(1-s({\mu}_{2}^{(k-1)}))}\right){x}_{1}^{(k)}$$

(17)

$$I({x}_{2}^{(k)})=\text{ln}s({x}_{2}^{(k)})+{x}_{2}^{(k)}({\mu}_{1}^{(k)}-1)-\frac{1}{2({\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega})}{({x}_{2}^{(k)}-{\mu}_{2}^{(k-1)})}^{2}$$

(18)

$$\begin{array}{l}I({x}_{3}^{(k)})=-\frac{1}{2}\text{ln}({\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {x}_{3}^{(k)}+\omega})-\frac{1}{2}\frac{{\sigma}_{2}^{(k)}+{({\mu}_{2}^{(k)}-{\mu}_{2}^{(k-1)})}^{2}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {x}_{3}^{(k)}+\omega}}\hfill \\ \text{\u221212(\sigma 3(k\u22121)+)(x3(k)\u2212\mu 3(k\u22121))2}\hfill \end{array}$$

(19)

Substituting these variational energies into Eq. 15, we obtain the posterior distribution of *x* under the mean-field approximation.

According to Eqs 15 and 17–19, $\widehat{q}({x}_{i}^{(k)})$ depends on ${u}^{(k)},\text{\lambda i(k\u22121),\lambda i(k\u22121),\lambda i(k)}$, and χ. In the next section, we show that it is possible to derive simple closed-form Markovian update equations of the form

$${\lambda}_{i}^{(k)}={f}_{i}({u}^{(k)},{\lambda}_{i}^{(k-1)},{\lambda}_{\backslash i}^{(k-1)},{\lambda}_{\backslash i}^{(k)},\chi )$$

(20)

Update equations of this form allow the agent to update its approximate posteriors over ${x}_{i}^{(k)}$ very efficiently and thus optimize its beliefs about the environment in real time. We now consider the detailed form of these equations and how they relate to classical heuristics from RL.

At the first level of the model, it is simple to determine *q*(*x*_{1}) since $\widehat{q}({x}_{1})=1/{\mathcal{Z}}_{1}\text{exp}(I({x}_{1}))$ is a Bernoulli distribution with parameter

$${\mu}_{1}^{(k)}=\widehat{q}({x}_{1}^{(k)}=1)={u}^{(k)}$$

(21)

and therefore already has the form required of *q*(*x*_{1}) by Eq. 12. We can thus take $q({x}_{1})=\widehat{q}({x}_{1})$ and have in Eq. 21 an update rule of the form of Eq. 20.

At the second level, $\widehat{q}({x}_{2})$ does not have the form required of *q*(*x*_{2}) by Eq. 13 since it is only approximately Gaussian. It is proportional to the exponential of *I*(*x*_{2}) and would only be Gaussian if *I*(*x*_{2}) were quadratic. The problem of finding a Gaussian approximation $q({x}_{2})\approx \widehat{q}({x}_{2})$ can therefore be reformulated as finding a quadratic approximation to *I*(*x*_{2}). The obvious way to achieve this is to expand *I*(*x*_{2}) in powers of *x*_{2} up to second order. The choice of expansion point, however, is not trivial. One possible choice is the mode or maximum of *I*(*x*_{2}), resulting in the frequently used *Laplace approximation* (Friston et al., 2007). This has the disadvantage that the maximum of *I*(*x*_{2}) is unknown and has to be found by numerical optimization methods, precluding a single-step analytical update rule of the form of Eq. 20 (this is no problem in a continuous time setting, where the mode of the variational energy can be updated continuously (cf. Friston, 2008)). However, for the discrete updates we seek, one can use the expectation ${\mu}_{2}^{(k-1)}$ as the expansion point for time *k*, when the agent receives input *u*^{(k)} and the expectation of ${x}_{2}^{(k)}$ is ${\mu}_{2}^{(k-1)}$. In terms of computational and mnemonic costs, this is the most economical choice since this value is known. Moreover, it yields analytical update equations which (i) bear structural resemblance to those used by RL models (see Structural Interpretation of the Update Equations) and (ii) can be computed very efficiently in a single step:

$${\sigma}_{2}^{(k)}=\frac{1}{1/{\widehat{\sigma}}_{2}^{(k)}+{\widehat{\sigma}}_{1}^{(k)}}$$

(22)

$${\mu}_{2}^{(k)}={\mu}_{2}^{(k-1)}+{\sigma}_{2}^{(k)}{\delta}_{1}^{(k)}$$

(23)

where, for clarity, we have used the definitions

$${\widehat{\mu}}_{1}^{(k)}\stackrel{\text{def}}{=}s({\mu}_{2}^{(k-1)})$$

(24)

$${\delta}_{1}^{(k)}\stackrel{\text{def}}{=}{\mu}_{1}^{(k)}-{\widehat{\mu}}_{1}^{(k)}$$

(25)

$${\widehat{\sigma}}_{1}^{(k)}\stackrel{\text{def}}{=}{\widehat{\mu}}_{1}^{(k-1)}(1-{\widehat{\mu}}_{1}^{(k-1)})$$

(26)

$${\widehat{\sigma}}_{2}^{(k)}\stackrel{\text{def}}{=}{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}$$

(27)

formulated in terms of precisions (inverse variances) ${\pi}_{2}^{(k)}\stackrel{\text{def}}{=}1/{\sigma}_{2}^{(k)}$, ${\widehat{\pi}}_{1}^{(k)}\stackrel{\text{def}}{=}1/{\widehat{\sigma}}_{1}^{(k)}$, ${\widehat{\pi}}_{2}^{(k)}\stackrel{\text{def}}{=}1/{\widehat{\sigma}}_{2}^{(k)}$, the variance update (Eq. 22) takes the simple form

$${\pi}_{2}^{(k)}={\widehat{\pi}}_{2}^{(k)}+\frac{1}{{\widehat{\pi}}_{1}^{(k)}}$$

(28)

in the context of the update equations, we use the hat notation to indicate “referring to prediction.” While the ${\widehat{\mu}}_{i}^{(k)}$ are the predictions before seeing the input *u*^{(k)}, the ${\widehat{\sigma}}_{i}^{(k)}$ and ${\widehat{\pi}}_{i}^{(k)}$ are the variances (i.e., uncertainties) and precisions of these predictions (see Structural Interpretation of the Update Equations).

The approach that produces these update equations is conceptually similar to a Gauss–Newton ascent on the variational energy that would, by iteration, produce the Laplace approximation (cf. Figure Figure3).3). The same approach can be taken at the third level, where we also have a Gaussian approximate posterior:

$${\pi}_{3}^{(k)}={\widehat{\pi}}_{3}^{(k)}+\frac{{\kappa}^{2}}{2}{w}_{2}^{(k)}({w}_{2}^{(k)}+{r}_{2}^{(k)}\text{\delta 2(k))}$$

(29)

$${\mu}_{3}^{(k)}={\mu}_{3}^{(k-1)}+{\sigma}_{3}^{(k)}\frac{\kappa}{2}{w}_{2}^{(k)}{\delta}_{2}^{(k)}$$

(30)

with ${\pi}_{3}^{(k)}\stackrel{\text{def}}{=}1/{\sigma}_{3}^{(k)}$ and

$${\widehat{\pi}}_{3}^{(k)}\stackrel{\text{def}}{=}\frac{1}{{\sigma}_{3}^{(k-1)}+}$$

(31)

$${w}_{2}^{(k)}\stackrel{\text{def}}{=}\frac{{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}$$

(32)

$${r}_{2}^{(k)}\stackrel{\text{def}}{=}\frac{{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}-{\sigma}_{2}^{(k-1)}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}$$

(33)

$${\delta}_{2}^{(k)}\stackrel{\text{def}}{=}\frac{{\delta}_{2}^{(k)}+{({\mu}_{2}^{(k)}-{\mu}_{2}^{(k-1)})}^{2}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}-1$$

(34)

The derivation of these update equations, based on this novel quadratic approximation to *I*(*x*_{2}) and *I*(*x*_{3}), is described in detail in the next section, and Figure Figure33 provides a graphical illustration of the ensuing procedure. The meaning of the terms that appear in the update equations and the overall structure of the updates will be discussed in detail in the Section “Structural Interpretation of the Update Equations.”

While knowing *I*(*x _{i}*) (with

$$\begin{array}{l}p({x}_{i}^{(k)}|{u}^{(\mathrm{1...}k)},\chi )\stackrel{\text{meanfield}}{\to}\widehat{q}({x}_{i}^{(k)})\text{exp}I({x}_{i}^{(k)})\hfill & \text{\u2192Gaussianq(xi(k))expI\u02dc(xi(k))}\hfill \end{array}$$

(35)

We denote by Ĩ the quadratic function obtained by expansion of *I* around ${x}_{i}^{(k)}={\mu}_{i}^{(k-1)}$ (see Figure Figure33 for a graphical summary). Its exponential has the Gaussian form required by Eqs 13 and 14 (where ${\tilde{\mathcal{Z}}}_{i}$ is a normalization constant):

$$q\left({x}_{i}^{\left(k\right)}\right)=\frac{1}{\sqrt{2\pi {\sigma}_{i}^{\left(k\right)}}}\mathrm{exp}\left(-\frac{{\left({x}_{i}^{\left(k\right)}-{\mu}_{i}^{\left(k\right)}\right)}^{2}}{2{\sigma}_{i}^{\left(k\right)}}\right)=\frac{1}{{\tilde{\mathcal{Z}}}_{i}}\mathrm{exp}\left(\tilde{I}\left({x}_{i}^{\left(k\right)}\right)\right)$$

(36)

This equation lets us find ${\mu}_{i}^{(k)}$ and ${\sigma}_{i}^{(k)}$ Taking the logarithm on both sides and then differentiating twice with respect to ${x}_{i}^{(k)}$ gives

$${\sigma}_{i}^{(k)}=-\frac{1}{{2}^{\tilde{I}}}$$

(37)

where ^{2} denotes the second derivative, which is constant for a quadratic function. Because ^{2}*I* and ^{2}Ĩ agree at the expansion point ${x}_{i}^{(k)}={\mu}_{i}^{(k-1)}$, we may write

$${\sigma}_{i}^{(k)}=-\frac{1}{{2}^{I}}$$

(38)

A somewhat different line of reasoning leads to ${\mu}_{i}^{(k)}$. Since ${\mu}_{i}^{(k)}$ is the argument of the maximum (argmax) of *q* (and exponentiation preserves the argmax) ${\mu}_{i}^{(k)}$ has to be the argmax of *Ĩ*. Starting at any point ${x}_{i}^{(k)}$ the exact argmax of a quadratic function can be found in one step by Newton's method:

$${\mu}_{i}^{(k)}=\text{arg maxI\u02dc(xi(k))=xi(k)\u2212I\u02dc(xi(k))2I\u02dc(xi(k))}$$

(39)

If we choose ${x}_{i}^{(k)}$ to be the expansion point ${\mu}_{i}^{(k-1)}$, we have agreement of *I* with *Ĩ* up to the second derivative at this point and may therefore write

$${\mu}_{i}^{(k)}={\mu}_{i}^{(k-1)}-\frac{I({\mu}_{i}^{(k-1)}){2}^{I}={\mu}_{i}^{(k-1)}+{\sigma}_{i}^{(k)}I({\mu}_{i}^{(k-1)})}{}$$

(40)

Plugging *I*(*x*_{2}) from Eq. 18 into Eqs 38 and 40 now yields the parameter update equations (Eqs 22 and 23) of the form Eq. 20 for the second level. For the third level (or indeed any higher level with an approximately Gaussian posterior distribution we may wish to include), the same procedure gives Eqs 29 and 30. Note that this method of obtaining closed-form Markovian parameter update equations can readily be applied to multidimensional *x _{i}'s* with approximately multivariate Gaussian posteriors by reinterpreting

As we have seen, variational inversion of our model, using a new quadratic approximation to the variational energies, provides a set of simple trial-by-trial update rules for the sufficient statistics λ* _{i}*={

Crucially, the update equations for *μ*_{2} and *μ*_{3} have a form that is familiar from RL models such as Rescorla–Wagner learning (Figure (Figure4).4). The general structure of RL models can be summarized as:

$${\text{prediction}}^{(k)}={\text{prediction}}^{(k-1)}+\text{learningrate}\times \text{predictionerror}$$

As we explain in detail below, this same structure appears in Eqs 23 and 30–displayed here in their full forms:

$$\begin{array}{l}\underset{{\text{prediction}}^{(k)}}{\underbrace{{\mu}_{2}^{(k)}}}=\underset{{\text{prediction}}^{(k-1)}}{\underbrace{{\mu}_{2}^{(k-1)}}}+\underset{\text{learningrate}}{\underbrace{{\sigma}_{2}^{(k)}}}\underset{\text{predictionerror}}{\underbrace{\left({\mu}_{1}^{(k)}-s\left({\mu}_{2}^{(k-1)}\right)\right)}}\\ \underset{{\text{prediction}}^{(k)}}{\underbrace{{\mu}_{3}^{(k)}}}=\underset{{\text{prediction}}^{(k-1)}}{\underbrace{{\mu}_{3}^{(k-1)}}}+\underset{\text{learningrate}}{\underbrace{{\sigma}_{3}^{(k)}\frac{\kappa}{2}\frac{{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{k{m}_{3}^{(k-1)}+v}}}}\\ \text{}\times \underset{\text{predictionerror}}{\underbrace{\left(\frac{{\sigma}_{2}^{(k)}+{\left({\mu}_{2}^{(k)}-{\mu}_{2}^{(k-1)}\right)}^{2}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}-1\right)}}\end{array}$$

The term ${\delta}_{1}^{(k)}={\mu}_{1}^{(k)}-s({\mu}_{2}^{(k-1)})$ in Eq. 23 corresponds to the prediction error at the first level. This prediction error is the difference between the expectation ${\mu}_{1}^{(k)}$ of *x*_{1} having observed input *u*^{(k)} and the prediction ${\widehat{\mu}}_{1}^{(k)}=s({\mu}_{2}^{(k-1)})$ before receiving *u*^{(k)}; i.e., the softmax transformation of the expectation of *x*_{2} before seeing *u*^{(k)}. Furthermore, ${\sigma}_{2}^{(k)}$ in Eq. 23 can be interpreted as the equivalent of a (time-varying) learning rate in RL models (cf. Preuschoff and Bossaerts, 2007). Since σ_{2} represents the width of *x*_{2}’s posterior and thus the degree of our uncertainty about *x*_{2}, it makes sense that updates in *μ*_{2} are proportional to this estimate of posterior uncertainty: the less confident the agent is about what it knows, the greater the influence of new information should be.

According to its update Eq. 22, *μ*_{2} always remains positive since it contains only positive terms. Crucially, *σ*_{2}, through ${\widehat{\sigma}}_{2}$, depends on the log-volatility estimate *μ*_{3}, from the third level of the model. For vanishing volatility, i.e., ${\text{e}}^{\kappa {\mu}_{3}+\omega}\approx 0$
*σ*_{2} can only decrease from trial to trial. This corresponds to the case in which the agent believes that *x*_{2} is fixed; the information from every trial then has the same weight and new information can only shrink *σ*_{2}. On the other hand, even with ${\text{e}}^{\kappa {\mu}_{3}+\omega}\approx 0$, *σ*_{2} has a lower bound: when *σ*_{2} approaches zero, the denominator of Eq. 22 approaches 1/*σ*_{2} from above, leading to ever smaller decreases in *σ*_{2}. This means that after a long train of inputs, the agent still learns from new input, even when it infers that the environment is stable.

The precision formulation (cf. Eq. 28)

$${\pi}_{2}^{(k)}={\widehat{\pi}}_{2}^{(k)}+\frac{1}{{\widehat{\pi}}_{1}^{(k)}}=\frac{1}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}+{\widehat{\sigma}}_{1}^{(k)}$$

(41)

illustrates that three forms of uncertainty influence the posterior variance ${\sigma}_{2}^{(k)}$: the informational $({\sigma}_{2}^{(k-1)})$ and the environmental $({\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega})$ uncertainty at the second level (see discussion in the next paragraph), and the uncertainty ${\widehat{\sigma}}_{1}^{(k)}$ of the prediction at the first level. While environmental uncertainty at the second level thus decreases precision ${\pi}_{2}^{(k)}$ relative to its previous value ${\pi}_{2}^{(k-1)}=1/{\sigma}_{2}^{(k-1)}$, predictive uncertainty at the first level counteracts that decrease, i.e., it keeps the learning rate ${\sigma}_{2}^{(k)}$ smaller than it would otherwise be. This makes sense because prediction error should mean less when predictions are more uncertain.

The update rule (Eq. 30) for *μ*_{3} has a similar structure to that of *μ*_{2} (Eq. 23) and can also be interpreted in terms of RL. Although perhaps not obvious at first glance, ${\delta}_{2}^{(k)}$ (Eq. 34) represents prediction error. It is positive if the updates at the second level (of *μ*_{2} and *σ*_{2}) in response to input *u*^{(k)} indicate that the agent was underestimating *x*_{3}. Conversely, it is negative if the agent was overestimating *x*_{3}. This can be seen by noting that the uncertainty about *x*_{2} has two sources: *informational*, i.e., the lack of knowledge about *x*_{2} (represented by *σ*_{2}), and *environmental*, i.e., the volatility of the environment (represented by ${\text{e}}^{\kappa {\mu}_{3}+\omega}$). Before receiving input *u*^{(k)} the total uncertainty is ${\widehat{\sigma}}_{2}^{(k)}={\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}$. After receiving the input, the updated total uncertainty is ${\sigma}_{2}^{(k)}+{({\mu}_{2}^{(k)}-{\mu}_{2}^{(k-1)})}^{2}$ where *σ*_{2} has been updated according to Eq. 22 and ${\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}$ has been replaced by the squared update of *μ*_{2}. If the total uncertainty is greater after seeing *u*^{(k)}, the fraction in ${\delta}_{2}^{(k)}$ is greater than one, and *μ*_{3} increases. Conversely, if seeing *u*^{(k)} reduces total uncertainty, *μ*_{3} decreases. (Since *x*_{3} is on a logarithmic scale with respect to *x*_{2}, the ratio and not the difference of quantities referring to *x*_{2} is relevant for the prediction error in *x*_{3}). It is important to note that we did not construct the update equations with any of these properties in mind. It is simply a reflection of Bayes optimality that emerges on applying our variational update method.

The term corresponding to the learning rate of *μ*_{3} is

$${\sigma}_{3}^{(k)}\frac{\kappa}{2}{w}_{2}^{(k)}={\sigma}_{3}^{(k)}\frac{\kappa}{2}\frac{{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}$$

(42)

As at the second level, this is proportional to the variance *σ*_{3}of the posterior. But here, the learning rate is also is proportional to the parameter κ and a *weighting factor*
${w}_{2}^{(k)}$ for prediction error. κ determines the form of the Gaussian random walk at the second level and couples the third level to the second (cf. Eqs 4 and 30); ${w}_{2}^{(k)}$ is a measure of the (environmental) volatility ${\text{e}}^{\kappa {\mu}_{3}+\omega}$ of *x*_{2} relative to its (informational) conditional uncertainty, *σ*_{2}. It is bounded between 0 and 1 and approaches 0 as ${\text{e}}^{\kappa {\mu}_{3}+\omega}$ becomes negligibly small relative to *σ*_{2}; conversely, it approaches 1 as *σ*_{2} becomes negligibly small relative to ${\text{e}}^{\kappa {\mu}_{3}+\omega}$. This means that lack of knowledge about *x*_{2} (i.e., conditional uncertainty *σ*_{2}) suppresses updates of *μ*_{3} by reducing the learning rate, reflecting the fact that prediction errors in *x*_{2} are only informative if the agent is confident about its predictions of *x*_{2}. As with the prediction error term, this weighting factor emerged from our variational approximation.

The precision update (Eq. 29) at third level also has an interpretable form. In addition to ${\delta}_{2}^{(k)}$ and ${w}_{2}^{(k)}$, we now also have the term ${r}_{2}^{(k)}$ (Eq. 33), which is the *relative difference* of environmental and informational uncertainty (i.e., relative to their sum). Note that it is a simple affine function of the weighting factor ${w}_{2}^{(k)}$

$${r}_{2}^{(k)}=2{w}_{2}^{(k)}-1$$

(43)

As at the second level, the precision update is the sum of the precision ${\widehat{\pi}}_{3}^{(k)}$ of the prediction, reflecting the informational and environmental uncertainty at the third level, and the term

$$\frac{{\kappa}^{2}}{2}{w}_{2}^{(k)}({w}_{2}^{(k)}+{r}_{2}^{(k)}{\delta}_{2}^{(k)})$$

(44)

Proportionality to *κ*^{2} reflects the fact that stronger coupling between the second and third levels leads to higher posterior precision (i.e., less posterior uncertainty) at the third level, while proportionality to *w*_{2} depresses precision at the third level when informational uncertainty at the second level is high relative to environmental uncertainty; the latter also applies to the first summand in the brackets. The second summand *r*_{2}δ_{2} means that, when the agent regards environmental uncertainty at the second level as relatively high (*r*_{2}>0), volatility is held back from rising further if δ_{2}>0 by way of a decrease in the learning rate (which is proportional to the inverse precision), but conversely pushed to fall if δ_{2}<0. If, however, environmental uncertainty is relatively low (*r*_{2}<0), the opposite applies: positive volatility prediction errors increase the learning rate, allowing the environmental uncertainty to rise more easily, while negative prediction errors decrease the learning rate. The term *r*_{2}δ_{2} thus exerts a stabilizing influence on the estimate of *μ*_{3}.

This automatic integration of all the information relevant to a situation is typical of Bayesian methods and brings to mind a remark made by Jaynes (2003, p. 517) in a different context: “This is still another example where Bayes’ theorem detects a genuinely complicated situation and automatically corrects for it, but in such a slick, efficient way that one is [at first, we would say] unaware of what is happening.” In the next section we use simulations to illustrate the nature of this inference and demonstrate some of its more important properties.

Here, we present several simulations to illustrate the behavior of the update equations under different values of the parameters , *ω*, and *κ*. Figure Figure55 depicts a “reference” scenario, which will serve as the basis for subsequent variations. Figures Figures66–8 display the effects of selectively changing one of the parameters , *ω*, and *κ*, leading to distinctly different types of inference.

The reference scenario in Figure Figure55 (and Figure Figure9,9, top row) illustrates some basic features of the model and its update rules. For this reference, we chose the following parameter values: =0.5, *ω*=−2.2, and *κ*=1.4. Overall, the agent is exposed to 320 sensory outcomes (stimuli) that are administered in three stages. In a first stage, it is exposed to 100 trials where the probability that *x*_{1}=1 is 0.5. The posterior expectation of *x*_{1} accordingly fluctuates around 0.5 and that of *x*_{2} around 0; the expected volatility remains relatively stable. There then follows a second period of 120 trials with higher volatility, where the probability that *x*_{1}=1 alternates between 0.9 and 0.1 every 20 trials. After each change, the estimate of *x*_{1} reliably approaches the true value within about 20 trials. In accordance with the changes in probability, the expected outcome tendency *x*_{2} now fluctuates more widely around zero. At the third level, the expected log-volatility *x*_{3} shows a tendency to rise throughout this period, displaying upward jumps whenever the probability of an outcome changes (and thus *x*_{2} experiences sudden updates). As would be anticipated, the expected log-volatility declines during periods of stable outcome probability. In a third and final period, the first 100 trials are repeated in exactly the same order. Note how owing to the higher estimate of volatility (i.e., greater ${\text{e}}^{\kappa {\mu}_{3}+\omega}$), the learning rate has increased, now causing the same sequence of inputs to have a greater effect on *μ*_{2} than during the first stage of the simulation. As expected, a more volatile environment leads to a higher learning rate.

One may wonder why, in the third stage of the simulation, the expected log-volatility *μ*_{3} continues to rise even after the true *x*_{2} has returned to a stable value of 0 (corresponding to *p*(*x*_{1}=1)=0.5; see the fine black line in Figure Figure5).5). This is because a series of three *x*_{1}=1 outcomes, followed by three *x*_{1}=0 could just as well reflect a stable *p*(*x*_{1}=1)=0.5 or a jump from *p*(*x*_{1}=1)=1 to *p*(*x*_{1}=1)=0 after the first three trials. Depending on the particular choice of parameters , *ω*, and *κ*, the agent shows qualitatively different updating behavior: Under the parameters in the reference scenario, it has a strong tendency to increase its posterior expectation of volatility in response to unexpected stimuli. For other parameterizations (as in the scenarios described below), this is not the case. Importantly, this ambiguity disappears once the model is inverted by fitting it to behavioral data (see Discussion).

The nature of the simulation in Figure Figure55 is not only determined by the choice of values for the parameters , *ω*, and *κ*, but also by initial values for *μ*_{2}, *σ*_{2}, *μ*_{3}, and *σ*_{3} (the initial value of *μ*_{1} is *s*(*μ*_{2})). Any change in the initial value ${\mu}_{3}^{(0)}$ of *μ*_{3} can be neutralized by corresponding changes in *κ* and *ω*. We may therefore assume ${\mu}_{3}^{(0)}=0$ without loss of generality but remembering that *κ* and *ω* are only unique relative to this choice. However, the initial values of *μ*_{2}, *σ*_{2}, and *σ*_{3} are, in principle, not neutral. They are nevertheless of little consequence in practice since, when chosen reasonably, they let the time series of hidden states *x*^{(k)} quickly converge to values that do not depend on the choice of initial value to any appreciable extent. In the simulations of Figures Figures55–8, we used ${\mu}_{2}^{(0)}=0$, ${\sigma}_{2}^{(0)}=1$, and ${\sigma}_{3}^{(0)}=1$.

If we reduce (the log-volatility of *x*_{3}) from 0.5 in the reference scenario to 0.05, we find an agent which is overly confident about its prior estimate of environmental volatility and expects to see little change (Figure (Figure6;6; Figure Figure9,9, second row). This leads to a greatly diminished learning rate for *x*_{3}, while learning in *x*_{2} is not directly affected. There is, however, an indirect effect on *x*_{2} in that the learning rate at the second level during the third period is no longer noticeably increased by the preceding period of higher volatility. In other words, this agent shows superficial (low-level) adaptability but has higher-level beliefs that remain impervious to new information.

Figure Figure77 (and Figure Figure9,9, third row) illustrate the effect of reducing *ω*, the absolute (i.e., independent of *x*_{3}) component of log-volatility, to −4. The multiplicative scaling exp(*ω*) of volatility is thus reduced to a sixth of that in the reference scenario. This leads to a low learning rate for *x*_{2}, which in turn leads to little learning in *x*_{3}, since the agent can only infer changes in *x*_{3} from changes in *x*_{2} (cf. Eq. 30). This corresponds to an agent who pays little attention to new information, effectively filtering it out at an early stage.

The coupling between *x*_{2} and *x*_{3} can be diminished by reducing the value of *κ*, the relative (i.e., *x*_{3}-dependent) scaling of volatility in *x*_{2} (*κ*=0.2 in Figure Figure8;8; Figure Figure9,9, bottom row). This impedes the flow of information up the hierarchy of levels in such a way that the agent's belief about *x*_{3} is effectively insulated from the effects of prediction error in *x*_{2} (cf. Eq. 30). This leads to less learning about *x*_{3} and to a much larger posterior variance *σ*_{3} than in any of the above scenarios (see Figure Figure9,9, right panel). As with a reduced (Figure (Figure6),6), learning about *x*_{2} itself is not directly affected, except that, in the second stage of the simulation, higher volatility remains without effect on the learning rate of *x*_{2} in the third stage. This time, however, this effect is not caused by overconfidence about *x*_{3} (due to small ) as in the above scenario. Instead, it obtains despite uncertainty about *x*_{3} (large *σ*_{3}), which would normally be expected to lead to greater learning because of the dependency of the learning rate on *σ*_{3}. This paradoxical effect can be understood by examining Eqs 30 and 29, where smaller *κ* exerts opposite direct and indirect effects on the learning rate for *μ*_{3}. Indirectly, the learning rate is increased, in that smaller *κ* increases *σ*_{3}. But this is dominated by the dependency of the learning rate on *κ*, which leads to a decrease in learning for smaller *κ*. This is quite an important property of the model: it makes it possible to have low learning rates in a highly volatile environment. This scenario describes an agent which is keen to learn but fails because the levels of its model are too weakly coupled for information to be passed efficiently up the hierarchy. In other words, the agent's low-level adaptability is accompanied by uncertainty about higher-level variables (i.e., volatility), leading to inflexibility. In anthropomorphic terms, one might imagine a person who displays rigid behavior because he/she remains uncertain about how volatile the world is (e.g., in anxiety disorders).

The simulations described above switch between two probability regimes: *p*(*x*_{1}=1)=0.5 and *p*(*x*_{1}=1)=0.9 or 0.1. The stimulus distributions under these two regimes have different variances (or risk). Figure Figure1010 shows an additional simulation run where risk is constant, i.e., *p*(*x*_{1}=1)=0.85 or 0.15, throughout the entire simulation. One recovers the same effects as in the reference scenario. We now consider generalizations of the generative model that relax some of the simplifying assumptions about sensory mappings and outcomes we made during its exposition above.

The model can readily accommodate perceptual uncertainty at the first level. This pertains to the mapping from stimulus category *x*_{1} to sensory input *u*. To allow for perceptual uncertainty, for example when the sensory inputs are ambiguous or affected by noise, we replace the deterministic relation in Eq. 1 by a stochastic one (cf. Daunizeau et al., 2010b):

$$p(u|{x}_{1})=\mathcal{N}{(u;{\eta}_{a},\alpha )}^{{x}_{1}}\mathcal{N}{(u;{\eta}_{b},\alpha )}^{1-{x}_{1}}$$

(45)

Here, the input *u* is no longer binary but a real number whose distribution is a mixture of Gaussians. If *x*_{1}=1, the probability of *u* is normally distributed with constant variance a around a constant value *η _{a}*, corresponding to the most likely sensation if

$$\begin{array}{l}I({x}_{1}^{(k)})={(-\frac{{({u}^{(k)}-{\eta}_{a})}^{2}}{2\alpha}+\text{ln}s({\mu}_{2}^{(k-1)}))}^{{x}_{1}^{(k)}}\hfill \\ \text{\xd7(\u2212(u(k)\u2212\eta b)22\alpha +lns(1\u2212\mu 2(k\u22121)))1\u2212x1(k)}\hfill \end{array}$$

(46)

With Eq. 15, we find the update rule for *μ*_{1}:

$$\begin{array}{l}{\mu}_{1}^{(k)}=\widehat{q}({x}_{1}^{(k)}=1)\hfill \\ =\frac{\text{exp}(-\frac{{({u}^{(k)}-{\eta}_{a})}^{2}}{2\alpha})s({\mu}_{2}^{(k-1)})\text{exp}(-\frac{{({u}^{(k)}-{\eta}_{a})}^{2}}{2\alpha})s({\mu}_{2}^{(k-1)})+\text{exp}(-\frac{{({u}^{(k)}-{\eta}_{b})}^{2}}{2\alpha})(1-s({\mu}_{2}^{(k-1)}))}{}\hfill \end{array}$$

(47)

If *u*^{(k)}≈*η _{a}*, one sees that ${\mu}_{1}^{(k)}\approx 1$ regardless of ${\mu}_{2}^{(k)}$. Likewise, ${\mu}_{1}^{(k)}\approx 0$ if

A generative model that comprises a hierarchy of Gaussian random walks can be applied to many problems of learning and inference. Here, we provide a final example, where the bottom-level state being inferred is not binary but continuous (i.e., a real number). Our example concerns the exchange rate between the U.S. Dollar (USD) and the Swiss Franc (CHF) during the first 180 trading days of the year 2010 (source: http://www.oanda.com). In this example, the agent represented by our model can be seen as an individual market observer (e.g., a currency trader), with the model describing how he “perceives” the relative value (and volatility) of USD. To maintain notational continuity with the preceding sections, we call this value *x*_{2} even though it occupies the lowest level of the hierarchy. In other words, the input *u* is generated directly from *x*_{2} (without passing through a binary state *x*_{1}). The data *u* are taken to be the closing USD–CHF exchange rates of each trading day with a likelihood model

$$p(u|{x}_{2})=\mathcal{N}(u;{x}_{2},\alpha )$$

(48)

where a is the constant variance with which the input *u* is normally distributed around the true value *x*_{2}. This can be regarded as a measure of uncertainty (i.e., how uncertain the trader is about his “perception” of USD value relative to CHF). On top of this input level, we can now add as many coupled random walks as we please. For our example, we can describe the hidden states in higher levels with Eqs 4 and 5. By the method introduced above, we obtain

$$I({x}_{2}^{(k)})=-\frac{1}{2\alpha}{({x}_{2}^{(k)}-{u}^{(k)})}^{2}-\frac{1}{2({\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega})}{({x}_{2}^{(k)}-{\mu}_{2}^{(k-1)})}^{2}$$

(49)

$${\sigma}_{2}^{(k)}=\frac{\alpha ({\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega})\alpha +{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}{}$$

(50)

$$\begin{array}{c}{\mu}_{2}^{(k)}={\mu}_{2}^{(k-1)}+{\sigma}_{2}^{(k)}\frac{1}{\alpha}({u}^{(k)}-{\mu}_{2}^{(k-1)})\\ ={\mu}_{2}^{(k-1)}+\frac{{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+v}}{\alpha +{\sigma}_{2}^{(k-1)}+{\text{e}}^{\kappa {\mu}_{3}^{(k-1)}+\omega}}({u}^{(k)}-{\mu}_{2}^{(k-1)})\end{array}$$

(51)

For *α*=0 (no perceptual uncertainty), the last two equations reduce, as they should, to ${\sigma}_{2}^{(k)}=0,\text{\mu 2(k)=u(k)}$. Note also that since *I*(*x*_{2}) is already quadratic here, no further approximation to the mean-field approximation is needed, and *μ*_{2} and *σ*_{2} are the exact moments of the posterior under the mean-field approximation. Because the higher levels remain the same, the update equations for *μ*_{3} and *σ*_{3} are as given in Eqs 29 and 30.

Scenarios with different parameter values for the USD–CHF example are presented in Figures Figures1111–14. These scenarios can be thought of as corresponding to different individual traders who receive the same market data but process them differently. The reference scenario (Figure (Figure11)11) is based on the parameter values *κ*=1, *ω*=−12, =0.3, and *α*=2·10^{−5}. This parameterization conveys a small amount of perceptual uncertainty that leads to minor but visible deviations of *μ*_{2} from *u*. The updates to *μ*_{2} are conservative in the sense that they consider prior information along with new input. Note also that *μ*_{3} rises whenever the prediction error about *x*_{2} is large, that is when the green dots denoting *u* are outside the range ${\mu}_{2}\pm \sqrt{{\sigma}_{2}}$ indicated by the red lines. Conversely, *μ*_{3} falls when predictions of *x*_{2} are more accurate. In the next scenario (Figure (Figure12),12), the value of *α* is further reduced to 10^{−6}. This scenario thus shows an agent who is effectively without perceptual uncertainty. As prescribed by the update equations above, *μ*_{2} now follows *u* with great accuracy and *μ*_{3} tracks the amount of change in *x*_{2}. In Figure Figure13,13, perceptual uncertainty is increased by two orders of magnitude (*α*=10^{−4}). Here, the agent adapts more slowly to changes in the exchange rate since it cannot be sure whether prediction error is due to a change in the true value of *x*_{2} or to misperception. The final scenario in Figure Figure1414 shows an agent with the same perceptual uncertainty as in the reference scenario but a prior belief that the environment is not very volatile, i.e., is reduced from 0.3 to 0.01. Smaller values of smooth the trajectory of *μ*_{3} in a similar way that perceptual uncertainty smoothes the trajectory of *μ*_{2}.

This example again emphasizes the fact that Bayes-optimal behavior can manifest in many diverse forms. The different behaviors emitted by the agents above are all optimal under their implicit prior beliefs encoded by the parameters that control the evolution of high-level hidden states. Clearly, it would be possible to optimize these parameters using the same variational techniques we have considered for the hidden states. This would involve optimizing the free energy bound on the evidence for each agent's model (prior beliefs) integrated over time (i.e., learning). Alternatively, one could optimize performance by selecting those agents with prior beliefs (about the parameters) that had the best free energy (made the most accurate inferences over time). Colloquially, this would be the difference between training an expert to predict financial markets and simply hiring experts whose priors were the closest to the true values. We will deal with these issues of model inversion and selection in forthcoming work (Mathys et al., in preparation). We close this article with a discussion of the neurobiology behind variations in priors and the neurochemical basis of differences in the underlying parameters of the generative model.

In this article, we have introduced a generic hierarchical Bayesian framework that describes inference under uncertainty; for example, due to environmental volatility or perceptual uncertainty. The model assumes that the states evolve as Gaussian random walks at all but the first level, where their volatility (i.e., conditional variance of the state given the previous state) is determined by the next highest level. This coupling across levels is controlled by parameters, whose values may differ across subjects. In contrast to “ideal” Bayesian learning models, which prescribe a fixed process for any agent, this allows for the representation of inter-individual differences in behavior and how it is influenced by uncertainty. This variation is cast in terms of prior beliefs about the parameters coupling hierarchical levels in the generative model.

A major goal of our work was to eschew the complicated integrals in exact Bayesian inference and instead derive analytical update equations with algorithmic efficiency and biological plausibility. For this purpose, we used an approximate (variational) Bayesian approach, under a mean-field assumption and a novel approximation to the posterior energy function. The resulting single-step, trial-by-trial update equations have several important properties:

- (i) They have an analytical form and are extremely efficient, allowing for real-time inference.
- (ii) They are biologically plausible in that the mathematical operations required for calculating the updates are fairly basic and could be performed by single neurons (London and Hausser, 2005; Herz et al., 2006).
- (iii) Their structure is remarkably similar to update equations from standard RL models; this enables an interpretation that places RL heuristics, such as learning rate or prediction error, in a principled (Bayesian) framework.
- (iv) The model parameters determine processes, such as precision-weighting of prediction errors, which play a key role in current theories of normal and pathological learning and may relate to specific neuromodulatory mechanisms in the brain (see below).
- (v) They can accommodate states of either discrete or continuous nature and can deal with deterministic and probabilistic mappings between environmental causes and perceptual consequences (i.e., situations with and without perceptual uncertainty).

Crucially, the closed-form update equations do not depend on the details of the model but only on its hierarchical structure and the assumptions on which the mean-field and quadratic approximation to the posteriors rest. Our method of deriving the update equations may thus be adopted for the inversion of a large class of models. Above, we demonstrated this anecdotally by providing update equations for two extensions of the original model, which accounted for sensory states of a continuous (rather than discrete) nature and perceptual uncertainty, respectively.

As alternatives to our variational scheme, one could deal with the complicated integrals of Bayesian inference by sampling methods or avoid them altogether and use simpler RL schemes. We did not pursue these options because we wanted to take a principled approach to individual learning under uncertainty; i.e., one that rests on the inversion of a full Bayesian generative model. Furthermore, we wanted to avoid sampling approximations because of the computational burden they impose. Although it is conceivable that neuronal populations could implement sampling methods, it is not clear how exactly they would do that and at what temporal and energetic cost (Yang and Shadlen, 2007; Beck et al., 2008; Deneve, 2008).

We would like to emphasize that the examples of update equations derived here can serve as the building blocks for those of more complicated models. For example, if we have more than two categories at the first level, this can be accommodated by additional random walks at the second and subsequent levels; at those levels, Eqs 38 and 40 have a straightforward interpretation in *n* dimensions. Inference using our update scheme is thus possible on *n*-categorical discrete states, *n*-dimensional unbounded continuous states, and (by logarithmic or logistic transformation of variables) *n*-dimensional bounded continuous states.

One specific problem that has been addressed with Bayesian methods in the recent past concerns online inference of “changepoints,” i.e., sudden changes in the statistical structure of the environment (Corrado et al., 2005; Fearnhead and Liu, 2007; Krugel et al., 2009; Nassar et al., 2010; Wilson et al., 2010). Our generative model based on hierarchically coupled random walks describes belief updating without representing such discrete changepoints explicitly. As illustrated by the simulations in Figures Figures55–8, this does not diminish its ability to deal with volatile environments. See also previous studies where similar models were applied to data generated by changepoint models (Behrens et al., 2007, 2008) or where RL-type update models were equipped with an adjustable learning rate in order to deal with sudden changes in the environment (Krugel et al., 2009; Nassar et al., 2010). One such sudden change in the environment that is nicely picked up in the application of our model to the empirical exchange rate data is the outbreak of the Greek financial crisis in spring 2010. The sudden realization of the financial markets that Greece was insolvent led to a flight into the USD which is reflected by a sharply increasing value of the Dollar against the Swiss Franc visible in the lower panel of Figure Figure11.11. Because this sudden rise of the Dollar is unexpected, it immediately leads to a jump in the agent's belief about the volatility of its environment, as is clearly visible in the upper panel of Figures Figures1111 or or13.13. In this manner, a sudden event akin to a changepoint is detected without representations of changepoints being an explicit component of the model.

Clearly, our approach is not the first that has tried to derive tractable update equations from the full Bayesian formulation of learning. Although not described in this way in the original work, even the famous Kalman filter can be interpreted as a Bayesian scheme with RL-like update properties but is restricted to relatively simple (non-hierarchical) learning processes in stable environments (Kalman, 1960). Notably, none of the previous Bayesian learning models we know leads to analytical one-step update equations without resorting to additional assumptions that are specifically tailored for the update equations (Nassar et al., 2010) or the learning rate (Krugel et al., 2009). In contrast, in our scheme, the update equations and their critical components, such as learning rate or prediction error, emerge naturally by inverting a full Bayesian generative model of arbitrary hierarchical depth under a generic mean-field reduction. That is, once we have specified the nature of our approximate inversion, the update equations are fully defined and do not require any further assumptions. This distinguishes our framework conceptually and mathematically from any previously suggested approach to Bayesian learning we are aware of.

Our update scheme also has its limitations. The most important of these is that it depends on the variational energies being approximately quadratic. If they are not, the approximate posterior implied by our update equations might bear little resemblance (e.g., in terms of Kullback–Leibler divergence) to the true posterior. Specifically, the update fails if the curvature of the variational energy at the expansion point is negative (which implies that the conditional variance is negative; see Eq. 38). According to the update Eq. 22, *σ*_{2} can never become negative; *σ*_{3}=1/π_{3} however could become negative according to Eq. 29. However, the simulations and application to empirical (exchange rate) data in this paper suggest that this is not a problem in practice. We are currently investigating the properties of our scheme more systematically in a series of theoretical and empirical analyses that will be published in future work. In particular, we will compare our variational inversion scheme against other methods such as numerical integration and sampling-based methods.

In addition to the derivation of these update equations, we have provided simulations that demonstrate model behavior under different parameterizations (priors). These simulations confirmed the computational efficiency of our approach: the simulations in Figures Figures55–8 (with 320 trials) each take about 5 ms on a standard laptop computer. Furthermore, they demonstrate that changes in any of the parameters lead to plausible changes in the evolution of the states (i.e., as predicted from the structure of the model) and that each parameter produces distinctly different behavior.

Maladaptive behavior, owing to inappropriate learning and decision-making, is at the heart of most psychiatric diseases, and our framework may be particularly useful for modeling the underlying mechanisms. There are two complementary approaches one might consider: phenomenological and neurophysiological. To illustrate a phenomenological approach, we will consider extreme settings of the parameters in terms of psychopathology. In the example of Figures Figures55–8, variations in the parameters can explain a spectrum of different types of inference, some of which may be interpreted as aberrant or even pathological. Given the scenario in Figure Figure8,8, we could adopt an anthropomorphic interpretation of the agent and interpret underconfidence about estimates of environmental volatility (i.e., high *σ*_{3}) as a possible cause of anxiety. In other words, knowing that the world is changing quickly is frightening enough, but being uncertain about the extent of this change may be even more upsetting. Anxiety of this sort is often observed prior to (or in the early phase of) psychotic episodes (Häfner et al., 1998). One way to reduce anxiety (that is, to reduce the effects of high s_{3} due to abnormally low *κ*), would be to reduce , leading to a scenario akin to that in Figure Figure6.6. This, however, would induce a rigid high-level belief that is impervious to prediction error from the lower level. Rigid high-level priors of this sort that provide inappropriate predictions for lower levels may provide a metaphor for delusions and hallucinations that constitute the positive symptoms of schizophrenia. In contrast, negative symptoms could be related to the scenario in Figure Figure7:7: here, a reduction in *ω* renders the agent completely passive, such that new information is barely taken in and only weakly processed. Notably, in these simple and anecdotal simulations, we chose some parameter settings that lead to superficially similar behavior (e.g., Figures Figures66 and and8).8). While this indicates some degree of interdependence among the parameters, this does not mean that the parameters are non-identifiable. Informally, one can intuit this by noting the obvious differences expressed in the evolution of higher-level states of the model; these will be expressed in different behavioral predictions, given a suitably chosen sequence of stimuli. When fitting the model to empirical data, once can test for parameter identifiability more formally using a sensitivity analysis or, equivalently but more conveniently, their posterior covariance. The identifiability of our model parameters will be examined systematically in forthcoming work (Mathys et al., in preparation).

From a neurophysiological perspective, it has been proposed that dopamine might not encode the prediction error on value (Schultz et al., 1997) but instead the value of prediction error, i.e., the precision-weighting of prediction errors (Friston, 2009). In our model, this process is represented, at the second level, by the parameters k and v. It is apparent from the definition Eq. 27 and the update Eq. 22 that these parameters influence the precision of the prediction on the next trial, the precision of the posterior belief, and the learning rate. If dopaminergic midbrain activity encodes the conditional (posterior) precision of beliefs, this dopaminergic activity should be reflected by estimates of *κ* and *ω*, obtained from behavioral, fMRI, or electrophysiological data. This hypothesis can be tested using neuropharmacological experiments. In short, by harvesting subject-specific parameter estimates for group analyses of physiological measurements, hierarchical generative models (of the sort considered in this work) could be used to test hypotheses about the relations between computational and physiological processes. We are currently pursuing this approach in ongoing research. Alternatively, one can also use our model for analyses at the subject level: the sequence of inferred hidden states, as represented by their sufficient statistics *μ* and *σ*, can be used as predictor variables in analyses of fMRI, EEG, or behavioral data to shed light on the neurophysiological correlates of inference and learning (cf. den Ouden et al., 2010).

The distinction between hidden states, which vary in time and are the dynamic components of the agent's model of the world, and parameters, which are time-invariant and encode stable subject-specific learning styles, is a key component of our model. One might compare this to classical RL models where value estimates (states) are updated dynamically while the learning rate is an invariant parameter. In our case, however, the (implicit) learning rate is dynamic and results from an interaction between states and parameters: the latter determine how higher-level states influence lower-level ones. This effect of the static parameters on dynamic cross-level coupling can be seen directly from the update equations above (e.g., Eqs 23 and 27) and is illustrated in Figures Figures55 and and1010 where the learning rate visibly changes while the parameters are fixed. In other words, subject-specific learning mechanisms, represented by cross-level coupling in our model, have both dynamic (higher-level states) and static (parameters) components in our model.

One could, of course, consider alternative formulations of our model in which individual learning mechanisms, determined by the coupling across levels, are encoded entirely by states. There is a simple reason why we did not pursue this alternative. Clearly, modeling dynamic aspects of learning, such as rapid updating of learning rates, does require a representation involving states (see above). On the other hand, within an individual agent's brain, the physiological mechanism underlying this coupling must obey some general principles that have been shaped both by (life-long) experience and genetic background (cf. RL depends on individual genotype (Frank et al., 2007, 2009; Krugel et al., 2009)). Such stable subject-specific learning mechanisms could be represented in two ways. One could choose a relatively deep hierarchical model with high-level states that change very slowly. Alternatively, these mechanisms could be represented by time-invariant parameters. We chose the latter option simply because it provides a more concise and interpretable summary of subject-specific learning mechanisms. For example, when quantifying individual differences in computational learning mechanisms as a function of individual differences in physiology (e.g., pharmacological treatment) or genetics, it is not only statistically easier to deal with time-invariant parameters (i.e., a single number per subject) rather than temporal trajectories of states, but the results are also more readily interpretable.

It should be emphasized that the idea of fitting learning models to subject-specific data and using the ensuing individual parameter estimates for assessing inter-individual variability is not new and has been pursued by many previous studies (e.g., Steyvers and Brown, 2006; Frank et al., 2007, 2009; Krugel et al., 2009). This, however, is less straightforward with those “ideal” Bayesian models that have no free parameters; in this case, parameters can only be taken from adjunct models in which ideal Bayesian models are often embedded (e.g., observation or decision models; cf. Brodersen et al., 2008). The novelty of our approach is that we transform, by variational approximation, an ideal Bayesian learner into a near-optimal scheme in which parameters represent individual learning traits as an integral part of Bayesian learning. These parameters shape the ensuing update equations which are analytical and have an RL-like structure. By combining the principled nature of Bayesian approaches and the practical ease of RL models, we hope that our approach will facilitate future empirical investigations of individual variability.

The main goal of this paper was to introduce the mathematical basis of our approach and illustrate its functionality. Clearly, the simulations shown in this paper are anecdotal and cannot fully demonstrate or establish the practical utility of our approach. In particular, we have not yet demonstrated how our model can be inverted (fitted), given empirical measurements. This requires extending the present approach with a response model that connects states from the Bayesian learning model to measurable responses of a neurophysiological or behavioral sort (cf. Daunizeau et al., 2010a,b). We will present this extension in future work and demonstrate its use for the analysis of behavioral and neuroimaging data.

In summary, we have introduced a novel and generic framework for approximate Bayesian inference with computationally efficient and interpretable closed-form update equations. Simulations show that our approach is applicable to a range of situations beyond classical RL, including inductive inference on discrete and continuous states and situations with perceptual ambiguity. Crucially, our approach accommodates inter-individual differences, in terms of prior beliefs about key model parameters, and quantifies their computational effects: Some of these parameters may map to neurophysiological (neuromodulatory) mechanisms that have been implicated in the neurobiology of learning and psychopathology. As such, it may be a useful framework for modeling individual differences in behavior and to formally characterize behavioral stereotypes and pathophysiologically distinct subgroups in psychiatric spectrum diseases (Stephan et al., 2009).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This work was supported by the NCCR “Neural Plasticity and Repair” (CM, KES), the University Research Priority Program “Foundations of Human Social Behavior” at the University of Zurich (KES), the Neurochoice project of the Swiss Systems Biology initiative SystemsX.ch (JD, KES), and the Wellcome Trust (KJF).

- Beal M. J. (2003). Variational Algorithms for Approximate Bayesian Inference. Ph.D. thesis, University College London; Available at: http://www.cse.buffalo.edu/faculty/mbeal/papers/beal03.pdf
- Beck J., Ma W., Kiani R., Hanks T., Churchland A., Roitman J., Shadlen M., Latham P. E., Pouget A. (2008). Probabilistic population codes for Bayesian decision making. Neuron 60, 1142–115210.1016/j.neuron.2008.09.021 [PMC free article] [PubMed] [Cross Ref]
- Behrens T. E. J., Hunt L. T., Woolrich M. W., Rushworth M. F. S. (2008). Associative learning of social value. Nature 456, 245–24910.1038/nature07538 [PMC free article] [PubMed] [Cross Ref]
- Behrens T. E. J., Woolrich M. W., Walton M. E., Rushworth M. F. S. (2007). Learning the value of information in an uncertain world. Nat. Neurosci. 10, 1214–122110.1038/nn1954 [PubMed] [Cross Ref]
- Bresciani J., Dammeier F., Ernst M. O. (2006). Vision and touch are automatically integrated for the perception of sequences of events. J. Vis. 6, 554–56410.1167/6.5.2 [PubMed] [Cross Ref]
- Brodersen K. H., Penny W. D., Harrison L. M., Daunizeau J., Ruff C. C., Duzel E., Friston K. J., Stephan K. E. (2008). Integrated Bayesian models of learning and decision making for saccadic eye movements. Neural Netw. 21, 1247–126010.1016/j.neunet.2008.08.007 [PMC free article] [PubMed] [Cross Ref]
- Corrado G. S., Sugrue L. P., Sebastian Seung H., Newsome W. T. (2005). Linear-nonlinear-poisson models of primate choice dynamics. J. Exp. Anal. Behav. 84, 581–61710.1901/jeab.2005.23-05 [PMC free article] [PubMed] [Cross Ref]
- Cox R. T. (1946). Probability, frequency and reasonable expectation. Am. J. Phys. 14, 1–1310.1119/1.1990764 [Cross Ref]
- Daunizeau J., den Ouden H. E. M., Pessiglione M., Kiebel S. J., Stephan K. E., Friston K. J. (2010a). Observing the observer (I): meta-Bayesian models of learning and decision-making. PLoS ONE 5, e15554.10.1371/journal.pone.0015554 [PMC free article] [PubMed] [Cross Ref]
- Daunizeau J., den Ouden H. E. M., Pessiglione M., Kiebel S. J., Friston K. J., Stephan K. E. (2010b). Observing the observer (II): deciding when to decide. PLoS ONE 5, e15555.10.1371/journal.pone.0015555 [PMC free article] [PubMed] [Cross Ref]
- Daw N. D., O'Doherty J. P., Dayan P., Seymour B., Dolan R. J. (2006). Cortical substrates for exploratory decisions in humans. Nature 441, 876–87910.1038/nature04766 [PMC free article] [PubMed] [Cross Ref]
- Dayan P., Huys Q. J. (2009). Serotonin in affective control. Annu. Rev. Neurosci. 32, 95–12610.1146/annurev.neuro.051508.135607 [PubMed] [Cross Ref]
- Dayan P., Niv Y. (2008). Reinforcement learning: the good, the bad and the ugly. Curr. Opin. Neurobiol. 18, 185–19610.1016/j.conb.2008.08.003 [PubMed] [Cross Ref]
- den Ouden H. E. M., Daunizeau J., Roiser J., Friston K. J., Stephan K. E. (2010). Striatal prediction error modulates cortical coupling. J. Neurosci. 30, 3210–321910.1523/JNEUROSCI.4458-09.2010 [PMC free article] [PubMed] [Cross Ref]
- Deneve S. (2008). Bayesian spiking neurons II: learning. Neural. Comput. 20, 118–14510.1162/neco.2008.20.1.118 [PubMed] [Cross Ref]
- Doya K. (2008). Modulators of decision making. Nat. Neurosci. 11, 410–41610.1038/nn2077 [PubMed] [Cross Ref]
- Fearnhead P., Liu Z. (2007). On-line inference for multiple changepoint problems. J. R. Stat. Soc. Series B Stat. Methodol. 69, 589–60510.1111/j.1467-9868.2007.00601.x [Cross Ref]
- Frank M. J. (2008). Schizophrenia: a computational reinforcement learning perspective. Schizophr. Bull. 34, 1008–101110.1093/schbul/sbn123 [PMC free article] [PubMed] [Cross Ref]
- Frank M. J., Doll B. B., Oas-Terpstra J., Moreno F. (2009). Prefrontal and striatal dopaminergic genes predict individual differences in exploration and exploitation. Nat. Neurosci. 12, 1062–106810.1038/nn.2342 [PMC free article] [PubMed] [Cross Ref]
- Frank M. J., Moustafa A. A., Haughey H. M., Curran T., Hutchison K. E. (2007). Genetic triple dissociation reveals multiple roles for dopamine in reinforcement learning. Proc. Natl. Acad. Sci. U.S.A. 104, 16311–1631610.1073/pnas.0706111104 [PubMed] [Cross Ref]
- Friston K. (2008). Hierarchical models in the brain. PLoS Comput. Biol. 4, e1000211.10.1371/journal.pcbi.1000211 [PMC free article] [PubMed] [Cross Ref]
- Friston K. (2009). The free-energy principle: a rough guide to the brain? Trends Cogn. Sci. (Regul. Ed.) 13, 293–301 [PubMed]
- Friston K., Mattout J., Trujillo-Barreto N., Ashburner J., Penny W. (2007). Variational free energy and the Laplace approximation. Neuroimage 34, 220–23410.1016/j.neuroimage.2006.08.035 [PubMed] [Cross Ref]
- Friston K. J., Stephan K. E. (2007). Free-energy and the brain. Synthese 159, 417–45810.1007/s11229-007-9237-y [PMC free article] [PubMed] [Cross Ref]
- Geisler W. S., Diehl R. L. (2002). Bayesian natural selection and the evolution of perceptual systems. Philos. Trans. R. Soc. B Biol. Sci. 357, 419–448 [PMC free article] [PubMed]
- Gershman S. J., Niv Y. (2010). Learning latent structure: carving nature at its joints. Curr. Opin. Neurobiol. 20, 251–25610.1016/j.conb.2010.02.008 [PMC free article] [PubMed] [Cross Ref]
- Gluck M. A., Shohamy D., Myers C. (2002). How do people solve the weather prediction task? individual variability in strategies for probabilistic category learning. Learn. Mem. 9, 408–41810.1101/lm.45202 [PubMed] [Cross Ref]
- Gu Q. (2002). Neuromodulatory transmitter systems in the cortex and their role in cortical plasticity. Neuroscience 111, 815–83510.1016/S0306-4522(02)00026-X [PubMed] [Cross Ref]
- Häfner H., Maurer K., Löffler W., an der Heiden W., Munk-Jørgensen P., Hambrecht M., Riecher-Rössler A. (1998). The ABC schizophrenia study: a preliminary overview of the results. Soc. Psychiatry Psychiatr. Epidemiol. 33, 380–38610.1007/s001270050069 [PubMed] [Cross Ref]
- Herz A. V. M., Gollisch T., Machens C. K., Jaeger D. (2006). Modeling single-neuron dynamics and computations: a balance of detail and abstraction. Science 314, 80–8510.1126/science.1127240 [PubMed] [Cross Ref]
- Jaynes E. T. (1957). Information theory and statistical mechanics. Phys. Rev. 106, 620.10.1103/PhysRev.106.620 [Cross Ref]
- Jaynes E. T. (2003). Probability Theory: The Logic of Science. Cambridge, UK: Cambridge University Press; 10.1017/CBO9780511790423 [Cross Ref]
- Kalman R.E. (1960). A new approach to linear filtering and prediction problems. J.Basic Eng. 82, 35–45
- Kording K. P., Wolpert D. M. (2004). Bayesian integration in sensorimotor learning. Nature 427, 244–24710.1038/nature02169 [PubMed] [Cross Ref]
- Krugel L. K., Biele G., Mohr P. N. C., Li S., Heekeren H. R. (2009). Genetic variation in dopaminergic neuromodulation influences the ability to rapidly and flexibly adapt decisions. Proc. Natl. Acad. Sci. U.S.A. 106, 17951–1795610.1073/pnas.0905191106 [PubMed] [Cross Ref]
- Laplace P. (1774). Mémoire sur la probabilité des causes par les évènemens. Mém. Acad. Roy. Sci. 6, 621–656
- Laplace P. (1812). Théorie Analytique des Probabilités. Paris: Courcier Imprimeur
- London M., Hausser M. (2005). Dendritic computation. Annu. Rev. Neurosci. 28, 503–53210.1146/annurev.neuro.28.061604.135703 [PubMed] [Cross Ref]
- Montague P. R., Hyman S. E., Cohen J. D. (2004). Computational roles for dopamine in behavioural control. Nature 431, 760–76710.1038/nature03015 [PubMed] [Cross Ref]
- Murray G. K., Corlett P. R., Clark L., Pessiglione M., Blackwell A. D., Honey G., Jones P. B., Bullmore E. T., Robbins T. W., Fletcher P. C. (2007). Substantia nigra/ventral tegmental reward prediction error disruption in psychosis. Mol. Psychiatry 13, 267–27610.1038/sj.mp.4002058 [PMC free article] [PubMed] [Cross Ref]
- Nassar M. R., Wilson R. C., Heasly B., Gold J. I. (2010). An approximately Bayesian delta-rule model explains the dynamics of belief updating in a changing environment. J. Neurosci. 30, 12366–1237810.1523/JNEUROSCI.0822-10.2010 [PMC free article] [PubMed] [Cross Ref]
- O'Doherty J., Dayan P., Schultz J., Deichmann R., Friston K., Dolan R. J. (2004). Dissociable roles of ventral and dorsal striatum in instrumental conditioning. Science 304, 452–45410.1126/science.1094285 [PubMed] [Cross Ref]
- Orbán G., Fiser J., Aslin R. N., Lengyel M. (2008). Bayesian learning of visual chunks by human observers. Proc. Natl. Acad. Sci. U.S.A. 105, 2745–275010.1073/pnas.0708424105 [PubMed] [Cross Ref]
- Pessiglione M., Seymour B., Flandin G., Dolan R. J., Frith C. D. (2006). Dopamine-dependent prediction errors underpin reward-seeking behavior in humans. Nature 442, 1042–104510.1038/nature05051 [PMC free article] [PubMed] [Cross Ref]
- Preuschoff K., Bossaerts P. (2007). Adding prediction risk to the theory of reward learning. Ann. N. Y. Acad. Sci. 1104, 135–14610.1196/annals.1390.005 [PubMed] [Cross Ref]
- Rescorla R. A., Wagner A. R. (1972). “A theory of Pavlovian conditioning: Variations in the effectiveness of reinforcement,” in Classical Conditioning II: Current Research and Theory, eds Black A. H., Prokasy W. F., editors. (New York: Appleton-Century-Crofts; ), 64–99
- Schultz W., Dayan P., Montague P. R. (1997). A neural substrate of prediction and reward. Science 275, 1593–159910.1126/science.275.5306.1593 [PubMed] [Cross Ref]
- Smith A., Li M., Becker S., Kapur S. (2006). Dopamine, prediction error and associative learning: a model-based account. Network 17, 61–8410.1080/09548980500361624 [PubMed] [Cross Ref]
- Stephan K. E., Friston K. J., Frith C. D. (2009). Dysconnection in schizophrenia: from abnormal synaptic plasticity to failures of self-monitoring. Schizophr. Bull. 35, 509–52710.1093/schbul/sbn176 [PMC free article] [PubMed] [Cross Ref]
- Steyvers M., Brown S. (2006). Prediction and change detection. Adv. Neural Inf. Process Syst. 18, 1281–1288
- Steyvers M., Lee M. D., Wagenmakers E. (2009). A Bayesian analysis of human decision-making on bandit problems. J. Math. Psychol. 53, 168–17910.1016/j.jmp.2008.11.002 [Cross Ref]
- Sutton R. S., Barto A. G. (1998). Reinforcement Learning: An Introduction. Cambridge, MA: MIT Press
- Thiel C. M., Huston J. P., Schwarting R. K. (1998). Hippocampal acetylcholine and habituation learning. Neuroscience 85, 1253–126210.1016/S0306-4522(98)00030-X [PubMed] [Cross Ref]
- Wilson R. C., Nassar M. R., Gold J. I. (2010). Bayesian online learning of the hazard rate in change-point problems. Neural. Comput. 22, 2452–247610.1162/NECO_a_00007 [PMC free article] [PubMed] [Cross Ref]
- Xu F., Tenenbaum J. B. (2007). Sensitivity to sampling in Bayesian word learning. Dev. Sci. 10, 288–29710.1111/j.1467-7687.2007.00590.x [PubMed] [Cross Ref]
- Yang T., Shadlen M. N. (2007). Probabilistic reasoning by neurons. Nature 447, 1075–108010.1038/nature05852 [PubMed] [Cross Ref]
- Yu A. J., Dayan P. (2005). Uncertainty, neuromodulation, and attention. Neuron 46, 681–69210.1016/j.neuron.2005.04.026 [PubMed] [Cross Ref]
- Yuille A., Kersten D. (2006). Vision as Bayesian inference: analysis by synthesis? Trends Cogn. Sci. (Regul. Ed.), 10, 301–308 [PubMed]

Articles from Frontiers in Human Neuroscience are provided here courtesy of **Frontiers Media SA**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |