---
title: Recursively Computable Approximate Sufficient Statistics
url: https://www.emergentmind.com/topics/recursively-computable-approximate-sufficient-statistics-rcass
type: topic
---

# Recursively Computable Approximate Sufficient Statistics

Recursively Computable Approximate Sufficient Statistics (RCASS) are data summaries that render an approximate likelihood, posterior, regret bound, or state-evolution description a fixed function of a finite-dimensional state, while that state can be updated recursively as new observations, segments, or iterations arrive. Across the cited literature, RCASS appear as information-theoretically selected summaries in approximate Bayesian computation (ABC), polynomial approximate sufficient statistics for generalized linear models, additive sufficient statistics for Burkholder-based online learning, cross-correlation integrands for stochastic gravitational-wave background searches, code-length-optimal approximate sufficient statistics for parametric families, and scalar state variables for sufficient statistic memory approximate message passing [1106.6281], [1709.09216], [1803.07617], [2012.00907], [1612.02542], [2206.11674].

## 1. Definition and formal structure

At the most classical level, sufficiency means that a statistic $S$ retains all information in the data $X$ about a parameter $\theta$. The Fisher–Neyman factorization expresses this as
$$
p(x\mid \theta)=h(x)\,g(S(x),\theta),
$$
and an equivalent posterior formulation is
$$
p(\theta\mid x)=p(\theta\mid S(x)).
$$
In information-theoretic form, sufficiency is equivalent to
$$
I(\theta;S)=I(\theta;X),
$$
or, equivalently, $I(\theta;X\mid S)=0$ [1106.6281].

Approximate sufficiency relaxes this equality by minimizing the information lost when $S(X)$ replaces $X$. In the ABC setting, the loss for parameter inference is
$$
L_{\mathrm{param}}(S)=\mathbb{E}_{X}\!\left[D_{\mathrm{KL}}\!\big(p(\theta\mid X)\,\|\,p(\theta\mid S(X))\big)\right],
$$
and for model selection
$$
L_{\mathrm{model}}(S)=\mathbb{E}_{X}\!\left[D_{\mathrm{KL}}\!\big(p(M\mid X)\,\|\,p(M\mid S(X))\big)\right].
$$
A statistic is approximately sufficient when the relevant loss is small [1106.6281].

The “recursively computable” part refers to update rules that avoid storing the full history. In additive settings this takes the form
$$
S_t=S_{t-1}+T(x_t,\hat y_t,\delta_t), \qquad S_0=0,
$$
while in more general recursive settings one writes
$$
S_t=T(x_t,\hat y_t,\delta_t,S_{t-1}).
$$
The central operational idea is that both inference and computation are carried by the evolving state $S_t$, not by the raw data sequence [1803.07617].

This suggests that RCASS is best understood as a cross-domain design principle: one searches for a low-dimensional state that is approximately sufficient for the target task and closed under streaming, distributed, or iterative updates.

| Setting | RCASS state | Recursive mechanism |
|---|---|---|
| ABC | selected summary subset | add statistics that most reduce KL loss |
| PASS-GLM | $t_\kappa=\sum_i (y_i x_i)^\kappa$ or $T_k=\sum_i x_i^{\otimes k}$ | additive updates |
| Online learning | $\zeta_t=\sum_{j\le t} T(x_j,\hat y_j,\delta_j)$ | additive updates |
| SGWB searches | frequency integrand and inverse variance | segmentwise weighted sums |
| Parametric coding | quantized estimator or expectation-parameter index | update estimate, then re-quantize |
| SS-MAMP | $S_t=(v_t^\gamma,v_t^\phi)$ | state-evolution recursion |

## 2. Information-theoretic RCASS in approximate Bayesian computation

In ABC, likelihoods are intractable but simulation from the model is available. The posterior based on a summary statistic $S(x)$ is
$$
p_\varepsilon(\theta\mid S(x)) \propto \pi(\theta)\int K_\varepsilon\!\big(\rho(S(x),S(y))\big)\,p(y\mid\theta)\,dy,
$$
with the paper focusing on the special case $K_\varepsilon(u)=\mathbf{1}(u\le \varepsilon)$. The methodological problem is that comparing summaries rather than full data generally loses information unless the summaries are sufficient [1106.6281].

