Learning in Volatile Environments with the Bayes Factor Surprise

TL;DR

This paper introduces the Bayes Factor Surprise as a novel measure to improve adaptive learning in volatile environments, leading to new algorithms that outperform existing methods in parameter estimation.

Contribution

It presents a formalization of surprise as a Bayes Factor, derives three new surprise-based algorithms with simple updates, and demonstrates their effectiveness in non-stationary settings.

Findings

Surprise-modulated algorithms outperform traditional methods in parameter estimation.

The proposed algorithms have constant complexity and simple update rules.

Bayes Factor Surprise can be distinguished from Shannon Surprise in physiological predictions.

Abstract

Surprise-based learning allows agents to rapidly adapt to non-stationary stochastic environments characterized by sudden changes. We show that exact Bayesian inference in a hierarchical model gives rise to a surprise-modulated trade-off between forgetting old observations and integrating them with the new ones. The modulation depends on a probability ratio, which we call "Bayes Factor Surprise", that tests the prior belief against the current belief. We demonstrate that in several existing approximate algorithms the Bayes Factor Surprise modulates the rate of adaptation to new observations. We derive three novel surprised-based algorithms, one in the family of particle filters, one in the family of variational learning, and the other in the family of message passing, that have constant scaling in observation sequence length and particularly simple update dynamics for any distribution in the exponential family. Empirical results show that these surprise-based algorithms estimate parameters better than alternative approximate approaches and reach levels of performance comparable to computationally more expensive algorithms. The Bayes Factor Surprise is related to but different from Shannon Surprise. In two hypothetical experiments, we make testable predictions for physiological indicators that dissociate the Bayes Factor Surprise from Shannon Surprise. The theoretical insight of casting various approaches as surprise-based learning, as well as the proposed online algorithms, may be applied to the analysis of animal and human behavior, and to reinforcement learning in non-stationary environments.

Tables

Table 1. Table 1: Experimental Hypotheses and Predictions 1. Δ ​ M ¯ 1 ​ ( θ ^ , δ , C ) Δ subscript ¯ 𝑀 1 ^ 𝜃 𝛿 𝐶 \Delta\bar{M}_{1}(\hat{\theta},\delta,C) stands for M ¯ 1 ( θ ^ , δ , s = + 1 , C ) − M ¯ 1 ( θ ^ , δ , s = − 1 , C ) \bar{M}_{1}(\hat{\theta},\delta,s=+1,C)-\bar{M}_{1}(\hat{\theta},\delta,s=-1,C)

Hypothesis	Prediction
The indicator reflects ${\text{S}}_{\mathrm{BF}}$	and
The indicator reflects ${\text{S}}_{\mathrm{Sh}}$	$\Delta {\overline{M}}_1(\hat{\theta },\delta ,C)>0$ and $\frac{\partial \Delta {\overline{M}}_1(\hat{\theta },\delta ,C)}{\partial \delta }>0$
The prior is not used for inference	$\Delta {\overline{M}}_1(\hat{\theta },\delta ,C)=0$

Table 2. Table 2: Experimental Hypotheses and Predictions 2.

Hypothesis	Prediction
The indicator reflects ${\text{S}}_{\mathrm{BF}}$	$\frac{\partial {\overline{M}}_2\left(p\right)}{\partial p}=0$
The indicator reflects ${\text{S}}_{\mathrm{Sh}}$

Equations

$$\text{P}(c_{1:T},{\theta }_{1:T},y_{1:T})=\text{P}(c_1)\text{P}({\theta }_1)\text{P}(y_1|{\theta }_1)\prod _{t=2}^T\text{P}(c_t)\text{P}({\theta }_t|c_t,{\theta }_{t-1})\text{P}(y_t|{\theta }_t),$$

$$\text{P}(c_t)=\text{Bernoulli}(c_t;p_c),$$

$$\displaystyle\textbf{P}(\theta_{t}|c_{t},\theta_{t-1})=\left\{\begin{array}[]{lr}\delta(\theta_{t}-\theta_{t-1})&\text{if}\quad c_{t}=0\,,\\ \pi^{(0)}(\theta_{t})&\text{if}\quad c_{t}=1\,,\end{array}\right.$$

$$\text{P}(y_t|{\theta }_t)=P_Y(y_t|{\theta }_t).$$

$${\pi }^{(t+1)}(\theta )=\frac{P_Y(y_{t+1}|\theta )\text{P}({\Theta }_{t+1}=\theta |y_{1:t})}{\text{P}(y_{t+1}|y_{1:t})}.$$

$$P(y;{\pi }^{(t^{\prime })})=\int P_Y(y|\theta ){\pi }^{(t^{\prime })}(\theta )d\theta .$$

$${\text{S}}_{\mathrm{BF}}(y_{t+1};{\pi }^{(t)})=\frac{P(y_{t+1};{\pi }^{(0)})}{P(y_{t+1};{\pi }^{(t)})}.$$

$${\pi }_B^{(t+1)}(\theta )=\frac{P_Y(y_{t+1}|\theta ){\pi }^{(t)}(\theta )}{P(y_{t+1};{\pi }^{(t)})}.$$

$$\gamma (\text{S},m)=\frac{m\text{S}}{1+m\text{S}},$$

$$\pi^{(t+1)}(\theta)=\Big{(}1-\gamma\Big{(}\textbf{S}_{\mathrm{BF}}^{(t+1)},\frac{p_{c}}{1-p_{c}}\Big{)}\Big{)}\pi^{(t+1)}_{B}(\theta)+\gamma\Big{(}\textbf{S}_{\mathrm{BF}}^{(t+1)},\frac{p_{c}}{1-p_{c}}\Big{)}P(\theta|y_{t+1})\,,$$

$$P(\theta |y_{t+1})=\frac{P_Y(y_{t+1}|\theta ){\pi }^{(0)}(\theta )}{P(y_{t+1};{\pi }^{(0)})}$$

$${\pi }^{\text{new}}(\theta )=(1-\gamma ){\pi }^{\text{integration}}(\theta |y^{\text{new}},{\pi }^{\text{old}})+\gamma {\pi }^{\text{reset}}(\theta |y^{\text{new}},{\pi }^{(0)}).$$

$$\text{log}\big{(}\hat{\pi}^{(t+1)}(\theta)\big{)}=(1-\gamma_{t+1})\text{ log}\big{(}\hat{\pi}_{B}^{(t+1)}(\theta)\big{)}+\gamma_{t+1}\text{ log}\big{(}P(\theta|y_{t+1})\big{)}+\text{Const.}\,,$$

$${\hat{\pi }}^{(t+1)}(\theta )$$

$$\displaystyle\textbf{D}_{KL}\big{[}q(\theta)||P(\theta|y_{t+1})\big{]}\leq B_{t+1}\,,$$

$${\pi }^{(t)}(\theta )$$

$$w_{t+1}^{(t)}={\gamma }_{t+1},$$

$$w_{t+1}^{(k)}=(1-{\gamma }_{t+1})w_{B,t+1}^{(k)}=(1-{\gamma }_{t+1})\frac{P(y_{t+1};{\pi }_k^{(t)})}{P(y_{t+1};{\pi }^{(t)})}{w}_t^{(k)}.$$

$${\pi }^{(t)}(\theta )$$

$$\displaystyle=\mathbb{E}_{\textbf{P}(C_{1:t}|y_{1:t})}\big{[}\textbf{P}(\Theta_{t}=\theta|C_{1:t},y_{1:t})\big{]}.$$

$${\hat{\pi }}^{(t)}(\theta )=\sum _{i=1}^N{w}_t^{(i)}{\hat{\pi }}_i^{(t)}(\theta )=\sum _{i=1}^N{w}_t^{(i)}\text{P}({\Theta }_t=\theta |c_{1:t}^{(i)},y_{1:t}),$$

$$w_{t+1}^{(i)}=(1-{\gamma }_{t+1})w_{B,t+1}^{(i)}+{\gamma }_{t+1}{w}_t^{(i)},$$

$$w_{B,t+1}^{(i)}=\frac{P(y_{t+1};{\hat{\pi }}_i^{(t)})}{P(y_{t+1};{\hat{\pi }}^{(t)})}{w}_t^{(i)},$$

$$\displaystyle Q(c_{t+1}^{(i)}=1|c_{1:t}^{(i)},y_{1:t+1})=\gamma\Big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\hat{\pi}_{i}^{(t)}),\frac{p_{c}}{1-p_{c}}\Big{)}\,.$$

$${\text{S}}_{\mathrm{Sh}}(y_{t+1};{\pi }^{(t)})$$

$$\displaystyle\gamma_{t+1}=p_{c}\text{exp}\Big{(}\Delta\textbf{S}_{\mathrm{Sh}}(y_{t+1};\pi^{(t)},\pi^{(0)})\Big{)},$$

$$\text{ where }\Delta {\text{S}}_{\mathrm{Sh}}(y_{t+1};{\pi }^{(t)},{\pi }^{(0)})={\text{S}}_{\mathrm{Sh}}(y_{t+1};{\pi }^{(t)})-{\text{S}}_{\mathrm{Sh}}(y_{t+1};{\pi }^{(0)}),$$

$$\textbf{S}_{\mathrm{BF}}(y_{t+1};\pi^{(t)})=\frac{P(y_{t+1};\pi^{(0)})}{P(y_{t+1};\pi^{(t)})}=\frac{\mathbb{E}_{\pi^{(0)}}\big{[}P_{Y}(y_{t+1}|\Theta)\big{]}}{\mathbb{E}_{\pi^{(t)}}\big{[}P_{Y}(y_{t+1}|\Theta)\big{]}}.$$

$${\text{S}}_{\mathrm{BF}}$$

$${\pi }^{(t+1)}(\theta )\equiv \text{P}({\Theta }_{t+1}=\theta |y_{1:t+1}).$$

$${\pi }^{(t+1)}(\theta )=\frac{P_Y(y_{t+1}|\theta )\text{P}({\Theta }_{t+1}=\theta |y_{1:t})}{\text{P}(y_{t+1}|y_{1:t})}.$$

$$\text{P}({\Theta }_{t+1}=\theta |y_{1:t})$$

Full text

Learning in Volatile Environments with the

Bayes Factor Surprise111To appear in Neural Computation.

Vasiliki Liakoni†,‡,∗, Alireza Modirshanechi†,‡,∗, Wulfram Gerstner†**, and Johanni Brea†

*†***École Polytechnique Fédérale de Lausanne (EPFL), School of Computer and Communication Sciences and School of Life Sciences, Lausanne, Switzerland

*‡*Equal Contribution

*∗*Corresponding authors:

vasiliki.liakoni@epfl.ch (VL) and alireza.modirshanechi@epfl.ch (AM)

Abstract

Surprise-based learning allows agents to rapidly adapt to non-stationary stochastic environments characterized by sudden changes. We show that exact Bayesian inference in a hierarchical model gives rise to a surprise-modulated trade-off between forgetting old observations and integrating them with the new ones. The modulation depends on a probability ratio, which we call “Bayes Factor Surprise”, that tests the prior belief against the current belief. We demonstrate that in several existing approximate algorithms the Bayes Factor Surprise modulates the rate of adaptation to new observations. We derive three novel surprised-based algorithms, one in the family of particle filters, one in the family of variational learning, and the other in the family of message passing, that have constant scaling in observation sequence length and particularly simple update dynamics for any distribution in the exponential family. Empirical results show that these surprise-based algorithms estimate parameters better than alternative approximate approaches and reach levels of performance comparable to computationally more expensive algorithms. The Bayes Factor Surprise is related to but different from Shannon Surprise. In two hypothetical experiments, we make testable predictions for physiological indicators that dissociate the Bayes Factor Surprise from Shannon Surprise. The theoretical insight of casting various approaches as surprise-based learning, as well as the proposed online algorithms, may be applied to the analysis of animal and human behavior, and to reinforcement learning in non-stationary environments.

1 Introduction

Animals, humans, and similarly reinforcement learning agents may safely assume that the world is stochastic and stationary during some intervals of time interrupted by change points. The position of leafs on a tree, a stock market index, or the time it takes to travel from A to B in a crowded city is often well captured by stationary stochastic processes for extended periods of time. Then sudden changes may happen, such that the distribution of leaf positions becomes different due to a storm, the stock market index is affected by the enforcement of a new law, or a blocked road causes additional traffic jams. The violation of an agent’s expectation caused by such sudden changes is perceived by the agent as surprise, which can be seen as a measure of how much the agent’s current belief differs from reality.

Surprise, with its physiological manifestations in pupil dilation (Preuschoff \BOthers., \APACyear2011; Nassar \BOthers., \APACyear2012) and EEG signals (Modirshanechi \BOthers., \APACyear2019; Ostwald \BOthers., \APACyear2012; Mars \BOthers., \APACyear2008), is believed to modulate learning, potentially through the release of specific neurotransmitters (Yu \BBA Dayan, \APACyear2005; Gerstner \BOthers., \APACyear2018), so as to allow animals and humans to adapt quickly to sudden changes. The quick adaptation to novel situations has been demonstrated in a variety of learning experiments (Nassar \BOthers., \APACyear2012, \APACyear2010; Behrens \BOthers., \APACyear2007; Yu \BBA Dayan, \APACyear2005; Glaze \BOthers., \APACyear2015; Heilbron \BBA Meyniel, \APACyear2019). The bulk of computational work on surprise-based learning can be separated into two groups. Studies in the field of computational neuroscience have focused on biological plausibility with little emphasis on the accuracy of learning (Nassar \BOthers., \APACyear2012; Yu \BBA Dayan, \APACyear2005; Nassar \BOthers., \APACyear2010; Faraji \BOthers., \APACyear2018; Friston \BOthers., \APACyear2017; Schwartenbeck \BOthers., \APACyear2013; Friston, \APACyear2010; Behrens \BOthers., \APACyear2007; Bogacz, \APACyear2017; Ryali \BOthers., \APACyear2018), whereas exact and approximate Bayesian online methods (Adams \BBA MacKay, \APACyear2007; Fearnhead \BBA Liu, \APACyear2007) for change point detection and parameter estimation have been developed without any focus on biological plausibility (Aminikhanghahi \BBA Cook, \APACyear2017; Wilson \BOthers., \APACyear2010; Cummings \BOthers., \APACyear2018; Lin \BOthers., \APACyear2017; Masegosa \BOthers., \APACyear2017).

In this work, we take a top-down approach to surprise-based learning. We start with a generative model of change points similar to the one that has been the starting point of multiple experiments (Nassar \BOthers., \APACyear2012, \APACyear2010; Behrens \BOthers., \APACyear2007; Yu \BBA Dayan, \APACyear2005; Glaze \BOthers., \APACyear2015; Heilbron \BBA Meyniel, \APACyear2019; Findling \BOthers., \APACyear2019). We demonstrate that Bayesian inference on such a generative model can be interpreted as modulation of learning by surprise; we show that this modulation leads to a natural definition of surprise which is different, but closely related to Shannon Surprise (Shannon, \APACyear1948). Moreover, we derive three novel approximate online algorithms with update rules that inherit the surprise-modulated adaptation rate of exact Bayesian inference. The overall goal of the present study is to give a Bayesian interpretation for surprise-based learning in the brain, and to find approximate methods that are computationally efficient and biologically plausible while maintaining the learning accuracy at a high level. As a by-product, our approach provides theoretical insights on commonalities and differences among existing surprise-based and approximate Bayesian approaches. Importantly, our approach makes specific experimental predictions.

In the Results section, we first introduce the generative model, and then we present our surprise-based interpretation of Bayesian inference and our three approximate algorithms. Next, we use simulations to compare our algorithms with existing ones on two different tasks inspired by and closely related to real experiments (Nassar \BOthers., \APACyear2010, \APACyear2012; Behrens \BOthers., \APACyear2007; Mars \BOthers., \APACyear2008; Ostwald \BOthers., \APACyear2012). At the end of the Results section, we formalize two experimentally testable predictions of our theory and illustrate them with simulations. A brief review of related studies as well as a few directions for further work are supplied in the Discussion section.

2 Results

In order to study learning in an environment that exhibits occasional and abrupt changes, we consider a hierarchical generative model (Fig. 1A) in discrete time, similar to existing model environments (Nassar \BOthers., \APACyear2012, \APACyear2010; Behrens \BOthers., \APACyear2007; Yu \BBA Dayan, \APACyear2005). At each time point $t$, the observation $Y_{t}=y$ comes from a distribution with the time-invariant likelihood $P_{Y}(y|\theta)$ parameterized by $\Theta_{t}=\theta$, where both $y$ and $\theta$ can be multi-dimensional. In general, we indicate random variables by capital letters, and values by small letters. Whenever there is no risk of ambiguity, we drop the explicit notation of random variables to simplify notation. Abrupt changes of the environment correspond to sudden changes of the parameter $\theta_{t}$. At every time $t$, there is a change probability $p_{c}\in(0,1)$ for the parameter $\theta_{t}$ to be drawn from its prior distribution $\pi^{(0)}$ independently of its previous value, and a probability $1-p_{c}$ to stay the same as $\theta_{t-1}$. A change at time $t$ is specified by the event $C_{t}=1$; otherwise $C_{t}=0$. Therefore, the generative model can be formally defined, for any $T\geq 1$, as a joint probability distribution over $\Theta_{1:T}\equiv(\Theta_{1},\ldots,\Theta_{T})$, $C_{1:T}$, and $Y_{1:T}$ as

