Infer-19 — Analyse de survie / fiabilite bayesienne : inferer le temps jusqu’a un événement

Corpus bayesien Infer.NET. Ce notebook (19e du corpus) ouvre la famille des modèles de duree (time-to-event) : la quantite centrale n’est plus une moyenne ou un compte, mais une fonction de survie S(t) = P(T > t) — probabilite qu’un composant, un patient ou un client survive au-dela de l’instant t. Il complete la famille « sequences dans le temps » (Infer-14 HMM, Infer-17 Kalman, Infer-18 Change-Point) par le cas ou l’observation n’est pas un etat recurrent mais un delai unique avant événement.

Plan. (1) Pourquoi un modèle dedie aux durees. (2) Le modèle exponentiel (cas conjugue). (3) Fonction de survie predictive en forme fermee. (4) Le modèle de Weibull via une astuce de transformee. (5) Sélection de la forme k. (6) Trois exercices (censure, regression AFT, comparaison de modèles).

1. Motivation : pourquoi pas une gaussienne ?

Un ingenieur fiabilite observe les durees de vie (en heures) de N composants identiques soumis a un test. La grandeur qui l’interesse est : « quelle probabilite qu’un composant fonctionne encore après 1500 h ? » — c’est-a-dire S(1500) = P(T > 1500).

Modeliser T par une gaussienne est inadequat pour trois raisons :

  1. Le temps est positif. Une gaussienne attribue une masse non nulle aux durees negatives.
  2. La distribution est asymetrique a droite. Quelques composants vivent très longtemps (longue queue), la moyenne depasse la mediane.
  3. Ce qui importe est la queue, pas le centre : S(t) dans la region ou peu de composants sont encore en vie.

Deux lois canoniques repondent a ces contraintes :

  • Exponentielle T ~ Exp(lambda) : taux de defaillance constant dans le temps (memoire, usure nulle). Un seul paramètre lambda > 0 (le taux).
  • Weibull T ~ Weibull(k, eta) : taux de defaillance qui varie — k < 1 (mortalite infantile, le composant se rodit), k = 1 (retombe sur l’exponentielle), k > 1 (usure, le risque croit avec l’age).

Le defi d’inference : estimer lambda, ou le couple (k, eta), a partir de durees observees, puis propager l’incertitude jusqu’a S(t) — avec un intervalle de credibilite, pas un chiffre unique. C’est exactement ce que le calcul bayesien via Infer.NET fournit.

#r "nuget: Microsoft.ML.Probabilistic, 0.4.2504.701"
#r "nuget: Microsoft.ML.Probabilistic.Compiler, 0.4.2504.701"
// restore Infer.NET -- isole dans sa propre cellule (convention de la serie)
Installed Packages
  • Microsoft.ML.Probabilistic, 0.4.2504.701
  • Microsoft.ML.Probabilistic.Compiler, 0.4.2504.701

Configuration des espaces de noms Infer.NET. On retrouve le moteur EP (Expectation Propagation) et les helpers. On ajoute une graine fixee pour que les durees synthetiques soient reproductibles (le test porte sur l’inference, pas sur le tirage).

using Microsoft.ML.Probabilistic;
using Microsoft.ML.Probabilistic.Distributions;
using Microsoft.ML.Probabilistic.Models;
using Microsoft.ML.Probabilistic.Algorithms;
using Range = Microsoft.ML.Probabilistic.Models.Range;   // desambiguise vs System.Range
using System;
using System.Linq;

var rand = new Random(42);

// Helper : tirage uniforme et exponentiel reproductibles.
double Uniform() => rand.NextDouble();                                   // U ~ Uniform(0,1)
double SampleExponential(double rate) => -Math.Log(Uniform()) / rate;     // T ~ Exp(rate)
double SampleWeibull(double shape, double scale)                         // T ~ Weibull(k, eta)
    => scale * Math.Pow(-Math.Log(Uniform()), 1.0 / shape);

Console.WriteLine("Infer.NET charge. Helpers Uniform / SampleExponential / SampleWeibull definis.");
Infer.NET charge. Helpers Uniform / SampleExponential / SampleWeibull definis.

2. Un jeu de données synthetique (verite connue)

On simule les durees de vie de N = 60 composants dont le vrai taux de defaillance est lambda_true = 1/1000 (duree moyenne caractéristique 1000 h). Connaissant la verite, on pourra verifier que le postieur la recouvre. On garde également un jeu Weibull avec usure (k = 1.8) pour la section 4.

const int N = 60;
const double lambdaTrue = 1.0 / 1000.0;   // duree moyenne 1000 h

double[] lifetimesExp = Enumerable.Range(0, N)
    .Select(_ => SampleExponential(lambdaTrue)).ToArray();

// Second jeu : Weibull avec usure (k=1.8), meme duree caracteristique (eta=1000).
const double kTrue = 1.8, etaTrue = 1000.0;
double[] lifetimesWei = Enumerable.Range(0, N)
    .Select(_ => SampleWeibull(kTrue, etaTrue)).ToArray();