The RCASS construction in this setting is explicitly recursive. Starting from a candidate pool $S=\{s_1,\dots,s_K\}$, one either performs an exhaustive search over all subsets or, more practically, a greedy forward construction. The greedy parameter-inference procedure initializes $S_0=\varnothing$, chooses the first statistic to maximize mutual information with $\theta$, and thereafter adds the statistic that maximizes the KL change
$$
KL\!\big(\hat p(\theta\mid S_{t-1}\cup\{s_j\})\,\|\,\hat p(\theta\mid S_{t-1})\big).
$$
The recursion stops when the marginal gain falls below a tolerance $\delta$, and an optional backward elimination stage re-tests previously selected summaries for redundancy. The paper also notes an optional penalized objective,
$$
\min_S\, L_{\mathrm{param}}(S)+\lambda\,\mathrm{dim}(S),
$$
to favor parsimony [1106.6281].

For model selection, the construction is more restrictive because parameter sufficiency within each model does not by itself imply sufficiency for discriminating models. The joint-space criterion is
$$
I(M,\theta_1,\dots,\theta_q;X\mid S)=0,
$$
with decomposition
$$
I(M,\theta_1,\dots,\theta_q;X\mid S)
=
I(M;X\mid \theta_1,\dots,\theta_q,S)+\sum_{i=1}^q I(\theta_i;X\mid S).
$$
The practical consequence is a two-stage RCASS workflow: first construct minimal sufficient sets for parameters within each model, then add further statistics that materially change the model posterior or Bayes factors conditional on those parameter summaries [1106.6281].

The empirical demonstrations show why this distinction matters. In the two-normal example, the mean $\bar y$ is selected universally for parameter inference, while model selection requires adding $S^2$, so that $\{\bar y,S^2\}$ is sufficient for the joint space. In the coalescent-model experiments, selected summaries vary across datasets and models, with $S4$ frequently selected and the random statistic $S11$ rarely chosen. In the random-walk models, $S5$ is often selected across models, while $S2$ is preferred particularly for biased walks. These cases exemplify data-dependent approximate sufficiency rather than a fixed universal summary set [1106.6281].

The framework also makes the main limitations explicit. Greedy selection is order-dependent, high-dimensional summary spaces induce the ABC curse of dimensionality, multivariate KL estimation can be noisy, and summaries sufficient for parameter inference may remain insufficient for model choice. The prescribed mitigations are backward elimination, stochastic or heuristic search, scale-aware discrepancies, sensitivity analysis over $\varepsilon$, and explicit enforcement of joint-space sufficiency [1106.6281].

## 3. Polynomial and distributed RCASS for Bayesian generalized linear models

PASS-GLM instantiates RCASS by replacing the non-linear term in a generalized linear model log-likelihood with a polynomial approximation. For canonical GLMs,
$$
L(\beta)=\sum_{i=1}^n \big[y_i x_i^\top \beta-b(x_i^\top \beta)+c(y_i)\big],
$$
and PASS-GLM replaces $b(x_i^\top \beta)$ by
$$
q_m(x_i^\top \beta)=\sum_{k=0}^m a_k (x_i^\top \beta)^k.
$$
Because
$$
(x_i^\top \beta)^k=\langle x_i^{\otimes k},\beta^{\otimes k}\rangle,
$$
the approximate log-likelihood depends on the data only through sums of tensor powers of the covariates, yielding approximate sufficient statistics that are additive across data points, updated recursively in streaming, aggregated exactly by summation in distributed settings, and introduce no additional approximation at aggregation time [1709.09216].

In the logistic parameterization with $s_i=y_i x_i^\top \beta$, the paper uses a polynomial approximation
$$
f_M(s)=\sum_{m=0}^{M} c_m \psi_m(s)
$$
in an orthogonal basis, specifically Chebyshev polynomials on $[-R,R]$. The coefficients satisfy
$$
c_m=\int \psi_m(s)\,\mathrm{logit}(s)\,d\pi(s),
$$
and the monomial coefficients are obtained through the basis expansion $\psi_m(s)=\sum_{j=0}^m \alpha_{m,j}s^j$. The approximate likelihood then takes the exponential-family-like form
$$
\tilde L(\beta)=\sum_{0\le |\kappa|\le M} a(\kappa,m,M)\, t_\kappa\, \beta^\kappa,
$$
where
$$
t_\kappa=\sum_{i=1}^n (y_i x_i)^\kappa.
$$
These $t_\kappa$ are the RCASS in the precise sense that the whole approximation is a fixed function of $\beta$ and the accumulated statistics [1709.09216].