[TABLE]

where $\textbf{P}(\theta_{1})=\pi^{(0)}(\theta_{1})$, $\textbf{P}(c_{1})=\delta(c_{1}-1)$, and

[TABLE]

P stands for either probability density function (for the continuous variables) or probability mass function (for the discrete variables), and $\delta$ is the Dirac or Kronecker delta distribution, respectively.

Given a sequence of observations $y_{1:t}$, the agent’s belief $\pi^{(t)}(\theta)$ about the parameter $\theta$ at time $t$ is defined as the posterior probability distribution $\textbf{P}(\Theta_{t}=\theta|y_{1:t})$. In the online learning setting studied here, the agent’s goal is to update the belief $\pi^{(t)}(\theta)$ to the new belief $\pi^{(t+1)}(\theta)$, or an approximation thereof, upon observing $y_{t+1}$.

A simplified real-world example of such an environment is illustrated in Fig. 1B. Imagine that every day a friend of yours meets you at the coffee shop, starting after work from her office (Fig. 1B left). To do so, she needs to cross a river via a bridge. The time of arrival of your friend (i.e. $y_{t}$) exhibits some variability, due to various sources of stochasticity (e.g. traffic and your friend’s daily workload), but it has a stable average over time (i.e. $\theta_{t}$). However, if a new bridge is opened, your friend arrives earlier, since she no longer has to take detour (Fig. 1B right). The moment of opening the new bridge is indicated by $c_{t+1}=1$ in our framework, and the sudden change in the average arrival time of your friend by a sudden change from $\theta_{t}$ to $\theta_{t+1}$. Even without any explicit discussion with your friend about this situation and only by observing her actual arrival time, you can notice the abrupt change and hence adapt your schedule to the new situation.

2.1 Online Bayesian inference modulated by surprise

According to the definition of the hierarchical generative model (Fig. 1A and Eq. 1 to Eq. 6), the value $y_{t+1}$ of the observation at time $t+1$ depends only on the parameters $\theta_{t+1}$, and is (given $\theta_{t+1}$) independent of earlier observations and earlier parameter values. We exploit this Markovian property and update, using Bayes’ rule, the belief $\pi^{(t)}(\theta)\equiv\textbf{P}(\Theta_{t}=\theta|y_{1:t})$ at time $t$ to the new belief at time $t+1$

[TABLE]

So far, Eq. 7 remains rather abstract. The aim of this section is to rewrite it in the form of a surprise-modulated recursive update. The first term in the numerator of Eq. 7 is the likelihood of the current observation given the parameter $\Theta_{t+1}=\theta$, and the second term is the agent’s estimated probability distribution of $\Theta_{t+1}$ before observing $y_{t+1}$. Because there is always the possibility of an abrupt change, the second term is not the agent’s previous belief $\pi^{(t)}$, but $\textbf{P}(\Theta_{t+1}=\theta|y_{1:t})=(1-p_{c})\pi^{(t)}(\theta)+p_{c}\pi^{(0)}(\theta)$. As a result, it is possible to find a recursive formula for updating the belief. For the derivation of this recursive rule, we define the following terms.

Definition 1**.**

The probability or density (for discrete and continuous variables respectively) of observing $y$ with a belief $\pi^{(t^{\prime})}$ is denoted as

[TABLE]

Note that if $\pi$ is the exact Bayesian belief defined as above in Eq. 7, then $P(y;\pi^{(t^{\prime})})=\textbf{P}(Y_{t^{\prime}+1}=y|y_{1:t^{\prime}},c_{t^{\prime}+1}=0)$. In Section 2.2 we will use also $P(y;\hat{\pi}^{(t^{\prime})})$ for an arbitrary $\hat{\pi}^{(t^{\prime})}$. Two particularly interesting cases of Eq. 8 are $P(y_{t+1};\pi^{(t)})$, i.e. the probability of a new observation $y_{t+1}$ with the current belief $\pi^{(t)}$, and $P(y_{t+1};\pi^{(0)})$, i.e. the probability of a new observation $y_{t+1}$ with the prior belief $\pi^{(0)}$.

Definition 2**.**

The “Bayes Factor Surprise” $\textbf{S}_{\mathrm{BF}}$ of the observation $y_{t+1}$ is defined as the ratio of the probability of observing $y_{t+1}$ given $c_{t+1}=1$ (i.e. when there is a change), to the probability of observing $y_{t+1}$ given $c_{t+1}=0$ (i.e. when there is no change), i.e.

[TABLE]

This definition of surprise measures how much more probable the current observation is under the naive prior $\pi^{(0)}$ relative to the current belief $\pi^{(t)}$ (see the Discussion section for further interpretation). This probability ratio is the Bayes factor (Kass \BBA Raftery, \APACyear1995; Efron \BBA Hastie, \APACyear2016) that tests the prior belief $\pi^{(0)}$ against the current belief $\pi^{(t)}$. We emphasize that our definition of surprise is not arbitrary, but essential in order to write the exact inference in Eq. 7 on the generative model in the compact recursive form indicated in the Proposition that follows. Moreover, as we show later, this term can be identified in multiple learning algorithms (among them Nassar \BOthers. (\APACyear2010, \APACyear2012)), but it has never been interpreted as a surprise measure. In the following sections we establish the generality of this computational mechanism and identify it as a common feature of many learning algorithms.

Definition 3**.**

Under the assumption of no change $c_{t+1}=0$, and using the most recent belief $\pi^{(t)}$ as prior, the exact Bayesian update for $\pi^{(t+1)}$ is denoted as

[TABLE]

$\pi_{B}^{(t+1)}(\theta)$ describes the incorporation of the new information into the current belief via Bayesian updating.

Definition 4**.**

The “Surprise-Modulated Adaptation Rate” is a function $\gamma:\mathbb{R}^{+}\times\mathbb{R}^{+}\to[0,1]$ specified as

[TABLE]

where $\text{S}\geq 0$ is a surprise value, and $m\geq 0$ is a parameter controlling the effect of surprise on learning.

Using the above definitions and Eq. 7, we have for the generative model of Fig. 1A and Eq. 1 to Eq. 6 the following Proposition.

Proposition**.**

Exact Bayesian inference on the generative model is equivalent to the recursive update rule

[TABLE]

where $\textbf{S}_{\mathrm{BF}}^{(t+1)}=\textbf{S}_{\mathrm{BF}}(y_{t+1};\pi^{(t)})$ is the Bayes Factor Surprise and

[TABLE]

is the posterior if we take $y_{t+1}$ as the only observation.

The proposition indicates that the exact Bayesian inference on the generative model discussed above (Fig. 1) leads to an explicit trade-off between (i) integrating a new observation $y^{\text{new}}$ (corresponding to $y_{t+1}$) with the old belief $\pi^{\text{old}}$ (corresponding to $\pi^{(t)}$) into a distribution $\pi^{\text{integration}}$ (corresponding to $\pi_{B}^{(t+1)}$) and (ii) forgetting the past observations, so as to restart with the belief $\pi^{\text{reset}}$ (corresponding to $P(\theta|y_{t+1})$) which relies only on the new observation and the prior $\pi^{(0)}$

[TABLE]

This trade-off is governed by a surprise-modulated adaptation rate $\gamma(\text{S},m)\in[0,1]$, where $\text{S}=\textbf{S}_{\mathrm{BF}}\geq 0$ (corresponding to the Bayes Factor Surprise) can be interpreted as the surprise of the most recent observation, and $m=\frac{p_{c}}{1-p_{c}}\geq 0$ is a parameter controlling the effect of surprise on learning. Because the parameter of modulation $m$ is equal to $\frac{p_{c}}{1-p_{c}}$, for a fixed value of surprise $S$, the adaptation rate $\gamma$ is an increasing function of $p_{c}$. Therefore, in more volatile environments, the same value of surprise $S$ leads to a higher adaptation rate than in a less volatile environment; in the case of $p_{c}\to 1$, any surprise value leads to full forgetting, i.e. $\gamma=1$.

As a conclusion, our first main result is that a split as in Eq. 14 with a weighting factor (“adaptation rate” $\gamma$) as in Eq. 11 is exact and always possible for the class of environments defined by our hierarchical generative model. This surprise-modulation gives rise to specific testable experimental predictions discussed later.

2.2 Approximate algorithms modulated by surprise

Despite the simplicity of the recursive formula in Eq. 12, the updated belief $\pi^{(t+1)}$ is generally not in the same family of distributions as the previous belief $\pi^{(t)}$, e.g. the result of averaging two normal distributions is not a normal distribution. Hence it is in general impossible to find a simple and exact update rule for e.g. some sufficient statistic. As a consequence, the memory demands for $\pi^{(t+1)}$ scale linearly in time, and updating $\pi^{(t+1)}$ using $\pi^{(t)}$ needs $\mathcal{O}(t)$ operations. In the following sections, we investigate three approximations (Algo. 1-3) that have simple update rules and finite memory demands, so that the updated belief remains tractable over a long sequence of observations.

As our second main result, we show that all three novel approximate algorithms inherit the surprise-modulated adaptation rate from the exact Bayesian approach, i.e. Eq. 11 and Eq. 14. The first algorithm adapts an earlier algorithm of surprise minimization learning (SMiLe, Faraji \BOthers. (\APACyear2018)) to variational learning. We refer to our novel algorithm as Variational SMiLe and abbreviate it by VarSMiLe (see Algo. 1). The second algorithm is based on message passing (Adams \BBA MacKay, \APACyear2007) restricted to a finite number of messages $N$ . We refer to this algorithm as MP$N$ (see Algo. 2). The third algorithm uses the ideas of particle filtering (Gordon \BOthers., \APACyear1993) for an efficient approximation for our hierarchical generative model. We refer to our approximate algorithm as Particle Filtering with $N$ particles and abbreviate it by pf$N$ (see Algo. 3). All algorithms are computationally efficient, have finite memory demands and are biologically plausible; Particle Filtering has possible neuronal implementations (Kutschireiter \BOthers., \APACyear2017; Shi \BBA Griffiths, \APACyear2009; Huang \BBA Rao, \APACyear2014; Legenstein \BBA Maass, \APACyear2014), MP$N$ can be seen as a greedy version of pf$N$ without sampling, and Variational SMiLe may be implemented by neo-Hebbian (Lisman \BOthers., \APACyear2011; Gerstner \BOthers., \APACyear2018) update rules. Simulation results show that the performance of the three approximate algorithms is comparable to and more robust across environments than other state-of-the-art approximations.

2.2.1 Variational SMiLe Rule (Algo. 1)

A simple heuristic approximation to keep the updated belief in the same family as the previous beliefs consists in applying the weighted averaging of the exact Bayesian update rule (Eq. 12) to the logarithm of the beliefs rather than the beliefs themselves, i.e.

[TABLE]

where \gamma_{t+1}=\gamma\big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\hat{\pi}^{(t)}),m\big{)} is given by Eq. 11 with a free parameter $m>0$ which can be tuned to each environment. By doing so, we still have the explicit trade-off between two terms as in Eq. 12, but in the logarithms; yet an advantageous consequence of averaging over logarithms is that, if the likelihood function $P_{Y}$ is in the exponential family, and if the initial belief $\pi^{(0)}$ is its conjugate prior, then $\hat{\pi}^{(t+1)}$ and $\pi^{(0)}$ are members of the same family. In this particular case, we arrive at a simple update rule for the parameters of $\hat{\pi}^{(t+1)}$ (see Algorithm 1 for pseudocode and Methods for details). As it is common in variational approaches (Beal, \APACyear2003), the price of this simplicity is that, except for the trivial cases of $p_{c}=0$ and $p_{c}=1$, there is no evidence other than simulations that the update rule of Eq. 15 will end up at an approximate belief close to the exact Bayesian belief.

One way to interpret the update rule of Eq. 15 is to rewrite it as the solution of a constraint optimization problem. The new belief $\hat{\pi}^{(t+1)}$ is a variational approximation of the Bayesian update $\hat{\pi}_{B}^{(t+1)}$ (see Methods)

[TABLE]

with a family of functions $q(\theta)$ constrained by the Kullback-Leibler divergence

[TABLE]

where the bound B_{t+1}\in\big{[}0,\textbf{D}_{KL}[\hat{\pi}_{B}^{(t+1)}(\theta)||P(\theta|y_{t+1})]\big{]} is a decreasing function of the Bayes Factor surprise $\textbf{S}_{\mathrm{BF}}(y_{t+1};\hat{\pi}^{(t)})$ (see Methods for proof), and $P(\theta|y_{t+1})$ is given by Eq. 13.

Because of the similarity of the constraint optimization problem in Eq. 16 and Eq. 17 to the Surprise Minimization Learning rule “SMiLe” (Faraji \BOthers., \APACyear2018), we call this algorithm “Variational Surprise Minimization Learning” rule, or in short “Variational SMiLe” rule. The differences between SMiLe and Variational SMiLe are discussed in the Methods section.

Our variational method, and particularly its surprise-modulated adaptation rate, is complementary to earlier studies (Özkan \BOthers., \APACyear2013; Masegosa \BOthers., \APACyear2017) in machine learning which assumed different generative models and used additional assumptions and different approaches for deriving the learning rule.

2.2.2 Message-Passing $N$ (Algo. 2)

For a hierarchical generative model similar to ours, a message passing algorithm has been used to perform exact Bayesian inference (Adams \BBA MacKay, \APACyear2007), where the algorithm’s memory demands and computational complexity scale linearly in time $t$. In this section, we first explain the idea of the message passing algorithm of Adams \BBA MacKay (\APACyear2007) and its relation to our Proposition. We then present our approximate version of this algorithm which has a constant (in time) computational complexity and memory demands.

The history of change points up to time $t$ is a binary sequence, e.g. $c_{1:t}=\{1,0,0,1,0,1,1\}$, where the value $1$ indicates a change in the corresponding time step. Following the idea of Adams \BBA MacKay (\APACyear2007), we define the random variable $R_{t}=\min\{n\in\mathbb{N}:C_{t-n+1}=1\}$ in order to describe the time since the last change point, which takes values between 1 to $t$. We can write the exact Bayesian expression for $\pi^{(t)}(\theta)$ by marginalizing $\textbf{P}(\Theta_{t}=\theta,r_{t}|y_{1:t})$ over the $t$ possible values of $r_{t}$ in the following way

[TABLE]

For consistency with Algorithm 3 (i.e. Particle Filtering), we call each term in the sum of Eq. 18 a “particle”, and denote as $\pi^{(t)}_{k}(\theta)=\textbf{P}(\Theta_{t}=\theta|R_{t}=t-k,y_{1:t})$ the belief of the particle corresponding to $R_{t}=t-k$, and $w^{(k)}_{t}=\textbf{P}(R_{t}=t-k|y_{1:t})$ its corresponding weight at time $t$, i.e.

[TABLE]

For each particle, the term $\pi^{(t)}_{k}(\theta)$ is simple to compute, because when $r_{t}$ is known, inference depends only on the observations after the last change point. Therefore, the goal of online inference is to find an update rule for the evolution of the weights $w^{(k)}_{t}$ over time.

We can apply the exact update rule of our Proposition (Eq. 12) to the belief expressed in the form of Eq. 19. Upon each observation of a new sample $y_{t+1}$, a new particle is generated and added to the set of particles, corresponding to $P(\theta|y_{t+1})$ (i.e. $\pi^{reset}$), modelling the possibility of a change point occurring at $t+1$. According to the proposition, the weight of the new particle (i.e. $k=t$) is equal to

[TABLE]

where \gamma_{t+1}=\gamma\big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\pi^{(t)}),\frac{p_{c}}{1-p_{c}}\big{)} (cf. Eq. 11). The other $t$ particles coming from $\pi^{(t)}$ corresponds to $\pi_{B}^{(t+1)}$ (i.e. $\pi^{integration}$) in the proposition. The update rule (see Methods for derivation) for the weights of these particles (i.e. $0\leq k\leq t-1$) is

[TABLE]

So far, we used the idea of (Adams \BBA MacKay, \APACyear2007) to write the belief as in Eq. 19 and used our proposition to arrive at the surprise-modulated update rules in LABEL:Eq:MP_w_newparticle and LABEL:Eq:MP_w_oldparticles.