Console.WriteLine($"Durees Exponentiel : N={N}, moyenne observee={lifetimesExp.Average():F1} h (attendu ~1000)");
Console.WriteLine($"Durees Weibull      : N={N}, moyenne observee={lifetimesWei.Average():F1} h (k=1.8 -> queue plus courte)");
Durees Exponentiel : N=60, moyenne observee=1262,9 h (attendu ~1000)
Durees Weibull      : N=60, moyenne observee=1041,0 h (k=1.8 -> queue plus courte)

3. Le modèle exponentiel : le cas conjugue

Le modèle le plus simple : T_i ~ Exp(lambda) pour chaque composant, avec un a priori Gamma sur le taux lambda. Le prior Gamma est conjugue a la vraisemblance exponentielle : EP renvoie donc un postieur exact, lui-même Gamma.

Lien de conjugaison. Si lambda ~ Gamma(a0, b0) (paramètre forme a0, paramètre taux b0, donc moyenne a0/b0) et T_i ~ Exp(lambda) pour i = 1..N, alors le postieur est Gamma(a0 + N, b0 + somme(T_i)). On retrouve ainsi, a la main, ce que le moteur calcule : chaque observation ajoute 1 a la forme et sa duree T_i au taux. Le prior faible Gamma(0.001, 0.001) laisse les données parler.

// Modele exponentiel conjugue sur le jeu de durees exponentielles.
Range item = new Range(N);
var lambda = Variable.GammaFromShapeAndRate(0.001, 0.001).Named("lambda");   // prior vague sur le taux
var T = Variable.Array<double>(item).Named("T");
T[item] = Variable.GammaFromShapeAndRate(1.0, lambda).ForEach(item);   // shape=1 => Exponentielle(lambda)
T.ObservedValue = lifetimesExp;

var engine = new InferenceEngine { Algorithm = new ExpectationPropagation() };
Gamma lambdaPost = engine.Infer<Gamma>(lambda);

Console.WriteLine("=== Posterieur du taux lambda (modele exponentiel) ===");
Console.WriteLine($"  Forme A   = {lambdaPost.Shape:F3}");
Console.WriteLine($"  Taux  B   = {lambdaPost.Rate:F3}");
Console.WriteLine($"  Moyenne   = {lambdaPost.GetMean():F6}  (vrai lambda = {lambdaTrue:F6})");
Console.WriteLine($"  Ecart-type= {Math.Sqrt(lambdaPost.GetVariance()):F6}");
Compiling model...done.
=== Posterieur du taux lambda (modele exponentiel) ===
  Forme A   = 60,001
  Taux  B   = 75774,766
  Moyenne   = 0,000792  (vrai lambda = 0,001000)
  Ecart-type= 0,000102

Lecture du postérieur : la conjugaison se vérifie chiffres en main

Le postérieur Gamma(A, B) renvoyé par EP coincide avec la formule de conjugaison :

  • A = 60,001 = 0,001 + N — la forme est l’a priori plus le nombre d’observations, au centieme pres.
  • B = 75 774,77 = 0,001 + Σ T_i — le taux est l’a priori plus la somme des durees (or Σ T_i / 60 = 1262,9 h, exactement la moyenne observee).

La moyenne postérieure A/B = 7,92e-4 est identique a l’estimateur du maximum de vraisemblance 1 / moyenne(T_i), parce que le prior Gamma(0,001 ; 0,001) est negligible devant 60 données. Elle est legerement inferieure au vrai lambda = 1e-3 non par biais du modèle, mais parce que cet echantillon a tire une moyenne de 1262,9 h (queue haute) au lieu de 1000 : avec N = 60, l’intervalle de confiance sur la moyenne est large, et le postérieur le reflete fidelement (ecart-type ~1,0e-4, soit un coefficient de variation de 13 %). C’est précisément l’intérêt du bayesien : on recupere non un chiffre, mais une distribution dont on propagera l’incertitude jusqu’a la fonction de survie.

4. Fonction de survie predictive : une forme fermee

La question de l’ingenieur — « probabilite de survivre au-dela de 1500 h » — se traduit par la survie predictive : on moyenne la survie exp(-lambda * t) sur le postieur de lambda.

Forme fermee. Si lambda ~ Gamma(A, B), alors par transformee de Laplace :

S(t) = E_lambda[ exp(-lambda * t) ] = ( B / (B + t) )^A

Pas besoin d’echantillonner : pour chaque t, on branche les moments du postieur Gamma et on obtient S(t) exactement, avec son incertitude qui decroit quand le postieur se concentre. On trace S(t) sur t ∈ [0, 3000] et on le compare a la courbe empirique (Kaplan-Meier simplifiee, ici sans censure : fraction d’observations superieures a t).

#load "SvgChartHelper.cs"
// Survie predictive (forme fermee) + comparaison empirique.
double A = lambdaPost.Shape, B = lambdaPost.Rate;
double Sbayes(double t) => Math.Pow(B / (B + t), A);
double Semp(double t, double[] data) => data.Count(d => d > t) / (double)data.Length;

