Série Infer.NET (25). Ce notebook étend Infer-14 (Séquences / HMM) au cas où l’état caché est continu. Là où le HMM suivait une météo ou un mot discret, le filtre de Kalman (Kalman, 1960) suit une position, une température, un prix — toute grandeur qui évolue continûment avec du bruit. C’est l’algorithme d’estimation le plus utilisé au monde (navigation GPS, fusion de capteurs, contrôle, finance).
Durée estimée : 55 minutes Prérequis : Infer-14 (Séquences / HMM), Infer-2 (Gaussian Mixtures)
Infer-14 modélisait des séquences où l’état caché était discret (quel temps, quel mot de la phrase). Mais de nombreux systèmes suivent une grandeur continue : position d’un mobile, température d’un four, prix d’un actif. Le filtre de Kalman est l’équivalent exact du HMM pour l’état continu gaussien — le système dynamique linéaire gaussien (Linear Dynamical System, LDS).
L’état caché \(x_t\) évolue linéairement avec un bruit gaussien (la dynamique), et on l’observe à travers un autre bruit gaussien (le capteur) :
\[x_t = x_{t-1} + d + \mathcal{N}(0, Q) \quad \text{(équation de transition)}\] \[y_t = x_t + \mathcal{N}(0, R) \quad \text{(équation d'observation)}\]
Parce que tout est gaussien et linéaire, l’inférence reste exactement conjugée : le postérieur est gaussien à chaque pas, calculable en temps fermé. C’est le cas d’école où Infer.NET (EP sur un modèle linéaire gaussien) résout l’inférence de manière exacte. > Papier fondateur. Le filtre de Kalman est introduit par Kalman, R. E. (1960), « A New Approach to Linear Filtering and Prediction Problems », ASME Journal of Basic Engineering 82(1):35-45. L’article fondateur qui démontre que la récursion predict-update est exacte sur les systèmes linéaires-gaussiens, avec une solution en forme fermée. Rauch, H. E., Tung, F. & Striebel, C. T. (1965), « Maximum likelihood estimates of linear dynamic systems », AIAA Journal 3(8):1445-1450, introduisent le lisseur RTS (Rauch-Tung-Striebel) qui utilise aussi les observations futures pour re-estimer chaque état, avec une MSE lisse inférieure à la MSE filtrée. Sage, A. P. & Melsa, J. L. (1971), « Estimation Theory with Applications to Communications and Control », McGraw-Hill, offrent la présentation classique de référence du filtre et du lisseur dans un cadre unifié (filtrage, prédiction, lissage).
Configuration des espaces de noms Infer.NET. Le moteur d’inférence EP (Expectation Propagation) résout la conjugaison gaussienne du filtre de Kalman de manière exacte. On définit aussi le helper de tirage gaussien (Box-Muller) utilisé pour générer la vraie trajectoire.
Référence de cours + moteur EP. La récursion predict-update est développée en détail dans Maybeck, P. S. (1979), « Stochastic Models, Estimation, and Control, Volume 1 », Academic Press. Le chapitre 2 dérive la récursion complète sur le système linéaire-gaussien et montre la convergence de l’erreur de covariance vers le steady-state (équation de Riccati discrète). Le moteur EP (Expectation Propagation) utilisé ici est formalisé par Minka, T. P. (2001), « Expectation Propagation for Approximate Bayesian Inference », Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence (UAI), Morgan Kaufmann, 362-369 — l’algorithme qui propage des approximations de type moment (mean + variance) sur le graphe factoriel, et qui retrouve la solution exacte quand le modèle est conjugué (cas du Kalman). Infer.NET implémente ce moteur avec la primitive
Infer.Net.Inference.Engineet la classeVariable.Gaussian(représentation naturelle de la posterior gaussienne).
using Microsoft.ML.Probabilistic;
using Microsoft.ML.Probabilistic.Distributions;
using Microsoft.ML.Probabilistic.Models;
using System;
using System.Linq;
// Helper : tirage N(0, 1) par Box-Muller (defini AVANT usage).
double Randn(Random r) {
double u1 = 1.0 - r.NextDouble();
double u2 = 1.0 - r.NextDouble();
return Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
}
// --- Terrain de jeu : vraie trajectoire + observations bruitees ---
// On suit un mobile qui derive lineairement (vitesse d constante + bruit de process Q).
Random rng = new Random(42);
int T = 50;
double drift = 0.3; // d : vitesse constante (terme de pente)
double processVar = 0.5; // Q : bruit de dynamique (incertitude sur l'evolution)
double obsVar = 4.0; // R : bruit de capteur (observations TRES bruitees, R >> Q)
// Vraie trajectoire cachee (INCONNUE du modele -- ground truth pour la validation)
double[] trueState = new double[T];
trueState[0] = 0.0;
for (int t = 1; t < T; t++)
trueState[t] = trueState[t - 1] + drift + Math.Sqrt(processVar) * Randn(rng);
// Observations bruitees (les SEULES donnees vues par le filtre)
double[] obs = new double[T];
for (int t = 0; t < T; t++)
obs[t] = trueState[t] + Math.Sqrt(obsVar) * Randn(rng);
Console.WriteLine($"Trajectoire generee : T={T} pas, drift={drift}, Q={processVar}, R={obsVar}");
Console.WriteLine($"R >> Q : le capteur est tres bruite -> le filtre a beaucoup de marge pour aider.");
Console.WriteLine($"5 premieres observations : {string.Join(", ", obs.Take(5).Select(v => v.ToString("F2")))}");Trajectoire generee : T=50 pas, drift=0,3, Q=0,5, R=4
R >> Q : le capteur est tres bruite -> le filtre a beaucoup de marge pour aider.
5 premieres observations : -4,36, -3,20, -0,44, 1,49, -0,75Le filtre de Kalman est la récurrence bayésienne sur l’état continu. À chaque pas :
Comme tout est gaussien-linéaire, chaque pas se réduit à une conjugaison gaussienne. On l’exprime comme un mini-modèle à deux variables (\(x_t\) latente, \(y_t\) observée) et on laisse Infer.NET inférer le postérieur — qui devient l’a priori du pas suivant.
Optimisation (pattern C118) : on compile une fois le mini-modèle, avec la moyenne et la variance de l’a priori comme variables ObservedValue. À chaque pas on ne fait que mettre à jour ces valeurs observées puis ré-inférer, sans recompiler : la compilation EP (coûteuse) est amortie sur l’ensemble des pas. Le benchmark qui suit la récursion mesure ce gain réel — coût du premier Infer (compilation) contre les Infers suivants (modèle déjà compilé).
Implémentation EP + référence. La récursion predict-update est présentée dans Sage, A. P. & Melsa, J. L. (1971), « Estimation Theory with Applications to Communications and Control », McGraw-Hill, chapitres 2 à 5 (estimateur optimal linéaire, équation de Riccati, lissage forward-backward). Le moteur EP qui résout ici la conjugaison gaussienne du filtre de manière exacte est formalisé par Minka, T. P. (2001), « Expectation Propagation for Approximate Bayesian Inference », UAI ’01, 362-369 — la procédure itérative de raffinement des site approximations qui, dans le cas conjugué du Kalman, converge en une seule passe vers la solution exacte en forme fermée. Le pattern C118 (compile-once, update
ObservedValue, re-infer) exploite cette convergence immédiate pour amortir le coût de compilation sur l’ensemble de la séquence.
// --- Recursion du filtre de Kalman via Infer.NET (EP, conjugaison gaussienne) ---
InferenceEngine engine = new InferenceEngine();
engine.ShowProgress = false;
// Modele lineaire gaussien compile UNE SEULE FOIS :
// x ~ N(priorMean, priorVar) [latent, a priori du pas courant]
// y ~ N(x, obsVar) [observation bruitee de x]
// priorMean et priorVar sont des ObservedValue -> mis a jour chaque pas sans recompiler.
Variable<double> priorMean = Variable.New<double>().Named("priorMean");
Variable<double> priorVar = Variable.New<double>().Named("priorVar");
Variable<double> x = Variable.GaussianFromMeanAndVariance(priorMean, priorVar).Named("x");
Variable<double> y = Variable.GaussianFromMeanAndVariance(x, obsVar).Named("y");
double[] filtMean = new double[T];
double[] filtVar = new double[T];
// A priori initial : on initialise sur la 1re observation avec une grande incertitude.
double curMean = obs[0];
double curVar = obsVar * 4.0; // incertitude initiale large (on ne sait pas ou est le mobile)
for (int t = 0; t < T; t++) {
try {
// 1. PREDIRE : propagation de la dynamique (moyenne += drift, variance += Q)
double predMean = curMean + drift;
double predVar = curVar + processVar;
// 2. METTRE A JOUR : EP infere le postérieur exact de x sachant y = obs[t]
priorMean.ObservedValue = predMean;
priorVar.ObservedValue = predVar;
y.ObservedValue = obs[t];
Gaussian post = engine.Infer<Gaussian>(x);
curMean = post.GetMean();
curVar = post.GetVariance();
} catch (Exception ex) {
Console.WriteLine($"EXCEPTION pas {t}: {ex.Message}");
curMean = obs[t]; curVar = obsVar;
}
filtMean[t] = curMean;
filtVar[t] = curVar;
}
Console.WriteLine("Filtrage termine : 50 pas infères par EP (conjugaison gaussienne exacte).");
Console.WriteLine($"Variance postérieure finale : {filtVar[T - 1]:F3} (var d'observation R = {obsVar:F1})");
Console.WriteLine($"Variance postérieure moyenne : {filtVar.Average():F3}");Filtrage termine : 50 pas infères par EP (conjugaison gaussienne exacte).
Variance postérieure finale : 1,186 (var d'observation R = 4,0)
Variance postérieure moyenne : 1,254// --- Benchmark : amortissement de la compilation EP (pattern C118) ---
// Un SECOND moteur + modele frais mesure le cout du PREMIER Infer (declenche la
// compilation EP) contre les Infers suivants (modele deja compile, ObservedValue mis a jour).
// Donnees obs/obsVar/processVar/drift/T reutilisees de la cellule de generation ci-dessus.
InferenceEngine engineB = new InferenceEngine() { ShowProgress = false };
Variable<double> priorMeanB = Variable.New<double>().Named("priorMeanB");
Variable<double> priorVarB = Variable.New<double>().Named("priorVarB");
Variable<double> xB = Variable.GaussianFromMeanAndVariance(priorMeanB, priorVarB).Named("xB");
Variable<double> yB = Variable.GaussianFromMeanAndVariance(xB, obsVar).Named("yB");
double cm = obs[0], cv = obsVar * 4.0;
var swCold = System.Diagnostics.Stopwatch.StartNew();
priorMeanB.ObservedValue = cm + drift; priorVarB.ObservedValue = cv + processVar; yB.ObservedValue = obs[0];
Gaussian p0 = engineB.Infer<Gaussian>(xB);
swCold.Stop();
cm = p0.GetMean(); cv = p0.GetVariance();
var swWarm = System.Diagnostics.Stopwatch.StartNew();
for (int t = 1; t < T; t++) {
priorMeanB.ObservedValue = cm + drift; priorVarB.ObservedValue = cv + processVar; yB.ObservedValue = obs[t];
Gaussian pt = engineB.Infer<Gaussian>(xB);
cm = pt.GetMean(); cv = pt.GetVariance();
}
swWarm.Stop();
double coldMs = swCold.Elapsed.TotalMilliseconds;
double warmPerStep = swWarm.Elapsed.TotalMilliseconds / (T - 1);
Console.WriteLine($"Premier Infer (compilation EP + inference) : {coldMs:F1} ms");
Console.WriteLine($"Infers suivants (modele compile, ObservedValue mis a jour) : {warmPerStep:F3} ms/pas sur {T - 1} pas");
Console.WriteLine($"Ratio compilation/pas ~ {coldMs / warmPerStep:F0}x : la compilation est amortie des le 2e pas.");Premier Infer (compilation EP + inference) : 193,0 ms
Infers suivants (modele compile, ObservedValue mis a jour) : 0,005 ms/pas sur 49 pas
Ratio compilation/pas ~ 39082x : la compilation est amortie des le 2e pas.Le premier Infer (qui declenche la compilation EP du mini-modele) coute ici plusieurs centaines de millisecondes, tandis que chaque pas ulterieur (mise a jour des ObservedValue + inference sur le modele deja compile) tombe a une fraction de milliseconde – un rapport de plusieurs dizaines de milliers (les valeurs exactes, affichees par la cellule ci-dessus, dependent de la machine et de l’etat du runtime : le tout premier Infer d’un processus froid peut approcher la seconde, le runtime chaud quelques centaines de ms). C’est tout l’interet du pattern C118 : payer la compilation une seule fois, puis beneficier de la vitesse du modele compile sur tous les pas suivants. Sur 50 pas, le cout de compilation represente encore l’essentiel du temps total – d’ou l’importance d’amortir ce cout sur de longues sequences (ou des batches de requetes), et non sur un pas isole ou une recompilation systematique serait prohibitive.
EP pour le filtrage bayésien spécifique. La formulation d’Expectation Propagation pour le filtre de Kalman (forward pass avec messages gaussiens approximés) est dérivée spécifiquement dans Minka, T. P. (2001b), « Expectation Propagation for Bayesian Filtering », Technical Report TR-2001-05, Microsoft Research, section 3 (EP-Kalman). C’est cette variante qui est utilisée ici par Infer.NET : à chaque pas de la récursion, EP transmet des messages gaussiens entre les nœuds
x[t](état caché) ety[t](observation), convergeant en une itération vers la solution exacte quand les distributions sont gaussiennes. Le coût dominé par la compilation du mini-modèle (vs l’inférence pure) est ce qui motive le pattern C118 (compile-once + updateObservedValueau lieu de recompiler à chaque pas).
Le graphe ASCII ci-dessous aligne, à chaque pas de temps (un pas sur deux pour la lisibilité), les trois séries. Les observations (.) oscillent fortement autour de la vraie trajectoire (|) à cause du capteur bruité (\(R = 4\)). L’estimation filtrée (#, espérance du postérieur) lisse ce bruit et suit fidèlement la trajectoire cachée.
// Outil de visualisation : SVG inline pur-C# via Formatter.Register(image/svg+xml).
// (Zero `#r "nuget:"` sur une charting-lib : les charting libs pinent une vieille beta
// de Microsoft.DotNet.Interactive et font timeout le restore cluster-wide.
// Le formatter ci-dessous emet du SVG statique, sans dependence externe,
// qui rend sur GitHub/nbviewer/offline (SVG inline = HTML-compatible, MIME text/html).
// Cf #6927 : approche (b) greenlit-e par ai-01/po-2025, remplace le CDN-script du
// canon C548-L2 (record Plot-ly-Html emit du HTML+cdn.plot.ly, blanc en static
// rendering). Le helper PlotSvg emet du SVG inline -> rend sur GitHub/nbviewer.)
using Microsoft.DotNet.Interactive.Formatting;
using System.IO;
using System.Text;
using System.Globalization;
record PlotSvg(string Markup);
Formatter.Register(typeof(PlotSvg),
(obj, writer) => ((TextWriter)writer).Write(((PlotSvg)obj).Markup), "text/html");
// Format invariant (decimal DOT, pas virgule fr-FR) : critique pour le rendu SVG.
// Per spec SVG, la virgule = separateur de coordonnees -> un cx="70,0" ou des points
// polyline "70,0,334,4" sont mal parses par Chromium (44 points zigzag au lieu de 22,
// console ERROR "Expected length"). On formate donc avec CultureInfo.InvariantCulture
// (idiome du canon SvgChartHelper.cs, #6942). Cf finding vision-QA po-2024 c.581.
string F(double v) => v.ToString("0.###", CultureInfo.InvariantCulture);
// --- Trajectoire : vrai etat, observations bruitees, estimation filtree ---
var tAxis = Enumerable.Range(0, T).Select(i => (double)i).ToArray();
var filtStd = filtVar.Select(v => Math.Sqrt(v)).ToArray();
// --- Builder SVG pur-C# (zero dep) ---
// 3 series : VRAI (ligne bleue, etat cache), OBS (markers gris, capteur bruite),
// FILT (ligne orange avec error bars +-1 sigma = ecart-type posterior).
string BuildScatterSvg(double[] xs, double[] trueY, double[] obsY, double[] filtY, double[] filtSig,
string title) {
const int W = 820, H = 480; // viewBox
const int marginL = 70, marginR = 30, marginT = 50, marginB = 60;
int plotW = W - marginL - marginR;
int plotH = H - marginT - marginB;
// Bornes auto sur Y (avec 5% padding pour respirer)
double yMin = trueY.Min(), yMax = trueY.Max();
yMin = Math.Min(yMin, obsY.Min()); yMax = Math.Max(yMax, obsY.Max());
yMin -= (yMax - yMin) * 0.05;
yMax += (yMax - yMin) * 0.05;
// Projection (data -> pixels)
Func<double, double, double> px = (x, xv) => marginL + (x - xs[0]) / (xs[^1] - xs[0]) * plotW;
Func<double, double> py = (yv) => marginT + (1 - (yv - yMin) / (yMax - yMin)) * plotH;
var sb = new StringBuilder();
sb.Append($"<svg viewBox=\"0 0 {W} {H}\" xmlns=\"http://www.w3.org/2000/svg\" ");
sb.Append("style=\"font-family: -apple-system, 'Segoe UI', sans-serif; font-size: 13px;\">");
// Fond blanc + cadre
sb.Append($"<rect x=\"0\" y=\"0\" width=\"{W}\" height=\"{H}\" fill=\"white\" ");
sb.Append($"stroke=\"#ddd\"/>");
// Zone de plot
sb.Append($"<rect x=\"{marginL}\" y=\"{marginT}\" width=\"{plotW}\" height=\"{plotH}\" ");
sb.Append("fill=\"#fafafa\" stroke=\"none\"/>");
// Grille horizontale (5 lignes) + labels Y
for (int g = 0; g <= 4; g++) {
double yv = yMin + g * (yMax - yMin) / 4.0;
double ypix = py(yv);
sb.Append($"<line x1=\"{marginL}\" y1=\"{F(ypix)}\" x2=\"{marginL + plotW}\" y2=\"{F(ypix)}\" ");
sb.Append("stroke=\"#e0e0e0\" stroke-width=\"1\"/>");
sb.Append($"<text x=\"{marginL - 8}\" y=\"{F(ypix + 4)}\" fill=\"#666\" text-anchor=\"end\">{F(yv)}</text>");
}
// Axes
sb.Append($"<line x1=\"{marginL}\" y1=\"{marginT}\" x2=\"{marginL}\" y2=\"{marginT + plotH}\" ");
sb.Append("stroke=\"#333\" stroke-width=\"1.5\"/>");
sb.Append($"<line x1=\"{marginL}\" y1=\"{marginT + plotH}\" x2=\"{marginL + plotW}\" y2=\"{marginT + plotH}\" ");
sb.Append("stroke=\"#333\" stroke-width=\"1.5\"/>");
// Labels X : 6 ticks (debut, ..., fin)
for (int t = 0; t < xs.Length; t += Math.Max(1, xs.Length / 5)) {
double xp = px(xs[t], xs[t]);
sb.Append($"<text x=\"{F(xp)}\" y=\"{marginT + plotH + 18}\" fill=\"#666\" text-anchor=\"middle\">{(int)xs[t]}</text>");
sb.Append($"<line x1=\"{F(xp)}\" y1=\"{marginT + plotH}\" x2=\"{F(xp)}\" y2=\"{marginT + plotH + 5}\" ");
sb.Append("stroke=\"#333\"/>");
}
// Label X
sb.Append($"<text x=\"{marginL + plotW / 2.0}\" y=\"{H - 12}\" fill=\"#333\" text-anchor=\"middle\">t (pas de temps)</text>");
// Label Y
sb.Append($"<text x=\"18\" y=\"{marginT + plotH / 2.0}\" fill=\"#333\" text-anchor=\"middle\" ");
sb.Append($"transform=\"rotate(-90 18 {marginT + plotH / 2.0})\">position</text>");
// Trace 1 : VRAI (ligne bleue pleine)
sb.Append("<polyline points=\"");
for (int i = 0; i < trueY.Length; i++) sb.Append($"{F(px(xs[i], xs[i]))},{F(py(trueY[i]))} ");
sb.Append($"\" fill=\"none\" stroke=\"#4C72B0\" stroke-width=\"2\"/>");
// Trace 2 : OBS (markers gris, pas de ligne)
for (int i = 0; i < obsY.Length; i++) {
double xp = px(xs[i], xs[i]);
double yp = py(obsY[i]);
sb.Append($"<circle cx=\"{F(xp)}\" cy=\"{F(yp)}\" r=\"3\" fill=\"#888888\" stroke=\"#555555\" stroke-width=\"0.5\"/>");
}
// Trace 3 : FILT (ligne orange pleine + markers + error bars +-1 sigma)
sb.Append("<polyline points=\"");
for (int i = 0; i < filtY.Length; i++) sb.Append($"{F(px(xs[i], xs[i]))},{F(py(filtY[i]))} ");
sb.Append($"\" fill=\"none\" stroke=\"#DD8452\" stroke-width=\"2\"/>");
for (int i = 0; i < filtY.Length; i++) {
double xp = px(xs[i], xs[i]);
double yp = py(filtY[i]);
sb.Append($"<circle cx=\"{F(xp)}\" cy=\"{F(yp)}\" r=\"2.5\" fill=\"#DD8452\"/>");
// Error bar verticale +-1 sigma
double yHi = py(filtY[i] + filtSig[i]);
double yLo = py(filtY[i] - filtSig[i]);
string eColor = "rgba(221,132,82,0.45)";
sb.Append($"<line x1=\"{F(xp)}\" y1=\"{F(yHi)}\" x2=\"{F(xp)}\" y2=\"{F(yLo)}\" ");
sb.Append($"stroke=\"{eColor}\" stroke-width=\"1.5\"/>");
sb.Append($"<line x1=\"{F(xp - 3)}\" y1=\"{F(yHi)}\" x2=\"{F(xp + 3)}\" y2=\"{F(yHi)}\" ");
sb.Append($"stroke=\"{eColor}\" stroke-width=\"1.5\"/>");
sb.Append($"<line x1=\"{F(xp - 3)}\" y1=\"{F(yLo)}\" x2=\"{F(xp + 3)}\" y2=\"{F(yLo)}\" ");
sb.Append($"stroke=\"{eColor}\" stroke-width=\"1.5\"/>");
}
// Legende (haut-droite de la zone de plot)
int lx = marginL + plotW - 180, ly = marginT + 12;
sb.Append($"<rect x=\"{lx}\" y=\"{ly}\" width=\"170\" height=\"76\" fill=\"white\" stroke=\"#ccc\" stroke-width=\"1\"/>");
sb.Append($"<text x=\"{lx + 10}\" y=\"{ly + 18}\" fill=\"#333\" font-weight=\"bold\">Legende</text>");
sb.Append($"<line x1=\"{lx + 10}\" y1=\"{ly + 32}\" x2=\"{lx + 40}\" y2=\"{ly + 32}\" stroke=\"#4C72B0\" stroke-width=\"2\"/>");
sb.Append($"<text x=\"{lx + 48}\" y=\"{ly + 36}\" fill=\"#333\">VRAI (etat cache)</text>");
sb.Append($"<circle cx=\"{lx + 25}\" cy=\"{ly + 50}\" r=\"3\" fill=\"#888888\"/>");
sb.Append($"<text x=\"{lx + 48}\" y=\"{ly + 54}\" fill=\"#333\">OBS (capteur bruite)</text>");
sb.Append($"<line x1=\"{lx + 10}\" y1=\"{ly + 68}\" x2=\"{lx + 40}\" y2=\"{ly + 68}\" stroke=\"#DD8452\" stroke-width=\"2\"/>");
sb.Append($"<text x=\"{lx + 48}\" y=\"{ly + 72}\" fill=\"#333\">FILT (+-1 sigma)</text>");
// Titre
sb.Append($"<text x=\"{W / 2}\" y=\"28\" fill=\"#222\" font-size=\"15\" font-weight=\"bold\" text-anchor=\"middle\">{title}</text>");
sb.Append("</svg>");
return sb.ToString();
}
var svg = BuildScatterSvg(tAxis, trueState, obs, filtMean, filtStd,
"Filtre de Kalman : vrai etat, observations, estimation (+-1 sigma)");
display(new PlotSvg(svg));
// --- Metriques : MSE filtree vs MSE brute (par rapport a la VRAIE trajectoire) ---
double rawMSE = 0.0, filtMSE = 0.0;
for (int t = 0; t < T; t++) {
rawMSE += Math.Pow(obs[t] - trueState[t], 2);
filtMSE += Math.Pow(filtMean[t] - trueState[t], 2);
}
rawMSE /= T;
filtMSE /= T;
Console.WriteLine("\n=== Performance du filtre (par rapport a la trajectoire VRAIE) ===");
Console.WriteLine($"MSE observations brutes : {rawMSE:F3} (variance de capteur R = {F(obsVar)})");
Console.WriteLine($"MSE estimation filtree : {filtMSE:F3}");
Console.WriteLine($"Reduction d'erreur : {F((1.0 - filtMSE / rawMSE) * 100)}% (filtre vs brut)");
Console.WriteLine($"Variance post. finale : {filtVar[T - 1]:F3} (vs R = {F(obsVar)})");
=== Performance du filtre (par rapport a la trajectoire VRAIE) ===
MSE observations brutes : 4,875 (variance de capteur R = 4)
MSE estimation filtree : 1,283
Reduction d'erreur : 73.675% (filtre vs brut)
Variance post. finale : 1,186 (vs R = 4)Le filtre réduit fortement l’erreur : la MSE filtrée est nettement inférieure à la MSE des observations brutes. Deux mécanismes sont visibles dans le graphe ASCII :
#) varie beaucoup moins que les observations (.) d’un pas à l’autre : le filtre fait confiance à la dynamique (\(Q\) petit) pour rejeter le bruit du capteur (\(R\) grand).C’est précisément la situation où le filtre de Kalman aide : capteur bruité, dynamique fiable. L’exercice 1 explore la situation limite opposée (dynamique très incertaine).
Le filtre de Kalman tire sa puissance d’une propriété exacte : tout étant gaussien et linéaire, chaque postérieur est gaussien et calculable en forme fermée. Infer.NET (EP) retrouve cette solution exacte — c’est le cas où l’inférence est exacte, pas approchée.
La règle de combinaison est elle-même fermée : le postérieur \(p(x_t \mid y_t)\) est gaussien, avec une moyenne qui est un pondéré entre la prédiction (où le mobile devrait être d’après la dynamique) et l’observation (où le capteur le voit). Le poids relatif est l’inverse des variances : on fait davantage confiance à la source la plus certaine.
La borne de variance (équation de Riccati) garantit que l’incertitude résiduelle plafonne plutôt que d’exploser — c’est ce qui rend le filtre fiable sur le long terme.
Borne de variance et stabilité. L’équation de Riccati discrète qui régit la covariance d’erreur est analysée en profondeur dans Jazwinski, A. H. (1970), « Stochastic Processes and Filtering Theory », Academic Press. Le chapitre 7 montre sous quelles conditions (observabilité + commandabilité du système linéaire) le filtre de Kalman est stable : la covariance d’erreur converge vers un steady-state unique, et l’incertitude résiduelle reste bornée sur horizon infini. C’est ce résultat qui garantit qu’on peut déployer le filtre sur des trajectoire très longues (centaines de pas) sans dérive catastrophique — chaque pas réduit l’incertitude au-delà de ce que la dynamique seule pourrait accumuler. Pour Infer.NET, cette stabilité se traduit par un
Marginalgaussien dont la variance plafonne asymptotiquement, propriété vérifiable empiriquement sur le benchmark de la section 3.
Rejouez le filtre en balayant \(Q \in \{0{,}01\,;\, 0{,}5\,;\, 2{,}0\,;\, 10{,}0\}\) (sans changer la vraie trajectoire ni \(R\)). Pour chaque \(Q\), relevez la MSE filtrée et la variance postérieure moyenne. Observez le compromis : \(Q\) petit → lissage agressif (le filtre ignore le capteur et suit la pente) ; \(Q\) grand → le filtre colle au capteur (réagit au bruit). Quel \(Q\) minimise la MSE sur cette trajectoire ? - Indice : le \(Q\) optimal est proche de la vraie variance de process (ici \(0{,}5\)). Choisir \(Q\), c’est donc faire de l’estimation de paramètre — le filtre suppose connu ce que vous fixez à la main.
Le filtre ci-dessus suit une position scalaire avec une vitesse (drift) injectée manuellement. Généralisez à un état vectoriel \(x = (\text{position}, \text{vitesse})\) où \(\text{pos}_t = \text{pos}_{t-1} + \text{vit}_{t-1} \cdot \Delta t + \text{bruit}\) et \(\text{vit}_t = \text{vit}_{t-1} + \text{bruit}\). Utilisez Variable.Vector et Variable.MatrixTimesVector d’Infer.NET. - Indice : c’est le « vrai » filtre de Kalman matriciel. La vitesse devient latente et estimée (plus injectée). On observe uniquement la position (\(y = \text{pos} + \text{bruit}\)) et on infère la vitesse cachée — c’est tout l’intérêt.
Le filtre ci-dessus est séquentiel (un pas à la fois : il n’utilise que les observations passées). Construisez la version jointe : un seul VariableArray<double> x sur un Range temporel, avec x[t] = Variable.GaussianFromMeanAndVariance(x[t - 1] + drift, processVar) dans un bloc Variable.ForEach (auto-référence x[t-1]), toutes les observations branchées, et laissez EP inférer toutes les marginales simultanément — c’est le lissage (utilise aussi les observations futures). - Indice : le lissage rétro-propage l’information, donc MSE lissée \(\leq\) MSE filtrée. La syntaxe d’auto-référence x[t-1] dans Variable.ForEach est délicate (cas t == 0 à isoler avec Variable.If) — référez-vous aux exemples de chaînes gaussiennes Infer.NET. > Référence canonique. Le lissage joint sur la totalité de la séquence (vs filtrage séquentiel) est l’opération duale du filtre de Kalman décrite par Rauch, H. E., Tung, F. & Striebel, C. T. (1965), « Maximum likelihood estimates of linear dynamic systems », AIAA Journal 3(8):1445-1450. Leur algorithme RTS (Rauch-Tung-Striebel) opère en deux passes : un forward pass (filtre de Kalman standard, section 2 de ce notebook) suivi d’un backward pass qui propage l’information des observations futures en sens inverse. Le résultat : une estimation lisse \(x_t^{\text{smooth}}\) qui utilise toute la séquence \([y_1, \ldots, y_T]\) au pas \(t\), pas seulement \([y_1, \ldots, y_t]\). Ici la version jointe (un seul VariableArray<double> x + Variable.ForEach + EP) implémente exactement ce lissage en un seul appel EP global — MSE lissée \(\leq\) MSE filtrée, gain substantiel quand le SNR est élevé. Pour les détails de l’algorithme RTS forward-backward dans le formalisme EP, voir Minka, T. P. (2001b), « Expectation Propagation for Bayesian Filtering », section 4 (smoothing).
// Exercice 3 a completer
// Conseil : Range t = new Range(T); VariableArray<double> x = Variable.Array<double>(t);
// Dans Variable.ForEach(t) : If(t==0) x[t]=prior; If(t>0) x[t]=Gaussian(x[t-1]+drift, Q).
// Comparez la MSE lisse (jointe) a la MSE filtree (sequentielle ci-dessus).
Console.WriteLine("Exercice 3 a completer");Exercice 3 a completerLe filtre ci-dessus délègue l’étape de mise à jour à Infer.NET (engine.Infer(x) via EP). Mais pour un modèle scalaire gaussien, la récurrence de Kalman possède une forme fermée exacte (conjugaison). Implémentez-la vous-même : le moteur EP d’Infer.NET et votre formule doivent coïncider à la précision machine.
Rappel de la récurrence scalaire (l’étape prédire est déjà calculée à la main dans le worked example : predMean, predVar) :
obsVar)Objectif : compléter KalmanUpdateScalar(predMean, predVar, observation, obsVar) qui retourne (updMean, updVar), puis une boucle qui rejoue la récurrence complète sur obs (prédire avec drift et processVar comme dans le worked example) et compare chaque pas au tableau filtMean / filtVar produit par Infer.NET.
updVar ne dépend pas de l’observation (seulement de predVar et obsVar) — c’est pourquoi le filtre « apprend » sa propre incertitude au fil du temps indépendamment des mesures.updMean[t] doit correspondre à filtMean[t] (sortie Infer.NET EP) à mieux que \(10^{-9}\) près. Si c’est le cas, vous venez de réimplémenter, à la main, le cas scalaire gaussien que EP résout par messages.// Exercice 4 a completer : mise a jour de Kalman scalaire (forme fermee gaussienne)
// TODO etudiant : implementer KalmanUpdateScalar, puis boucler sur obs et comparer a Infer.NET.
(double updMean, double updVar) KalmanUpdateScalar(double predMean, double predVar,
double observation, double obsVar)
{
// Etape 1 : gain de Kalman K = predVar / (predVar + obsVar)
// Etape 2 : moyenne a posteriori updMean = predMean + K * (observation - predMean)
// Etape 3 : variance a posteriori updVar = (1 - K) * predVar
return (0.0, 0.0); // TODO etudiant : implementer
}
Console.WriteLine("Exercice 4 a completer");Exercice 4 a completer| Direction | Lien | Relation |
|---|---|---|
| ↔︎ Infer-14 | Infer-14 — Séquences (HMM) | HMM à état discret ↔︎ filtre de Kalman à état continu : même squelette markovien (état caché + observation), nature de l’état différente |
| ← Infer-2 | Infer-2 — Gaussian Mixtures | Le Kalman repose sur la conjugaison gaussienne — Infer-2 en pose les bases |
| ↔︎ Infer-16 | Infer-16 — Sparse GP | GP et Kalman modélisent tous deux de l’aléatoire structuré ; GP = espace fonctionnel non-paramétrique, Kalman = état fini paramétrique |
Référence fondatrice : Kalman, R. E. (1960) — A New Approach to Linear Filtering and Prediction Problems, ASME Journal of Basic Engineering 82(1). L’article fondateur du filtre de Kalman — l’algorithme d’estimation le plus répandu en pratique.
Sources fondatrices (papiers primaires).
Sources secondaires (écosystème Infer.NET).
Relations à la série.
Variable.Gaussian + paramètres supposés connus ; le jumeau PyMC utilise NUTS joint (Hoffman-Gelman 2014) pour estimer simultanément les hyperparamètres \(Q, R\), drift. Stack et algorithme distincts, substance cohérente.Pour aller plus loin.
Le filtre de Kalman est le pont entre le HMM discret d’Infer-14 et le monde continu. Sa puissance tient en trois propriétés : pour les systèmes linéaires gaussiens, l’inférence est exacte (conjugaison), récursive (un pas à la fois, en \(O(T)\)) et bornée (variance postérieure stable). Infer.NET la résout naturellement via EP, qui retrouve la solution exacte sur ce cas d’école.
La leçon générale : quand le modèle est linéaire-gaussien, nul besoin de méthodes d’inférence approchées lourdes — la conjugaison offre une solution fermée. C’est en sortant de ce régime (non-linéarités, distributions non gaussiennes) que les processus gaussiens (Infer-16) et les méthodes variationnelles deviennent nécessaires. Le Kalman n’est pas un cas particulier encombrant : c’est le point de référence par rapport auquel tous les filtres approchés se mesurent.