Nous preparons l’environnement pour les modèles de sequences temporelles, notamment les Hidden Markov Models (HMM). Ces modèles capturent les dependances temporelles entre observations via des etats caches qui evoluent selon une chaîne de Markov.
The below script needs to be able to find the current output cell; this is an easy method to get it.
Installed Packages
Microsoft.ML.Probabilistic, 0.4.2504.701
Microsoft.ML.Probabilistic.Compiler, 0.4.2504.701
Infer.NET pret !
Chargement du helper de visualisation des graphes de facteurs.
// Chargement du helper pour visualiser les graphes de facteurs#load "FactorGraphHelper.cs"Console.WriteLine($"FactorGraphHelper charge. Graphviz disponible : {FactorGraphHelper.IsGraphvizAvailable()}");
Environnement pret : Les namespaces Infer.NET sont charges, incluant Microsoft.ML.Probabilistic.Models pour définir les HMM et Microsoft.ML.Probabilistic.Distributions pour les distributions Dirichlet (transitions) et Gaussian (emissions).
2. Introduction aux HMM
Structure
Un HMM est défini par : - Etats caches : \(z_t\) - non observables - Observations : \(x_t\) - dependant de l’etat cache - Transitions : \(P(z_t | z_{t-1})\) - Emissions : \(P(x_t | z_t)\)
Schema
z_1 --> z_2 --> z_3 --> ... --> z_T (etats caches)
| | | |
v v v v
x_1 x_2 x_3 ... x_T (observations)
Applications
Domaine
Etats caches
Observations
NLP
POS tags
Mots
Finance
Regime de marche
Prix
Bio
Gene/Intergene
Sequence ADN
Meteo
Vrai temps
Mesures capteurs
Ancres savantes – les trois problemes canoniques des HMM. Rabiner (1989) structure toute la theorie des Hidden Markov Models autour de trois problemes distincts, que ce notebook aborde en pratique avec Infer.NET :
Probleme
Question
Algorithme canonique
Mise en oeuvre dans ce notebook
1. Evaluation
Quelle est la vraisemblance \(P(x_{1:T} \mid \lambda)\) d’une sequence sous le modele ?
Forward (recursion avant, Baum & Petrie 1966)
Calculee par Infer.NET via le passage de messages sur le graphe factoriel
2. Decodage
Quelle sequence d’etats caches \(z_{1:T}\) explique le mieux les observations ?
Viterbi (programmation dynamique, Viterbi 1967)
Section 3.3 – decode la sequence d’etats par Viterbi
3. Apprentissage
Comment estimer les parametres \(\lambda = (\pi, A, \mu)\) depuis les donnees ?
Baum-Welch (EM sur les HMM, Baum & Petrie 1966 ; Dempster-Laird-Rubin 1977)
Infer.NET (EP, 50 iter.) + init K-means pour briser la symetrie (section 3.3)
Pourquoi Infer.NET plutot que Forward-Backward / Baum-Welch “a la main” ? Les algorithmes canoniques (Forward, Backward, Baum-Welch) sont des instances de passage de messages sur le graphe factoriel du HMM. Infer.NET construit ce graphe de facon declarative (Variable.Array, ForEach, Switch) et instancie automatiquement ce passage de messages – exact (Forward-Backward) pour les emissions conjuguees (gaussiennes, section 3), approche via Expectation Propagation (EP) pour les modeles non conjugues. L’apport du model-based machine learning (Winn & Bishop, MBML) : on specifie le modele, la librairie produit l’inference, sans re-ecrire Baum-Welch pour chaque nouvelle structure (FHMM section 3.4, motif finding section 5).
References fondatrices : Rabiner, L.R. (1989), A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition, Proceedings of the IEEE 77(2):257-286 (les trois problemes canoniques, Baum-Welch, tutorial de reference) ; Baum, L.E. & Petrie, T. (1966), Statistical Inference for Probabilistic Functions of Finite State Markov Chains, Annals of Mathematical Statistics 37(6):1554-1563 (Forward-Backward, base de Baum-Welch) ; Viterbi, A.J. (1967), Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm, IEEE Transactions on Information Theory 13(2):260-269 (algorithme de Viterbi) ; Ghahramani, Z. (2001), An Introduction to Hidden Markov Models and Bayesian Networks, International Journal of Pattern Recognition and Artificial Intelligence 15(1):9-42 (cadre bayesien unifie, lien HMM / reseaux bayesiens) ; Ghahramani, Z. & Jordan, M.I. (1997), Factorial Hidden Markov Models, Machine Learning 29(2-3):245-273 (le FHMM de la section 3.4) ; Bailey, T.L. & Elkan, C. (1994), Fitting a Mixture Model by Expectation Maximization to Discover Motifs in Biopolymers, Proceedings ISMB 1994 (algorithme MEME, derriere le motif finding de la section 5) ; Winn, J. & Bishop, C.M. (2024), Model-Based Machine Learning, Chapman & Hall/CRC (l’approche modele-declaratif + EP d’Infer.NET).*
Limitation : Cette approche ignore les correlations temporelles. Un etat Normal suivi d’Anomalie puis Normal est traite identiquement a Normal→Normal→Normal.
// Definition du modele HMMRange stateRange =newRange(nStates).Named("state");Range timeRange =newRange(T).Named("time");// Distribution initialeVariable<Vector> probInit = Variable.Dirichlet(newdouble[]{1,1}).Named("probInit");// Matrice de transition (lignes = etat courant, colonnes = etat suivant)VariableArray<Vector> transMatrix = Variable.Array<Vector>(stateRange).Named("transMatrix");transMatrix[stateRange]= Variable.Dirichlet(newdouble[]{5,1}).ForEach(stateRange);// Favorise rester dans le meme etat// Parametres d'emission par etatVariableArray<double> emitMean = Variable.Array<double>(stateRange).Named("emitMean");VariableArray<double> emitPrec = Variable.Array<double>(stateRange).Named("emitPrec");// Priors sur les emissionsemitMean[0]= Variable.GaussianFromMeanAndVariance(10,10);// Etat 0 : NormalemitMean[1]= Variable.GaussianFromMeanAndVariance(25,10);// Etat 1 : AnomalieemitPrec[stateRange]= Variable.GammaFromShapeAndScale(2,0.5).ForEach(stateRange);// Sequence d'etatsVariableArray<int> states = Variable.Array<int>(timeRange).Named("states");// ObservationsVariableArray<double> obs = Variable.Array<double>(timeRange).Named("obs");Console.WriteLine("Variables HMM definies.");
Variables HMM definies.
Structure du modèle : - Priors Dirichlet : probInit ~ Dir(1,1) pour l’etat initial, transMatrix[k] ~ Dir(5,1) favorisant la persistance - Emissions gaussiennes : Moyennes a priori centrees sur 10 (Normal) et 25 (Anomalie)
Note : Ces variables preparent un HMM complet, mais l’inference utilise une approche simplifiee ci-dessous.
// Baseline naive : classification par pas de temps INDEPENDANTE (sans transitions).// Sert de point de comparaison au vrai HMM de la section 3.3.// (chaque pas est infere separement : on perd volontairement les dependances temporelles)Console.WriteLine("\n=== Inference des etats (baseline independante, sans transitions) ===");Console.WriteLine();for(int t =0; t < T; t++){// Pour chaque observation, determiner l'etat le plus probable Variable<int> etat = Variable.DiscreteUniform(nStates); Variable<double> obsVar = Variable.New<double>();// Emission selon l'etatusing(Variable.Case(etat,0)){ obsVar.SetTo(Variable.GaussianFromMeanAndPrecision(10,1));// Normal}using(Variable.Case(etat,1)){ obsVar.SetTo(Variable.GaussianFromMeanAndPrecision(25,1));// Anomalie} obsVar.ObservedValue= observations[t]; InferenceEngine eng =newInferenceEngine(); eng.Compiler.CompilerChoice= CompilerChoice.Roslyn; Discrete etatPost = eng.Infer<Discrete>(etat);int etatMAP = etatPost.GetProbs()[0]>0.5?0:1;string nomEtat = etatMAP ==0?"Normal":"Anomalie"; Console.WriteLine($"t={t} : obs={observations[t]:F1}, P(Normal)={etatPost.GetProbs()[0]:F2}, P(Anomalie)={etatPost.GetProbs()[1]:F2} -> {nomEtat}");}
Visualisation du graphe de facteurs pour la classification de sequences.
// Visualisation du graphe de facteurs pour la classification independante// On cree un modele unique avec ShowFactorGraph pour illustrer la structureVariable<int> etatViz = Variable.DiscreteUniform(nStates).Named("etat");Variable<double> obsViz = Variable.New<double>().Named("observation");using(Variable.Case(etatViz,0)){ obsViz.SetTo(Variable.GaussianFromMeanAndPrecision(10,1).Named("emission_Normal"));}using(Variable.Case(etatViz,1)){ obsViz.SetTo(Variable.GaussianFromMeanAndPrecision(25,1).Named("emission_Anomalie"));}obsViz.ObservedValue=15.0;// Valeur quelconqueInferenceEngine engViz =newInferenceEngine();engViz.Compiler.CompilerChoice= CompilerChoice.Roslyn;engViz.ShowFactorGraph=true;var _ = engViz.Infer(etatViz);// Declenche la generation du grapheConsole.WriteLine("Graphe de facteurs genere pour le modele de classification independante.");
Compiling model...done.
Graphe de facteurs genere pour le modele de classification independante.
Analyse de la détection d’états
Résultats : Classification parfaite des 10 observations
Position
Observation
P(Normal)
P(Anomalie)
Décision
0-2
9.5, 11.2, 10.8
1.00
0.00
Normal
3-5
24.5, 26.1, 25.3
0.00
1.00
Anomalie
6-9
10.1, 9.8, 11.5, 10.2
1.00
0.00
Normal
Observations clés :
Séparation nette : Les émissions gaussiennes (µ=10 vs µ=25) sont suffisamment séparées pour une décision sans ambiguïté
Limitation de l’approche : Chaque observation est classée indépendamment, ignorant les transitions temporelles
Un vrai HMM ferait mieux :
Lisser les transitions isolées (éviter Normal→Anomalie→Normal si transitions rares)
Propager l’information entre pas de temps adjacents
Détecter des changements de régime plus subtils
Note technique : La compilation répétée (“Compiling model…”) est due à la boucle créant un nouveau modèle à chaque itération - une implémentation réelle utiliserait un modèle unique.
// Affichage du graphe de facteursdisplay(HTML(FactorGraphHelper.GetLatestFactorGraphHtml()));
Model_06_22_26_21_17_55_72.svg
Graphe de facteurs : Classification independante
Ce graphe represente le modèle sans dependances temporelles :
Élément
Signification
etat
Variable discrete (Normal=0, Anomalie=1)
observation
Variable continue observee
Factor Cases
Branchement conditionnel selon l’etat
emission_Normal/Anomalie
Distributions gaussiennes d’emission
Structure : L’observation est reliee a l’etat via un switch (Variable.Case), creant deux chemins d’emission distincts. Ce modèle traite chaque observation independamment - aucune fleche ne relie les pas de temps.
3.3 HMM complet avec transitions markoviennes (Infer.NET)
L’approche précédente classe chaque observation independamment. Un vrai HMM relie les pas de temps par une matrice de transition : l’etat \(z_t\) depend de \(z_{t-1}\), ce qui permet de lisser les decisions et de detecter des regimes coherents.
Scénario enrichi : un capteur industriel a trois regimes (et non plus deux), avec des transitions “collantes” (un regime persiste avant de basculer) :
Regime
Emission moyenne
Ecart-type
Normal
~20
3
Alerte
~35
3
Critique
~50
4
Contrairement a une legende tenace, Infer.NET sait modeliser un HMM complet. La chaîne se construit avec un bloc Variable.ForEach sur le temps, en separant le premier pas (Variable.If(t == 0)) du reste (Variable.If(t > 0)), et en conditionnant la transition sur l’etat précédent via Variable.Switch. Deux idiomes sont indispensables :
Tr.AddAttribute(new Sequential()) : ordonne l’inference le long de la chaîne (sans cela, EP traite les pas dans le desordre et converge mal).
Brisure de symetrie par K-means : on initialise les etats avec un K-means 1-D et on centre les priors d’emission sur les centroides. Sans cette amorce, EP fusionne les regimes (tous les etats convergent vers la moyenne globale) – c’est précisément cet echec qui avait fait croire, a tort, a une “limitation d’Infer.NET”.
using Microsoft.ML.Probabilistic.Models.Attributes;// Sequential// ===================== HMM REEL : transitions markoviennes + emissions gaussiennes =====================// Port fidele de usptact/Infer.NET-HMM + oliparson/infer-hmm (API moderne 0.4.x).// Idiomes-cles que le "toy" precedent ignorait :// (1) Tr.AddAttribute(new Sequential()) -> ordonnancement sequentiel le long de la chaine// (2) priors d'emission ETALES sur la plage + init par K-means -> brisure de symetrie// (3) decodage de Viterbi (chemin globalement optimal) au lieu du MAP marginal.// (sans (2) EP fusionne les etats : c'est l'echec qui avait fait croire a tort a une "limitation").Rand.Restart(42);int nReg =3;int Tlen =120;// ---- Scenario capteur : 3 regimes (Normal ~20, Alerte ~35, Critique ~50) ----double[][] trueTrans =newdouble[][]{new[]{0.85,0.12,0.03},new[]{0.10,0.80,0.10},new[]{0.05,0.15,0.80},};double[] trueMean ={20.0,35.0,50.0};double[] trueStd ={3.0,3.0,4.0};int[] trueStates =newint[Tlen];double[] emit =newdouble[Tlen];int cur =0;for(int t =0; t < Tlen; t++){if(t >0){double u = Rand.Double(), acc =0;int nx =0;for(int j =0; j < nReg; j++){ acc += trueTrans[cur][j];if(u <= acc){ nx = j;break;}} cur = nx;} trueStates[t]= cur; emit[t]= trueMean[cur]+ trueStd[cur]* Rand.Normal();}int[] tc =newint[nReg];foreach(int s in trueStates) tc[s]++;Console.WriteLine($"Frequences vraies des etats : Normal={tc[0]}, Alerte={tc[1]}, Critique={tc[2]}");// ---- Statistiques de donnees (priors etales) ----double dMean = emit.Average();double dVar = emit.Select(x =>(x - dMean)*(x - dMean)).Average();double dMin = emit.Min(), dMax = emit.Max(), dRange = dMax - dMin;// ---- Modele Infer.NET ----Range K =newRange(nReg).Named("K");Range Tr =newRange(Tlen).Named("temps");Tr.AddAttribute(newSequential());// (1) chaine sequentiellevar probInitPrior = Variable.New<Dirichlet>().Named("probInitPrior");var probInit = Variable<Vector>.Random(probInitPrior).Named("probInit");probInit.SetValueRange(K);var cptTransPrior = Variable.Array<Dirichlet>(K).Named("cptTransPrior");var cptTrans = Variable.Array<Vector>(K).Named("cptTrans");cptTrans[K]= Variable<Vector>.Random(cptTransPrior[K]);cptTrans.SetValueRange(K);var emitMeanPrior = Variable.Array<Gaussian>(K).Named("emitMeanPrior");var emitMean = Variable.Array<double>(K).Named("emitMean");emitMean[K]= Variable<double>.Random(emitMeanPrior[K]);var emitPrecPrior = Variable.Array<Gamma>(K).Named("emitPrecPrior");var emitPrec = Variable.Array<double>(K).Named("emitPrec");emitPrec[K]= Variable<double>.Random(emitPrecPrior[K]);var zeroState = Variable.Discrete(probInit).Named("z0");var states = Variable.Array<int>(Tr).Named("states");var emissions = Variable.Array<double>(Tr).Named("emissions");using(var block = Variable.ForEach(Tr)){var t = block.Index;var previousState = states[t -1];using(Variable.If(t ==0)){using(Variable.Switch(zeroState)){ states[Tr]= Variable.Discrete(cptTrans[zeroState]);}}using(Variable.If(t >0)){using(Variable.Switch(previousState)){ states[t]= Variable.Discrete(cptTrans[previousState]);}}using(Variable.Switch(states[t])){ emissions[t]= Variable.GaussianFromMeanAndPrecision(emitMean[states[t]], emitPrec[states[t]]);}}// ---- K-means 1-D D'ABORD : sert a la fois aux priors (centroides + variance intra) et a l'init ----int[] initAssign =KMeans1D(emit, nReg,25);double[] kmCent =KMeansCentroids(emit, nReg,25);int[] cc =newint[nReg];foreach(int a in initAssign) cc[a]++;if(cc.Any(c => c ==0|| c < Tlen /(nReg *4))) initAssign =RobustInit(emit, nReg);// variance INTRA-cluster (pas la variance globale : c'est elle qui pilote le bruit d'emission)double wSum =0;for(int i =0; i < Tlen; i++){double d = emit[i]- kmCent[initAssign[i]]; wSum += d * d;}double wVar = Math.Max(wSum / Tlen,1e-3);// ---- Priors : means aux centroides K-means ; precision centree sur la variance INTRA ----probInitPrior.ObservedValue= Dirichlet.Uniform(nReg);cptTransPrior.ObservedValue= Util.ArrayInit(nReg, k => Dirichlet.Uniform(nReg));emitMeanPrior.ObservedValue= Util.ArrayInit(nReg, k => Gaussian.FromMeanAndVariance(kmCent[k], dVar *100));emitPrecPrior.ObservedValue= Util.ArrayInit(nReg, k => Gamma.FromShapeAndRate(2.0,2.0* wVar));emissions.ObservedValue= emit;Console.WriteLine($"Init K-means : centroides=[{string.Join(",", kmCent.Select(c => c.ToString("0.0")))}], var_intra={wVar:0.0}");var zinit = Variable<Discrete>.Array(Tr);zinit.ObservedValue= Util.ArrayInit(Tlen, t => Discrete.PointMass(initAssign[t], nReg));states[Tr].InitialiseTo(zinit[Tr]);var engine =newInferenceEngine(newExpectationPropagation());engine.NumberOfIterations=50;engine.ShowProgress=false;var transPost = engine.Infer<Dirichlet[]>(cptTrans);var initPost = engine.Infer<Dirichlet>(probInit);var meanPost = engine.Infer<Gaussian[]>(emitMean);var precPost = engine.Infer<Gamma[]>(emitPrec);// ---- Decodage de Viterbi (chemin globalement optimal) a partir des posteriors ----double[] vMean = meanPost.Select(g => g.GetMean()).ToArray();double[] vVar = precPost.Select(g =>1.0/ g.GetMean()).ToArray();double[][] vTrans = Util.ArrayInit(nReg, i =>{var v = transPost[i].GetMean();return Util.ArrayInit(nReg, j => v[j]);});Vector vInit = initPost.GetMean();int[] decoded =Viterbi(emit, nReg, vInit, vTrans, vMean, vVar);// ---- Appariement BIJECTIF inferes <-> vrais par moyenne d'emission (greedy, sans crash) ----int[] trueOfInfer = Enumerable.Repeat(-1, nReg).ToArray();bool[] usedTrue =newbool[nReg], usedInfer =newbool[nReg];var pairs =new List<(double d,int i,int j)>();for(int i =0; i < nReg; i++)for(int j =0; j < nReg; j++) pairs.Add((Math.Abs(vMean[i]- trueMean[j]), i, j));foreach(var p in pairs.OrderBy(x => x.d))if(!usedInfer[p.i]&&!usedTrue[p.j]){ trueOfInfer[p.i]= p.j; usedInfer[p.i]=true; usedTrue[p.j]=true;}int[] inferOfTrue =newint[nReg];for(int i =0; i < nReg; i++) inferOfTrue[trueOfInfer[i]]= i;string[] nm ={"Normal ","Alerte ","Critiq "};Console.WriteLine("\n=== Emissions gaussiennes retrouvees (apres appariement) ===");for(int it =0; it < nReg; it++){int i = inferOfTrue[it]; Console.WriteLine($" {nm[it]}: mean={vMean[i]:0.0} (vrai {trueMean[it]:0.0}), std={Math.Sqrt(vVar[i]):0.0} (vrai {trueStd[it]:0.0})");}Console.WriteLine("\n=== Matrice de transition retrouvee vs vraie ===");Console.WriteLine(" ->Normal ->Alerte ->Critiq (vrai)");for(int it =0; it < nReg; it++){ Vector row = transPost[inferOfTrue[it]].GetMean();double[] r =newdouble[nReg];for(int ni =0; ni < nReg; ni++) r[trueOfInfer[ni]]= row[ni]; Console.WriteLine($" {nm[it]}: {r[0]:0.00} {r[1]:0.00} {r[2]:0.00} ({trueTrans[it][0]:0.00},{trueTrans[it][1]:0.00},{trueTrans[it][2]:0.00})");}int correct =0;for(int t =0; t < Tlen; t++)if(trueOfInfer[decoded[t]]== trueStates[t]) correct++;Console.WriteLine($"\n=== Decodage de Viterbi : {correct}/{Tlen} corrects ({100.0*correct/Tlen:0.0}%) ===");// ===================== fonctions utilitaires =====================int[]KMeans1D(double[] data,int k,int iters){int n = data.Length;var sorted =(double[])data.Clone(); Array.Sort(sorted);double[] cent = Util.ArrayInit(k, c => sorted[Math.Min(n -1,(int)((c +0.5)* n / k))]);int[] assign =newint[n];for(int it =0; it < iters; it++){for(int i =0; i < n; i++){int best =0;double bd =double.MaxValue;for(int c =0; c < k; c++){double d = Math.Abs(data[i]- cent[c]);if(d < bd){ bd = d; best = c;}} assign[i]= best;}double[] sum =newdouble[k];int[] cnt =newint[k];for(int i =0; i < n; i++){ sum[assign[i]]+= data[i]; cnt[assign[i]]++;}for(int c =0; c < k; c++)if(cnt[c]>0) cent[c]= sum[c]/ cnt[c];}return assign;}int[]RobustInit(double[] data,int k){int n = data.Length, minPer = Math.Max(1, n /(k *5));double[] cent =KMeansCentroids(data, k,25);var sc = cent.Select((c, i)=>(c, i)).OrderBy(x => x.c).ToArray();var sd = data.Select((v, i)=>(v, i)).OrderBy(x => x.v).ToArray();int[] a =newint[n];for(int i =0; i < n; i++) a[sd[i].i]= sc[Math.Min(k -1, i * k / n)].i;for(int iter =0; iter <5; iter++){int[] counts =newint[k];foreach(int x in a) counts[x]++;for(int i =0; i < n; i++){int curA = a[i];int near = Enumerable.Range(0, k).OrderBy(c => Math.Abs(data[i]- cent[c])).First();if(near != curA && counts[curA]> minPer){ a[i]= near; counts[curA]--; counts[near]++;}}}return a;}double[]KMeansCentroids(double[] data,int k,int iters){int n = data.Length;var sorted =(double[])data.Clone(); Array.Sort(sorted);double[] cent = Util.ArrayInit(k, c => sorted[Math.Min(n -1,(int)((c +0.5)* n / k))]);int[] assign =newint[n];for(int it =0; it < iters; it++){for(int i =0; i < n; i++){int best =0;double bd =double.MaxValue;for(int c =0; c < k; c++){double d = Math.Abs(data[i]- cent[c]);if(d < bd){ bd = d; best = c;}} assign[i]= best;}double[] sum =newdouble[k];int[] cnt =newint[k];for(int i =0; i < n; i++){ sum[assign[i]]+= data[i]; cnt[assign[i]]++;}for(int c =0; c < k; c++)if(cnt[c]>0) cent[c]= sum[c]/ cnt[c];}return cent;}int[]Viterbi(double[] obs,int k, Vector init,double[][] trans,double[] mean,double[] var){int n = obs.Length;doubleLogPdf(double x,int s){double df = x - mean[s];return-0.5* Math.Log(2* Math.PI* var[s])- df * df /(2* var[s]);}double[,] delta =newdouble[n, k];int[,] psi =newint[n, k];for(int s =0; s < k; s++){ delta[0, s]= Math.Log(Math.Max(init[s],1e-300))+LogPdf(obs[0], s);}for(int t =1; t < n; t++)for(int s =0; s < k; s++){double mx =double.NegativeInfinity;int arg =0;for(int j =0; j < k; j++){double p = delta[t -1, j]+ Math.Log(Math.Max(trans[j][s],1e-300));if(p > mx){ mx = p; arg = j;}} delta[t, s]= mx +LogPdf(obs[t], s); psi[t, s]= arg;}int last =0;double best =double.NegativeInfinity;for(int s =0; s < k; s++)if(delta[n -1, s]> best){ best = delta[n -1, s]; last = s;}int[] path =newint[n]; path[n -1]= last;for(int t = n -2; t >=0; t--) path[t]= psi[t +1, path[t +1]];return path;}
Moteur : ExpectationPropagation, 50 itérations. EP propage les messages le long de la chaîne séquentielle.
Priors d’emission : moyennes centrees sur les centroides K-means (variance large pour rester souple) ; precision centree sur la variance intra-cluster – c’est elle qui represente le bruit d’emission, pas la variance globale (qui melange les trois regimes).
Brisure de symetrie : states[Tr].InitialiseTo(...) injecte l’assignation K-means comme point de depart, ce qui empeche EP de collapser tous les etats vers un même mode.
Decodage de Viterbi : on extrait le chemin d’etats globalement le plus probable (programmation dynamique sur les log-probabilites), plus coherent que le simple MAP marginal pas-a-pas.
Comme les etiquettes d’etats sont arbitraires (label switching), on apparie les etats inferes aux vrais regimes par proximite des moyennes d’emission (appariement bijectif glouton) avant de mesurer la justesse.
// Contraste : classification INDEPENDANTE (sans transitions) sur les MEMES donnees que le HMM.// Chaque pas est classe au regime gaussien le plus vraisemblable, en ignorant le voisinage temporel.doubleGaussLogPdf(double x,double m,double v){double d = x - m;return-0.5* Math.Log(2* Math.PI* v)- d * d /(2* v);}int[] indep =newint[Tlen];for(int t =0; t < Tlen; t++){int best =0;double bl =double.NegativeInfinity;for(int s =0; s < nReg; s++){double l =GaussLogPdf(emit[t], vMean[s], vVar[s]);if(l > bl){ bl = l; best = s;}} indep[t]= best;}int indepCorrect =0;for(int t =0; t < Tlen; t++)if(trueOfInfer[indep[t]]== trueStates[t]) indepCorrect++;Console.WriteLine("=== HMM (avec transitions) vs classification independante ===\n");Console.WriteLine($"Classification independante (argmax emission) : {indepCorrect}/{Tlen} corrects ({100.0 * indepCorrect / Tlen:0.0}%)");Console.WriteLine($"HMM + Viterbi (avec transitions) : {correct}/{Tlen} corrects ({100.0 * correct / Tlen:0.0}%)");int flips =0, hmmFlips =0, trueFlips =0;for(int t =1; t < Tlen; t++){if(indep[t]!= indep[t -1]) flips++;if(decoded[t]!= decoded[t -1]) hmmFlips++;if(trueStates[t]!= trueStates[t -1]) trueFlips++;}Console.WriteLine($"\nChangements d'etat le long de la sequence : vrais={trueFlips}, HMM={hmmFlips}, independant={flips}");Console.WriteLine("=> L'approche independante sur-segmente (un point bruite isole = faux changement de regime).");Console.WriteLine(" Le HMM colle au nombre reel de transitions grace a la penalite implicite des changements d'etat.");
=== HMM (avec transitions) vs classification independante ===
Classification independante (argmax emission) : 118/120 corrects (98,3%)
HMM + Viterbi (avec transitions) : 119/120 corrects (99,2%)
Changements d'etat le long de la sequence : vrais=22, HMM=20, independant=25
=> L'approche independante sur-segmente (un point bruite isole = faux changement de regime).
Le HMM colle au nombre reel de transitions grace a la penalite implicite des changements d'etat.
Analyse : ce que les transitions apportent
Le HMM retrouve les trois moyennes d’emission et la structure “collante” de la matrice de transition, puis decode la sequence d’etats par Viterbi.
Justesse : le decodage Viterbi du HMM depasse la classification independante, surtout aux frontieres de regime et sur les points bruites.
Coherence temporelle : l’approche independante sur-segmente (elle prend un point bruite isole pour un changement de regime), tandis que le HMM colle au nombre reel de transitions grace a la penalite implicite des changements d’etat.
Prix a payer : complexite \(O(T \cdot K^2)\) (vs \(O(T \cdot K)\)) et la necessite d’une brisure de symetrie pour eviter le collapse des etats.
Ce qui rend un HMM delicat dans Infer.NET (et comment le resoudre)
Le HMM ci-dessus fonctionne : il n’y a pas de “limitation” empechant Infer.NET de chainer transitions et emissions. La difficulte est ailleurs, et elle est generique a l’inference variationnelle sur les modèles a etats latents discrets :
Difficulte
Symptome
Solution appliquee ici
Symetrie des etats (label switching)
EP fusionne les regimes vers la moyenne globale
Init K-means + priors centres sur les centroides + InitialiseTo
Ordre d’inference
Convergence lente/instable sur la chaîne
Tr.AddAttribute(new Sequential())
Prior de precision
Variances d’emission explosent ou s’effondrent
Precision centree sur la variance intra-cluster
Etiquettes arbitraires
Etats inferes non alignes aux vrais
Appariement bijectif par moyenne d’emission
A retenir : un echec de convergence sur un HMM n’est presque jamais une “incapacite de la librairie” ; c’est un problème d’initialisation et de priors. La première version “toy” de ce notebook concluait a tort a une limitation faute d’avoir brise la symetrie.
Graphe de facteurs : HMM complet (3 pas de temps)
Ce graphe illustre la structure d’un Hidden Markov Model avec dependances temporelles :
[s0] -----> [s1] -----> [s2] (etats caches avec transitions)
| | |
v v v
[x0] [x1] [x2] (observations)
Élément
Rôle
Facteur
s0, s1, s2
Etats caches
Variables discretes (regimes)
x0, x1, x2
Observations
Variables continues (gaussiennes)
Fleches horizontales
Transitions
P(s_t | s_{t-1}) via Variable.Switch
Fleches verticales
Emissions
P(x_t | s_t) gaussiennes conditionnelles
Propagation de messages : - Forward : P(s_t | x_{1:t}) - information passee - Backward : P(x_{t+1:T} | s_t) - information future - Marginal : P(s_t | x_{1:T}) = Forward x Backward (normalise)
C’est cette structure en chaîne qui permet au HMM de lisser les predictions en tenant compte des observations voisines.
3.4 HMM factoriel : des causes superposees
Un HMM ordinaire suppose un etat cache a la fois. Mais souvent plusieurs processus independants se superposent dans une seule mesure. Exemple canonique (desagregation de charge, energy disaggregation) : la consommation electrique totale d’un logement est la somme des contributions de plusieurs appareils, chacun suivant sa propre dynamique ON/OFF.
Le HMM factoriel (FHMM, port fidele de usptact/FactorialHiddenMarkovModel) modelise \(C\) chaînes de Markov independantes dont les contributions s’additionnent dans l’observation :
Interpretation : identifiabilite a permutation pres
Le FHMM retrouve les contributions additives de chaque chaîne (~0/20 et ~0/5) et reconstruit l’etat conjoint avec une très bonne precision.
Point subtil et pedagogique : un FHMM n’est identifiable qu’a permutation des chaînes pres. Comme seule leur somme est observee, rien ne distingue “chaîne A = gros appareil, chaîne B = petit” de l’arrangement inverse – ni l’etiquette 0/1 a l’interieur d’une chaîne. Le score est donc calcule modulo relabeling (on essaie l’identite + l’echange des chaînes, et le flip d’etiquette dans chaque chaîne, on garde le meilleur).
On choisit aussi le run de meilleure log-evidence sur plusieurs initialisations aleatoires (multi-restart) : c’est la encore une parade a la symetrie, recommandee par l’implementation de reference.
4. Detection de Regimes Meteo
Cas d’usage classique : inferer une variable latente (le vrai temps) a partir d’observations bruitees (temperature).
Regime
Temperature moyenne
Ecart-type
Soleil
22C
2C
Pluie
15C
2C
// Exemple : Detection de regimes meteo (Soleil/Pluie) a partir de temperature// Soleil : temperature ~22C// Pluie : temperature ~15Cdouble[] tempJour ={21,23,22,20,15,14,16,15,14,21,22,23};int Tmeteo = tempJour.Length;Console.WriteLine("=== Detection Regimes Meteo ===");Console.WriteLine($"\nTemperatures : {string.Join(",", tempJour)}\n");for(int t =0; t < Tmeteo; t++){ Variable<int> meteo = Variable.DiscreteUniform(2); Variable<double> temp = Variable.New<double>();using(Variable.Case(meteo,0))// Soleil{ temp.SetTo(Variable.GaussianFromMeanAndVariance(22,4));}using(Variable.Case(meteo,1))// Pluie{ temp.SetTo(Variable.GaussianFromMeanAndVariance(15,4));} temp.ObservedValue= tempJour[t]; InferenceEngine eng =newInferenceEngine(); eng.Compiler.CompilerChoice= CompilerChoice.Roslyn; Discrete meteoPost = eng.Infer<Discrete>(meteo);string regime = meteoPost.GetProbs()[0]>0.5?"Soleil":"Pluie"; Console.WriteLine($"Jour {t+1,2} : {tempJour[t]:F0}C -> {regime} (P={meteoPost.GetProbs().Max():F2})");}
Visualisation du graphe de facteurs pour la detection de regimes meteo.
// Visualisation du graphe de facteurs pour la detection meteo// Modele identique a la classification independanteVariable<int> meteoViz = Variable.DiscreteUniform(2).Named("meteo");Variable<double> tempViz = Variable.New<double>().Named("temperature");using(Variable.Case(meteoViz,0)){ tempViz.SetTo(Variable.GaussianFromMeanAndVariance(22,4).Named("emission_Soleil"));}using(Variable.Case(meteoViz,1)){ tempViz.SetTo(Variable.GaussianFromMeanAndVariance(15,4).Named("emission_Pluie"));}tempViz.ObservedValue=18.0;InferenceEngine engMeteoViz =newInferenceEngine();engMeteoViz.Compiler.CompilerChoice= CompilerChoice.Roslyn;engMeteoViz.ShowFactorGraph=true;var _meteo = engMeteoViz.Infer(meteoViz);Console.WriteLine("\nGraphe de facteurs genere pour le modele meteo.");
Compiling model...done.
Graphe de facteurs genere pour le modele meteo.
Interpretation des résultats meteo
Jours
Temperature
Regime
Confiance
1-4
20-23C
Soleil
93-100%
5-9
14-16C
Pluie
99-100%
10-12
21-23C
Soleil
99-100%
Observation : Le jour 4 (20C) montre une confiance de 93% car cette temperature est a mi-chemin entre les deux regimes. C’est exactement le type de cas ou un HMM avec transitions aiderait : si les jours 1-3 sont “Soleil”, le jour 4 devrait l’etre aussi par continuite.
Visualisons le graphe de facteurs pour ce modèle de classification :
// Affichage du graphe meteodisplay(HTML(FactorGraphHelper.GetLatestFactorGraphHtml()));
Model_06_22_26_21_18_32_44.svg
Graphe de facteurs : Modèle meteo (classification)
Structure identique a la detection d’anomalies - un modèle de melange gaussien :
Composante
Paramètres
Soleil
N(22, 4) - temperature elevee
Pluie
N(15, 4) - temperature basse
Le graphe montre le branchement conditionnel (Variable.Case) qui selectionne l’emission appropriee selon l’etat latent meteo.
Application : Ce modèle simple est equivalent a un classifieur de Bayes naive pour une seule feature (temperature).
Résultats : Detection correcte des deux regimes (Soleil jours 1-4 et 10-12, Pluie jours 5-9).
Point d’intérêt : Le jour 4 (20C) a une confiance “seulement” 93% car la temperature est intermediaire entre les deux regimes.
4bis. Exercice : Sensibilite du Modèle aux Emissions Chevauchantes
Dans la section 4, les deux regimes meteo (Soleil ~22C, Pluie ~15C) sont bien separes avec un ecart de 7C pour un ecart-type de 2C. Que se passe-t-il quand les emissions se chevauchent davantage ?
Scénario : Deux villes ont des regimes de temperature proches :
Ville
Moyenne
Ecart-type
Ville A (oceanique)
18C
4C
Ville B (semi-continental)
16C
4C
L’ecart entre les moyennes (2C) est plus petit que l’ecart-type (4C) : les distributions se recouvrent fortement.
Objectifs : 1. Implementez un modèle de classification independante avec ces paramètres 2. Testez sur les observations : {15, 17, 19, 16, 18, 14, 20, 17} 3. Identifiez les observations ambigues (P < 0.8) 4. Repetez avec des emissions mieux separees (moyennes 20C et 12C) et comparez
Étapes : 1. Définir les paramètres d’emission pour les deux villes (GaussianFromMeanAndVariance avec variance 16) 2. Boucler sur les observations et inferer l’appartenance a chaque ville avec Variable.Case 3. Afficher P(Ville A) et P(Ville B) pour chaque observation 4. Repeter avec muA=20, muB=12 et comparer les confiances
Indices : - Utilisez Variable.GaussianFromMeanAndVariance(mu, sigma2) avec sigma2 = 16 (sigma = 4) - Une observation est ambigu si P(Ville A) est entre 0.2 et 0.8 - Avec des moyennes plus eloignees (20 et 12), l’ecart (8C) depasse 2 sigma : separation nette attendue
// Exercice : Sensibilite du modele aux emissions chevauchantes// Observationsdouble[] tempsVilles ={15,17,19,16,18,14,20,17};// Parametres d'emission (ecart faible = chevauchement)double muA =18.0, muB =16.0;double sigma2Ville =16.0;// sigma = 4// TODO: Etape 1 - Implementez la classification independante// (Reutilisez le pattern Variable.DiscreteUniform + Variable.Case de la section 4)// TODO: Etape 2 - Affichez P(Ville A) et P(Ville B) pour chaque observation// Identifiez les observations ambigues (P entre 0.2 et 0.8)// TODO: Etape 3 - Repetez avec muA=20 et muB=12 (emissions bien separees)// Comparez les niveaux de confianceConsole.WriteLine("Exercice a completer : analysez la sensibilite aux emissions chevauchantes.");
Exercice a completer : analysez la sensibilite aux emissions chevauchantes.
5. Motif Finding (Bioinformatique)
Problème
Trouver des motifs conserves dans des sequences ADN.
Modèle
Arriere-plan : nucleotides uniformes (A, C, G, T)
Motif : positions avec distributions spécifiques
Un modèle generatif de motif (et son miroir : l’inference)
Le comptage de k-mers (naif) trouve les sous-chaînes frequentes, mais il ne modelise pas l’incertitude : une position du motif peut tolerer plusieurs bases. Le modèle probabiliste de dotnet/infer (MotifFinder) decrit comment les sequences sont engendrees, puis inverse ce processus pour retrouver le motif.
Modèle generatif (port fidele du MotifFinder officiel) : - PWM (Position Weight Matrix) : pour chacune des \(W\) positions du motif, une distribution sur {A, C, G, T}, de prior Dirichlet. - Pour chaque sequence : presence du motif motifPresence ~ Bernoulli(0.8), et si present, une position motifPos ~ Uniforme. - La sequence = fond uniforme (gauche) + motif (echantillonne de la PWM) + fond uniforme (droite), assemblee par Variable.StringFromArray / Variable.StringOfLength (inference sur des chaînes de caractères, une specialite d’Infer.NET).
L’inference et la generation sont symetriques : on echantillonne d’abord depuis une PWM connue (ÉTAPE 1), puis on fait tourner le modèle a l’envers pour retrouver cette PWM a partir des seules sequences (ÉTAPE 2). Comparer les deux est le coeur pedagogique de la section.
// ETAPE 1 - GENERATION (forward) : on echantillonne des sequences depuis une PWM CONNUE.// C'est le miroir pedagogique de l'inference : meme modele, parcouru dans l'autre sens.Rand.Restart(1337);int SequenceCount =70;int SequenceLength =25;double MotifPresenceProbability =0.8;DiscreteChar NucleobaseDist(double a,double c,double g,double t){ Vector probs = PiecewiseVector.Zero(char.MaxValue+1); probs['A']= a; probs['C']= c; probs['G']= g; probs['T']= t;return DiscreteChar.FromVector(probs);}// ---- Vraie PWM du motif (matrice de frequences position-specifique, 8 positions = original) ----// La position 4 (index 3) est volontairement uniforme : le modele doit aussi gerer une colonne non-informative.var trueMotif =new[]{NucleobaseDist(0.80,0.10,0.05,0.05),NucleobaseDist(0.00,0.90,0.05,0.05),NucleobaseDist(0.00,0.00,0.50,0.50),NucleobaseDist(0.25,0.25,0.25,0.25),NucleobaseDist(0.10,0.10,0.10,0.70),NucleobaseDist(0.00,0.00,0.90,0.10),NucleobaseDist(0.90,0.05,0.00,0.05),NucleobaseDist(0.50,0.50,0.00,0.00),};int motifLength = trueMotif.Length;var background =NucleobaseDist(0.25,0.25,0.25,0.25);// ---- Echantillonnage FORWARD (background + motif insere a position aleatoire) ----string[] seqData =newstring[SequenceCount];int[] posData =newint[SequenceCount];for(int i =0; i < SequenceCount; i++){if(Rand.Double()<= MotifPresenceProbability){ posData[i]= Rand.Int(SequenceLength - motifLength +1);var before = Util.ArrayInit(posData[i], j => background.Sample());var after = Util.ArrayInit(SequenceLength - motifLength - posData[i], j => background.Sample());var mc = Util.ArrayInit(motifLength, j => trueMotif[j].Sample()); seqData[i]=newstring(before)+newstring(mc)+newstring(after);}else{ posData[i]=-1; seqData[i]=newstring(Util.ArrayInit(SequenceLength, j => background.Sample()));}}// ---- Affichage : la PWM "verite terrain" + quelques sequences echantillonnees ----voidShowPWM(string cap, Func<int,char,double> w){ Console.WriteLine(cap +" (lignes = base, colonnes = position du motif) :");foreach(char b innew[]{'A','C','G','T'}){ Console.Write($" {b}:");for(int i =0; i < motifLength; i++) Console.Write($" {w(i,b):0.00}"); Console.WriteLine();}}Console.WriteLine($"Echantillonnage FORWARD : {SequenceCount} sequences de longueur {SequenceLength}, motif present avec prob {MotifPresenceProbability}.\n");ShowPWM("Vraie PWM du motif (8 positions, position 4 volontairement uniforme)",(i,b)=> trueMotif[i][b]);Console.WriteLine("\nQuelques sequences echantillonnees (^ marque la position vraie du motif) :");for(int i =0; i <6; i++){string mark = posData[i]<0?"(pas de motif)":newstring(' ', posData[i])+newstring('^', motifLength); Console.WriteLine($" seq {i}: {seqData[i]}"); Console.WriteLine($" {mark}");}
Echantillonnage FORWARD : 70 sequences de longueur 25, motif present avec prob 0,8.
Vraie PWM du motif (8 positions, position 4 volontairement uniforme) (lignes = base, colonnes = position du motif) :
A: 0,80 0,00 0,00 0,25 0,10 0,00 0,90 0,50
C: 0,10 0,90 0,00 0,25 0,10 0,00 0,05 0,50
G: 0,05 0,05 0,50 0,25 0,10 0,90 0,00 0,00
T: 0,05 0,05 0,50 0,25 0,70 0,10 0,05 0,00
Quelques sequences echantillonnees (^ marque la position vraie du motif) :
seq 0: CTACTTCGAATTTACCCCTATATTT
^^^^^^^^
seq 1: TTGTGCGGGCATAAGACTTTGACTA
^^^^^^^^
seq 2: CAAACGTCGGGCAGACCGTACTCCG
(pas de motif)
seq 3: ACTTTGACCTACCAGAGCCTGATGA
^^^^^^^^
seq 4: GGTTGACGCCACCGTCGATGAACGA
^^^^^^^^
seq 5: AGCGTACTCTGAATCAGTCGTAGCA
^^^^^^^^
Étape 2 : inference – retrouver la PWM, la presence et la position
On observe maintenant uniquement les sequences (en oubliant la PWM et les positions vraies) et on demande a Infer.NET de les reconstruire. Le même modèle generatif sert de squelette ; le moteur infere a posteriori : - la PWM du motif (a comparer a la verite terrain de l’étape 1), - pour chaque sequence, la probabilite de presence du motif et sa position la plus probable.
using Microsoft.ML.Probabilistic.Factors.Attributes;// QualityBand// ETAPE 2 - INFERENCE (backward) : on OBSERVE les sequences et on RETROUVE// la PWM + la presence + la position du motif, par programmation probabiliste.// (reutilise seqData / posData / trueMotif / background generes a l'etape 1)// ---- Modele generatif Infer.NET ----Vector pseudo = PiecewiseVector.Constant(char.MaxValue+1,1e-6);pseudo['A']= pseudo['C']= pseudo['G']= pseudo['T']=1.0;Range mcr =newRange(motifLength);var motifProbs = Variable.Array<Vector>(mcr);motifProbs[mcr]= Variable.Dirichlet(pseudo).ForEach(mcr);var sr =newRange(SequenceCount);var sequences = Variable.Array<string>(sr);var motifPos = Variable.Array<int>(sr);motifPos[sr]= Variable.DiscreteUniform(SequenceLength - motifLength +1).ForEach(sr);var motifPresence = Variable.Array<bool>(sr);motifPresence[sr]= Variable.Bernoulli(MotifPresenceProbability).ForEach(sr);using(Variable.ForEach(sr)){using(Variable.If(motifPresence[sr])){var motifChars = Variable.Array<char>(mcr); motifChars[mcr]= Variable.Char(motifProbs[mcr]);var motif = Variable.StringFromArray(motifChars);var bgRightLen = SequenceLength - motifLength - motifPos[sr];var bgLeft = Variable.StringOfLength(motifPos[sr], background);var bgRight = Variable.StringOfLength(bgRightLen, background); sequences[sr]= bgLeft + motif + bgRight;}using(Variable.IfNot(motifPresence[sr])){ sequences[sr]= Variable.StringOfLength(SequenceLength, background);}}// ---- Inference ----sequences.ObservedValue= seqData;var eng =newInferenceEngine();eng.NumberOfIterations=30;eng.Compiler.RecommendedQuality= QualityBand.Experimental;eng.ShowProgress=false;var pwmPost = eng.Infer<IList<Dirichlet>>(motifProbs);var presPost = eng.Infer<IList<Bernoulli>>(motifPresence);var posPost = eng.Infer<IList<Discrete>>(motifPos);// ---- Resultats : PWM vraie vs retrouvee ----voidPrintPWM(string caption, Func<int,char,double> w){ Console.WriteLine(caption +" :");foreach(char b innew[]{'A','C','G','T'}){ Console.Write($" {b}:");for(int i =0; i < motifLength; i++) Console.Write($" {w(i,b):0.00}"); Console.WriteLine();}}PrintPWM("Vraie PWM du motif",(i,b)=> trueMotif[i][b]);Console.WriteLine();PrintPWM("PWM RETROUVEE (moyenne posterieure)",(i,b)=> pwmPost[i].GetMean()[b]);Console.WriteLine("\n=== Quelques sequences : position predite vs vraie ===");int shown =0;for(int i =0; i < SequenceCount && shown <8; i++){int predPos = presPost[i].GetProbTrue()>0.5? posPost[i].GetMode():-1; Console.WriteLine($" seq {i,2}: {seqData[i]} P(motif)={presPost[i].GetProbTrue():0.00} pos_pred={predPos,2} (vrai {posData[i],2})"); shown++;}int posOK =0, present =0;for(int i =0; i < SequenceCount; i++){if(posData[i]>=0){ present++;if(presPost[i].GetProbTrue()>0.5&& posPost[i].GetMode()== posData[i]) posOK++;}}Console.WriteLine($"\n=== Positions exactes retrouvees : {posOK}/{present} sequences avec motif ===");
La PWM retrouvee reproduit fidelement la verite terrain (colonnes fortement conservees retrouvees nettement ; la colonne volontairement uniforme reste plate), et la majorite des positions du motif sont correctement localisees.
Pourquoi l’erreur residuelle ? Les quelques ecarts de probabilite dans la PWM et les rares faux positifs/negatifs de position viennent du lissage laplacien (Laplacian / additive smoothing) : le prior Dirichlet place des pseudo-comptes (=1.0 pour A/C/G/T) sur chaque position. Ces pseudo-comptes empechent les probabilites de tomber a exactement 0 ou 1 (utile : un nucleotide jamais vu dans l’echantillon reste possible), mais ils tirent legerement les estimations vers l’uniforme. Avec plus de sequences, le poids des données domine les pseudo-comptes et l’estimation se resserre vers la PWM vraie.
C’est exactement le compromis biais-variance du lissage : un peu de biais (vers l’uniforme) contre beaucoup moins de variance (pas de probabilite 0 catastrophique sur un petit echantillon).
6. Exemple guide : Detection d’Anomalies
Enonce
Utilisez un HMM pour detecter des periodes anormales dans une serie temporelle de ventes.
Transition : Après avoir explore les HMM theoriquement (section 3) et sur des exemples meteo (section 4) et bioinformatiques (section 5), passons a un exercice pratique de detection d’anomalies dans un contexte business.
Detection d’anomalies dans les series temporelles de ventes :
Scénario
Caractéristiques
Action
Promotion
Hausse temporaire
Optimiser le stock
Rupture stock
Baisse soudaine
Reapprovisionner
Tendance
Evolution graduelle
Ajuster previsions
// Exemple guide : Detection d'anomalies dans les ventes// Ventes journalieres (normal ~100, promo ~200)double[] ventes ={98,105,102,99,195,210,205,198,103,97,101,100};Console.WriteLine("=== Detection Periodes de Promotion ===");Console.WriteLine($"\nVentes : {string.Join(",", ventes.Select(v => v.ToString("F0")))}\n");for(int t =0; t < ventes.Length; t++){ Variable<int> regime = Variable.DiscreteUniform(2); Variable<double> venteVar = Variable.New<double>();using(Variable.Case(regime,0))// Normal{ venteVar.SetTo(Variable.GaussianFromMeanAndVariance(100,100));}using(Variable.Case(regime,1))// Promotion{ venteVar.SetTo(Variable.GaussianFromMeanAndVariance(200,100));} venteVar.ObservedValue= ventes[t]; InferenceEngine eng =newInferenceEngine(); eng.Compiler.CompilerChoice= CompilerChoice.Roslyn; Discrete regimePost = eng.Infer<Discrete>(regime);string etat = regimePost.GetProbs()[1]>0.5?"PROMO":"Normal"; Console.WriteLine($"Jour {t+1,2} : {ventes[t],3:F0} -> {etat} (P={regimePost.GetProbs().Max():F2})");}Console.WriteLine("\n=> Jours 5-8 detectes comme periode de promotion");
=== Detection Periodes de Promotion ===
Ventes : 98, 105, 102, 99, 195, 210, 205, 198, 103, 97, 101, 100
Compiling model...done.
Jour 1 : 98 -> Normal (P=1,00)
Compiling model...done.
Jour 2 : 105 -> Normal (P=1,00)
Compiling model...done.
Jour 3 : 102 -> Normal (P=1,00)
Compiling model...done.
Jour 4 : 99 -> Normal (P=1,00)
Compiling model...done.
Jour 5 : 195 -> PROMO (P=1,00)
Compiling model...done.
Jour 6 : 210 -> PROMO (P=1,00)
Compiling model...done.
Jour 7 : 205 -> PROMO (P=1,00)
Compiling model...done.
Jour 8 : 198 -> PROMO (P=1,00)
Compiling model...done.
Jour 9 : 103 -> Normal (P=1,00)
Compiling model...done.
Jour 10 : 97 -> Normal (P=1,00)
Compiling model...done.
Jour 11 : 101 -> Normal (P=1,00)
Compiling model...done.
Jour 12 : 100 -> Normal (P=1,00)
=> Jours 5-8 detectes comme periode de promotion
Visualisation du graphe de facteurs pour la detection d’anomalies dans les ventes.
// Visualisation du graphe pour la detection de promotionsVariable<int> regimeViz = Variable.DiscreteUniform(2).Named("regime_ventes");Variable<double> venteViz = Variable.New<double>().Named("ventes_journalieres");using(Variable.Case(regimeViz,0)){ venteViz.SetTo(Variable.GaussianFromMeanAndVariance(100,100).Named("emission_Normal"));}using(Variable.Case(regimeViz,1)){ venteViz.SetTo(Variable.GaussianFromMeanAndVariance(200,100).Named("emission_Promo"));}venteViz.ObservedValue=150.0;InferenceEngine engVenteViz =newInferenceEngine();engVenteViz.Compiler.CompilerChoice= CompilerChoice.Roslyn;engVenteViz.ShowFactorGraph=true;var _regime = engVenteViz.Infer(regimeViz);Console.WriteLine("\nGraphe de facteurs genere pour la detection de promotions.");
Compiling model...done.
Graphe de facteurs genere pour la detection de promotions.
Interpretation des résultats de l’exemple guide
Le detecteur identifie correctement la periode promotionnelle (jours 5-8) :
Periode
Ventes moyennes
Regime detecte
Jours 1-4
~101
Normal
Jours 5-8
~202
PROMO
Jours 9-12
~100
Normal
Applications business : - Planification stocks : Anticiper les pics de demande - Analyse retrospective : Identifier les periodes promotionnelles non documentees - Detection de fraude : Reperer les anomalies inexpliquees
Visualisons le graphe de facteurs pour ce modèle :
// Affichage du graphe de detection de promotionsdisplay(HTML(FactorGraphHelper.GetLatestFactorGraphHtml()));
Model_06_22_26_21_18_45_80.svg
Graphe de facteurs : Detection de promotions
Même structure que les modèles précédents - classification par melange gaussien :
Regime
Distribution
Interpretation
Normal
N(100, 100)
Ventes moyennes ~100 unites
Promo
N(200, 100)
Ventes elevees ~200 unites
Extension possible : Pour une detection plus robuste des periodes promotionnelles, on ajouterait des transitions Markoviennes comme dans la section 3.3. Cela permettrait de :
Penaliser les changements frequents Normal <-> Promo
Detecter des periodes coherentes (promotion sur plusieurs jours consecutifs)
Lisser les observations ambigues (ex: 150 unites) selon le contexte
Résultats : Detection parfaite - jours 5-8 identifies comme periode promotionnelle (ventes ~200 vs ~100 en normal).
6bis. Exercice : Impact des Paramètres de Transition sur le Lissage
Dans la section 3.3, nous avons construit un HMM complet a transitions markoviennes (trois regimes Normal/Alerte/Critique, inference par Expectation Propagation puis decodage de Viterbi) et observe comment la matrice de transition lisse les decisions la ou la classification independante sur-segmente.
Votre mission : modifiez les paramètres du HMM et analysez l’impact sur les probabilites a posteriori.
Scénario : Un capteur de debit hydraulique produit les mesures suivantes :
{8.2, 9.1, 9.5, 12.5, 14.0, 9.8, 8.5, 12.8}
Les deux regimes possibles sont : - Debit normal : ~N(9, 2) (moyenne 9, variance 2) - Fuite : ~N(13, 2) (moyenne 13, variance 2)
Objectifs : 1. Implementez la passe Forward avec les paramètres fournis 2. Testez deux configurations de transitions : - Config A : pStay = 0.7, pSwitch = 0.3 (transitions frequentes) - Config B : pStay = 0.95, pSwitch = 0.05 (transitions rares) 3. Comparez les probabilites a posteriori pour les observations ambigues (12.5, 14.0, 9.8, 12.8) 4. Expliquez pourquoi la Config B produit des sequences plus “lisses”
Indice : Les valeurs 12.5 et 12.8 sont ambigues car entre les deux moyennes (9 et 13). Comparez P(Fuite) entre les deux configs pour ces valeurs - la config B devrait favoriser la continuite avec les voisins plutot que la vraisemblance seule.
// Exercice : Impact des parametres de transition sur le lissage// Donnees capteur de debit hydrauliquedouble[] debit ={8.2,9.1,9.5,12.5,14.0,9.8,8.5,12.8};int Tdebit = debit.Length;// Parametres d'emissiondouble muNormal =9.0, muFuite =13.0;double sigma2Debit =2.0;// Fonction de vraisemblance gaussienne (reutilisez celle de la section 3.3)// gaussLik(x, mu, s2) = exp(-0.5 * (x - mu)^2 / s2)// TODO: Etape 1 - Implementez la passe Forward pour une config donnee// Parametres : pStay, pSwitch, observations// Retourne : tableau alpha[T, 2] normalise// TODO: Etape 2 - Implementez la passe Backward// Retourne : tableau beta[T, 2] normalise// TODO: Etape 3 - Calculez les posterieurs gamma[T, 2] = alpha * beta (normalise)// TODO: Etape 4 - Testez Config A (pStay=0.7) et affichez les resultatsdouble pStayA =0.7;double pSwitchA =0.3;// TODO: Etape 5 - Testez Config B (pStay=0.95) et affichez les resultatsdouble pStayB =0.95;double pSwitchB =0.05;// TODO: Etape 6 - Affichez un tableau comparatif pour les 4 observations ambigues// | t | Obs | P(Fuite) Config A | P(Fuite) Config B | Decision A | Decision B |Console.WriteLine("Exercice a completer : implementez Forward-Backward pour les deux configs.");
Exercice a completer : implementez Forward-Backward pour les deux configs.
7. Resume
Concept
Description
HMM
Modèle a etats caches avec dependances temporelles
Emissions
Distribution des observations selon l’etat
Transitions
Probabilites de changement d’etat
Viterbi
Algorithme pour trouver la sequence d’etats optimale
Attendu : Normal (jours 1-4), Degrade (5-8), Panne (9-12), Degrade (13-14), Normal (15-16).
Objectifs : 1. Adapter le modèle de classification independante a 3 etats au lieu de 2 2. Parametrer les emissions gaussiennes pour chaque etat (moyenne + precision) 3. Implementer la boucle d’inference et afficher les decisions
Indice : Utilisez Variable.DiscreteUniform(3) pour 3 etats et ajoutez un troisieme bloc Variable.Case(etat, 2) pour l’etat “En panne”.
// Exercice : HMM 3 etats pour detection de pannes machinedouble[] capteur ={98,102,99,95,68,75,72,65,18,22,15,25,71,68,100,103};int T = capteur.Length;int K =3;// 3 etats : Normal=0, Degrade=1, En panne=2// Parametres d'emission par etat (moyennes et precisions)double[] moyennesEtat ={100.0,70.0,20.0};double[] precisionsEtat ={1.0/25.0,1.0/225.0,1.0/625.0};// sigma = 5, 15, 25// Matrice de transition suggeree (auto-persistance elevee):// Normal : [0.90, 0.09, 0.01]// Degrade : [0.10, 0.80, 0.10]// En panne : [0.01, 0.19, 0.80]Vector[] transMatrix =new Vector[]{ Vector.FromArray(0.90,0.09,0.01), Vector.FromArray(0.10,0.80,0.10), Vector.FromArray(0.01,0.19,0.80)};// TODO: Creer le moteur d'inference// TODO: Construire le HMM 3 etats// (Reutilisez la structure HMM de l'exemple guide en changeant nStates de 2 a 3)// TODO: Observer les donnees capteur et inferer les etats// TODO: Afficher les probabilites d'etat pour les jours clesConsole.WriteLine("Exercice a completer : adaptez le modele de classification a 3 etats.");
Exercice a completer : adaptez le modele de classification a 3 etats.
Conclusion
Ce notebook a explore les Hidden Markov Models avec Infer.NET : etats caches, emissions gaussiennes, transitions markoviennes, decodage de Viterbi, HMM factoriel et motif finding sur chaînes de caractères.
Concept
Point cle
HMM
Etats caches lies par transitions, observations via emissions
Classification independante
Chaque observation classee seule, sans contexte temporel (baseline)
HMM complet (Infer.NET)
ForEach + If(t==0)/If(t>0) + Switch ; possible, a condition de briser la symetrie
Brisure de symetrie
Init K-means + priors centres + precision intra-cluster ; sinon EP fusionne les etats
Chaînes independantes dont les contributions s’additionnent ; identifiable a permutation pres
Motif finding
Modèle generatif PWM + inference sur chaînes (StringFromArray), lissage laplacien
Distribution
Usage
Discrete
Etats caches (regimes, presence de motif)
Gaussian
Emissions conditionnelles aux etats
Dirichlet
Priors sur transitions et PWM (pseudo-comptes = lissage laplacien)
A retenir : Infer.NET sait chainer transitions et emissions (HMM, FHMM, motif finding). Les difficultes rencontrees ne sont pas des “limitations de la librairie” mais des problemes classiques d’initialisation, de priors et d’identifiabilite (label switching, permutation des chaînes) – qui se resolvent par K-means, priors informatifs et appariement.
Comment le modèle est infere
ExpectationPropagation, 50 itérations. EP propage les messages le long de la chaîne séquentielle.states[Tr].InitialiseTo(...)injecte l’assignation K-means comme point de depart, ce qui empeche EP de collapser tous les etats vers un même mode.Comme les etiquettes d’etats sont arbitraires (label switching), on apparie les etats inferes aux vrais regimes par proximite des moyennes d’emission (appariement bijectif glouton) avant de mesurer la justesse.