// Tableau : S(t) bayesienne vs empirique a des horizons choisis.
Console.WriteLine("=== Fonction de survie S(t) : bayesienne vs empirique ===");
Console.WriteLine($"{"t (h)",8}  {"S_bayes",8}  {"S_empir",8}");
foreach (var t in new[] { 0.0, 250, 500, 1000, 1500, 2000, 3000 })
    Console.WriteLine($"{t,8:F0}  {Sbayes(t),8:F3}  {Semp(t, lifetimesExp),8:F3}");

// Courbe continue S(t) : trace Plotly comparant le modele bayesien aux donnees empiriques.
int NP = 80;
var tGrid = Enumerable.Range(0, NP + 1).Select(j => j * 3000.0 / NP).ToArray();
SvgChartHelper.Overlay(
    "Fonction de survie S(t) : bayesienne vs empirique",
    "t (heures)", "S(t)",
    new[] {
        new SvgSeries("Bayesienne (modele exponentiel)",
            tGrid, tGrid.Select(Sbayes).ToArray(), TraceStyle.Line, "#2a6dba"),
        new SvgSeries("Empirique (donnees)",
            tGrid, tGrid.Select(t => Semp(t, lifetimesExp)).ToArray(),
            TraceStyle.Markers, "#e8743b"),
    })
=== Fonction de survie S(t) : bayesienne vs empirique ===
   t (h)   S_bayes   S_empir
       0     1,000     1,000
     250     0,821     0,900
     500     0,674     0,700
    1000     0,455     0,450
    1500     0,308     0,317
    2000     0,209     0,183
    3000     0,097     0,100
Fonction de survie S(t) : bayesienne vs empirique0.0430.2960.5490.8011.0540750150022503000t (heures)S(t)Bayesienne (modele exponentiel)Empirique (donnees)

Lecture : la question de l’ingenieur a maintenant une reponse quantifiee

La question « probabilite qu’un composant survive au-dela de 1500 h » admet la reponse S(1500) ≈ 0,31 : environ un composant sur trois est encore en vie a 1500 h. Deux points :

  • Accord bayesien / empirique. Les deux courbes se suivent (ecart maximal ~0,08 a t = 250) mais de nature différente : la bayesienne est parametrique et lisse (elle croit au modèle exponentiel et projette sa forme), l’empirique est en escalier (fraction brute des durees superieures a t). Leur proximite valide que l’exponentielle decrit bien ce jeu.
  • Forme fermee exacte. S(t) = (B/(B+t))^A n’est pas une approximation Monte-Carlo : c’est la transformee de Laplace du postieur Gamma, calculee en un seul coup. On peut donc donner un intervalle de credibilite sur S(t) en tirant des lambda du postieur — l’incertitude sur le taux se propage gratuitement jusqu’a la queue, la ou l’enjeu fiabilite se trouve.

5. Le modèle de Weibull : une astuce de transformee

L’exponentielle suppose un taux constant (ni rodage, ni usure). La loi de Weibull generalise : T ~ Weibull(k, eta) avec fonction de survie S(t) = exp(-(t/eta)^k).

Origine. La loi de Weibull a ete proposee par Waloddi Weibull (Weibull, 1951, J. Applied Mechanics 18(3):293-297) pour decrire la resistance de materiaux. Sa richesse tient a la forme k : k < 1 modelise un defaut initial (mortalite precoce), k = 1 redonne l’exponentielle (taux constant), k > 1 une usure accumulee – c’est ce meme k que la section 6 retrouvera par balayage.

Inferer directement la forme k par EP est delicat (la vraisemblance de Weibull n’est conjuguee a aucun prior usuel). Mais si l’on fixe k, une transformee ramene le problème au cas conjugue :

U_i = T_i^k  ==>  U_i ~ Exp(rate = eta^{-k})

On infere donc le taux transforme r = eta^{-k} avec un prior Gamma (conjugue), puis on remonte a eta = r^{-1/k}. La survie predictive devient S(t) = (B / (B + t^k))^A (même forme fermee, avec t remplace par t^k). La section suivante choisira le meilleur k par balayage.

// Weibull a forme k fixee : transformee U = T^k ~ Exp(r), r = eta^{-k}.
double InferWeibullRate(double k, double[] data, out Gamma rPost)
{
    Range it = new Range(data.Length);
    var r = Variable.GammaFromShapeAndRate(0.001, 0.001).Named("r_" + k.ToString("F1"));
    var U = Variable.Array<double>(it).Named("U_" + k.ToString("F1"));
    U[it] = Variable.GammaFromShapeAndRate(1.0, r).ForEach(it);   // shape=1 => Exponentielle(r)
    U.ObservedValue = data.Select(t => Math.Pow(t, k)).ToArray();
    var eng = new InferenceEngine { Algorithm = new ExpectationPropagation() };
    rPost = eng.Infer<Gamma>(r);
    return rPost.GetMean();
}

double kFixed = 1.8;   // (on retrouvera ce k par balayage a la section 6)
double rMean = InferWeibullRate(kFixed, lifetimesWei, out Gamma rPost);
double etaMean = Math.Pow(rMean, -1.0 / kFixed);