The computational complexity and memory requirements of the complete message passing algorithm increase linearly with time $t$. To deal with this issue and to have a constant computation and memory demands over time, we implemented a message passing algorithm of the form of Eq. 19 to LABEL:Eq:MP_w_oldparticles, but with a fixed number $N$ of particles, chosen as those with the highest weights $w^{(k)}_{t}$. Therefore, our second algorithm adds a new approximation step to the full message passing algorithm of Adams \BBA MacKay (\APACyear2007): Whenever $t>N$, after adding the new particle with the weight as in LABEL:Eq:MP_w_newparticle and updating the previous weights as in LABEL:Eq:MP_w_oldparticles, we discard the particle with the smallest weight (i.e. set its weight equal to 0), and renormalize the weights. By doing so, we always keep the number of particles with non-zero weights equal to $N$. Note that, for $t\leq N$, our algorithm is exact, and identical to the message passing algorithm of (Adams \BBA MacKay, \APACyear2007). We call our modification of the message passing algorithm of Adams \BBA MacKay (\APACyear2007) “Message Passing $N$” and abbreviate it by “MP$N$”.

To deal with the computational complexity and memory requirements, one may alternatively keep only the particles with weights greater than a cut-off threshold (Adams \BBA MacKay, \APACyear2007). However, such a constant cut-off leads to a varying number (smaller or equal to $t$) of particles in time. Our approximation MP$N$ can therefore be seen as a variation of the thresholding algorithm in Adams \BBA MacKay (\APACyear2007) with fixed number of particles $N$, and hence a variable cut-off threshold. The work of Fearnhead \BBA Liu (\APACyear2007) follows the same principle as Adams \BBA MacKay (\APACyear2007), but employs stratified resampling to eliminate particles with negligible weights, in order to reduce the total number of particles. Their resampling algorithm involves solving a complicated non-linear equation at each time step, which makes it unsuitable for a biological implementation. In addition, we experienced that in some cases, the small errors introduced in the resampling step of the algorithm of Fearnhead \BBA Liu (\APACyear2007) accumulated and led to a worse performance than our MP$N$ algorithm which simply keeps the $N$ particles with the highest weight at each time step.

For the case where the likelihood function $P_{Y}(y|\theta)$ is in the exponential family and $\pi^{(0)}$ is its conjugate prior, the resulting algorithm of MP$N$ has a simple update rule for the belief parameters (see Algorithm 2 and Methods for details). For the sake of comparison, we also implemented in our simulations the full message passing algorithm of Adams \BBA MacKay (\APACyear2007) with an almost zero cut-off (machine precision), which we consider as our benchmark “Exact Bayes”, as well as the stratified optimal resampling algorithm of Fearnhead \BBA Liu (\APACyear2007), called “SOR$N$”.

2.2.3 Particle Filtering (Algo. 3)

Eq. 19 demonstrates that the exact Bayesian belief $\pi^{(t)}$ can be expressed as a sum of two factors, i.e. as the marginalization of $\textbf{P}(\Theta_{t}=\theta,r_{t}|y_{1:t})$ over the time since the last change point $r_{t}$. Equivalently, one can compute the exact Bayesian belief as the marginalization of $\textbf{P}(\Theta_{t}=\theta,c_{1:t}|y_{1:t})$ over the history of change points $c_{1:t}$, i.e.

[TABLE]

The idea of our third algorithm is to approximate this expectation by particle filtering, i.e. sequential Monte Carlo sampling (Gordon \BOthers., \APACyear1993; Doucet \BOthers., \APACyear2000) from $\textbf{P}(C_{1:t}|y_{1:t})$.

We then approximate $\pi^{(t)}$ by

[TABLE]

where $\{c_{1:t}^{(i)}\}_{i=1}^{N}$ is a set of $N$ realizations (or samples) of $c_{1:t}$ (i.e. $N$ particles) drawn from a proposal distribution $Q(c_{1:t}|y_{1:t})$, $\{w_{t}^{(i)}\}_{i=1}^{N}$ are their corresponding weights at time $t$, and $\hat{\pi}_{i}^{(t)}(\theta)=\textbf{P}(\Theta_{t}=\theta|c_{1:t}^{(i)},y_{1:t})$ is the approximate belief corresponding to particle $i$.

Upon observing $y_{t+1}$, the update procedure for the approximate belief $\hat{\pi}^{(t+1)}$ of Eq. 23 includes two steps: (i) updating the weights, and (ii) sampling the new hidden state $c_{t+1}$ for each particle. The two steps are coupled together through the choice of the proposal distribution $Q$, for which, we choose the optimal proposal function (Doucet \BOthers., \APACyear2000) (see Methods). As a result, given this choice of proposal function, we show (see Methods) that the first step amounts to

[TABLE]

where \gamma_{t+1}=\gamma\big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\hat{\pi}^{(t)}),m\big{)} with $m=\frac{p_{c}}{1-p_{c}}$ (cf. Eq. 11), and $\{w_{B,t+1}^{(i)}\}_{i=1}^{N}$ are the weights corresponding to the Bayesian update $\hat{\pi}_{B}^{(t+1)}$ of Eq. 10 (see Methods). In the second step, we update each particle’s history of change points by going from the sequence $\{c_{1:t}^{(i)}\}_{i=1}^{N}$ to $\{c_{1:t+1}^{(i)}\}_{i=1}^{N}$, for which we always keep the old sequence up to time $t$, and for each particle $i$, we add a new element $c_{t+1}^{(i)}\in\{0,1\}$ representing no change $c_{1:t+1}^{(i)}=[c_{1:t}^{(i)},0]$ or change $c_{1:t+1}^{(i)}=[c_{1:t}^{(i)},1]$. Note, however, that it is not needed to keep the whole sequences $c_{1:t+1}^{(i)}$ in memory, but instead one can use $c_{t+1}^{(i)}$ to update $\hat{\pi}_{i}^{(t)}$ to $\hat{\pi}_{i}^{(t+1)}$. We sample the new element $c^{(i)}_{t+1}$ from the optimal proposal distribution $Q(c_{t+1}^{(i)}|c_{1:t}^{(i)},y_{1:t+1})$ (Doucet \BOthers., \APACyear2000), which is given by (see Methods)

[TABLE]

Interestingly, the above formulas entail the same surprise modulation and the same trade-off as proposed by the Proposition Eq. 12. For the weight update, there is a trade-off between an exact Bayesian update and keeping the value of the previous time step, controlled by a adaptation rate modulated exactly in the same way as in Eq. 12. Note that in contrast to Eq. 12, the trade-off for the particles’ weights is not between forgetting and integrating, but between maintaining the previous knowledge and integrating. However, the change probability (Eq. 25) for sampling is equal to the adaptation rate and is an increasing function of surprise. As a result, although the weights are updated less for surprising events, a higher surprise causes a higher probability for change, indicated by $c_{t+1}^{(i)}=1$, which implies forgetting, because for a particle $i$ with $c_{t+1}^{(i)}=1$, the associated belief $\hat{\pi}_{i}^{(t+1)}=\textbf{P}(\Theta_{t+1}=\theta|c_{t+1}^{(i)}=1,c_{1:t}^{(i)},y_{1:t+1})$ is equal to $\textbf{P}(\Theta_{t+1}=\theta|c_{t+1}^{(i)}=1,y_{t+1})=P(\theta|y_{t+1})$ (see Fig. 1A), which is equivalent to a reset of the belief as in Eq. 14. In other words, while in MP$N$ and the exact Bayesian inference in Proposition Eq. 12, the trade-off between integration and reset is accomplished by adding at each time step a new particle with weight $\gamma_{t+1}$, in Particle Filtering, it is accomplished via sampling. As a conclusion, the above formulas are essentially the same as the update rules of MP$N$ (c.f. LABEL:Eq:MP_w_oldparticles and LABEL:Eq:MP_w_newparticle) and have the same spirit as the recursive update of the Proposition Eq. 12.

Equations 23 and LABEL:Eq:wB can be applied to the case where the likelihood function $P_{Y}(y|\theta)$ is in the exponential family and $\pi^{(0)}$ is its conjugate prior. The resulting algorithm (Algorithm 3) has a particularly simple update rule for the belief parameters (see Methods for details).

The theory of particle filter methods is well established (Gordon \BOthers., \APACyear1993; Doucet \BOthers., \APACyear2000; Särkkä, \APACyear2013). Particle filters in simpler (Brown \BBA Steyvers, \APACyear2009) or more complex (Findling \BOthers., \APACyear2019) forms have also been employed to explain human behaviour. Here we derived a simple particle filter for the general case of generative models of Eq. 2, Eq. 5, and Eq. 6. Our main contribution is to show that the use of the optimal proposal distribution in this particle filter leads to a surprise-based update scheme.

2.2.4 Surprise-modulation as a framework for other algorithms

Other existing algorithms (Adams \BBA MacKay, \APACyear2007; Fearnhead \BBA Liu, \APACyear2007; Nassar \BOthers., \APACyear2010, \APACyear2012; Faraji \BOthers., \APACyear2018) can also be formulated in the surprise-modulation framework of Eq. 11 and Eq. 14 (see Methods). Moreover, in order to allow for a transparent discussion and for fair comparisons in simulations, we extended the algorithms of Nassar \BOthers. (\APACyear2010, \APACyear2012) to a more general setting. Here we give a brief summary of the algorithms we considered. A detailed analysis is provided in subsection “Surprise-modulation as a framework for other algorithms” in the Methods.

The algorithms of Nassar \BOthers. (\APACyear2010, \APACyear2012) were originally designed for a Gaussian estimation task (see Simulations for details of the task) with a broad uniform prior. We extended them to the more general case of Gaussian tasks with Gaussian priors, and we call our extended versions Nas10*∗* and Nas12*∗* for Nassar \BOthers. (\APACyear2010) and Nassar \BOthers. (\APACyear2012) respectively (for a performance comparison between our extended algorithms and their original versions see Supplementary Fig. A1 and Supplementary Fig. A2). Both algorithms have the same surprise-modulation as in our Proposition (Eq. 12). There are multiple interpretations of the approaches of Nas10*∗* and Nas12*∗* and links to other algorithms. One such link we identify is in relation to Particle Filtering with a single particle (pf$1$). More specifically, one can show that pf$1$ behaves in expectation similar to Nas10*∗* and Nas12*∗* (see Methods and Supplementary Material).

To summarize, the algorithms Exact Bayes and SOR$N$ come from the field of change point detection, and whereas the former has high memory demands, the latter has the same memory demands as our algorithms pf$N$ and MP$N$. The algorithms Nas10*∗, Nas12∗*, and SMiLe, on the other hand, come from the human learning literature and are more biologically oriented.

2.3 Simulations

With the goal of gaining a better understanding of different approximate algorithms, we evaluated the departure of their performance from the exact Bayesian algorithm in terms of mean squared error (MSE) of $\Theta_{t}$ (see Methods), on two tasks inspired by and closely related to real experiments (Nassar \BOthers., \APACyear2010, \APACyear2012; Behrens \BOthers., \APACyear2007; Mars \BOthers., \APACyear2008; Ostwald \BOthers., \APACyear2012; Maheu \BOthers., \APACyear2019): a Gaussian and a Categorical estimation task.

We compared our three novel algorithms VarSMiLe, Particle Filtering (pf$N$, where $N$ is the number of particles), and Message Passing with finite number of particles $N$ (MP$N$) to the online exact Bayesian Message Passing algorithm (Adams \BBA MacKay, \APACyear2007) (Exact Bayes), which yields the optimal solution with $\hat{\Theta}_{t}=\hat{\Theta}_{t}^{\text{Opt}}$. Furthermore, we included in the comparison the stratified optimal resampling algorithm (Fearnhead \BBA Liu, \APACyear2007) (SOR$N$, where $N$ is the number of particles), our variant of Nassar \BOthers. (\APACyear2010) (Nas10*∗) and of Nassar \BOthers. (\APACyear2012) (Nas12∗*), the Surprise-Minimization Learning algorithm of Faraji \BOthers. (\APACyear2018) (SMiLe), as well as a simple Leaky Integrator (Leaky - see Methods).

2.3.1 Gaussian estimation task

The task is a generalized version of the experiment of Nassar \BOthers. (\APACyear2010, \APACyear2012). The goal of the agent is to estimate the mean $\theta_{t}=\mu_{t}$ of observed samples, which are drawn from a Gaussian distribution with known variance $\sigma^{2}$, i.e. $y_{t+1}|\mu_{t+1}\sim\mathcal{N}(\mu_{t+1},\sigma^{2})$. The mean $\mu_{t+1}$ is itself drawn from a Gaussian distribution $\mu_{t+1}\sim\mathcal{N}(0,1)$ whenever the environment changes. In other words, the task is a special case of the generative model of Eq. 2, Eq. 5, and Eq. 6, with $\pi^{(0)}(\mu_{t})=\mathcal{N}(\mu_{t};0,1)$ and $P_{Y}(y_{t}|\mu_{t})=\mathcal{N}(y_{t};\mu_{t},\sigma^{2})$. An example of the task can be seen in Fig. 2A.

We simulated the task for all combinations of $\sigma\in\{0.1,0.5,1,2,5\}$ and $p_{c}\in\{0.1,0.05,0.01,\allowbreak 0.005,\allowbreak 0.001,0.0001\}$. For each combination of $\sigma$ and $p_{c}$, we first tuned the free parameter of each algorithm, i.e. $m$ for SMiLe and Variational SMiLe, the leak parameter for the Leaky Integrator, and the $p_{c}$ of Nas10*∗* and Nas12*∗*, by minimizing the MSE on three random initializations of the task. For the Particle Filter (pf$N$), the Exact Bayes, the MP$N$, and the SOR$N$ we empirically checked that the true $p_{c}$ of the environment was indeed the value that gave the best performance, and we used this value for the simulations. We evaluated the performance of the algorithms on ten different random task instances of $10^{5}$ steps each for $p_{c}\in\{0.1,0.05,0.01,0.005\}$ and $10^{6}$ steps each for $p_{c}\in\{0.001,0.0001\}$ (in order to sample more change points). Note that the parameter $\sigma$ is not estimated and its actual value is used by all algorithms except the Leaky Integrator.

In Fig. 2B we show the $\textbf{MSE}[\hat{\Theta}_{t}|R_{t}=n]$ in estimating the parameter after $n$ steps since the last change point, for each algorithm, computed over multiple changes, for two exemplar task settings. The Particle Filter with 20 particles (pf20), the VarSMiLe and the Nas12*∗* have an overall performance very close to that of the Exact Bayes algorithm (i.e. $\textbf{MSE}[\hat{\Theta}_{t}^{\text{Opt}}|R_{t}=n]$), with much lower memory requirements. VarSMiLe sometimes slightly outperforms the other two early after an environmental change (Fig. 2B, right), but shows slightly higher error values at later phases. The MP$N$ algorithm is the closest one to the optimal solution (i.e. Exact Bayes) for low $\sigma$ (Fig. 2B, left), but its performance is much worse for the case of high $\sigma$ and low $p_{c}$ (Fig. 2B, right). For the Stratisfied Optimal Resampling (SOR$N$) we observe a counter-intuitive behavior in the regime of low $\sigma$; the inclusion of more particles leads to worse performance (Fig. 2B, left). At higher $\sigma$ levels the performance of SOR20 is close to optimal and better than the MP20 in later time-steps. This may be due to the fact that the MP$N$ discards particles in a deterministic and greedy way (i.e. the one with the lowest weight), whereas for the SOR$N$ there is a component of randomness in the process of particle elimination, which may be important for environments with higher stochasticity.

For the Leaky Integrator we observe a trade-off between good performance in the transient phase and the stationary phase; a fixed leak value cannot fulfill both requirements. The SMiLe rule, by construction, never narrows its belief $\hat{\pi}(\theta)$ below some minimal value, which allows it to have a low error immediately after a change, but leads later to high errors. Its performance deteriorates for higher $\sigma$ (Fig. 2B, right). The Nas10*∗* performs well for low, but not for higher values of $\sigma$. Despite the fact that a Particle Filter with 1 particle (pf1) is in expectation similar to Nas10*∗* and Nas12*∗* (see Methods), it performs worse than these two algorithms on trial-by-trial measures. Still, it performs better than the MP1 and identically to the SOR1.