The recursive update rule is purely additive. When a new datum $(x_{n+1},y_{n+1})$ arrives,
$$
t_\kappa^{(n+1)}=t_\kappa^{(n)}+(y_{n+1}x_{n+1})^\kappa.
$$
For quadratic logistic regression, the RCASS reduce to
$$
T_0=n,\qquad T_1=\sum_{i=1}^n y_i x_i,\qquad T_2=\sum_{i=1}^n x_i x_i^\top,
$$
with updates
$$
T_0\leftarrow T_0+1,\quad T_1\leftarrow T_1+yx,\quad T_2\leftarrow T_2+xx^\top.
$$
With a Gaussian prior and $M=2$, the approximate posterior is Gaussian with closed-form mean and covariance [1709.09216].

Theoretical guarantees are attached to the polynomial approximation itself. Chebyshev approximation error on $[-R,R]$ is exponentially small in the degree $M$, and the paper gives bounds for MAP error, approximate posterior quality in Wasserstein-$1$ distance, and posterior mean and uncertainty estimates. For logistic regression, the reported approximation obeys
$$
\sup_{t\in[-4,4]} |f_2(t)-\mathrm{logit}(t)|<0.069,
$$
with the empirical observation that for many datasets at least $98\%$ of the inner products $s_i=y_i x_i^\top \beta$ fall in $[-4,4]$. The paper also states a pathology for degrees $M=4k$, $k\ge 1$, and recommends $M=2+4k$ [1709.09216].

The empirical results show how these structural claims translate into computation. PASS-LR2 on ChemReact, Webspam, CovType, and CodRNA was reported as approximately $10\times$ faster than SGD and $100$–$1000\times$ faster than Laplace, while remaining competitive in posterior mean, posterior variance, and test log-likelihood; CodRNA is identified as a known outlier because many $s_i$ fall outside $[-4,4]$. On Criteo, with $40$ million points and $20{,}000$ dimensions after random projection, streaming PASS-LR2 achieved negative test log-likelihood $0.045$ versus $0.07$ for SGD, slightly worse AUC, and distributed near-linear speedups with $K$ cores giving approximately $K/2$ speedup from a baseline of $3$ minutes at $1$ core on $6$ million points [1709.09216].

## 4. Additive RCASS in online learning and the Burkholder method

In online learning, RCASS arise when regret can be upper bounded by a function of cumulative sufficient statistics rather than the full data sequence. The paper defines a sufficient statistic pair $(T,V)$ by the inequality
$$
\sum_{t=1}^n \ell(\hat y_t,y_t)-\phi(x_1,y_1,\ldots,x_n,y_n)
\le
V\!\left(\sum_{t=1}^n T(x_t,\hat y_t,\partial \ell(\hat y_t,y_t))\right).
$$
This representation is “approximate” when the desired regret expression is upper bounded by $V$ rather than exactly equal to it. The recursively computable state is
$$
\zeta_t=\sum_{j=1}^t T(x_j,\hat y_j,\delta_j)
=
\zeta_{t-1}+T(x_t,\hat y_t,\delta_t),
\qquad \delta_t=\partial \ell(\hat y_t,y_t),
$$
so the learner need only keep $\zeta_t$ in memory [1803.07617].

The Burkholder method gives the structural condition under which such RCASS support an actual online algorithm. The key object is a function $U:\mathcal T\to\mathbb R$ satisfying $U(0)\le 0$, $U(\tau)\ge V(\tau)$, and the restricted concavity condition
$$
\mathbb E_{\alpha\sim p}\big[U(\tau+T(z,\alpha))\big]\le U(\tau)
$$
for mean-zero $\alpha\in[-L,L]$. In the convex case it suffices to verify a Rademacher version. Existence of such a Burkholder function is equivalent to the needed martingale inequality, and the resulting prediction strategy evaluates $U$ only on the cumulative sufficient statistic [1803.07617].

