Building a Gaussian hierarchical state-space model¶
We unroll a linear-Gaussian state-space model into affine gates that reuse the same fixed transition parameters at inference time, condition on observations, and query the posterior for an event-time score. To check the work, a dense NumPy recursion plus a Schur complement independently verify the prior propagation and the conditioning algebra.
In this tutorial, we build a linear-Gaussian state-space model in Torx and read an event-time score from its exact posterior.
A linear-Gaussian state-space model turns a multi-channel feature stream into a smoothed latent trajectory. Here, the example is a small brain-computer-interface decoder. Six observation channels stand in for continuous spike features, and a 3-dimensional latent state carries the underlying drive.
Each latent state and each observation is a pmode: a continuous Torx site valued in $\mathbb{R}^N$ and tracked by its mean and covariance.
An SSM is time-translation equivariant, so the same dynamics act between every pair of steps. Torx represents that by explicitly unrolling the sequence into per-step gates while reusing the same transition theta for each transition gate.
The full program has three parts: the shared-parameter transition gates, an initial gate, and per-step emission gates. All three are AffineGaussianGate instances.
By the end, you'll be able to:
- build the model as a stack of
AffineGaussianGateinstances: one initial gate, one gate per transition, and one per emission, - condition on the observations with
AffineGaussianSimulator.conditionto get the exact smoothed posterior, the same posterior that Kalman filtering followed by Rauch-Tung-Striebel smoothing produces, and - read an event-time score from the posterior drive and compare it with a dense NumPy implementation of the LGSSM prior recursion plus conditioning formula.
Setup¶
The setup code is folded away; it only wires imports, paths, and plot styling:
from pathlib import Path
import sys
import jax.numpy as jnp
import numpy as np
ROOT = Path.cwd()
if not (ROOT / "helpers").exists() and (ROOT.parent / "helpers").exists():
ROOT = ROOT.parent
HELPER_DIR = ROOT / "helpers"
sys.path.insert(0, str(HELPER_DIR))
from _notebook_paths import figure_dir
from _notebook_style import apply_notebook_style, make_savefig
import _plots_fields as P_fld
import _plots_schematics as P_sch
from torx.psc import AffineGaussianGate, HybridPCircuit
apply_notebook_style()
FIGURE_DIR = figure_dir(ROOT)
SEED = 23
rng = np.random.default_rng(SEED)
savefig = make_savefig(FIGURE_DIR)
def _det_log_var(dim):
"""Log-variance of a noise-free (deterministic) affine channel."""
return jnp.full((dim,), -jnp.inf)
The model¶
A 3-dimensional latent state (latent_dim = 3) carries the hidden drive.
Six observation channels (num_channels = 6) stand in for continuous features that would arrive after binning and variance stabilization.
The model is linear-Gaussian:
$$ z_0 \sim \mathcal{N}(m_0,\,P_0),\qquad z_t = \underbrace{A\,z_{t-1}}_{\vphantom{\big|}\text{drift}} + \underbrace{\varepsilon_t}_{\vphantom{\big|}\text{process noise}}, $$ $$ x_t = C\,z_t + \eta_t. $$
The initial state draws from the prior $\mathcal{N}(m_0,P_0)$. The process noise is $\varepsilon_t \sim \mathcal{N}(0,Q)$, the measurement noise is $\eta_t \sim \mathcal{N}(0,R)$, and $C$ is the observation map.
The transition sub-product uses the same transition parameters at every adjacent pair:
$$ G_{\mathrm{trans}}(\theta)=\underbrace{\prod_{t=1}^{T-1}}_{\vphantom{\big|}\text{$T-1$ transitions}} \underbrace{K_\theta^{(t-1,\,t)}}_{\vphantom{\big|}\text{same theta reused}} . $$
Each factor $K_\theta^{(t-1,\,t)}$ is the transition kernel from latent site $t-1$ to latent site $t$, and every factor reuses the same theta.
This sharing is inference-time reuse of the same fixed theta values: the built circuit is explicitly unrolled into one gate per transition and one gate per emission, and each gate aliases the same theta dictionary. There is no compact scan primitive and no trainable parameter tying involved. Under a normal JAX optimizer the repeated gates would receive separate gradients, so training would require an explicit gate-to-theta reduction. The initial gate carries the prior $\mathcal{N}(m_0,P_0)$, and the emission gates carry the observation map $C$ and measurement noise $R$.
One length-$T$ sequence is therefore one initial gate, $T-1$ transition gates that reuse the same transition theta, and $T$ emission gates that reuse the same emission theta.
The exact-inference construction below uses log_var = -inf on deterministic passthrough coordinates. That is exact for conditioning, but not differentiable through those coordinates. Clamp to a large finite-negative value before training with gradients.
The next cell sets the dimensions that define the latent and observation sites, plus the fixed intent_axis used later to score the drive.
T = 32
latent_dim = 3
num_channels = 6
latent_sites = list(range(T))
obs_sites = list(range(T, 2 * T))
intent_axis = np.array([1.0, 0.35, -0.25])
The transition matrix A specifies the latent dynamics from one time step to the next.
A = np.array(
[
[0.92, 0.08, 0.00],
[-0.04, 0.88, 0.06],
[0.02, -0.08, 0.90],
]
)
The observation matrix C maps the latent state into six continuous channels.
C = np.array(
[
[1.15, 0.10, -0.15],
[0.85, -0.25, 0.05],
[-0.35, 0.95, 0.20],
[0.20, 0.45, 0.85],
[-0.10, -0.55, 1.05],
[0.55, 0.15, -0.65],
]
)
The prior mean m0 and the diagonal covariance parameters set the initial uncertainty (P0_diag), process noise (Q_diag), and measurement noise (R_diag).
m0 = np.array([-0.45, 0.15, 0.10])
P0_diag = np.array([0.35, 0.25, 0.20])
Q_diag = np.array([0.055, 0.045, 0.040])
R_diag = np.array([0.16, 0.18, 0.20, 0.17, 0.19, 0.16])
simulate_lgssm draws one latent trajectory and its noisy observations with plain NumPy.
def simulate_lgssm():
z = np.zeros((T, latent_dim))
x = np.zeros((T, num_channels))
# Draw the prior state and its first noisy observation.
z[0] = rng.multivariate_normal(m0, np.diag(P0_diag))
x[0] = C @ z[0] + rng.normal(scale=np.sqrt(R_diag))
for t in range(1, T):
# Reuse the same dynamics and noise scale at each later step.
z[t] = A @ z[t - 1] + rng.normal(scale=np.sqrt(Q_diag))
x[t] = C @ z[t] + rng.normal(scale=np.sqrt(R_diag))
return z, x
A synthetic trial supplies observations for conditioning and hidden true_event labels for evaluation.
true_latent, observations = simulate_lgssm()
true_drive = true_latent @ intent_axis
true_event_threshold = 0.0
true_event = true_drive > true_event_threshold
true_event_rate = float(true_event.mean())
majority_baseline_acc = max(true_event_rate, 1.0 - true_event_rate)
print(f"T={T}, latent_dim={latent_dim}, num_channels={num_channels}")
print(f"true events active: {true_event.sum()} / {T} (rate = {true_event_rate:.1%})")
print(f"majority-class baseline accuracy: {majority_baseline_acc:.1%}")
T=32, latent_dim=3, num_channels=6 true events active: 12 / 32 (rate = 37.5%) majority-class baseline accuracy: 62.5%
A shape check catches simulation mistakes before the circuit uses the data.
np.testing.assert_equal(true_latent.shape, (T, latent_dim))
np.testing.assert_equal(observations.shape, (T, num_channels))
The printed true_event_rate gives the majority_baseline_acc used later for the detector comparison.
Building the circuit¶
Each time index has one latent pmode and one observation pmode.
The transition gates implement the kernel $K_\theta^{(t-1,\,t)}$ above; the initial and emission gates carry the prior and the observation model:
- Initial gate ($t=0$): encodes the prior $z_0 \sim \mathcal{N}(m_0, P_0)$ as a single-site
AffineGaussianGate. - Transition gate ($t=1\ldots T-1$): uses a block matrix that copies $z_{t-1}$ and writes $z_t = A z_{t-1}$, with process noise $Q$ on the new site.
- Emission gate ($t=0\ldots T-1$): uses a block matrix that copies $z_t$ and writes $x_t = C z_t$, with measurement noise $R$ on the observation site.
Passthrough coordinates use log_var = -inf (via _det_log_var, a one-line helper defined in the hidden setup), so they are deterministic copies. The block structure carries the model, and the gate list repeats it across time by aliasing the same fixed theta, the inference-time reuse described above.
The matrices transition_matrix and emission_matrix make the deterministic passthrough coordinates explicit.
# First block copies z_{t-1}; second block writes the next latent state.
transition_matrix = np.block(
[
[np.eye(latent_dim), np.zeros((latent_dim, latent_dim))],
[A, np.zeros((latent_dim, latent_dim))],
]
)
# First block keeps z_t; second block writes the predicted observation.
emission_matrix = np.block(
[
[np.eye(latent_dim), np.zeros((latent_dim, num_channels))],
[C, np.zeros((num_channels, num_channels))],
]
)
The initial gate initial_gate writes the Gaussian prior onto the first latent site.
initial_gate = AffineGaussianGate(
sites=[latent_sites[0]],
dims=(latent_dim,),
)
initial_theta = {
"A": jnp.zeros((latent_dim, latent_dim)),
"b": jnp.asarray(m0),
"log_var": jnp.log(jnp.asarray(P0_diag)),
}
The transition gates in transition_gates unroll adjacent time steps while reusing the same transition parameter dictionary.
transition_gates = [
AffineGaussianGate(
sites=[latent_sites[t - 1], latent_sites[t]],
dims=(latent_dim, latent_dim),
)
for t in range(1, T)
]
transition_theta = {
"A": jnp.asarray(transition_matrix),
"b": jnp.zeros(2 * latent_dim),
"log_var": jnp.concatenate(
[
_det_log_var(latent_dim),
jnp.log(jnp.asarray(Q_diag)),
]
),
}
transition_thetas = [transition_theta for _ in range(1, T)]
The emission gates in emission_gates attach one observed pmode to each latent time step.
emission_gates = [
AffineGaussianGate(
sites=[latent_sites[t], obs_sites[t]],
dims=(latent_dim, num_channels),
)
for t in range(T)
]
emission_theta = {
"A": jnp.asarray(emission_matrix),
"b": jnp.zeros(latent_dim + num_channels),
"log_var": jnp.concatenate(
[
_det_log_var(latent_dim),
jnp.log(jnp.asarray(R_diag)),
]
),
}
emission_thetas = [emission_theta for _ in range(T)]
HybridPCircuit combines the gate groups into the full unrolled circuit, with one theta per gate collected in a matching list.
circuit = HybridPCircuit([initial_gate, *transition_gates, *emission_gates])
# Parameters live outside the circuit, one theta per gate in gate order.
thetas = [initial_theta, *transition_thetas, *emission_thetas]
np.testing.assert_equal(len(thetas), len(circuit.gates))
print(f"{len(circuit.gates)} AffineGaussianGate instances")
print(f" 1 initial + {T - 1} transition + {T} emission = {1 + (T - 1) + T}")
64 AffineGaussianGate instances 1 initial + 31 transition + 32 emission = 64
A dimension check confirms the latent and observed site layout in circuit.
np.testing.assert_equal(circuit.continuous_dims, (latent_dim,) * T + (num_channels,) * T)
The printed count confirms that one initial gate, 31 transition gates, and 32 emission gates represent the model.
One slice of the circuit contains a transition gate followed by its emission gate. plot_ssm_slice_circuit draws that local pattern.
fig_c = P_sch.plot_ssm_slice_circuit([("Affine", [0, 1]), ("Affine", [1, 2])])
savefig(fig_c, "11_ssm_slice_circuit")
The slice shows the repeated kernel used to unroll the state-space model over time.
Exact posterior¶
For a linear-Gaussian model, conditioning the joint Gaussian gives the smoothed posterior in closed form, so we do not sample.
AffineGaussianSimulator.condition propagates the joint prior through all observation sites at once and returns the smoothed posterior over the latent sites.
We compile the circuit for AffineGaussianSimulator, then condition on the observed pmodes in observed and query the latent pmodes in latent_sites.
observed = {site: jnp.asarray(observations[t]) for t, site in enumerate(obs_sites)}
from torx.psc import AffineGaussianSimulator
affine_sim = AffineGaussianSimulator()
compiled = affine_sim.build_circuit(circuit, thetas)
# flat initial continuous state over all latent + observed sites
initial_continuous = jnp.zeros(T * latent_dim + T * num_channels)
# Condition on observed sites, then return moments only for the latent sites.
posterior = affine_sim.condition(
compiled,
observations=observed,
initial_continuous=initial_continuous,
query_sites=latent_sites,
)
posterior_mean = np.asarray(posterior.mean).reshape(T, latent_dim)
posterior_cov = np.asarray(posterior.covariance)
A dense NumPy baseline independently propagates the LGSSM moments from (A, C, Q, R, m0, P0) and then applies the Schur-complement conditioning formula. Torx must match both the dense prior and dense posterior to within the relative and absolute tolerances asserted below.
def dense_lgssm_joint_moments(A, C, m0, P0_diag, Q_diag, R_diag, T):
A = np.asarray(A, dtype=float)
C = np.asarray(C, dtype=float)
m0 = np.asarray(m0, dtype=float)
Q = np.diag(np.asarray(Q_diag, dtype=float))
R = np.diag(np.asarray(R_diag, dtype=float))
latent_dim = A.shape[0]
num_channels = C.shape[0]
z_mean = np.zeros((T, latent_dim))
z_cov = np.zeros((T, T, latent_dim, latent_dim))
z_mean[0] = m0
z_cov[0, 0] = np.diag(np.asarray(P0_diag, dtype=float))
for t in range(1, T):
z_mean[t] = A @ z_mean[t - 1]
for s in range(t):
cov_ts = A @ z_cov[t - 1, s]
z_cov[t, s] = cov_ts
z_cov[s, t] = cov_ts.T
z_cov[t, t] = A @ z_cov[t - 1, t - 1] @ A.T + Q
total_dim = T * latent_dim + T * num_channels
mean = np.zeros(total_dim)
cov = np.zeros((total_dim, total_dim))
obs_mean = z_mean @ C.T
mean[: T * latent_dim] = z_mean.reshape(-1)
mean[T * latent_dim :] = obs_mean.reshape(-1)
def z_slice(t):
return slice(t * latent_dim, (t + 1) * latent_dim)
def x_slice(t):
start = T * latent_dim + t * num_channels
return slice(start, start + num_channels)
for t in range(T):
for s in range(T):
cov_zz = z_cov[t, s]
cov[z_slice(t), z_slice(s)] = cov_zz
cov[z_slice(t), x_slice(s)] = cov_zz @ C.T
cov[x_slice(t), z_slice(s)] = C @ cov_zz
cov_xx = C @ cov_zz @ C.T
if t == s:
cov_xx = cov_xx + R
cov[x_slice(t), x_slice(s)] = cov_xx
return mean, cov
First, propagate the Torx prior through every site and check it against the dense joint moments.
prior = affine_sim.propagate(compiled, initial_continuous)
dense_prior_mean, dense_prior_cov = dense_lgssm_joint_moments(
A, C, m0, P0_diag, Q_diag, R_diag, T
)
np.testing.assert_equal(prior.sites, tuple(latent_sites + obs_sites))
np.testing.assert_allclose(prior.mean, dense_prior_mean, rtol=1e-5, atol=3e-5)
np.testing.assert_allclose(prior.covariance, dense_prior_cov, rtol=1e-5, atol=3e-5)
Then condition the dense joint on the observations. The Schur complement gives the reference posterior mean and covariance for the Torx result to match.
latent_total = T * latent_dim
query_idx = np.arange(latent_total)
obs_idx = np.arange(latent_total, latent_total + T * num_channels)
y = observations.reshape(-1)
# condition returns query_sites in latent_sites order, matching query_idx.
Soo = dense_prior_cov[np.ix_(obs_idx, obs_idx)]
Sqo = dense_prior_cov[np.ix_(query_idx, obs_idx)]
Sqq = dense_prior_cov[np.ix_(query_idx, query_idx)]
dense_mean = dense_prior_mean[query_idx] + Sqo @ np.linalg.solve(
Soo, y - dense_prior_mean[obs_idx]
)
dense_cov = Sqq - Sqo @ np.linalg.solve(Soo, Sqo.T)
np.testing.assert_allclose(posterior.mean, dense_mean, rtol=1e-5, atol=3e-5)
np.testing.assert_allclose(posterior.covariance, dense_cov, rtol=1e-5, atol=3e-5)
np.testing.assert_equal(posterior_mean.shape, (T, latent_dim))
np.testing.assert_array_less(-1e-6, np.linalg.eigvalsh(posterior_cov))
A structural check bins the smoothed covariance into per-timestep blocks. Most of the mass should sit on the diagonal, with the nearest-neighbor lag dominating the off-diagonal terms.
block_norms = np.array(
[
[
np.linalg.norm(
posterior_cov[
t * latent_dim : (t + 1) * latent_dim,
s * latent_dim : (s + 1) * latent_dim,
]
)
for s in range(T)
]
for t in range(T)
]
)
lag_mean_norms = np.array(
[np.mean([block_norms[t, t + lag] for t in range(T - lag)]) for lag in range(T)]
)
offdiag_lag_means = lag_mean_norms[1:]
np.testing.assert_equal(int(np.argmax(offdiag_lag_means)), 0)
np.testing.assert_array_less(np.mean(lag_mean_norms[4:]), lag_mean_norms[1])
print(f"mean block norms by lag 0..4: {lag_mean_norms[:5].round(4).tolist()}")
mean block norms by lag 0..4: [0.056699998676776886, 0.027799999341368675, 0.014000000432133675, 0.0071000000461936, 0.003700000001117587]
The next cell computes the model log-evidence from the dense reference: the log-probability of the observed data under the model.
# Soo is an SPD observation covariance, so take its log-det from the Cholesky
# factor: the stable SPD path, and it raises LinAlgError if Soo is ever not
# positive definite instead of returning a quietly wrong number.
chol_Soo = np.linalg.cholesky(Soo)
logdet = float(2.0 * np.sum(np.log(np.diag(chol_Soo))))
y_centered = y - dense_prior_mean[obs_idx]
log_evidence = float(
-0.5
* (
len(y) * np.log(2 * np.pi)
+ logdet
+ y_centered @ np.linalg.solve(Soo, y_centered)
)
)
print("Torx matches the independent dense LGSSM baseline")
print(f"log-evidence: {log_evidence:.3f}")
Torx matches the independent dense LGSSM baseline log-evidence: -152.580
The printed log_evidence scores how probable the observed trial is under this model, and it is the number to compare across candidate models.
Detecting events on the posterior¶
The downstream quantity is the time at which the latent drive crosses into an active state.
Projecting posterior_mean onto intent_axis gives a scalar latent drive, and each step's per-site covariance marginal (its own spread, ignoring the other steps) gives its variance.
intent_axis is fixed in advance and is the same vector that defines the ground-truth labels, so the test measures how well the posterior recovers the drive along a direction that is given in advance.
A threshold at the median of posterior_drive turns that drive into a binary event predictor. The median is observable at deployment, so the hidden labels play no part in choosing it.
The posterior drive posterior_drive and its standard deviation posterior_drive_std come from the smoothed latent moments.
posterior_drive = posterior_mean @ intent_axis
posterior_drive_var = np.array(
[
# Project each smoothed covariance onto the same one-dimensional intent axis.
intent_axis @ np.asarray(posterior.site_moments(s)[1]) @ intent_axis
for s in latent_sites
]
)
posterior_drive_std = np.sqrt(np.maximum(posterior_drive_var, 0.0))
The detector thresholds the posterior drive at its median. The logistic transform is only a display score (the figure's y-axis labels it P(event)); the 0.5 decision is identical to posterior_drive > detection_threshold.
detection_threshold = float(np.median(posterior_drive))
# Center the logistic display score at the median; the decision is the same threshold on posterior_drive.
event_score = 1.0 / (1.0 + np.exp(-3.0 * (posterior_drive - detection_threshold)))
predicted_event = posterior_drive > detection_threshold
With predictions in hand, we score the detector against the hidden labels. Accuracy, a tie-aware AUC, and precision and recall summarize the fit, and the last check gates on the sanity thresholds.
accuracy = float(np.mean(predicted_event == true_event))
pos = posterior_drive[true_event]
neg = posterior_drive[~true_event]
# AUC and precision/recall are undefined without both classes present.
assert pos.size and neg.size, "metrics require both event and non-event steps"
# tie-aware AUC (Mann-Whitney): ties between a positive and negative get half credit
gt = float((pos[:, None] > neg[None, :]).mean())
eq = float((pos[:, None] == neg[None, :]).mean())
auc = gt + 0.5 * eq
tp = int(np.sum(predicted_event & true_event))
fp = int(np.sum(predicted_event & ~true_event))
fn = int(np.sum(~predicted_event & true_event))
# report undefined metrics as nan rather than masking an empty denominator
precision = tp / (tp + fp) if (tp + fp) else float("nan")
recall = tp / (tp + fn) if (tp + fn) else float("nan")
if accuracy < majority_baseline_acc + 0.05:
raise AssertionError(
f"detector accuracy {accuracy:.2%} barely beats majority "
f"({majority_baseline_acc:.2%}); posterior isn't informative enough."
)
if auc < 0.80:
raise AssertionError(f"posterior-drive AUC {auc:.2f} below 0.80")
print(f"detection threshold (median posterior): {detection_threshold:.4f}")
print(
f"detector accuracy: {accuracy:.1%} "
f"(vs majority baseline {majority_baseline_acc:.1%})"
)
print(f"AUC (threshold-independent): {auc:.3f}")
print(f"precision: {precision:.2f} recall: {recall:.2f}")
detection threshold (median posterior): -0.0654 detector accuracy: 81.2% (vs majority baseline 62.5%) AUC (threshold-independent): 0.921 precision: 0.69 recall: 0.92
The evaluation output reports accuracy, auc, precision, and recall for this single seeded trajectory. Treat these checks as a sanity gate for the notebook example; one seeded trajectory cannot benchmark detector performance.
The first plot shows the six channels seen by the multi-channel front end, with gray bands on time steps where the true drive is positive.
fig = P_fld.plot_observations(observations, true_event)
savefig(fig, "11_ssm_observations")
The gray bands identify the active event steps used only for evaluation.
The next plot compares true_drive with the Torx-smoothed posterior_drive both already projected onto intent_axis. The figure also shows the $\pm 2\sigma$ band and the dashed detection_threshold.
The band shows $\pm 2$ marginal posterior standard deviations for the projected drive at each step, read from the same per-site moments as in 10_pmode_gaussian_gates.ipynb, here swept over time.
fig = P_fld.plot_posterior_drive(
true_drive,
posterior_drive,
posterior_drive_std,
detection_threshold,
)
savefig(fig, "11_ssm_posterior_drive")
The posterior mean tracks the latent drive, and the uncertainty band stays tight in the interior. It widens near the start and end, where the smoother has less context.
The covariance plot shows the block structure induced by smoothing. Diagonal blocks describe same-step uncertainty. The nearest-neighbor band is the strongest off-diagonal coupling, and longer-lag block norms are smaller in the check above.
fig = P_fld.plot_posterior_covariance(posterior_cov)
savefig(fig, "11_ssm_covariance")
The covariance remains concentrated near the diagonal, as expected for a local dynamical model.
The final plot shows the event-time output used by a downstream consumer. The top panel shows the logistic display score from the smoothed posterior drive against the 0.5 decision boundary, and the bottom panel compares predicted_event with true_event.
time = np.arange(T)
fig = P_fld.plot_detector(time, event_score, true_event, predicted_event, accuracy)
savefig(fig, "11_ssm_detector")
The predicted events align with most active intervals and stay above the majority-class baseline.
Conclusion¶
We built a linear-Gaussian state-space model in Torx and read an event-time score from its exact posterior.
- The linear-Gaussian model is represented as 64
AffineGaussianGateinstances: one initial gate plus one gate per transition and per emission. - Transition and emission gates alias the same fixed
thetadictionaries at inference time, while the circuit itself is explicitly unrolled. AffineGaussianSimulator.conditionreturns the exact smoothed posterior, matching an independent dense NumPy prior recursion plus Schur-complement baseline to a relative tolerance of1e-5.- A median threshold on
posterior_drivebeats the majority-class accuracy baseline by at least five percentage points and exceeds the AUC sanity threshold (auc >= 0.80) for this single seeded check, as the detector cell prints. - The same block-matrix recipe, with
-inflog_varon passthrough coordinates for exact inference, expresses affine relations between sites. Clamp those passthrough variances before gradient training.
This tutorial applies the affine-Gaussian density layer of 10_pmode_gaussian_gates.ipynb to a time series. To continue, work through 12_langevin_graph_ising.ipynb.
References¶
- Kalman, R.E. 1960. A new approach to linear filtering and prediction problems. Trans. ASME J. Basic Eng. 82(1), 35-45.
- Rauch, H.E., Tung, F., Striebel, C.T. 1965. Maximum likelihood estimates of linear dynamic systems. AIAA J. 3(8), 1445-1450.