In Fig. 3A, we have plotted the average of $\textbf{MSE}[\hat{\Theta}_{t}^{\text{Opt}}]$ of the Exact Bayes algorithm over the whole simulation time for each of the considered $\sigma$ and $p_{c}$ levels. The difference between the other algorithms and this benchmark is called $\Delta\textbf{MSE}[\hat{\Theta}_{t}]$ (see Methods) and is plotted in Fig. 3C–F. All algorithms except for the SOR20 have lower average error values for low $\sigma$ and low $p_{c}$, than high $\sigma$ and high $p_{c}$. The Particle Filter pf20 and the Message Passing MP20 have the smallest difference from the optimal solution. The average error of MP20 is higher than that of pf20 for high $\sigma$ and low $p_{c}$, whereas pf20 is more robust across levels of environmental parameters. The worst case performance for pf20 is $\Delta\textbf{MSE}[\hat{\Theta}_{t}]=0.033$ for $\sigma=5$ and $p_{c}=0.0001$, and for SOR20 it is $\Delta\textbf{MSE}[\hat{\Theta}_{t}]=0.061$ for $\sigma=0.1$ and $p_{c}=0.1$. The difference between these two worst case scenarios is significant ($p\text{-value}=2.79\times 10^{-6}$, two-sample t-test, 10 random seeds for each algorithm). Next in performance is the algorithm Nas12*∗* and VarSMiLe. VarSMiLe exhibits its largest deviation from the optimal solution for high $\sigma$ and low $p_{c}$, but is still more resilient compared to the MP$N$ algorithms for this type of environments. Among the algorithms with only one unit of memory demands, i.e. pf1, MP1, SOR1, VarSMiLe, SMiLe, Leaky, Nas10*∗* and Nas12*∗, the winners are VarSMiLe and Nas12∗*. The SOR20 has low error overall, but unexpectedly high error for environmental settings that are presumably more relevant for biological agents (intervals of low stochasticity marked by abrupt changes). The simple Leaky Integrator performs well at low $\sigma$ and $p_{c}$ but deviates more from the optimal solution as these parameters increase (Fig. 3F). The SMiLe rule performs best at lower $\sigma$, i.e. in more deterministic environments.

A summary graph, where we collect the $\Delta\textbf{MSE}[\hat{\Theta}_{t}]$ across all levels of $\sigma$ and $p_{c}$, is shown in Fig. 4. We can see that pf20, Nas12*∗*, and VarSMiLe give the lowest worst case (lowest maximum value) $\Delta\textbf{MSE}[\hat{\Theta}_{t}]$ and are statistically better than the other 8 algorithms (the errorbars indicate the standard error of the mean across the ten random task instances).

2.3.2 Categorical estimation task

The task is inspired by the experiments of Behrens \BOthers. (\APACyear2007); Mars \BOthers. (\APACyear2008); Ostwald \BOthers. (\APACyear2012); Maheu \BOthers. (\APACyear2019). The goal of the agent is to estimate the occurrence probability of five possible states. Each observation $y_{t+1}\in\{1,...,5\}$ is drawn from a categorical distribution with parameters $\theta_{t+1}=\boldsymbol{p}_{t+1}$, i.e. $y_{t+1}|\boldsymbol{p}_{t+1}\sim\text{Cat}(y_{t+1};\boldsymbol{p}_{t+1})$. When there is a change $C_{t+1}=1$ in the environment, the parameters $\boldsymbol{p}_{t+1}$ are drawn from a Dirichlet distribution $\text{Dir}(s\cdot\boldsymbol{1})$, where $s\in(0,\infty)$ is the stochasticity parameter. In relation to the generative model of Eq. 2, Eq. 5, and Eq. 6 we, thus, have $\pi^{(0)}(\boldsymbol{p}_{t})=\text{Dir}(\boldsymbol{p}_{t};s\cdot\boldsymbol{1})$ and $P_{Y}(y_{t}|\boldsymbol{p}_{t})=\text{Cat}(y_{t};\boldsymbol{p}_{t})$. An illustration of this task is depicted in Fig. 5A.

We considered the combinations of stochasticity levels $s\in\{0.01,0.1,0.14,0.25,\allowbreak 1,\allowbreak 2,5\}$ and change probability levels $p_{c}\in\{0.1,0.05,0.01,0.005,0.001,0.0001\}$. The algorithms of Nassar \BOthers. (\APACyear2012, \APACyear2010) were specifically developed for a Gaussian estimation task and cannot be applied here. All other algorithms were first optimized for each combination of environmental parameters before an experiment starts, and then evaluated on ten different random task instances, for $10^{5}$ steps each for $p_{c}\in\{0.1,0.05,0.01,0.005\}$ and $10^{6}$ steps each for $p_{c}\in\{0.001,0.0001\}$. The parameter $s$ is not estimated and its actual value is used by all algorithms except the Leaky Integrator.

The Particle Filter pf20, the MP20 and the SOR20 have a performance closest to that of Exact Bayes, i.e. the optimal solution (Fig. 5B). VarSMiLe is the next in the ranking, with a behavior after a change similar to the Gaussian task. pf20 performs better for $s>2$ and MP20 performs better for $s\leq 2$ (Fig. 6). For this task the biologically less plausible SOR20 is the winner in performance and it behaves most consistently across environmental parameters. Its worst case performance is $\Delta\textbf{MSE}[\hat{\Theta}_{t}]=8.16\times 10^{-5}$ for $s=2$ and $p_{c}=0.01$, and the worst case performance for pf20 is $\Delta\textbf{MSE}[\hat{\Theta}_{t}]=0.0048$ for $s=0.25$ and $p_{c}=0.005$ ($p\text{-value}=1.148\times 10^{-12}$, two-sample t-test, 10 random seeds for each algorithm). For all the other algorithms, except for MP20, the highest deviations from the optimal solution are observed for medium stochasticity levels (Fig. 6B–F). When the environment is nearly deterministic (e.g. $s=0.001$ so that the parameter vectors $\boldsymbol{p}_{t}$ have almost all mass concentrated in one component), or highly stochastic (e.g. $s>1$ so that nearly uniform categorical distributions are likely to be sampled), these algorithms achieve higher performance, while the Particle Filter is the algorithm that is most resilient to extreme choices of the stochasticity parameter $s$. For VarSMiLe in particular, the lowest mean error is achieved for high $s$ and high $p_{c}$ or low $s$ and low $p_{c}$.

A summary graph, with the $\Delta\textbf{MSE}[\hat{\Theta}_{t}]$ across all levels of $s$ and $p_{c}$, can be seen in Fig. 7. The algorithms with the lowest “worst case” are SOR20 and pf20. The top-4 algorithms SOR20, pf20, MP20 and VarSMiLe are significantly better than the others (the errorbars indicate the standard error of the mean across the ten random task instances), whereas MP$1$ and SMiLe have the largest error with a maximum at 0.53.

2.3.3 Summary of simulation results

In summary, our simulation results of the two tasks collectively suggest that our Particle Filtering (pf$N$) and Message Passing (MP$N$) algorithms achieve a high level of performance, very close to the one of biologically less plausible algorithms with higher (Exact Bayes) and same (SOR$N$) memory demands. Moreover, their behavior is more consistent across tasks. Finally, among the algorithms with memory demands of one unit, VarSMiLe performs best.

2.3.4 Robustness against suboptimal parameter choice

In all algorithms we considered, the environment’s hyper-parameters are assumed to be known. We can distinguish between two types of hyper-parameters in our generative model:

1. the parameters of the likelihood function (e.g. $\sigma$ in the Gaussian task), and
2. the $p_{c}$ and the parameters of the conjugate prior (e.g. $s$ in the Categorical task). Hyper-parameters of the first type can be added to the parameter vector $\theta$ and be inferred with the same algorithm. However, learning the second type of hyper-parameters is not straightforward. By assuming that these hyper-parameters are learned more slowly than $\theta$, one can fine-tune them after each $n$ (e.g. 10) change points, while change points can be detected by looking at the particles for the Particle Filter and at the peaks of surprise values for VarSMiLe. Other approaches to hyper-parameter estimation can be found in George \BBA Doss (\APACyear2017); Liu \BBA West (\APACyear2001); Doucet \BBA Tadić (\APACyear2003); Wilson \BOthers. (\APACyear2010).

When the hyper-parameters are fixed, a mismatch between the assumed values and the true values is a possible source of errors. In this section, we investigate the robustness of the algorithms to a mismatch between the assumed and the actual probability of change points. To do so, we first tuned each algorithm’s parameter for an environment with a change probability $p_{c}$, and then tested the algorithms in environments with different change probabilities, while keeping the parameter fixed. For each new environment with a different change probability, we calculated the difference between the MSE of these fixed parameters and the optimal MSE, i.e. the resulting MSE for the case that the Exact Bayes’ parameter is tuned for the actual $p_{c}$.

More precisely, if we denote as $\textbf{MSE}[\hat{\Theta}_{t};p_{c}^{\prime},p_{c}]$ the MSE of an algorithm with parameters tuned for an environment with $p_{c}^{\prime}$, applied in an environment with $p_{c}$, we calculated the mean regret, defined as $\textbf{MSE}[\hat{\Theta}_{t};p_{c}^{\prime},p_{c}]-\textbf{MSE}[\hat{\Theta}_{t}^{\text{Opt}},p_{c}]$ over time; note that the second term is equal to $\textbf{MSE}[\hat{\Theta}_{t};p_{c},p_{c}]$ when the algorithm Exact Bayes is used for estimation. The lower the values and the flatter the curve of the mean regret, the better the performance and the robustness of the algorithm in the face of lacking knowledge of the environment. The slope of the curve indicates the degree of deviations of the performance as we move away from the optimally tuned parameter. We ran three random (and same for all algorithms) tasks initializations for each $p_{c}$ level.

In Fig. 8 we plot the mean regret for each algorithm for the Gaussian task for four pairs of $s$ and $p_{c}^{\prime}$ levels. For $\sigma=0.1$ and $p_{c}^{\prime}=0.04$ (Fig. 8A) the Exact Bayes and the MP20 show the highest robustness (smallest regret) and are closely followed by the pf20, VarSMiLe, and Nas12*∗* (note the regret’s small range of values). The lower the actual $p_{c}$, the higher the regret, but still the changes are very small. The curves for the SMiLe and the Leaky Integrator are also relatively flat, but the mean regret is much higher. The SOR$20$ is the least robust algorithm.

Similar observations can be made for $\sigma=0.1$ and $p_{c}^{\prime}=0.004$ (Fig. 8B). In this case, the performance of all algorithms deteriorates strongly when the actual $p_{c}$ is higher than the assumed one.

However, for $\sigma=5$ (Fig. 8C and Fig. 8D), the ranking of algorithms changes. The SOR$20$ is very robust for this level of stochasticity. The pf$20$ and MP$20$ perform similarly for $p_{c}=0.04$, but for lower $p_{c}^{\prime}$ the pf$20$ is more robust and the MP$20$ exhibits high fluctuations in its performance. The Nas12*∗* is quite robust at this $\sigma$ level. Overall for Exact Bayes, SOR$20$, pf$20$, VarSMiLe and Nas12*∗*, a mismatch of the assumed $p_{c}$ from the actual one does not deteriorate the performance dramatically for $\sigma=5$, $p_{c}^{\prime}=0.004$ (Fig. 8D). The SMiLe and the Leaky Integrator outperform the other algorithms for higher $p_{c}^{\prime}$ if $p_{c}<p_{c}^{\prime}$ (Fig. 8C). A potential reason is that the optimal behavior for the Leaky Integrator (according to the tuned parameters) is to constantly integrate new observations into its belief (i.e. to act like a Perfect Integrator) regardless of the $p_{c}^{\prime}$ level. This feature makes it blind to the $p_{c}$ and therefore very robust against the lack of knowledge of it (Fig. 8C).

In summary, most of the time, the mean regret for Exact Bayes and MP20 is less than the mean regret for pf20 and VarSMiLe. However, the variability in the mean regret for pf20 and VarSMiLe is smaller, and their curves are flatter across $p_{c}$ levels, which makes their performance more predictable. The results for the Categorical estimation task are similar to those of the Gaussian task, with the difference that the SOR$20$ is very robust for this case (Fig. 9).

2.4 Experimental prediction

It has been experimentally shown that some important behavioral and physiological indicators statistically correlate with a measure of surprise or a prediction error. Examples of such indicators are the pupil diameter (Preuschoff \BOthers., \APACyear2011; Nassar \BOthers., \APACyear2012; Joshi \BBA Gold, \APACyear2019), the amplitude of the P300, N400, and MMN components of EEG (Mars \BOthers., \APACyear2008; Ostwald \BOthers., \APACyear2012; Lieder \BOthers., \APACyear2013; Kopp \BBA Lange, \APACyear2013; Meyniel \BOthers., \APACyear2016; Modirshanechi \BOthers., \APACyear2019; Musiolek \BOthers., \APACyear2019), the amplitude of MEG in specific time windows (Maheu \BOthers., \APACyear2019), BOLD responses in fMRI (Konovalov \BBA Krajbich, \APACyear2018; Loued-Khenissi \BOthers., \APACyear2020), and reaction time (Huettel \BOthers., \APACyear2002; Meyniel \BOthers., \APACyear2016). The surprise measure is usually the negative log-probability of the observation, known as Shannon Surprise (Shannon, \APACyear1948), and denoted here as $\textbf{S}_{\mathrm{Sh}}$. However, as we show in this section, as long as there is an uninformative prior over observations, Shannon Surprise $\textbf{S}_{\mathrm{Sh}}$ is just an invertible function of our modulated adaptation rate $\gamma$ and hence an invertible function of the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$. Thus, based on the results of previous works (Nassar \BOthers., \APACyear2010, \APACyear2012; Meyniel \BOthers., \APACyear2016; Ostwald \BOthers., \APACyear2012; Modirshanechi \BOthers., \APACyear2019), that always used uninformative priors, one cannot determine whether the aforementioned physiological and behavioral indicators correlate with $\textbf{S}_{\mathrm{Sh}}$ or $\textbf{S}_{\mathrm{BF}}$.

In this section, we first investigate the theoretical differences between the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$ and Shannon Surprise $\textbf{S}_{\mathrm{Sh}}$. Then, based on their observed differences, we formulate two experimentally testable predictions, with a detailed experimental protocol. Our predictions make it possible to discriminate between the two measures of surprise, and to determine whether physiological or behavioral measurements are signatures of $\textbf{S}_{\mathrm{BF}}$ or of $\textbf{S}_{\mathrm{Sh}}$.

2.4.1 Theoretical difference between $\textbf{S}_{\mathrm{BF}}$ and $\textbf{S}_{\mathrm{Sh}}$

Shannon Surprise (Shannon, \APACyear1948) is defined as

[TABLE]

where for computing $\textbf{P}(y_{t+1}|y_{1:t})$, one should know the structure of the generative model. For the generative model of Fig. 1A, we find \textbf{S}_{\mathrm{Sh}}(y_{t+1};\pi^{(t)})=-\text{log}\Big{(}(1-p_{c})P(y_{t+1};\pi^{(t)})+p_{c}P(y_{t+1};\pi^{(0)})\Big{)}. While the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$ depends on a ratio between the probability of the new observation under the prior and the current beliefs, Shannon Surprise depends on a weighted sum of these probabilities. Interestingly, it is possible to express (see Methods for derivation) the adaptation rate $\gamma_{t+1}$ as a function of the “difference in Shannon Surprise”

[TABLE]

where \gamma_{t+1}=\gamma\big{(}\textbf{S}_{\mathrm{BF}}^{(t+1)},m\big{)} depends on the Bayes Factor Surprise and the saturation parameter $m$ (cf. Eq. 11). Equation LABEL:gamma_Shannon shows that the modulated adaptation rate is not just a function of Shannon Surprise upon observing $y_{t+1}$, but a function of the difference between the Shannon Surprise of this observation under the current and under the prior beliefs. In the next subsections, we exploit differences between $\textbf{S}_{\mathrm{BF}}$ and $\textbf{S}_{\mathrm{Sh}}$ to formulate our experimentally testable predictions.

2.4.2 Experimental protocol

Consider the variant of the Gaussian task of Nassar \BOthers. (\APACyear2010, \APACyear2012) which we used in our simulations, i.e. $P_{Y}(y|\theta)=\mathcal{N}(y;\theta,\sigma^{2})$ and $\pi^{(0)}(\theta)=\mathcal{N}(\theta;0,1)$. Human subjects are asked to predict the next observation $y_{t+1}$ given what they have observed so far, i.e. $y_{1:t}$. The experimental procedure is as follows:

Fix the hyper parameters $\sigma^{2}$ and $p_{c}$. 2. 2.

At each time $t$, show the observation $y_{t}$ (produced in the aforementioned way) to the subject, and measure a physiological or behavioral indicator $M_{t}$, e.g. pupil diameter (Nassar \BOthers., \APACyear2010, \APACyear2012). 3. 3.

At each time $t$, after observing $y_{t}$, ask the subject to predict the next observation $\hat{y}_{t+1}$ and their confidence $C_{t}$ about their prediction.

Note that the only difference between our task and the task of Nassar \BOthers. (\APACyear2010, \APACyear2012) is the choice of prior for $\theta$ (i.e. Gaussian instead of uniform). The assumption is that, according to the previous studies, there is a positive correlation between $M_{t}$ and a measure of surprise.

2.4.3 Prediction 1