Two explicit instantiations make the RCASS structure concrete. For parameter-free supervised learning with linear classes and a $\beta$-smooth norm, the sufficient statistics are
$$
T(x_t,\hat y_t,\delta_t)=(\delta_t \hat y_t,\delta_t x_t),
$$
so the running state is $(b_t,x_t^{\mathrm{sum}})=(\sum_{s\le t}\delta_s \hat y_s,\sum_{s\le t}\delta_s x_s)$, maintained in $O(d)$ time and space per step. The corresponding time-varying Burkholder family is
$$
U_t(b,x)
=
b+\frac{1}{\sqrt n}\exp\!\left(\frac{\|x\|^2}{2\beta t}+\frac{1}{2}\sum_{s=t+1}^n \frac{1}{s}\right)-1,
$$
which yields the comparator-dependent regret bound stated in the paper [1803.07617].

For matrix prediction, the sufficient statistics are
$$
T(X_t,\hat y_t,\delta_t)=\big(\delta_t \hat y_t,\ \delta_t \mathcal H(X_t),\ \mathcal M(X_t)\big),
$$
with
$$
\mathcal H(X)=
\begin{pmatrix}
0 & X\\
X^\top & 0
\end{pmatrix},
\qquad
\mathcal M(X)=\mathcal H(X)^2.
$$
The corresponding Burkholder function,
$$
U(a,H,M)=a+\frac{r}{\eta}\log\mathrm{tr}\exp\!\left(\eta H-\frac{1}{2}\eta^2L^2M\right)-\frac{c}{\eta},
$$
certifies a variance-adaptive regret bound involving the spectral norm of $\sum_t \mathcal M(X_t)$. The RCASS state is $(a_t,H_t,M_t)$, so storage and computation depend on the maintained aggregates rather than on $t$ [1803.07617].

A common misconception is that online sufficient statistics in this framework are merely bookkeeping devices. In fact, the paper shows that the existence of the sufficient statistic representation and the existence of a Burkholder function are jointly algorithmic and analytic: the same state both certifies the martingale inequality and drives the online strategy [1803.07617].

## 5. Cross-correlation RCASS in stochastic gravitational-wave background searches

For stochastic gravitational-wave background searches with ground-based interferometers, the approximate sufficient statistics are the frequency integrand of the cross-correlation statistic and its variance. The sufficiency is approximate because the analysis uses the weak-signal approximation,
$$
\frac{P_h^2(f)}{P_1(f)P_2(f)}\ll 1,
$$
and replaces the unknown auto-powers by measured estimates from neighboring segments, with a bias-correction factor $(1+2/N_{\rm avg})$ [2012.00907].

At coarse-grained frequency bin $f_\ell$, the segment-level cross-power estimator is
$$
\hat P_{hI}(f_\ell)
=
\frac{2}{T}\,\frac{1}{M}\sum_{k=\ell-M/2}^{\ell+M/2-1}
\mathrm{Re}\!\left(\tilde d_{1I}^*(f_k)\tilde d_{2I}(f_k)\right)\frac{1}{\gamma(f_\ell)},
$$
and in the $\Omega_{\rm gw}(f)$ parameterization the sufficient statistics are
$$
\frac{\hat\Omega_{\rm gw}(f_\ell)}{\bar\sigma_{\rm gw}^2(f_\ell)}
\quad\text{and}\quad
\bar\sigma_{\rm gw}^2(f_\ell).
$$
For a fixed spectral shape $\Omega_{\rm gw}(f)=\Omega_\alpha (f/f_{\rm ref})^\alpha$, they reduce further to the amplitude-level pair $(\hat\Omega_\alpha,\bar\sigma_\alpha^2)$ [2012.00907].

The reduced likelihood depends on the data only through these statistics. In the $\Omega_{\rm gw}$ representation,
$$
p(d\mid \{\bar P_{1I}\},\{\bar P_{2I}\},\Omega_{\rm gw})
=
\prod_\ell \text{const}\times
\exp\!\left[\frac{\hat\Omega_{\rm gw}^2(f_\ell)}{2\bar\sigma_{\rm gw}^2(f_\ell)}\right]
\exp\!\left[-\frac{(\hat\Omega_{\rm gw}(f_\ell)-\Omega_{\rm gw}(f_\ell))^2}{2\bar\sigma_{\rm gw}^2(f_\ell)}\right].
$$
This is precisely the factorization-theorem statement of approximate sufficiency under the reduced model [2012.00907].