Console.WriteLine("=== Modele Weibull (k=1.8 fixe) sur le jeu Weibull ===");
Console.WriteLine($"  r = eta^-k      = {rMean:E3}");
Console.WriteLine($"  eta (back-out)  = {etaMean:F1} h  (vrai eta = {etaTrue:F1})");
Console.WriteLine($"  Survie a 1000 h = {Math.Pow(rPost.Rate/(rPost.Rate + Math.Pow(1000,kFixed)), rPost.Shape):F3}");
Compiling model...done.
=== Modele Weibull (k=1.8 fixe) sur le jeu Weibull ===
  r = eta^-k      = 2,926E-006
  eta (back-out)  = 1186,7 h  (vrai eta = 1000,0)
  Survie a 1000 h = 0,482

Lecture : la transformee ramene Weibull au cas conjugue

L’astuce U = T^k ~ Exp(rate = eta^{-k}) fonctionne : EP infere un postérieur Gamma sur r sans aucune fragilite (pas d’inference directe sur la forme). On remonte eta = r^{-1/k} ≈ 1187 h, a comparer au vrai eta = 1000 h. La surestimation vient, comme pour l’exponentiel, de l’echantillon : sa moyenne (1041 h) est un peu haute, et avec k = 1,8 on a moyenne = eta · Γ(1 + 1/k) ≈ eta · 0,89, d’ou un eta implique par les données autour de 1170 h.

La survie estimee S(1000) ≈ 0,48 (contre e^{-1} ≈ 0,37 pour le vrai Weibull(1,8 ; 1000)) reflete cette même surestimation de eta. Ce n’est pas un defaut de la méthode mais du bruit d’echantillonnage a N = 60 ; l’apport pedagogique est ailleurs : toute la richesse de Weibull (usure, rodage) est accessible sans payer le cout d’une inference EP sur la forme, des lors que l’on fixe k (ou qu’on le choisit par balayage, section suivante).

6. Sélection de la forme k : balayage par log-vraisemblance predictive

Comment choisir k sans payer le cout d’une inference EP sur la forme ? On balaye k sur une grille et, pour chaque k, on calcule la log-vraisemblance predictive (leave-one-out approchee) du jeu observe : pour chaque duree T_i, la densite Weibull evaluee avec le postieur entraene sur les autres. Le k qui maximise cette score est le meilleur compromis biais/variance.

Sur le jeu exponentiel (taux constant), k ≈ 1 doit gagner ; sur le jeu Weibull k = 1.8, c’est k ≈ 1.8 qui doit gagner. C’est un test de coherence interne.

// Balayage de k : pour chaque k, log-vraisemblance predictive (LOO approchee).
double ShapeScore(double k, double[] data)
{
    // Posterieur du taux transforme r sur TOUT le jeu, densite Weibull evaluee en chaque point.
    // (approximation LOO : le postieur est peu sensible a un seul point quand N est grand.)
    double rMean = InferWeibullRate(k, data, out Gamma rPost);
    double A_ = rPost.Shape, B_ = rPost.Rate;
    // Densite de T_i sous Weibull(k, eta) avec eta = r^{-1/k}, r ~ Gamma(A,B).
    // On evalue la log-densite marginale approchee a la moyenne post.: r_hat = A/B.
    double rHat = A_ / B_;
    double eta = Math.Pow(rHat, -1.0 / k);
    double lp = 0.0;
    foreach (var t in data)
    {
        double z = t / eta;
        lp += Math.Log(k) - k * Math.Log(eta) + (k - 1) * Math.Log(t) - Math.Pow(z, k);
    }
    return lp;
}

Console.WriteLine("=== Selection de la forme k (log-vraisemblance predictive) ===");
Console.WriteLine($"{"k",5}  {"score Exp",12}  {"score Weibull",14}");
var ks = new[] { 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 2.5 };
double bestExp = double.NegativeInfinity, bestWei = double.NegativeInfinity;
double kExp = 1, kWei = 1;
foreach (var k in ks)
{
    double se = ShapeScore(k, lifetimesExp);
    double sw = ShapeScore(k, lifetimesWei);
    Console.WriteLine($"{k,5:F1}  {se,12:F1}  {sw,14:F1}");
    if (se > bestExp) { bestExp = se; kExp = k; }
    if (sw > bestWei) { bestWei = sw; kWei = k; }
}
Console.WriteLine($"\n  k* (jeu exponentiel) = {kExp:F1}  (attendu ~1.0)");
Console.WriteLine($"  k* (jeu Weibull)     = {kWei:F1}  (attendu ~1.8)");
=== Selection de la forme k (log-vraisemblance predictive) ===
    k     score Exp   score Weibull
Compiling model...done.
Compiling model...done.
  0,8        -494,0          -485,7
Compiling model...done.
Compiling model...done.
  1,0        -488,5          -476,9
Compiling model...done.
Compiling model...done.
  1,2        -486,8          -471,3
Compiling model...done.
Compiling model...done.
  1,5        -489,4          -467,1
Compiling model...done.
Compiling model...done.
  1,8        -496,7          -466,7
Compiling model...done.
Compiling model...done.
  2,0        -503,5          -468,1
Compiling model...done.
Compiling model...done.
  2,5        -526,1          -475,7

  k* (jeu exponentiel) = 1,2  (attendu ~1.0)
  k* (jeu Weibull)     = 1,8  (attendu ~1.8)