Based on the results of Nassar \BOthers. (\APACyear2010, \APACyear2012), in such a Gaussian task, the best fit for subjects’ prediction $\hat{y}_{t+1}$ is $\hat{\theta}_{t}$, and the confidence $C_{t}$ is a monotonic function of $\hat{\sigma}_{t}$. In order to formalize our experimental prediction, we define, at time $t$, the prediction error as $\delta_{t}=y_{t}-\hat{y}_{t}$ and the “sign bias” as $s_{t}=\text{sign}(\delta_{t}\hat{y}_{t})$. The variable $s_{t}$ is a crucial variable for our analysis. It shows whether the prediction $\hat{y}_{t}$ is an overestimation in absolute value ($s_{t}=+1$) or an underestimation in absolute value ($s_{t}=-1$). Fig. 10A shows a schematic for the case that both the current and prior beliefs are Gaussian distributions. The two observations indicated by dashed lines have same absolute error $|\delta_{t}|$, but differ in the sign bias $s$.

Given an absolute prediction value $\hat{y}>0$, an absolute prediction error $\delta>0$, a confidence value $C>0$, and a sign bias $s\in\{-1,1\}$, we can compute the average of $M_{t}$ over time for the time points with $|\hat{y}_{t}|\approx\hat{y}$, $|\delta_{t}|\approx\delta$, $C_{t}\approx C$, and $s_{t}=s$, which we denote as $\bar{M}_{1}(\hat{y},\delta,s,C)$ – the index 1 stands for experimental prediction 1. The approximation notation $\approx$ is used for continuous variables instead of equality, due to practical limitations, i.e. for obtaining adequate number of samples for averaging. Note that for our theoretical proofs we use equality, but in our simulation we include the practical limitations of a real experiment, and hence, use an approximation. The formal definitions can be found in Methods. It is worth noting that the quantity $\bar{M}_{1}(\hat{y},\delta,s,C)$ is model independent; its calculation does not require any assumption on the learning algorithm the subject may employ. Depending on whether the measurement $\bar{M}_{1}(\hat{y},\delta,s,C)$ reflects $\textbf{S}_{\mathrm{Sh}}$ or $\textbf{S}_{\mathrm{BF}}$, its relationship to the defined four variables (i.e. $\hat{y}$, $\delta$, $s$, $C$) is qualitatively and quantitatively different.

In order to prove and illustrate our prediction, let us consider each subject as an agent enabled with one of the learning algorithms that we discussed. Similar to above, given an absolute prediction $\hat{\theta}>0$ (corresponding to the subjects’ absolute prediction $\hat{y}$), an absolute prediction error $\delta>0$, a standard deviation $\sigma_{C}$ (corresponding to the subjects’ confidence value $C$), and a sign bias $s\in\{-1,1\}$, we can compute the average Shannon Surprise $\textbf{S}_{\mathrm{Sh}}(y_{t};\hat{\pi}^{(t-1)})$ and the average Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}(y_{t};\hat{\pi}^{(t-1)})$ over time, for the time points with $|\hat{\theta}_{t-1}|\approx\hat{\theta}$, $|\delta_{t}|\approx\delta$, $\hat{\sigma}_{t}\approx\sigma_{C}$, and $s_{t}=s$, which we denote as $\bar{\textbf{S}}_{\mathrm{Sh}}(\hat{\theta},\delta,s,\sigma_{C})$ and $\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s,\sigma_{C})$ respectively. We can show theoretically (see Methods) and in simulations (see Fig. 10B and Methods) that for any value of $\hat{\theta}$, $\delta$, and $\sigma_{C}$, we have $\bar{\textbf{S}}_{\mathrm{Sh}}(\hat{\theta},\delta,s=+1,\sigma_{C})>\bar{\textbf{S}}_{\mathrm{Sh}}(\hat{\theta},\delta,s=-1,\sigma_{C})$ for the Shannon Surprise, and exactly the opposite relation, i.e. $\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s=+1,\sigma_{C})<\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s=-1,\sigma_{C})$ for the Bayes Factor Surprise. Moreover, this effect increases with increasing $\delta$.

It should be noted that such an effect is due to the essential difference of $\textbf{S}_{\mathrm{Sh}}$ and $\textbf{S}_{\mathrm{BF}}$ in using the prior belief $\pi^{(0)}(\theta)$. Our experimental prediction is theoretically provable for the cases that each subject’s belief $\hat{\pi}^{(t)}$ is a Gaussian distribution, which is the case if they employ VarSMiLe, Nas10*∗, Nas12∗*, pf1, MP1, or Leaky Integrator as their learning rule (see Methods). For the cases that different learning rules (e.g. pf20) are used, where the posterior belief is a weighted sum of Gaussians, the theoretical analysis is more complicated, but our simulations show the same results (see Fig. 10B and Methods). Therefore, independent of the learning rule, we have the same experimental prediction on the manifestation of different surprise measures on physiological signals, such as pupil dilation. Our first experimental prediction can be summarized as a set of hypotheses shown in Table 1.

2.4.4 Prediction 2

Our second prediction follows the same experimental procedure as the one for the first prediction. The main difference is that for the second prediction we need to fit a model to the experimental data. Given one of the learning algorithms, the fitting procedure can be done by tuning the free parameters of the algorithm with the goal of minimizing the mean squared error between the model’s prediction $\hat{\theta}_{t}$ and a subject’s prediction $\hat{y}_{t+1}$ (similar to Nassar \BOthers. (\APACyear2010, \APACyear2012)) or with the goal of maximizing the likelihood of subject’s prediction $\hat{\pi}^{(t)}(\hat{y}_{t+1})$. Our prediction is independent of the learning algorithm, but in an actual experiment, we recommend to use model selection to find the model that fits the human data best.

Having a fitted model, we can compute the probabilities $P(y_{t+1};\hat{\pi}^{(t)})$ and $P(y_{t+1};\hat{\pi}^{(0)})$. For the case that these probabilities are equal, i.e. $P(y_{t+1};\hat{\pi}^{(t)})=P(y_{t+1};\hat{\pi}^{(0)})=p$, the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$ is equal to 1, independent of the value of $p$ (cf. Eq. 9). However, the Shannon Surprise $\textbf{S}_{\mathrm{Sh}}$ is equal to $-\log p$, and varies with $p$. Fig. 11A shows a schematic for the case that both current and prior beliefs are Gaussian distributions. Two cases for which we have $P(y_{t+1};\hat{\pi}^{(t)})=P(y_{t+1};\hat{\pi}^{(0)})=p$, for two different $p$ values, are marked by black dots at the intersections of the curves.

Given a probability $p>0$, we can compute the average of $M_{t}$ over time for the time points with $P(y_{t+1};\hat{\pi}^{(t)})\approx p$ and $P(y_{t+1};\hat{\pi}^{(0)})\approx p$, which we denote as $\bar{M}_{2}(p)$ – the index 2 stands for experimental prediction 2. Analogous to the first prediction, the approximation notation $\approx$ is used due to practical limitations. Then, if $\bar{M}_{2}(p)$ is independent of $p$, its behavior is consistent with $\textbf{S}_{\mathrm{BF}}$, whereas if it decreases by increasing $p$, it can be a signature of $\textbf{S}_{\mathrm{Sh}}$. Our second experimental prediction can be summarized as two hypotheses shown in Table 2. Note that in contrast to our first prediction, with the assumption that the standard deviation of the prior belief is fitted using the behavioral data, we do not consider the hypothesis that the prior is not used for inference, because this is indistinguishable from a very large variance of the prior belief.

In order to illustrate the possible results and the feasibility of the experiment, we ran a simulation and computed $\bar{\textbf{S}}_{\mathrm{BF}}(p)$ and $\bar{\textbf{S}}_{\mathrm{Sh}}(p)$ for the time points with $P(y_{t+1};\hat{\pi}^{(t)})\approx p$ and $P(y_{t+1};\hat{\pi}^{(0)})\approx p$ (see Methods for details). The results of the simulation are shown in Fig. 11B.

3 Discussion

We have shown that performing exact Bayesian inference on a generative world model naturally leads to a definition of surprise and a surprise-modulated adaptation rate. We have proposed three approximate algorithms (VarSMiLe, MP$N$, and pf$N$) for learning in non-stationary environments, which all exhibit the surprise-modulated adaptation rate of the exact Bayesian approach and are biologically plausible. Empirically we observed that our algorithms achieve levels of performance comparable to approximate Bayesian methods with higher memory demands (Adams \BBA MacKay, \APACyear2007), and are more resilient across different environments compared to methods with similar memory demands (Fearnhead \BBA Liu, \APACyear2007; Nassar \BOthers., \APACyear2010, \APACyear2012; Faraji \BOthers., \APACyear2018).

Learning in a volatile environment has been studied for a long time in the fields of Bayesian learning, neuroscience, and signal processing. In the following, we discuss the biological relevance of our work, and we briefly review some of the previously developed algorithms, with particular focus on the ones that have studied environments which can be modeled with a generative model similar to the one in Fig. 1. We then discuss further our results, and propose directions for future work on surprise-based learning.

3.1 Biological interpretation

Humans are able to quickly adapt to changes (Nassar \BOthers., \APACyear2012, \APACyear2010; Behrens \BOthers., \APACyear2007), but human behaviour is also often observed to be suboptimal, compared to the normative approach of exact Bayesian inference (Mathys \BOthers., \APACyear2011; Wilson \BOthers., \APACyear2013; Nassar \BOthers., \APACyear2010; Glaze \BOthers., \APACyear2015; Prat-Carrabin \BOthers., \APACyear2020). In general, biological agents have limited resources and possibly inaccurate assumptions about hyper-parameters, yielding sub-optimal behaviour, as we also see with our algorithms whose accuracies degrade with a sub-optimal choice of hyper-parameters. Performance also deteriorates with a decreasing number of particles in the sampling-based algorithms, which might be another possible explanation of suboptimal human behavior. Previously, Particle Filtering has been shown to explain the behavior of human subjects in changing environments: Daw \BBA Courville (\APACyear2008) use a single particle, (Brown \BBA Steyvers, \APACyear2009) use a simple heuristic form of particle filtering based on direct simulation, Findling \BOthers. (\APACyear2019) combine Particle Filtering with a noisy inference, and Prat-Carrabin \BOthers. (\APACyear2020) use it for a task with temporal structure.

At the level of neuronal implementation, we do not propose a specific suggestion. However, there are several hypotheses about neural implementations of related particle filters (Kutschireiter \BOthers., \APACyear2017; Shi \BBA Griffiths, \APACyear2009; Huang \BBA Rao, \APACyear2014; Legenstein \BBA Maass, \APACyear2014), on which, a neural model of pf$N$ and – its greedy version – MP$N$ could be based. In a similar spirit, the updating scheme of Variational SMiLe may be implemented in biological neural networks (for distributions in the exponential family).

Our theoretical framework for modulation of learning by the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$ is related to the body of literature on neo-Hebbian three-factor learning rules (Lisman \BOthers., \APACyear2011; Frémaux \BBA Gerstner, \APACyear2016; Gerstner \BOthers., \APACyear2018), where a third factor indicating reward or surprise enables or modulates a synaptic change or a belief update (Yu \BBA Dayan, \APACyear2005; Yu, \APACyear2012). We have shown how Bayesian or approximate Bayesian inference naturally leads to such a third factor that modulates learning via the surprise modulated adaptation rate $\gamma(\textbf{S}_{\mathrm{BF}},m)$. This may offer novel interpretations of behavioural and neurophysiological data, and help in understanding how three-factor learning computations may be implemented in the brain.

3.2 Related work

Exact Bayesian inference

As already described in the “Message-Passing $N$” section of the Results, for the generative model in Fig. 1, it is possible to find an exact online Bayesian update of the belief using a message passing algorithm (Adams \BBA MacKay, \APACyear2007). The space and time complexity of the algorithm increases linearly with $t$, which makes it unsuitable for an online learning setting. However, approximations like dropping messages below a certain threshold (Adams \BBA MacKay, \APACyear2007) or stratified resampling (Fearnhead \BBA Liu, \APACyear2007) allow to reduce the computational complexity. The former has a variable number of particles in time, and the latter needs solving a complicated non-linear equation at each time step in order to reduce the number of particles to $N$ (called SOR$N$ in the Results section).

Our message passing algorithm with finite number of particles (messages) $N$ (MP$N$, Algo. 3) is closely related to these algorithms and can be seen as a biologically more plausible variant of the other two. All three algorithms have the same update rules given by LABEL:Eq:MP_w_oldparticles and LABEL:Eq:MP_w_newparticle. Hence the algorithms of both Adams \BBA MacKay (\APACyear2007) and Fearnhead \BBA Liu (\APACyear2007) have the same surprise modulation as our MP$N$. The difference lies in their approaches to eliminate less “important” particles.

In the literature of switching state-space models (Barber, \APACyear2012), the generative models of the kind in Fig. 1 are known as “reset models”, and the message passing algorithm of Adams \BBA MacKay (\APACyear2007) is known to be the standard algorithm for inference over these models (Barber, \APACyear2012). See Barber (\APACyear2012, \APACyear2006); Ghahramani \BBA Hinton (\APACyear2000) for other variations of switching state-space models and examples of approximate inference over them.

Leaky integration and variations of delta-rules

In order to estimate some statistics, leaky integration of new observations is a particularly simple form of a trade-off between integrating and forgetting. After a transient phase, the update of a leaky integrator takes the form of a delta-rule that can be seen as an approximation of exact Bayesian updates (Meyniel \BOthers., \APACyear2016; Heilbron \BBA Meyniel, \APACyear2019; Yu \BBA Cohen, \APACyear2009; Ryali \BOthers., \APACyear2018). This update rule was found to be biologically plausible and consistent with human behavioral data (Meyniel \BOthers., \APACyear2016; Yu \BBA Cohen, \APACyear2009). However, Behrens \BOthers. (\APACyear2007) and Heilbron \BBA Meyniel (\APACyear2019) demonstrated that in some situations, the exact Bayesian model is significantly better than leaky integration in explaining human behavior. The inflexibility of leaky integration with a single, constant leak parameter can be overcome by a weighted combination of multiple leaky integrators (Wilson \BOthers., \APACyear2013), where the weights are updated in a similar fashion as in the exact online methods (Adams \BBA MacKay, \APACyear2007; Fearnhead \BBA Liu, \APACyear2007), or by considering an adaptive leak parameter (Nassar \BOthers., \APACyear2012, \APACyear2010). We have shown that the two algorithms of Nassar \BOthers. (\APACyear2010, \APACyear2012) can be generalized to Gaussian prior beliefs (Nas10*∗* and Nas12*∗*). Our results show that these algorithms also inherit the surprise-modulation of the exact Bayesian inference. Our surprise-dependent adaptation rate $\gamma$ can be interpreted as a surprise-modulated leak parameter.

Other approaches

Learning in the presence of abrupt changes has also been considered without explicit assumptions about the underlying generative model. One approach uses a surprise-modulated adaptation rate (Faraji \BOthers., \APACyear2018) similar to Eq. 11. The Surprise-Minimization Learning (SMiLe) algorithm of Faraji \BOthers. (\APACyear2018) has an updating rule similar to the one of VarSMiLe (Eq. 16 and Eq. 17). The adaptation rate modulation, however, is based on the Confidence Corrected Surprise (Faraji \BOthers., \APACyear2018) rather than the Bayes Factor Surprise, and the trade-off in its update rule is between resetting and staying with the latest belief rather than between resetting and integrating (see Methods).

Other approaches use different generative models, such as conditional sampling of the parameters also when there is a change (Yu \BBA Dayan, \APACyear2005; Glaze \BOthers., \APACyear2015), a deeper hierarchy without fixed change probability $p_{c}$ (Wilson \BOthers., \APACyear2010), or drift in the parameters (Mathys \BOthers., \APACyear2011; Gershman \BOthers., \APACyear2014). A recent work shows that inference on the generative model of Fig. 1 can explain human behavior well even when the true generative model of the environment is different and more complicated (Findling \BOthers., \APACyear2019). They develop a heuristic approach to add noise in the inference process of a Particle Filter. Their algorithm can be interpreted as a surprise-modulated Particle Filter, where the added noise scales with a measure of surprise (conceptually equivalent to Bayesian surprise (Storck \BOthers., \APACyear1995; Schmidhuber, \APACyear2010; Itti \BBA Baldi, \APACyear2006)). Moreover, another recent work (Prat-Carrabin \BOthers., \APACyear2020) shows that approximate sampling algorithms (like Particle Filtering) can explain human behavior better than their alternatives in tasks closely related to the generative model of Fig. 1. The signal processing literature provides further methods to address the problem of learning in non-stationary environments with abrupt changes; see Aminikhanghahi \BBA Cook (\APACyear2017) for a review, and Lin \BOthers. (\APACyear2017); Cummings \BOthers. (\APACyear2018); Özkan \BOthers. (\APACyear2013); Masegosa \BOthers. (\APACyear2017) for a few recent examples.

3.3 Surprise-modulation as a generic phenomenon