The RCASS property is explicit at the segment level. If
$$
C_{I\ell}=\frac{2}{T}\frac{1}{M}\sum_k
\mathrm{Re}\!\left(\tilde d_{1I}^*(f_k)\tilde d_{2I}(f_k)\right)
$$
and
$$
W_{I\ell}
=
\frac{2T\,\delta f}{\bar P_{1I}(f_\ell)\bar P_{2I}(f_\ell)}
\frac{1}{(1+2/N_{\rm avg})},
$$
then the per-frequency accumulated statistics are
$$
A_\ell=\frac{\hat\Omega_{\rm gw}(f_\ell)}{\bar\sigma_{\rm gw}^2(f_\ell)},
\qquad
V_\ell^{-1}=\frac{1}{\bar\sigma_{\rm gw}^2(f_\ell)},
$$
and a new segment updates them by simple addition. The same recursive accumulation holds for the amplitude-only summaries $(A_\alpha,V_\alpha^{-1})$ [2012.00907].

The principal significance of this reduction is methodological equivalence. The paper proves that LIGO–Virgo’s hybrid frequentist–Bayesian analysis—frequentist cross-correlation estimation followed by Bayesian parameter estimation on the integrands and variances—is approximately equivalent to a fully Bayesian analysis on the raw time-series under the stated approximations. The practical gain is extreme data-volume reduction: months-long time-frequency data are compressed to two frequency series, or further to a single amplitude estimate and variance for fixed spectral shape [2012.00907].

The limits of the approximation are also clear. Using the same analysis segment to estimate auto-powers induces bias; the weak-signal approximation breaks down for strong signals; and the treatment does not incorporate non-stationarity beyond the segment model, calibration uncertainties, or correlated non-gravitational-wave noise such as Schumann resonances [2012.00907].

## 6. Rate theory and code-length optimal approximate sufficient statistics

A different branch of the literature studies approximate sufficient statistics as compressed codes for parametric families. For $n$ i.i.d. samples from a $k$-nomial family with $d=k-1$ degrees of freedom, storing the exact sufficient statistic requires code length
$$
d\log n + O(1),
$$
because the number of types is $\Theta(n^d)$. The main result is that if small approximation error is allowed, the code length can be reduced to
$$
\frac{d}{2}\log n + O(1),
$$
and this rate is both achievable and, under the stated criteria, unavoidable [1612.02542].