Lecture : le balayage recupere la vraie forme et distingue les deux regimes

Le test de coherence reussit sur les deux jeux :

  • Jeu Weibull (k = 1,8) : le score est minimal en k = 1,8 (-466,7), c’est-a-dire exactement la valeur vraie. Le pic est net (k = 1,5 et k = 2,0 sont déjà moins bons). La méthode retrouve la signature d’usure.
  • Jeu exponentiel (k = 1) : le meilleur est k = 1,2 (-486,8), mais k = 1,0 (-488,5) est a moins de 2 unites : ex aequo dans le bruit d’echantillonnage. Leger tilt vers 1,2 parce que cet echantillon exponentiel n’etait pas un exponentiel parfait. On conclut donc raisonnablement k ≈ 1 (taux constant), ce qui est la bonne reponse.

Cette approche par grille evite l’inference EP directe sur la forme (notoirement instable) au prix de quelques compilations de modèle (7 ici, chacune instantanee). C’est le compromis operationnel : robustesse contre temps de calcul, dans l’esprit de la règle « realiser le moteur, pas le contourner ».

7. Censure a droite : le piege du modele naif, corrige par la vraisemblance (exemple execute)

En fiabilite reelle, beaucoup de composants survivent a la fin du test : on n’observe pas leur duree exacte T_i, seulement qu’elle depasse la duree du test c_i (censure administrative a c*, la meme pour tous). La contribution de vraisemblance d’un point censure n’est pas la densite f(c_i) mais la survie S(c_i) = exp(-lambda * c_i).

Nous executons la comparaison complete sur le jeu exponentiel de la section 3, censuré a c* = 800 h — environ 45 % de censures :

  • Modele naif : les durees enregistrees (censures figees a c*) traitees comme des evenements exacts. Chaque censure est comptee comme une mort precoce : le taux estime monte, la survie estimee s’effondre.
  • Modele censure : f(t_i) pour les evenements, S(c_i) pour les censures. La forme « temps latent contraint T_i > c* » (Variable.ConstrainPositive) n’est pas compilable dans Infer.NET — ni EP ni VMP n’ont d’operateur pour la difference T_i - c* d’une variable Gamma (verifie empiriquement : CompilationFailedException). On branche donc la vraisemblance par sa forme exacte en statistiques suffisantes : Gamma(n_obs, lambda) observe au temps total a risque, meme vraisemblance, entierement conjuguee. La ou PyMC-19 ecrit S(c_i) a la main via pm.Potential (le MCMC accepte des facteurs arbitraires), Infer.NET exige une structure conjuguee — une vraie difference de paradigme entre les deux moteurs.
  • Kaplan–Meier : reference non parametrique (Kaplan & Meier, 1958), l’escalier empirique adapte a la censure — la ou l’empirique naif de la section 4 est fausse par construction.

L’ancienne version de cette section etait un exercice non resolu ; le contraste naif/censure etant le message central de l’analyse de survie, il est desormais demontre execute, et les exercices (section 8) l’etendent.

// --- Censure administrative a c* = 800 h sur le jeu exponentiel (meme protocole que PyMC-19) ---
const double cStar = 800.0;
bool[] isEvent = lifetimesExp.Select(t => t <= cStar).ToArray();   // evenement observe avant la fin du test
double[] tRecorded = lifetimesExp.Select((t, i) => isEvent[i] ? t : cStar).ToArray();
int nObs = isEvent.Count(b => b), nCens = N - nObs;
double[] tEvents = lifetimesExp.Where(t => t <= cStar).ToArray();