Learning rate modulation similar to the one in Eq. 11 has been previously proposed in the neuroscience literature with either heuristic arguments (Faraji \BOthers., \APACyear2018) or with Bayesian arguments for a particular experimental task, e.g. when samples are drawn from a Gaussian distribution (Nassar \BOthers., \APACyear2010, \APACyear2012). The fact that the same form of modulation is at the heart of Bayesian inference for our relatively general generative model, that it is derived without any further assumptions, and is not a-priori defined is in our view an important contribution to the field of adaptive learning algorithms in computational neuroscience.

Furthermore, the results of our three approximate methods (Particle Filtering, Variational SMiLe, and Message Passing with fixed $N$ number of messages) as well as some previously developed ones (Adams \BBA MacKay, \APACyear2007; Fearnhead \BBA Liu, \APACyear2007; Nassar \BOthers., \APACyear2010, \APACyear2012) demonstrate that the surprise-based modulation of the learning rate is a generic phenomenon. Therefore, regardless of whether the brain uses Bayesian inference or an approximate algorithm (Mathys \BOthers., \APACyear2011; Friston, \APACyear2010; Nassar \BOthers., \APACyear2010, \APACyear2012; Gershman \BOthers., \APACyear2014; Bogacz, \APACyear2017, \APACyear2019; Gershman, \APACyear2019; Findling \BOthers., \APACyear2019; Prat-Carrabin \BOthers., \APACyear2020), the notion of Bayes Factor Surprise and the way it modulates learning (i.e. Eq. 14 and Eq. 11) look generic.

The generality of the way surprise should modulate learning depends on an agent’s inductive biases about its environment and is directly associated with the assumed generative model of the world. The generative model we considered in this work involves abrupt changes. However, one can think of other realistic examples, where an improbable observation does not indicate a persistent change, but a singular event or an outlier, similar to d’Acremont \BBA Bossaerts (\APACyear2016); Nassar \BOthers. (\APACyear2019). In such situations, the belief should not be changed and surprise should attenuate learning, rather than accelerate it. Interestingly, we can show that exact and approximate Bayesian inference on such a generative model naturally lead to a surprise-modulated adaptation rate $\gamma(\textbf{S}_{\mathrm{BF}},m)$, with the same definition of $\textbf{S}_{\mathrm{BF}}$, where the trade-off is not between integrating and resetting, but between integrating and ignoring the new observation 222Liakoni, V., Lehmann, M. P., Modirshanechi, A., Brea, J., Herzog, M. H., Gerstner, W., & Preuschoff, K. (in preparation). Dissociating brain regions encoding reward prediction error and surprise.. This extends previous work on such environments (d’Acremont \BBA Bossaerts, \APACyear2016; Nassar \BOthers., \APACyear2019) to a general setting and highlights the general principle of surprise-based modulation, given the prior knowledge on the structure of the environment.

An aspect that the generative model we considered does not capture is the potential return to a previous state of the environment, rather than a change to a completely new situation. If in our example of Fig. 1B, the bridge with the shortest path is temporarily closed for repairs, your friend would again have to take the longer detour, therefore, her arrival times will return to their previous values, i.e. increase. In such cases, an agent should infer whether the surprising observation stems from a new hidden state or from an old state stored in memory. Relevant generative models have been studied in Fox \BOthers. (\APACyear2011); Gershman \BOthers. (\APACyear2014, \APACyear2017) and are out of the scope of our present paper.

3.4 Bayes Factor Surprise as a novel measure of surprise

In view of a potential application in the neurosciences, a definition of surprise should exhibit two properties: (i) surprise should reflect how unexpected an event is, and, (ii) surprise should modulate learning. Surprising events indicate that our belief is far from the real world and suggest to update our model of the world, or, for large surprise, simply forget it. Forgetting is the same as returning to the prior belief. However, an observation $y_{t+1}$ can be unexpected under both the prior $\pi^{(0)}$ and the current beliefs $\pi^{(t)}$. In these situations, it is not obvious whether forgetting helps. Therefore, the modulation between forgetting or not should be based on a comparison between the probability of an event under the current belief $P(y_{t+1};\pi^{(t)})$ and its probability under the prior belief $P(y_{t+1};\pi^{(0)})$.

The definition of the Bayes Factor Surprise $\textbf{S}_{\mathrm{BF}}$ as the ratio of $P(y_{t+1};\pi^{(t)})$ and $P(y_{t+1};\pi^{(0)})$ exploits this insight. The Bayes Factor Surprise appears as a modulation factor in the recursive form of the exact Bayesian update rule for a hierarchical generative model of the environment. When two events are equally probable under the prior belief, the one which is less expected under the current belief is more surprising - satisfying the first property. At the same time, when two events are equally probable under the current belief, the one which is more expected under the prior belief is more surprising - signaling that forgetting may be beneficial.

$\textbf{S}_{\mathrm{BF}}$ can be written (using Eq. 8) in a more explicit way as

[TABLE]

Note that the definition by itself is independent of the specific form of the generative model. In other words, even in the cases where data is generated with another generative model (e.g. the real world), $\textbf{S}_{\mathrm{BF}}$ could be a candidate surprise measure in order to interpret brain activity or pupil dilation.

We formally discussed the connections between the Bayes Factor Surprise and Shannon Surprise (Shannon, \APACyear1948), and showed that they are closely linked. We showed that the modulated adaptation rate ($\gamma$) used in (approximate) Bayesian inference is a function of the difference between the Shannon Surprise under the current and the prior beliefs, but cannot be expressed solely by the Shannon Surprise under the current one. Our formal comparisons between these two different measures of surprise lead to specific experimentally testable predictions.

The Bayesian Surprise $\textbf{S}_{\mathrm{Ba}}$ (Storck \BOthers., \APACyear1995; Schmidhuber, \APACyear2010; Itti \BBA Baldi, \APACyear2006) and the Confidence Corrected Surprise $\textbf{S}_{\mathrm{CC}}$ (Faraji \BOthers., \APACyear2018) are two other measures of surprise in neuroscience. The learning modulation derived in our generative model cannot be expressed as a function of $\textbf{S}_{\mathrm{Ba}}$ and $\textbf{S}_{\mathrm{CC}}$. However, one can hypothesize that $\textbf{S}_{\mathrm{Ba}}$ is computed after the update of the belief to measure the information gain of the observed event, and is therefore not a good candidate for online learning modulation. The Confidence Corrected surprise $\textbf{S}_{\mathrm{CC}}$ takes into account the shape of the belief, and therefore includes the effects of confidence, but it does not consider any information about the prior belief. Hence, a result of $\bar{M}_{1}(\hat{\theta},\delta,s=+1,C)=\bar{M}_{1}(\hat{\theta},\delta,s=-1,C)$ in our first experimental prediction would be consistent with the corresponding behavioral or physiological indicator reflecting the $\textbf{S}_{\mathrm{CC}}$.

3.5 Difference in Shannon Surprise, an alternative perspective

Following our formal comparison in the “Experimental prediction” section, $\textbf{S}_{\mathrm{BF}}$ can be expressed as a deterministic function of the difference in Shannon Surprise as

[TABLE]

All of our theoretical results can be rewritten by replacing $\textbf{S}_{\mathrm{BF}}$ with this function of $\Delta\textbf{S}_{\mathrm{Sh}}$. Moreover, because there is a 1-to-1 mapping between $\textbf{S}_{\mathrm{BF}}$ and $\Delta\textbf{S}_{\mathrm{Sh}}$, from a systemic point of view, it is not possible to specify whether the brain computes the former or the latter by analysis of behavioral data and biological signals. This suggests an alternative interpretation of surprise-modulated learning as an approximation of Bayesian inference: What the brain computes and perceives as surprise or prediction error may be Shannon Surprise, but the modulating factor in a three-factor synaptic plasticity rule (Lisman \BOthers., \APACyear2011; Frémaux \BBA Gerstner, \APACyear2016; Gerstner \BOthers., \APACyear2018) may be implemented by comparing the Shannon Surprise values under the current and the prior beliefs.

3.6 Future directions

A natural continuation of our study is to test our experimental predictions in human behavior and physiological signals, in order to investigate which measures of surprise are used by the brain. Along a similar direction, our approximate learning algorithms can be evaluated on human behavioral data from experiments that use a similar generative model (Nassar \BOthers., \APACyear2012, \APACyear2010; Wilson \BOthers., \APACyear2013; Behrens \BOthers., \APACyear2007; Yu \BBA Dayan, \APACyear2005; Glaze \BOthers., \APACyear2015; Heilbron \BBA Meyniel, \APACyear2019) in order to assess if our proposed algorithms achieve similar or better performance in explaining data.

Finally, our methods can potentially be applied to model-based reinforcement learning in non-stationary environments. In recent years, there has been a growing interest in adaptive or continually learning agents in changing environments in the form of Continual learning and Meta-learning (Lomonaco \BOthers., \APACyear2019; Traoré \BOthers., \APACyear2019). Many Continual learning model-based approaches make use of some procedure to detect changes (Lomonaco \BOthers., \APACyear2019; Nagabandi \BOthers., \APACyear2018). Integrating $\textbf{S}_{\mathrm{BF}}$ and a learning rate $\gamma(\textbf{S}_{\mathrm{BF}})$ into a reinforcement learning agent would be an interesting future direction.

4 Methods

4.1 Proof of the proposition

By definition

[TABLE]

We exploit the Markov property of the generative model in Eq. 2, Eq. 5, and Eq. 6, condition on the fixed past $y_{1:t}$ and rewrite

[TABLE]

By marginalization over the hidden state $C_{t+1}$, the second factor in the numerator of Eq. 31 can be written as

[TABLE]

The denominator in Eq. 31 can be written as

[TABLE]

where we used the definition in Eq. 8. Using these two expanded forms, Eq. 31 can be rewritten

[TABLE]

We define $P(\theta|y_{t+1})$ as the posterior given a change in the environment as

[TABLE]

Then, we can write Eq. 34 as

[TABLE]

where $\pi^{(t+1)}_{B}(\theta)$ is defined in Eq. 10, and

[TABLE]

with $\textbf{S}_{\mathrm{BF}}$ defined in Eq. 9, and $\gamma(\text{S},m)$ defined in Eq. 11. Thus our calculation yields a specific choice of surprise ($\text{S}=\textbf{S}_{\mathrm{BF}}$) and a specific value for the saturation parameter $m=\frac{p_{c}}{1-p_{c}}$.

4.2 Derivation of the optimization-based formulation of VarSMiLe (Algo. 1)

To derive the optimization-based update rule for the Variational SMiLe rule and the relation of the bound $B_{t+1}$ with surprise, we used the same approach used in Faraji \BOthers. (\APACyear2018).

Derivation of the update rule

Consider the general form of the following variational optimization problem:

[TABLE]

where B\in\big{[}0,\textbf{D}_{KL}[p_{1}(\theta)||p_{2}(\theta)]\big{]}. On the extremes of $B$, we will have trivial solutions

[TABLE]

Note that the Kullback–Leibler divergence is a convex function with respect to its first argument, i.e. $q$ in our setting. Therefore, both the objective function and the constraints of the optimization problem in Eq. 38 are convex. For convenience, we assume that the parameter space for $\theta$ is discrete, but the final results can be generalized also to the continuous case with some considerations - see Beal (\APACyear2003) and Faraji \BOthers. (\APACyear2018). For the discrete setting, the optimization problem in Eq. 38 can be rewritten as

[TABLE]

For solving the mentioned problem, one should find a $q$ which satisfies the Karush–Kuhn–Tucker (KKT) conditions (Boyd \BBA Vandenberghe, \APACyear2004) for

[TABLE]

[TABLE]

where $\lambda$ and $\alpha$ are the parameters of the dual problem. Defining $\gamma=\frac{\lambda}{1+\lambda}$, and considering the partial derivative to be zero, we have

[TABLE]

where $\alpha$ is always specified in a way to have $\text{Const}(\alpha,\gamma)$ as the normalization factor

[TABLE]

According to the KKT conditions, $\lambda\geq 0$, and as a result $\gamma\in[0,1]$. Therefore, considering $p_{1}(\theta)=\hat{\pi}_{B}^{(t+1)}(\theta)$ and $p_{2}(\theta)=P(\theta|y_{t+1})$, the solution to the optimization problem of Eq. 16 and Eq. 17 is Eq. 15.

Proof of the claim that $B$ is a decreasing function of Surprise

According to the KKT conditions

[TABLE]

For the case that $\lambda\neq 0$ (i.e. $\gamma\neq 0$), we have $B$ as a function of $\gamma$

[TABLE]

Now, we show that the derivate of $B(\gamma)$ with respect to $\gamma$ is always non-positive. To do so, we first compute the derivative of $Z(\gamma)$ as

[TABLE]

and the derivate of $\mathbb{E}_{q^{*}}[h(\theta)]$ for an arbitrary $h(\theta)$ as

[TABLE]

Using the last three equations, we have

[TABLE]

which means that $B$ is a decreasing function of $\gamma$. Because $\gamma$ is an increasing function of surprise, $B$ is also a decreasing function of surprise.

4.3 Derivations of Message Passing $N$ (Algo. 2)

For the sake of clarity and coherence, we repeat here some steps performed in the Results section.

Following the idea of Adams \BBA MacKay (\APACyear2007) let us first define the random variable $R_{t}=\min\{n\in\mathbb{N}:C_{t-n+1}=1\}$. This is the time window from the last change point. Then the exact Bayesian form for $\pi^{(t)}(\theta)$ can be written as

[TABLE]

To have a formulation similar to the one of Particle Filtering we rewrite the belief as

[TABLE]

where $\pi^{(t)}_{k}(\theta)=\textbf{P}(\Theta_{t}=\theta|R_{t}=t-k,y_{1:t})$ is the term corresponding to $R_{t}=t-k$, and $w^{(k)}_{t}=\textbf{P}(R_{t}=t-k|y_{1:t})$ is its corresponding at time $t$.

To update the belief after observing $y_{t+1}$, one can use the exact Bayesian recursive formula Eq. 12, for which one needs to compute $\pi_{B}^{(t+1)}(\theta)$ as

[TABLE]

Using Bayes’ rule and the conditional independence of observations, we have

[TABLE]

and once again, by using the Bayes’ rule and the conditional independence of observations, we find

[TABLE]

This gives us

[TABLE]

and finally

[TABLE]

Using the recursive formula, the update rule for the weights for $0\leq k\leq t-1$ is

[TABLE]

and for the newly added particle $t$

[TABLE]

where \gamma_{t+1}=\gamma\big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\pi^{(t)}),m=\frac{p_{c}}{1-p_{c}}\big{)} of Eq. 11.

The MP$N$ algorithm uses Eq. 51, LABEL:Eq:SMP_w_oldparticles, and LABEL:Eq:SMP_w_newparticle for computing the belief for $t\leq N$ - which is same as the exact Bayesian inference. For $t>N$, it first updates the weights in the same fashion as LABEL:Eq:SMP_w_oldparticles and LABEL:Eq:SMP_w_newparticle, keeps the greatest $N$ weights, and sets the rest weights equal to 0. After normalizing the new weights, it uses Eq. 51 (but only over the particles with non-zero weights) to compute the belief $\hat{\pi}^{(t)}$. For the particular case of exponential family, see Algorithm 2 for the pseudocode.

4.4 Derivation of the weight update for Particle Filtering (Algo. 3)

We derive here the weight update for the particle filter. The difference in our formalism from a standard derivation (Särkkä, \APACyear2013) is the absence of the Markov property of conditional observations (i.e. $\textbf{P}(y_{t+1}|c_{1:t+1},y_{1:t})\neq\textbf{P}(y_{t+1}|c_{t+1})$). Our goal is to perform the approximation

[TABLE]

Given a proposal sampling distribution $Q$, for the weight of particle $i$ at time $t+1$ we have

[TABLE]

where the only assumption for the proposal distribution $Q$ is that the previous hidden states $c_{1:t}^{(i)}$ are independent of the next observation $y_{t+1}$, which allows to keep the previous samples $c_{1:t}^{(i)}$ when going from $c_{1:t}^{(i)}$ to $c_{1:t+1}^{(i)}$ and to write the update of the weights in a recursive way (Särkkä, \APACyear2013).

Notice that $w_{t}^{(i)}\propto\frac{\textbf{P}(c_{1:t}^{(i)}|y_{1:t})}{Q(c_{1:t}^{(i)}|y_{1:t})}$ are the weights calculated at the previous time step. Therefore

[TABLE]

For the choice of $Q$, we use the optimal proposal function in terms of variance of the weights (Doucet \BOthers., \APACyear2000)

[TABLE]

Using Bayes’ rule and Eq. 61 and Eq. 62, after a few steps of algebra, we have

[TABLE]