The paper distinguishes blind and visible encoders and studies errors under KL and TV criteria. Its summarized rate statements include
$$
R_v^{(2)}(\delta_2)=\frac d2 \quad \text{for all } \delta_2\in[0,\infty),
$$
$$
R_b^{(1)}(\delta_1')=R_v^{(1)}(\delta_1)=\frac d2 \quad \text{for all } \delta\in[0,2),
$$
and, under exponential-family structure,
$$
R_b^{(2)}(\delta_2')=\frac d2 \quad \text{for all } \delta_2'\in[0,\infty).
$$
The paper also proves strong converses: even if non-vanishing error is allowed, the pre-log factor cannot be reduced below $d/2$ [1612.02542].

The information-theoretic mechanism behind these rates is local Gaussian geometry. Clarke–Barron asymptotics for the Bayesian mixture give
$$
D(P_\theta^n\|M_\pi^n)
=
\frac d2 \log\!\frac{n}{2\pi}
+\frac 12 \log \det I(\theta)-\log \pi(\theta)+o(1),
$$
while local asymptotic normality yields
$$
D(P_\theta^n\|P_{\tilde\theta}^n)
\approx
\frac n2\, \delta^\top I(\theta)\delta
\quad\text{for } \tilde\theta=\theta+\delta.
$$
Quantization at resolution $\|\delta\|=\Theta(n^{-1/2})$ therefore leads naturally to a grid of size $\Theta(n^{d/2})$ and hence to code length $(d/2)\log n$ [1612.02542].

The RCASS viewpoint enters through recursive maintenance of low-dimensional estimators or sufficient statistics followed by periodic quantization. In exponential families one may maintain
$$
T_n=\sum_{t=1}^n T(X_t),\qquad \eta_n=\frac{T_n}{n},
$$
quantize the expectation parameter on a lattice of span proportional to $n^{-1/2}$, and decode using a conditional mixture over sufficient-statistic fibers. In more general parametric families, one may instead maintain a running estimator and quantize it in the Fisher metric. The paper’s synthesis describes geometric or doubling schedules for emitting the quantized index, keeping the cumulative code length at approximately $(d/2)\log n$ with $O(d)$ working memory [1612.02542].

A recurrent misunderstanding is that approximate sufficient statistics in this sense merely approximate a parameter. The actual object is stronger: a compressed index from which one reconstructs a full distribution on $X^n$ with controlled KL or TV error. The distinction is what makes the strong converse results nontrivial [1612.02542].

## 7. Sufficient-statistic memory AMP and asymptotic RCASS

In approximate message passing for large random linear systems,
$$
y=Ax+n,
$$
state evolution describes the dynamics but does not by itself ensure convergence. Sufficient Statistic Memory Approximate Message Passing (SS-MAMP) addresses this by imposing a sufficient statistic condition on a memory-AMP iteration and constructing it by damping from an arbitrary MAMP [2206.11674].

The general MAMP structure is
$$
r_t=\gamma_t(X_t)=\mathcal Q_t y+\sum_{i=1}^t \mathcal P_{t,i} x_i,
\qquad
x_{t+1}=\phi_t(R_t),
$$
with memory linear estimator $\gamma_t$ and memory nonlinear estimator $\phi_t$. Orthogonality conditions,
$$
\frac1N\langle g_t,x\rangle=0,\qquad
\frac1N\langle g_t,F_t\rangle=0,\qquad
\frac1N\langle f_{t+1},G_t\rangle=0,
$$
play the Onsager-like role that preserves Gaussianity and valid state evolution [2206.11674].

The sufficient statistic condition is
$$
\mathrm{var}\{x\mid R_t\}=\mathrm{var}\{x\mid r_t\},
\qquad
\mathrm{var}\{x\mid y,X_t\}=\mathrm{var}\{x\mid y,x_t\}.
$$
Under the paper’s assumptions, this is equivalent to the covariance matrices being L-banded:
$$
v_{i,j}^\gamma=v_{\max(i,j)}^\gamma,
\qquad
v_{i,j}^\phi=v_{\max(i,j)}^\phi.
$$
For Gaussian observation vectors, the paper proves that the last component is sufficient if and only if the last row and column of the covariance are constant, which is the same L-banded structure in finite-dimensional form [2206.11674].

Once L-bandedness holds, the matrix-valued state evolution collapses to a scalar recursion,
$$
v_t^\gamma=\gamma_{\rm SE}(v_t^\phi),\qquad
v_{t+1}^\phi=\phi_{\rm SE}(v_t^\gamma),\qquad
v_1^\phi=1.
$$
The RCASS is therefore
$$
S_t=(v_t^\gamma,v_t^\phi)\in\mathbb R_+^2,
$$
a finite-dimensional state that fully determines the asymptotic dynamics and is recursively computable. The diagonal variance sequences are monotonically nonincreasing and convergent, so the RCASS recursion itself converges [2206.11674].

The construction from an arbitrary MAMP to SS-MAMP uses optimal damping. If $\tilde R_t$ and $\tilde X_{t+1}$ are the raw outputs, the damped outputs are
$$
r_t=\tilde R_t\,\zeta_t^\gamma,\qquad
x_{t+1}=\tilde X_{t+1}\,\zeta_{t+1}^\phi,
$$
with
$$
\zeta_t^\gamma=\frac{(\tilde V_t^\gamma)^{-1}1}{1^\top(\tilde V_t^\gamma)^{-1}1},
\qquad
\zeta_{t+1}^\phi=\frac{(\tilde V_{t+1}^\phi)^{-1}1}{1^\top(\tilde V_{t+1}^\phi)^{-1}1}.
$$
This damping preserves orthogonality, yields L-banded covariance matrices, and ensures convergence while retaining correct state evolution [2206.11674].

This setting makes explicit an asymptotic version of RCASS. The sufficient statistic is not a classical finite-sample statistic for a parametric likelihood; it is a low-dimensional state that is sufficient in the state-evolution sense, with all asymptotic empirical laws parameterized by $(v_t^\gamma,v_t^\phi)$. The paper’s contribution is to show that, under the sufficient statistic condition, memory and damping can be organized so that convergence and state evolution coexist rather than conflict [2206.11674].

Source: https://www.emergentmind.com/topics/recursively-computable-approximate-sufficient-statistics-rcass