Console.WriteLine($"Censure administrative a c* = {cStar:F0} h sur N = {N} composants");
Console.WriteLine($"  Evenements observes : {nObs}   (durees exactes)");
Console.WriteLine($"  Censures a droite   : {nCens}   (T_i > {cStar:F0}, seule l'information 'survit a c*' est connue)");
Console.WriteLine($"  Fraction censuree   : {nCens / (double)N:P1}   (theorie : exp(-lambda c*) = {Math.Exp(-lambdaTrue * cStar):P1})");
Censure administrative a c* = 800 h sur N = 60 composants
  Evenements observes : 29   (durees exactes)
  Censures a droite   : 31   (T_i > 800, seule l'information 'survit a c*' est connue)
  Fraction censuree   : 51,7 %   (theorie : exp(-lambda c*) = 44,9 %)
// --- Modele NAIF : toutes les durees enregistrees traitees comme des evenements exacts ---
Range itemN = new Range(N);
var lambdaN = Variable.GammaFromShapeAndRate(0.001, 0.001).Named("lambdaN");
var Tn = Variable.Array<double>(itemN).Named("Tn");
Tn[itemN] = Variable.GammaFromShapeAndRate(1.0, lambdaN).ForEach(itemN);
Tn.ObservedValue = tRecorded;   // ERREUR du naif : les censures comptent comme des morts a c*

var engineN = new InferenceEngine(new ExpectationPropagation());
Gamma lambdaNPost = engineN.Infer<Gamma>(lambdaN);
double lambdaNaif = lambdaNPost.GetMean();

Console.WriteLine("=== Modele NAIF (censures traitees comme des evenements a c*) ===");
Console.WriteLine($"  lambda_naif   = {lambdaNaif:F6}   (vrai lambda = {lambdaTrue:F6}, biais de {lambdaNaif / lambdaTrue:F2}x)");
Console.WriteLine($"  S_naif(1500)  = {Math.Exp(-lambdaNaif * 1500):F3}   (vraie survie = {Math.Exp(-lambdaTrue * 1500):F3})");
Compiling model...done.
=== Modele NAIF (censures traitees comme des evenements a c*) ===
  lambda_naif   = 0,001617   (vrai lambda = 0,001000, biais de 1,62x)
  S_naif(1500)  = 0,088   (vraie survie = 0,223)
// --- Modele CENSURE : f(t_i) pour les evenements, S(c_i) pour les censures ---
// Vraisemblance complete : prod f(t_i) x prod S(c_i) = lambda^nObs * exp(-lambda * (somme t_obs + somme c_i)).
// Deux voies pour la brancher dans Infer.NET :
//  (a) forme contrainte : garder T_i ~ Exp(lambda) latent et contraindre T_i > c* --
//      EMPIRIQUEMENT non compilable (ni EP ni VMP : aucun operateur pour
//      Factor.Difference(Gamma, const), CompilationFailedException) ;
//  (b) forme exacte par statistiques suffisantes : Gamma(nObs, lambda) observe au temps total
//      a risque porte EXACTEMENT la meme vraisemblance lambda^nObs * exp(-lambda * total),
//      entierement conjuguee. C'est la voie retenue -- et c'est une lecon : la conjugaison
//      ne consomme jamais les donnees brutes, seulement leurs statistiques suffisantes.
double totalTempsARisque = tEvents.Sum() + nCens * cStar;   // somme t_obs + somme c_i
var lambdaC = Variable.GammaFromShapeAndRate(0.001, 0.001).Named("lambdaC");
var Tsum = Variable.GammaFromShapeAndRate((double)nObs, lambdaC).Named("Tsum");  // nObs evenements, temps cumule
Tsum.ObservedValue = totalTempsARisque;

var engineC = new InferenceEngine(new ExpectationPropagation());
Gamma lambdaCPost = engineC.Infer<Gamma>(lambdaC);
double lambdaCensure = lambdaCPost.GetMean();

// Verification conjuguee : posterieur exact = Gamma(a + nObs, b + somme(t_obs) + somme(c_i))
Console.WriteLine("=== Modele CENSURE (f(t_i) pour les evenements, S(c_i) pour les censures) ===");
Console.WriteLine($"  lambda_censure = {lambdaCensure:F6}   (vrai lambda = {lambdaTrue:F6})");
Console.WriteLine($"  Forme fermee conjuguee = {(0.001 + nObs) / (0.001 + totalTempsARisque):F6}  -> l'EP retrouve la conjugaison");
Console.WriteLine($"  S_censure(1500) = {Math.Exp(-lambdaCensure * 1500):F3}   (vraie survie = {Math.Exp(-lambdaTrue * 1500):F3})");
Console.WriteLine($"\n  Le terme S(c_i) ajoute {nCens * cStar:F0} h de temps-a-risque sans ajouter d'evenement :");
Console.WriteLine("  chaque composant censure a 'travaille' au-dela de c* sans qu'on puisse le compter comme mort.")
Compiling model...done.
=== Modele CENSURE (f(t_i) pour les evenements, S(c_i) pour les censures) ===
  lambda_censure = 0,000782   (vrai lambda = 0,001000)
  Forme fermee conjuguee = 0,000782  -> l'EP retrouve la conjugaison
  S_censure(1500) = 0,310   (vraie survie = 0,223)

  Le terme S(c_i) ajoute 24800 h de temps-a-risque sans ajouter d'evenement :
  chaque composant censure a 'travaille' au-dela de c* sans qu'on puisse le compter comme mort.
// --- Kaplan-Meier (non parametrique, sans modele) + confrontation graphique ---
(double[] kmT, double[] kmS) KaplanMeier(double[] tRec, bool[] evt)
{
    var order = Enumerable.Range(0, tRec.Length).OrderBy(i => tRec[i]).ToArray();
    double S = 1.0;
    var xs = new List<double> { 0.0 }; var ys = new List<double> { 1.0 };
    foreach (int i in order)
    {
        if (!evt[i]) continue;                       // seuls les evenements reduisent la survie
        int atRisk = order.Count(j => tRec[j] >= tRec[i]);
        S *= 1.0 - 1.0 / atRisk;
        xs.Add(tRec[i]); ys.Add(S);
    }
    return (xs.ToArray(), ys.ToArray());
}

var (kmT, kmS) = KaplanMeier(tRecorded, isEvent);
double SKm(double t) { double s = 1.0; for (int j = 0; j < kmT.Length; j++) if (kmT[j] <= t) s = kmS[j]; return s; }
double SMod(double lam, double t) => Math.Exp(-lam * t);

Console.WriteLine("=== S(t) : verite simulee vs naive vs censuree vs Kaplan-Meier ===");
Console.WriteLine($"{"t (h)",8}  {"verite",7}  {"naif",7}  {"censure",8}  {"K-M",6}");
foreach (var t in new[] { 0.0, 250, 500, 1000, 1500, 2000, 3000 })
    Console.WriteLine($"{t,8:F0}  {SMod(lambdaTrue, t),7:F3}  {SMod(lambdaNaif, t),7:F3}  {SMod(lambdaCensure, t),8:F3}  {SKm(t),6:F3}");
Console.WriteLine($"\nEmpirique NAIF (comme en section 4) : S(1500) = {tRecorded.Count(t => t > 1500) / (double)N:F3}"
                  + $"  <- fausse aussi : les censures a {cStar:F0} h y sont comptees mortes.");

int NPK = 80;
var tGridK = Enumerable.Range(0, NPK + 1).Select(j => j * 3000.0 / NPK).ToArray();
SvgChartHelper.Overlay(
    "Censure a droite : verite vs naive vs censuree vs Kaplan-Meier",
    "t (heures)", "S(t)",
    new[] {
        new SvgSeries("Verite simulee (lambda = 1/1000)", tGridK, tGridK.Select(t => SMod(lambdaTrue, t)).ToArray(), TraceStyle.Line, "#4daf4a"),
        new SvgSeries("Modele naif (censures = morts)", tGridK, tGridK.Select(t => SMod(lambdaNaif, t)).ToArray(), TraceStyle.Line, "#e41a1c"),
        new SvgSeries("Modele censure (S(c_i))", tGridK, tGridK.Select(t => SMod(lambdaCensure, t)).ToArray(), TraceStyle.Line, "#2a6dba"),
        new SvgSeries("Kaplan-Meier", kmT, kmS, TraceStyle.Markers, "#e8743b"),
    })
=== S(t) : verite simulee vs naive vs censuree vs Kaplan-Meier ===
   t (h)   verite     naif   censure     K-M
       0    1,000    1,000     1,000   1,000
     250    0,779    0,667     0,822   0,900
     500    0,607    0,445     0,676   0,700
    1000    0,368    0,198     0,458   0,517
    1500    0,223    0,088     0,310   0,517
    2000    0,135    0,039     0,209   0,517
    3000    0,050    0,008     0,096   0,517

Empirique NAIF (comme en section 4) : S(1500) = 0,000  <- fausse aussi : les censures a 800 h y sont comptees mortes.
Censure a droite : verite vs naive vs censuree vs Kaplan-Meier-0.0520.2260.5040.7821.060750150022503000t (heures)S(t)Verite simulee (lambda = 1/1000)Modele naif (censures = morts)Modele censure (S(c_i))Kaplan-Meier

Lecture : EP, contrainte de positivite et forme fermee

  • Le biais naif est un biais de comptage : les nCens censures enregistrees comme des morts a c* gonflent le taux d’un facteur exact N / nObs (naif/censure) — sans qu’aucun composant supplementaire soit mort, la fiabilite estimee s’effondre (S_naif(1500) contre la vraie survie dans la sortie ci-dessus).
  • Statistiques suffisantes = ce que la conjugaison consomme : le posterieur reste Gamma(a + nObs, b + somme(t_obs) + somme(c_i)) — chaque censure ajoute c* heures de temps a risque sans ajouter d’evenement, et l’EP retrouve la forme fermee chiffres en main. La contrainte « T_i > c* », non compilable dans Infer.NET (aucun operateur Difference pour une Gamma), et l’observation Gamma(nObs, lambda) au temps total portent la MEME vraisemblance : la sufficience est la passerelle entre les deux formes.
  • Kaplan–Meier s’arrete ou finit l’observation : l’escalier suit la verite tant que des evenements restent observables, puis devient plat au-dela de c* = 800 h — ce n’est pas une estimation de S(3000) : avec une censure administrative commune, plus aucune duree observee ne depasse c*, le risque est vide et l’estimateur non parametrique ne peut pas extrapoler. Le modele parametrique, lui, extrapole (au prix de son hypothese de loi) — la complementarite non-parametrique/parametrique.
  • Pont moteur : PyMC-19 mene la meme comparaison cote NUTS — la censure y entre par un pm.Potential(-lambda * c_i) ecrit a la main, ici par Variable.ConstrainPositive. Meme protocole, meme forme fermee, deux machines.

8. Exercices

Les trois exercices ci-dessous etendent le modele de duree au-dela de l’exemple execute de la section 7. Ils sont laisses a completer (stubs sans erreur) — le notebook s’execute de bout en bout meme non rempli.

Exercice 1 — Regression de temps accelere (AFT)

La duree depend d’une covariable (temperature, contrainte). Modèle AFT (Accelerated Failure Time) : le taux devient spécifique au composant, lambda_i = lambda0 * exp(beta * x_i) ou x_i est la covariable centree et beta le coefficient a inferer.

Indice : declarez lambda0 ~ Gamma et beta ~ Gaussian(0, grande), puis bouclez sur les composants avec Variable.Exponential(lambda0 * Variable.Exp(beta * x[i])). Inferer beta donne l’effet de la contrainte sur la duree de vie (signe et amplitude). Attention : le produit dans le taux rend le postieur non conjugue — EP le traite quand même (approximation).

// Exercice 1 : regression AFT (a completer)
// TODO etudiant : lambda_i = lambda0 * exp(beta * x_i), inferer beta.
Console.WriteLine("Exercice 1 a completer : regression de temps accelere (AFT).");
Exercice 1 a completer : regression de temps accelere (AFT).

Exercice 2 — Exponentiel vs Weibull : facteur de Bayes

Sur un même jeu de durees, comparer le modèle exponentiel (k = 1) au modèle Weibull (k libre) via un facteur de Bayes (rapport des vraisemblances marginales). Le verdict depend de la taille N et du vrai k.

Indice : la section 6 calcule déjà un score par k. Generalisez-le en vraisemblance marginale (intégrez sur le postieur de r, pas juste evaluez a la moyenne), puis formez le rapport BF = p(D | Weibull) / p(D | exponentiel). A partir de quel N le Weibull k = 1.8 est-il prefere de facon decisive (log BF > 5) ?

// Exercice 2 : facteur de Bayes exponentiel vs Weibull (a completer)
// TODO etudiant : vraisemblance marginale par k, rapport BF, seuil en N.
Console.WriteLine("Exercice 2 a completer : facteur de Bayes exponentiel vs Weibull.");
Exercice 2 a completer : facteur de Bayes exponentiel vs Weibull.

Exercice 3 — Sensibilite au taux de censure

L’exemple de la section 7 fixe c* = 800 h (~45 % de censures). Refaites tourner le trio naif / censure / Kaplan–Meier pour c* dans {600, 1000, 1400} h et rapportez le biais relatif du naif, S_naif(1500) / S_verite(1500), en fonction de la fraction censurée.

  • Indice : seul cStar change — encapsulez la construction du jeu et les deux modeles dans une methode ComparerCensure(double cStar) et bouclez. A c* = 1400 h, peu de censures subsistent et le biais naif devient petit : le phenomene est non lineaire en la fraction censurée.
  • Etape 1 : ecrire la boucle sur les trois valeurs de c*.
  • Etape 2 : rapporter fraction censurée et biais relatif pour chaque c*.
// Exercice 3 : sensibilite du biais naif au taux de censure (a completer)
// TODO etudiant : boucler sur cStar in {600, 1000, 1400}, recalculer lambda_naif / lambda_censure
// et le biais relatif S_naif(1500)/S_verite(1500) en fonction de la fraction censurée.
Console.WriteLine("Exercice 3 a completer : sensibilite du biais naif au taux de censure.");
Exercice 3 a completer : sensibilite du biais naif au taux de censure.

Conclusion

L’analyse de survie complete la famille des modèles temporels du corpus Infer :

Notebook Nature du temps Quantite inferee
Infer-14 (HMM) discret, recurrent trajectoire d’etats caches
Infer-17 (Kalman) continu, recurrent etat filtre pas-a-pas
Infer-18 (Change-Point) continu, a rupture instant d’un changement de regime
Infer-19 (Survie) continu, duree unique fonction de survie S(t)

A retenir.

  • Le temps jusqu’a un événement se modelise par une loi positive asymetrique (exponentielle, Weibull), pas une gaussienne — ce qui compte est la queue, i.e. S(t).
  • Le cas exponentiel est conjugue : prior Gamma sur le taux, postieur Gamma exact. La survie predictive se ramene a une forme fermée (B/(B+t))^A (transformee de Laplace du Gamma).
  • Le cas Weibull se ramene a l’exponentiel par transformee U = T^k des que k est fixe, evitant la fragilite de EP sur la forme. On choisit k par balayage (log-vraisemblance predictive), non par inference directe.
  • La censure, la regression AFT et la sélection par facteur de Bayes sont les prolongements naturels (exercices).

Pour aller plus loin. La censure a gauche/par intervalle, les modèles a risques concurrents (competing risks) et les modèles de fragilite partagee (frailty, analogue hiérarchique de Infer-12 applique a la survie) sont les extensions classiques au-dela de ce notebook.

References

Sources fondatrices (papiers primaires).

  • Weibull, W. (1951), « A Statistical Distribution Function of Wide Applicability », Journal of Applied Mechanics 18(3):293-297 — origine de la loi de Weibull (sections 5-6).

  • Kaplan, E. L. & Meier, P. (1958), « Nonparametric Estimation from Incomplete Observations », JASA 53(282):457-481 — estimateur de survie sous censure à droite (exercice 1).

  • Cox, D. R. (1972), « Regression Models and Life-Tables », JRSS Series B 34(2):187-202 — modèle à risques proportionnels, l’extension de référence pour la régression sur la survie (au-delà du present notebook).

  • Gelman A. et al., Bayesian Data Analysis (3e ed.), §2.6 (modèle exponentiel pour données de durées de vie).

  • Lawless J. F., Statistical Models and Methods for Lifetime Data — reference pour Weibull/AFT.

  • Infer.NET documentation : Variable.GammaFromShapeAndRate (prior sur le taux ; shape=1 donne l’exponentielle), Variable.Weibull.

  • Relation a la serie : Infer-12 (modèles hiérarchiques), Infer-18 (rupture — ou le taux change brusquement).

Retour au sommet