Using the definition in Eq. 8, we have $\textbf{P}(y_{t+1}|c_{1:t}^{(i)},y_{1:t},c^{(i)}_{t+1}=0)=P(y_{t+1};\hat{\pi}_{i}^{(t)})$ and $\textbf{P}(y_{t+1}|c_{1:t}^{(i)},y_{1:t},c^{(i)}_{t+1}=1)=P(y_{t+1};\pi^{(0)})$. Therefore, we have

[TABLE]

where $Z$ is the normalization factor

[TABLE]

where we have

[TABLE]

We now compute the weights corresponding to $\pi_{B}^{(t+1)}$ as defined in Eq. 10

[TABLE]

Combining Eq. 64, Eq. 65 and Eq. 67 we can then re-write the weight update rule as

[TABLE]

where \gamma_{t+1}=\gamma\Big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\hat{\pi}^{(t)}),\frac{p_{c}}{1-p_{c}}\Big{)} of Eq. 11.

At every time step $t+1$ we sample each particle’s hidden state $c_{t+1}$ from the proposal distribution. Using Eq. 62, we have

[TABLE]

We implemented the Sequential Importance Resampling algorithm (Gordon \BOthers., \APACyear1993; Doucet \BOthers., \APACyear2000), where the particles are resampled when their effective number falls below a threshold. The effective number of the particles is defined as (Doucet \BOthers., \APACyear2000; Särkkä, \APACyear2013)

[TABLE]

When $N_{\text{eff}}$ is below a critical threshold, the particles are resampled with replacement from the categorical distribution defined by their weights, and all their weights are set to $w_{t}^{(i)}=1/N$. We did not optimize the parameter $N_{\text{eff}}$, and following Doucet \BBA Johansen (\APACyear2009), we performed resampling when $N_{\text{eff}}\leq N/2$.

4.5 Surprise-modulation as a framework for other algorithms

4.5.1 SMiLe Rule

The Confidence Corrected Surprise (Faraji \BOthers., \APACyear2018) is

[TABLE]

where $\tilde{P}(\theta|y_{t+1})$ is the scaled likelihood defined as

[TABLE]

Note that this becomes equal to $P(\theta|y_{t+1})$ if the prior belief $\pi^{(0)}$ is a uniform distribution; cf. Eq. 13.

With the aim of minimizing the Confidence Corrected Surprise by updating the belief during time, Faraji \BOthers. (\APACyear2018) suggested an update rule solving the optimization problem

[TABLE]

where B_{t+1}\in\big{[}0,\textbf{D}_{KL}[P(\theta|y_{t+1})||\hat{\pi}^{(t)}(\theta)]\big{]} is an arbitary bound. The authors showed that the solution to this optimization problem is

[TABLE]

where $\gamma_{t+1}\in[0,1]$ is specified so that it satisfies the constraint in Eq. 73.

Although Eq. 74 looks very similar to Eq. 15, it signifies a trade-off between the latest belief $\hat{\pi}^{(t)}$ and the belief updated by only the most recent observation $\tilde{P}(\theta|y_{t+1})$, i.e. a trade-off between adherence to the current belief and reset. While SMiLe adheres to the current belief $\hat{\pi}^{(t)}$, Variational SMiLe integrates the new observation with the current belief to get $\hat{\pi}_{B}^{(t)}$, which leads to a trade-off similar to the one of the exact Bayesian inference (Eq. 14 and Eq. 12).

To modulate the learning rate by surprise, Faraji \BOthers. (\APACyear2018) considered the boundary $B_{t+1}$ as a function of the Confidence Corrected Surprise, i.e.

[TABLE]

where $m$ is a free parameter. Then, $\gamma_{t+1}$ is found by satisfying the constraint of the optimization problem in Eq. 73 using Eq. 74 and Eq. 75.

4.5.2 Nassar’s algorithm

For the particular case that observations are drawn from a Gaussian distribution with known variance and unknown mean, i.e. $y_{t+1}|\mu_{t+1}\sim\mathcal{N}(\mu_{t+1},\sigma^{2})$ and $\theta_{t}=\mu_{t}$, Nassar \BOthers. (\APACyear2012, \APACyear2010) considered the problem of estimating the expected $\mu_{t}$ and its variance rather than a probability distribution (i.e. belief) over it, implicitly assuming that the belief is always a Gaussian distribution. The algorithms of Nassar \BOthers. (\APACyear2012, \APACyear2010) were developed for the case that, whenever the environment changes, the mean $\mu_{t+1}$ is drawn from a uniform prior with a range of values much larger than the width of the Gaussian likelihood function. The authors showed that in this case, the expected $\mu_{t+1}$ (i.e. $\hat{\mu}_{t+1}$) estimated by the agent upon observing a new sample $y_{t+1}$ is

[TABLE]

with $\alpha_{t+1}$ the adaptive learning rate given by

[TABLE]

where $\hat{r_{t}}$ is the estimated time since the last change point (i.e. the estimated $R_{t}=\min\{n\in\mathbb{N}:C_{t-n+1}=1$) and $\Omega_{t+1}=\textbf{P}(c_{t+1}=1|y_{1:t+1})$ the probability of a change given the observation. Note that this quantity, i.e. the posterior change point probability, is the same as our adaptation rate $\gamma_{t+1}$ of Eq. 11.

In the next subsection, we extend their approach to a more general case where the prior is a Gaussian distribution with arbitrary variance, i.e. $\mu_{t+1}\sim\mathcal{N}(\mu_{0},\sigma_{0}^{2})$. We then discuss the relation of this method to Particle Filtering. A performance comparison between our extended algorithms Nas10*∗* and Nas12*∗* and their original versions Nas10 and Nas12 is depicted in Supplementary Fig. A1 and Supplementary Fig. A2.

4.5.3 Nas10*∗* and Nas12*∗* algorithms

Let us consider that $y_{1:t}$ are observed, the time since the last change point $r_{t}$ is known, and the agent’s current estimation of $\mu_{t}$ is $\hat{\mu}_{t}$. It can be shown (see Appendix for the derivation) that the expected $\mu_{t+1}$ (i.e. $\hat{\mu}_{t+1}$) upon observing the new sample $y_{t+1}$ is

[TABLE]

where $\rho=\frac{\sigma^{2}}{\sigma^{2}_{0}}$, $\mu_{0}$ is the mean of the prior distribution and $\gamma_{t+1}$ is the adaptation rate of Eq. 11.

We can see that the updated mean is a weighted average, with surprise-modulated weights, between integrating the new observation with the current mean $\hat{\mu}_{t}$ and integrating it with the prior mean $\mu_{0}$, in the same spirit as the other algorithms we considered here. Eq. 78 can also be seen as a surprise-modulated weighted sum of two delta rules: one including a prediction error between the new observation and the current mean $(y_{t+1}-\hat{\mu}_{t})$ and one including a prediction error between the observed sample and the prior mean $(y_{t+1}-\mu_{0})$.

In order to obtain a form similar to the one of Nassar \BOthers. (\APACyear2012, \APACyear2010), we can rewrite the above formula as

[TABLE]

where we have defined $\alpha_{t+1}=\frac{\rho+\gamma_{t+1}r_{t}+1}{\rho+r_{t}+1}$. Hence the update rule takes the form of a weighted average, with fixed weights, between two delta rules: one including a prediction error between the prior mean and the current mean $(\mu_{0}-\hat{\mu}_{t})$ and one including a prediction error between the observed sample and the current mean $(y_{t+1}-\hat{\mu}_{t})$, both with surprise-modulated learning rates.

In Nassar \BOthers. (\APACyear2012, \APACyear2010) the true new mean after a change point is drawn from a uniform distribution with a range of values much larger than the width of the Gaussian likelihood. Their derivations implicitly approximate the uniform distribution with a Gaussian distribution with $\sigma_{0}\gg\sigma$. Note that if $\sigma_{0}\gg\sigma$ then $\rho\rightarrow 0$, so that the first term of Eq. 79 disappears, and $\alpha_{t+1}=\frac{1+\gamma_{t+1}r_{t}}{1+r_{t}}$. This results in the delta-rule of the original algorithm in Eq. 76 and Eq. 77, with $\gamma_{t+1}=\Omega_{t+1}$.

All of the calculations so far were done by assuming that $r_{t}$ is known. However, for the case of a non-stationary regime with a history of change points, the time interval $r_{t}$ is not known. Nassar \BOthers. (\APACyear2012, \APACyear2010) used the expected time interval $\hat{r}_{t}$ as an estimate. We make a distinction here between Nassar \BOthers. (\APACyear2012) and Nassar \BOthers. (\APACyear2010):

In Nassar \BOthers. (\APACyear2010) $\hat{r}_{t}$ is calculated recursively on each trial in the same spirit as Eq. 12: $\hat{r}_{t+1}=(1-\gamma_{t+1})(\hat{r}_{t}+1)+\gamma_{t+1}$, i.e., at each time step, there is a probability $(1-\gamma_{t+1})$ that $\hat{r}_{t}$ increments by 1 and a probability $\gamma_{t+1}$ that it is reset to 1. So $\hat{r}_{t+1}$ is the weighted sum of these two outcomes. Hence, Eq. 79 combined with the expected time interval $\hat{r}_{t}$ constitutes a generalization of the update rule of Nassar \BOthers. (\APACyear2010) for the case of Gaussian prior $\mathcal{N}(\mu_{0},\sigma_{0}^{2})$. We call this algorithm Nas10*∗* (see Appendix for pseudocode).

In Nassar \BOthers. (\APACyear2012), the variance $\hat{\sigma}^{2}_{t+1}=\text{Var}[\mu_{t+1}|y_{1:t+1}]$ is estimated given $\hat{\mu}_{t}$, $\hat{r}_{t}$, and $\hat{\sigma}^{2}_{t}$. Based on this variance, $\hat{r}_{t+1}=\frac{\sigma^{2}}{\hat{\sigma}^{2}_{t+1}}-\frac{\sigma^{2}}{\sigma^{2}_{0}}$ is computed. The derivation of the recursive computation of $\hat{\sigma}^{2}_{t+1}$ for the case of Gaussian priors can be found in the Appendix. We call the combination of Eq. 79 with this way of computing the expected time interval $\hat{r}_{t}$ Nas12*∗* (see Appendix for pseudocode). These two versions of calculating $\hat{r}_{t}$ in Nassar \BOthers. (\APACyear2010) and Nassar \BOthers. (\APACyear2012) give different results, and we compare our algorithms with both Nas10*∗* and Nas12*∗* in our simulations. Note that, as discussed in the section “Online Bayesian inference modulated by surprise” of the Results, the posterior belief at time $t+1$ does not generally belong to the same family of distributions as the belief of time $t$. However, we therefore approximate for both algorithms the posterior belief $P(\theta|y_{1:t+1})$ by a Gaussian.

4.5.4 Nassar’s algorithm and Particle Filtering with one particle

In the case of Particle Filtering (cf. Eq. 25) with only one particle, at each time step we sample the particle’s hidden state with change probability $Q(c_{t+1}^{(1)}=1|c_{1:t}^{(1)},y_{1:t+1})=\gamma_{t+1}$, generating a posterior belief that takes two possible values with probability (according to the proposal distribution)

[TABLE]

So, in expectation, the updated belief will be

[TABLE]

If we apply Eq. 81 to $\hat{\mu}_{t+1}$, we find that $\mathbb{E}_{Q}[\hat{\mu}^{(t+1)}]$, is identical to the generalization of Nassar \BOthers. (\APACyear2010) (see Eq. 78).

Moreover, in Particle Filtering with a single particle, we sample the particle’s hidden state, which is equivalent to sampling the interval $\hat{R}_{t+1}$. Because $\hat{R}_{t+1}$ takes the value $\hat{r}_{t}+1$ with $(1-\gamma_{t+1})$ and the value of $1$ (=reset) with probability $\gamma_{t+1}$, the expected value of $\hat{R}_{t+1}$ is

[TABLE]

In other words, in Nassar \BOthers. (\APACyear2010), the belief is updated based on the expected $\hat{r}_{t}$, whereas in Particle Filtering with one particle, the belief is updated using the sampled $\hat{r}_{t}$.

In summary, the two methods will give different estimates on a trial-per-trial basis, but the same result in expectation. The pseudocode for Particle Filtering with one particle for the particular case of the Gaussian estimation task can be found in the Appendix.

4.6 Application to the exponential family

For our all three algorithms Variational SMiLe, Message Passing with fixed number $N$ of particles, and Particle Filtering, we derive compact update rules for $\hat{\pi}^{(t+1)}(\theta)$ when the likelihood function $P_{Y}(y|\theta)$ is in the exponential family and $\pi^{(0)}(\theta)$ is its conjugate prior. In that case, the likelihood function has the form

[TABLE]

where $\theta$ is the vector of natural parameters, $h(y)$ is a positive function, $\phi(y)$ is the vector of sufficient statistics, and $A(\theta)$ is the normalization factor. Then, the conjugate prior $\pi^{(0)}$ has the form

[TABLE]

where $\chi^{(0)}$ and $\nu^{(0)}$ are the distribution parameters, $\tilde{h}(\theta)$ is a positive function, and f\big{(}\chi^{(0)},\nu^{(0)}\big{)} is the normalization factor. For this setting and while \pi^{(t)}=\textbf{P}_{\pi}\big{(}\Theta=\theta;\chi^{(t)},\nu^{(t)}\big{)}, the “Bayes Factor Surprise” has the compact form

[TABLE]

The pseudocode for Variational SMiLe, MP$N$, and Particle Filtering can be seen in Algorithms 1, 2, and 3, respectively.

4.7 Simulation task

In this subsection, we first argue why the mean squared error is a proper measure for comparing different algorithms with each other, and then we explain the version of Leaky integrator which we used for simulations.

4.7.1 Mean squared error as an optimality measure

Consider the case that at each time point $t$, the goal of an agent is to have an estimation of the parameter $\Theta_{t}$ as a function of the observations $Y_{1:t}$, i.e. $\hat{\Theta}_{t}=f(Y_{1:t})$. The estimator which minimizes the mean squared error \textbf{MSE}[\hat{\Theta}_{t}]=\mathbb{E}_{\textbf{P}(Y_{1:t},\Theta_{t})}\big{[}(\hat{\Theta}_{t}-\Theta_{t})^{2}\big{]} is

[TABLE]

which is the expected value of $\Theta_{t}$ conditioned on the observations $Y_{1:t}$, or in other words under the Bayes-optimal current belief (see Papoulis \BBA Saunders (\APACyear1989) for a proof). The MSE for any other estimator $\hat{\Theta}_{t}$ can be written as (see below for the proof)

[TABLE]

This means that the MSE for any arbitrary estimator $\hat{\Theta}_{t}$ includes two terms: the optimal MSE and the mismatch of the actual estimator from the optimal estimator $\hat{\Theta}_{t}^{\text{Opt}}$. As a result, if the estimator we are interested in is the expected value of $\Theta_{t}$ under the approximate belief $\hat{\pi}^{(t)}$ computed by each of our algorithms (i.e. \hat{\Theta}_{t}^{\prime}=\mathbb{E}_{\hat{\pi}^{(t)}}\big{[}\Theta_{t}\big{]}), the second term in LABEL:eq:MSE_dMSE, i.e. the deviation from optimality, is a measure of how good the approximation is.

Proof for the algorithms without sampling:

Consider $\hat{\Theta}^{\text{Opt}}_{t}=f_{\text{Opt}}(Y_{1:t})$. Then, for any other arbitrary estimator $\hat{\Theta}_{t}=f(Y_{1:t})$ (except for the ones with sampling), we have

[TABLE]

The quadratic term in the last line can be expanded and written as

[TABLE]

The random variables in the expected value of the first line are not dependent on $\Theta_{t}$, so it can be computed over $Y_{1:t}$. The expected value of the second line is equal to $\textbf{MSE}[\hat{\Theta}^{\text{Opt}}_{t}]$. It can also be shown that the expected value of the third line is equal to 0, i.e.

[TABLE]

where in the last line we used the definition of the optimal estimator. All together, we have

[TABLE]

Proof for the algorithms with sampling:

For particle filtering (and any kind of estimator with sampling), the estimator is not a deterministic function of observations $Y_{1:t}$. Rather, the estimator is a function of observations as well as a set of random variables (samples) which are drawn from a distribution which is also a function of observations $Y_{1:t}$. In our case, the samples are the sequence of hidden states $C_{1:t}$. The estimator can be written as

[TABLE]

where $C_{1:t}^{(1:N)}$ are $N$ iid samples drawn from the proposal distribution $Q(C_{1:t}|Y_{1:t})$. MSE for this estimator should also be averaged over the samples, which leads to

[TABLE]

Similar to what we did before, the MSE for particle filtering can be written as

[TABLE]

which can be written in terms of bias and variance over samples as

[TABLE]

4.7.2 Leaky integration

Gaussian task: The goal is to have an estimation of the mean of the Gaussian distribution at each time $t$, denoted by $\hat{\theta}_{t}$. Given a leak parameter $\omega\in(0,1]$, the leaky integrator estimation is

[TABLE]

Categorical task: The goal is to have an estimation of the parameters of the categorical distribution at each time $t$, denoted by $\hat{\theta}_{t}=[\hat{\theta}_{i,t}]_{i=1}^{N}$ for the case that there are $N$ categories. Given a leak parameter $\omega\in(0,1]$, the leaky integrator estimation is

[TABLE]

where $\delta$ is the Kronecker delta function.

4.8 Derivation of the Formula Relating Shannon Surprise to the Modulated Learning Rate

Given the defined generative model, the Shannon surprise upon observing $y_{t+1}$ can be written as

[TABLE]

where \gamma_{t+1}=\gamma\Big{(}\textbf{S}_{\mathrm{BF}}(y_{t+1};\pi^{(t)}),m=\frac{p_{c}}{1-p_{c}}\Big{)} of Eq. 11. As a result, the modulated adaptation rate can be written as in LABEL:gamma_Shannon and the Bayes Factor Surprise as in Eq. 29.

4.9 Experimental predictions

4.9.1 Setting

Consider a Gaussian task where $Y_{t}$ can take values in $\mathbb{R}$. The likelihood function $P_{Y}(y|\theta)$ is defined as

[TABLE]

where $\sigma\in\mathbb{R}^{+}$ is the standard deviation, and $\theta\in\mathbb{R}$ is the mean of the distribution, i.e. the parameter of the likelihood. Whenever there is a change in the environment (with probability $p_{c}\in(0,1)$), the value $\theta$ is drawn from the prior distribution $\pi^{(0)}(\theta)=\mathcal{N}(\theta;0,1)$.

4.9.2 Theoretical proofs for prediction 1

For our theoretical derivations for our first prediction, we consider the specific but relatively mild assumption that the subjects’ belief $\pi^{(t)}$ at each time is a Gaussian distribution

[TABLE]

where $\hat{\theta}_{t}$ and $\hat{\sigma}_{t}$ are determined by the learning algorithm and the sequence of observations $y_{1:t}$. This is the case when the subjects use either VarSMiLe, Nas10*∗, Nas12∗*, pf1, MP1, or Leaky Integration as their learning rule. With such assumptions, the inferred probability distribution $P(y;\pi^{(t)})$ can be written as

[TABLE]

As mentioned in the Results section, we define, at time $t$, the prediction error as $\delta_{t+1}=y_{t+1}-\hat{\theta}_{t}$ and the “sign bias” as $s_{t+1}=\text{sign}(\delta_{t+1}\hat{\theta}_{t})$. Then, given an absolute prediction $\hat{\theta}>0$, an absolute prediction error $\delta>0$, a standard deviation $\sigma_{C}$, and a sign bias $s\in\{-1,1\}$, the average Bayes Factor Surprise is computed as

[TABLE]

It can easily be shown that the value $\textbf{S}_{\mathrm{BF}}(y_{t};\hat{\pi}^{(t-1)})$ is same for all $t\in\mathcal{T}$, and hence the average surprise is same as the surprise for each time point. For example, the average surprise for $s=+1$ is equal to

[TABLE]

Similar formulas can be computed for $\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s=-1,\sigma_{C})$. Then, the difference $\Delta\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,\sigma_{C})=\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s=+1,\sigma_{C})-\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,s=-1,\sigma_{C})$ can be computed as

[TABLE]

It can be shown that

[TABLE]

for all $\hat{\theta}>0$, $\delta>0$, and $\sigma_{C}>0$. The first inequality is trivial, and the proof for the second inequality is given below.

The average Shannon Surprise can be computed in a similar way. For example, for $s=+1$, we have

[TABLE]

and then the difference $\Delta\bar{\textbf{S}}_{Sh}(\hat{\theta},\delta,\sigma_{C})=\bar{\textbf{S}}_{Sh}(\hat{\theta},\delta,s=+1,\sigma_{C})-\bar{\textbf{S}}_{Sh}(\hat{\theta},\delta,s=-1,\sigma_{C})$ can be computed as

[TABLE]

where $m=\frac{p_{c}}{1-p_{c}}$. Then, using the results for the Bayes Factor Surprise, we have

[TABLE]

for all $\hat{\theta}>0$, $\delta>0$, and $\sigma_{C}>0$. See below for the proof of the second inequality.

Proof of the 2nd inequality for $\textbf{S}_{\mathrm{BF}}$:

Let us define the variables

[TABLE]

as well as the functions

[TABLE]

The following inequalities hold true

[TABLE]

Then, the derivative of $f(\delta)$ can be compute as

[TABLE]

Therefore we have $\frac{\partial}{\partial\delta}\Delta\bar{\textbf{S}}_{\mathrm{BF}}(\hat{\theta},\delta,\sigma_{C})=\frac{d}{d\delta}f(\delta)<0$.

Proof of the 2nd inequality for $\textbf{S}_{\mathrm{Sh}}$:

Using the functions we defined for the previous proof, and after computing the partial derivative of Eq. 107, we have

[TABLE]

The derivative of the last term can be written in terms of the derivates of $f_{1}$ and $f_{2}$, indicated by $f_{1}^{\prime}$ and $f_{2}^{\prime}$, respectively,

[TABLE]

The 1st term is always positive based on the proof for $\textbf{S}_{\mathrm{BF}}$. The 2nd term is also always positive, because

[TABLE]

As a result, we have $\frac{\partial}{\partial\delta}\Delta\bar{\textbf{S}}_{\mathrm{Sh}}(\hat{\theta},\delta,\sigma_{C})>0$.

4.9.3 Simulation procedure for prediction 1

In order to relax the main assumption of our theoretical proofs (i.e. the belief is always a Gaussian distribution), to include the practical difficulties of a real experiment (e.g. to use $|\delta_{t}|\approx\delta$ instead of $|\delta_{t}|=\delta$), and to have an estimation of the effect size, we also performed simulations for our first experimental prediction.

For each simulated subject, the procedure of our simulation was as follows:

We fixed the hyper parameters $\sigma^{2}$ and $p_{c}$ for producing samples. 2. 2.

We selected a learning algorithm (e.g. pf20) and fixed its corresponding tuned parameters (based on our simulations in the Results section). 3. 3.

We applied the learning algorithm over a sequence of observations $y_{1:T}$. Note that in a real experiment, this step can be done through a few episodes, which makes it possible to have a long sequence of observations, i.e. large $T$. 4. 4.

At each time $t$, we saved the values $y_{t}$, $\hat{\theta}_{t}$, $\hat{\sigma}_{t}$, $\delta_{t}$, $s_{t}$, $\textbf{S}_{\mathrm{Sh}}(y_{t};\hat{\pi}^{(t-1)})$, and $\textbf{S}_{\mathrm{BF}}(y_{t};\hat{\pi}^{(t-1)})$.

Then, given an absolute prediction $\hat{\theta}>0$, an absolute prediction error $\delta>0$, a standard deviation $\sigma_{C}>0$, and a sign bias $s\in\{-1,1\}$, we defined the set of time points

[TABLE]

where $||\hat{\theta}_{t-1}|-\hat{\theta}|<\Delta\theta$, $||\delta_{t}|-\delta|<\Delta\delta$, and $|\hat{\sigma}_{t}-\sigma_{C}|<\Delta\sigma_{C}$ are equivalent to $|\hat{\theta}_{t-1}|\approx\hat{\theta}$, $|\delta_{t}|\approx\delta$, and $\hat{\sigma}_{t}\approx\sigma_{C}$, respectively. $\Delta\theta$, $\Delta\delta$, and $\Delta\sigma_{C}$ are positive real values that should be determined based on practical limitations (mainly the length of the observation sequence $T$). We then computed the average surprise values as

[TABLE]

We repeated this procedure for $N$ different simulated subjects (with different random seeds). The average of $\bar{\textbf{S}}_{\mathrm{BF}}$ and $\bar{\textbf{S}}_{Sh}$ over $N=20$ subjects, for two learning algorithms (i.e. Nas12*∗* and pf20), and for $T=500$, $\hat{\theta}=1$, $\sigma_{C}=0.5$, $\Delta\theta=0.25$, $\Delta\delta=0.1$, and $\Delta\sigma_{C}=1$ is shown in Fig. 10B. The results are the same as what was predicted by our theoretical analysis.

4.9.4 Simulation procedure for prediction 2

For our second prediction, the theoretical proof is trivial. However, in order to have a setting similar to a real experiment (e.g. to use $P(y_{t+1};\hat{\pi}^{(t)})\approx p$ instead of $P(y_{t+1};\hat{\pi}^{(t)})=p$), and to have an estimation of the effect size, we used simulations also for our second experimental predictions.

We followed the same procedure as the one for the simulation of the first prediction. For each simulated subject, and at each time $t$, we saved the quantities $P(y_{t+1};\hat{\pi}^{(t)})$, $P(y_{t+1};\hat{\pi}^{(0)})$, $\textbf{S}_{\mathrm{Sh}}(y_{t};\hat{\pi}^{(t-1)})$, and $\textbf{S}_{\mathrm{BF}}(y_{t};\hat{\pi}^{(t-1)})$. Then, for a given a probability value $p>0$, we defined the set of time points

[TABLE]

where $|P(y_{t+1};\hat{\pi}^{(t)})-p|<\Delta p$ and $|P(y_{t+1};\hat{\pi}^{(0)})-p|<\Delta p$ are equivalent to $P(y_{t+1};\hat{\pi}^{(t)})\approx p$ and $P(y_{t+1};\hat{\pi}^{(0)})\approx p$, respectively. $\Delta p$ is a positive real value that should be determined based on practical limitations (mainly the length of the observation sequence $T$). We then computed the average surprise $\bar{\textbf{S}}_{\mathrm{BF}}(p)$ and $\bar{\textbf{S}}_{\mathrm{Sh}}(p)$ over $\mathcal{T}$ for each value of $p$. We repeated this procedure for $N$ different simulated subjects (with different random seeds). The average of $\bar{\textbf{S}}_{\mathrm{BF}}$ and $\bar{\textbf{S}}_{Sh}$ over $N=20$ subjects, for two learning algorithms (i.e. Nas12*∗* and pf20), and for $T=500$ and $\Delta p=0.0125$ is shown in Fig. 11B.

Acknowledgments

We thank both reviewers for their constructive and helpful comments. This research was supported by Swiss National Science Foundation No. 200020_184615) and by the European Union Horizon 2020 Framework Program under grant agreement No.785907 (Human Brain Project, SGA2).

Appendix

Modified algorithm of Nassar \BOthers. (\APACyear2012, \APACyear2010): Adaptation for Gaussian prior

Recursive update of the estimated mean for Gaussian prior

Let us first consider the case of a stationary regime (i.e. no change points) where observed samples are drawn from a Gaussian distribution with known variance, i.e. $y_{t+1}|\theta\sim\mathcal{N}(\theta,\sigma^{2})$, and the parameter $\theta$ is also drawn from a Gaussian distribution $\theta\sim\mathcal{N}(\mu_{0},\sigma_{0}^{2})$. After having observed samples $y_{1},...,y_{t+1}$, it can be shown that, using Bayes’ rule, the posterior distribution $P(\theta|y_{1:t+1})=\pi^{(t+1)}_{B}(\theta)$ is

[TABLE]

An estimate of $\theta$ is its expected value $\mathbb{E}(\theta|y_{1:t+1})=\mu_{B,t+1}$.

In a non-stationary regime where, after having observed $y_{1},...,y_{t}$ from the same hidden state, there is the possibility for a change point upon observing $y_{t+1}$, the posterior distribution is

[TABLE]

To facilitate notation in this subsection we denote $c_{t+1}=0$ as “stay” and $c_{t+1}=1$ as “change” so that

[TABLE]

Note that the above is equivalent to Bayesian recursive formula (Eq. 12) of the main text, where $\gamma_{t+1}$ is the adaptation rate we saw in Eq. 11 of the main text, and is essentially the probability to change given the new observation, i.e. $\textbf{P}(c_{t+1}=1|y_{1:t+1})$. In Nassar \BOthers. (\APACyear2010) this quantity is denoted as $\Omega_{t+1}$. Taking Eq. A1 into account we have

[TABLE]

where $r_{t}$ is the time interval of observations coming from the same hidden state, calculated at time $t$. Taking the expectation of Eq. A3 the estimated mean upon observing the new sample $y_{t+1}$ is

[TABLE]

where we dropped the subscript $t+1$ in $\gamma$ to simplify notations. We have

[TABLE]

Because \hat{\mu}_{t}=\frac{1}{\frac{1}{\sigma^{2}_{0}}+\frac{r_{t}}{\sigma^{2}}}\Big{(}\frac{\mu_{0}}{\sigma^{2}_{0}}+\frac{\sum_{i=t+1-r_{t}}^{t}y_{i}}{\sigma^{2}}\Big{)}, after a few lines of algebra we have

[TABLE]

We now define $\rho=\frac{\sigma^{2}}{\sigma^{2}_{0}}$ and find

[TABLE]

A rearrangement of the terms and inclusion of the dependency of $\gamma$ on time yields

[TABLE]

In order to obtain a form similar to the one of Nassar \BOthers. (\APACyear2012, \APACyear2010) we continue and we spell out the terms that include the quantities $\hat{\mu}_{t},\mu_{0}$ and $y_{t+1}$

[TABLE]

Using that $\frac{1}{\rho+r_{t}+1}=\frac{1}{\rho+1}-\frac{r_{t}}{(\rho+1)(\rho+r_{t}+1)}$ we have

[TABLE]

After a further step of algebra we arrive at

[TABLE]

If we define $1-\alpha=(1-\gamma)\frac{r_{t}}{\rho+r_{t}+1}\Rightarrow\alpha=1-(1-\gamma)\frac{r_{t}}{\rho+r_{t}+1}\Rightarrow\alpha=\frac{\rho+\gamma r_{t}+1}{\rho+r_{t}+1}$ and rearrange the terms, we have

[TABLE]

Adding back the dependency of $\gamma$ and $\alpha$ on time we finally have

[TABLE]

Recursive update of the the Estimated Variance for Gaussian Prior

In Nassar \BOthers. (\APACyear2012) the authors calculate first the variance $\hat{\sigma}^{2}_{t+1}=\text{Var}(\theta|y_{1:t+1})$ and based on this compute then $\hat{r}_{t+1}$. We derive here these calculations for the case of Gaussian prior. We remind once again that

[TABLE]

Then for the variance $\hat{\sigma}^{2}_{t+1}=\text{Var}(\theta|y_{1:t+1})$ we have

[TABLE]

where $\sigma^{2}_{B,t+1}=\frac{1}{\frac{1}{\sigma^{2}_{0}}+\frac{r_{t}+1}{\sigma^{2}}}$ and $\sigma^{2}_{change}=\frac{1}{\frac{1}{\sigma^{2}_{0}}+\frac{1}{\sigma^{2}}}$.

We have defined earlier $\rho=\frac{\sigma^{2}}{\sigma^{2}_{0}}$ so that

[TABLE]

Using, as before, that $\frac{1}{\rho+r_{t}+1}=\frac{1}{\rho+1}-\frac{r_{t}}{(\rho+1)(\rho+r_{t}+1)}$ we have

[TABLE]

We have defined earlier the learning rate $\alpha=1-(1-\gamma)\frac{r_{t}}{\rho+r_{t}+1}$, so we can write

[TABLE]

Note that \mu_{t}=\frac{1}{\frac{1}{\sigma^{2}_{0}}+\frac{r_{t}}{\sigma^{2}}}\Big{(}\frac{\mu_{0}}{\sigma^{2}_{0}}+\frac{\sum_{i=t+1-r_{t}}^{t}y_{i}}{\sigma^{2}}\Big{)} so for the calculation of the last term we have

[TABLE]

We now rearrange terms

[TABLE]

and finally we have

[TABLE]

Implementation of Nas10*∗, Nas12∗* and Particle Filtering with 1 particle for the Gaussian estimation task

We provide here the pseudocode for the algorithms Nas10*∗, Nas12∗* and Particle Filtering with 1 particle for the Gaussian estimation task. Observations are drawn from a Gaussian distribution with known variance and unknown mean, i.e. $y_{t+1}|\mu_{t+1}\sim\mathcal{N}(\mu_{t+1},\sigma^{2})$ and $\theta_{t}=\mu_{t}$. When there is a change, the parameter $\mu$ is also drawn from a Gaussian distribution $\mu\sim\mathcal{N}(\mu_{0},\sigma_{0}^{2})$. All three algorithms estimate the expected $\mu_{t+1}$ (i.e. $\hat{\mu}_{t+1}$) upon observing a new sample $y_{t+1}$.

After re-writing Eq. 85 we have

[TABLE]

Note that the pseudocode for pf1 provided here, is a translation of Algorithm 3 to the case of a single particle (where there are no weights to calculate) and for the Gaussian distribution as a particular instance of the exponential family.