19. Analyse de survie / fiabilite bayesienne : inferer le temps jusqu’a un événement (jumeau PyMC)

Serie Parite .NET <=> Python (#4956). Ce notebook est le jumeau Python de Infer-19 – Analyse de survie (Infer.NET, EP). Il reprend le même terrain de jeu (même germe, même N, mêmes vrais paramètres) et y ajoute deux value-add PyMC : (1) inferer directement la forme Weibull k par MCMC (qu’Infer.NET evite par transformee + balayage, l’EP sur la forme etant fragile) ; (2) sélection de modèle par LOO cross-validation (arviZ) au lieu d’un balayage manuel.

Duree estimee : 50 minutes Prerequis : PyMC-2 – Melanges gaussiens (conjugaison Gamma-Exponentielle)

Objectifs

  • Comprendre pourquoi le temps jusqu’a un événement se modelise par une loi positive asymetrique (exponentielle, Weibull), pas une gaussienne
  • Implementer le modèle exponentiel conjugue (prior Gamma sur le taux) et propager l’incertitude jusqu’a la fonction de survie \(S(t) = P(T > t)\)
  • Value-add PyMC : inferer directement la forme Weibull k par NUTS (Infer.NET l’evite)
  • Value-add PyMC : sélectionner le modèle (Exp vs Weibull) par LOO (arviZ), pas un balayage manuel

1. Motivation : pourquoi pas une gaussienne ?

Un ingenieur fiabilite observe les durees de vie (en heures) de N composants. 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 : (1) le temps est positif (une gaussienne donne de la masse aux durees negatives) ; (2) la distribution est asymetrique a droite (longue queue) ; (3) ce qui importe est la queue, pas le centre.

Deux lois canoniques : - Exponentielle \(T \sim \text{Exp}(\lambda)\) : taux de defaillance constant (usure nulle). - Weibull \(T \sim \text{Weibull}(k, \eta)\) : taux qui varie – \(k < 1\) (mortalite infantile), \(k = 1\) (retombe sur l’exponentielle), \(k > 1\) (usure, le risque croit avec l’age).

# --- Imports + donnees synthetiques (meme vrais parametres qu'Infer-19) ---
import numpy as np
import pymc as pm
import arviz as az
import warnings
warnings.filterwarnings("ignore", message=".*data structure.*")
warnings.filterwarnings("ignore", message="PyTensor could not link to a BLAS")  # advisory pytensor (#3436)
print(f"PyMC {pm.__version__}, ArviZ {az.__version__}")

# Meme N et memes vrais parametres qu'Infer-19 (germe 42 ; NB: numpy.default_rng != .NET Random,
# donc l'echantillon tiré differe, mais la verite sous-jacente et les conclusions sont identiques).
rng = np.random.default_rng(42)
N = 60
lambda_true = 1.0 / 1000.0      # duree moyenne caracteristique 1000 h
k_true, eta_true = 1.8, 1000.0  # Weibull avec usure

lifetimes_exp = rng.exponential(scale=1.0 / lambda_true, size=N)   # T ~ Exp(lambda)
lifetimes_wei = (eta_true * rng.weibull(k_true, size=N))           # T ~ Weibull(k, eta)

print(f"Durees Exponentiel : N={N}, moyenne observee={lifetimes_exp.mean():.1f} h (attendu ~1000)")
print(f"Durees Weibull      : N={N}, moyenne observee={lifetimes_wei.mean():.1f} h (k=1.8 -> queue plus courte)")
PyMC 6.0.1, ArviZ 1.1.0
Durees Exponentiel : N=60, moyenne observee=823.1 h (attendu ~1000)
Durees Weibull      : N=60, moyenne observee=918.9 h (k=1.8 -> queue plus courte)

2. Le modèle exponentiel : le cas conjugue

Le modèle le plus simple : \(T_i \sim \text{Exp}(\lambda)\) avec un a priori Gamma sur le taux \(\lambda\). Le prior Gamma est conjugue a la vraisemblance exponentielle : le postieur est donc lui-même Gamma, exact (comme le retrouve l’EP d’Infer.NET dans Infer-19).

\[\lambda \sim \text{Gamma}(a_0, b_0), \quad T_i \sim \text{Exp}(\lambda), \quad \lambda \mid T \sim \text{Gamma}(a_0 + N,\, b_0 + \textstyle\sum_i T_i)\]

PyMC echantillonne ce postieur par NUTS (ici il retrouve la solution conjuguee, c’est un cas d’école).

# --- Modele exponentiel conjugue sur le jeu de durees exponentielles ---
with pm.Model() as modele_exp:
    # Prior faible Gamma(0.001, 0.001) sur le taux lambda (laisse les donnees parler).
    # pm.Gamma parameterise par shape (alpha) et rate (beta) : moyenne = alpha/beta.
    lam = pm.Gamma("lam", alpha=0.001, beta=0.001)
    pm.Exponential("T", lam=lam, observed=lifetimes_exp)
    trace_exp = pm.sample(1000, tune=1000, chains=2, target_accept=0.9,
                          random_seed=42, progressbar=False)

# Log-vraisemblance calculee a part : requise pour la selection de modele par LOO
# (cellule suivante) -- les PyMC recents ne la stockent plus pendant le sampling.
pm.compute_log_likelihood(trace_exp, model=modele_exp)

post_lam = trace_exp.posterior["lam"]
print("=== Posterieur du taux lambda (modele exponentiel) ===")
print(f"  Moyenne   = {float(post_lam.mean()):.6f}  (vrai lambda = {lambda_true:.6f})")
print(f"  Ecart-type= {float(post_lam.std()):.6f}")
print(f"  Intervalle 94% : [{float(post_lam.quantile(0.03)):.6f}, {float(post_lam.quantile(0.97)):.6f}]")






=== Posterieur du taux lambda (modele exponentiel) ===
  Moyenne   = 0.001217  (vrai lambda = 0.001000)
  Ecart-type= 0.000157
  Intervalle 94% : [0.000939, 0.001524]








Lecture du postérieur


La moyenne postérieure est proche de l’estimateur du maximum de vraisemblance \(1/\overline{T}\) (parce que le prior faible est negligible devant 60 données). L’ecart-type reflete l’incertitude d’echantillonnage a \(N = 60\). C’est 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.








3. Fonction de survie predictive : propager l’incertitude jusqu’a la queue


La question de l’ingenieur – « probabilite de surviver au-dela de 1500 h » – se traduit par la survie predictive : on evalue \(S(t) = e^{-\lambda t}\) pour chaque echantillon \(\lambda\) du postérieur, puis on resume la distribution de \(S(t)\) (medianne + intervalle de credibilite).


Infer-19 (Infer.NET) exploitait une forme fermee exacte \(S(t) = (B/(B+t))^A\) (transformee de Laplace du postérieur Gamma) – economisant le Monte-Carlo. Ici, PyMC echantillonnant déjà le postérieur, on propage directement les tirages : plus general (marche pour tout modèle, pas seulement le conjugue), au prix d’un peu de Monte-Carlo.




# --- Fonction de survie predictive (propagation des tirages postérieurs) ---
ts = np.array([0, 250, 500, 1000, 1500, 2000, 3000.0])
lam_draws = post_lam.values.flatten()  # tous les tirages de lambda

# S(t) pour chaque tirage : matrice (n_draws, n_t)
S_draws = np.exp(-np.outer(lam_draws, ts))
S_med = np.median(S_draws, axis=0)
S_lo = np.percentile(S_draws, 3, axis=0)
S_hi = np.percentile(S_draws, 97, axis=0)
S_emp = np.array([(lifetimes_exp > t).mean() for t in ts])

print("=== Fonction de survie S(t) : bayesienne vs empirique ===")
print(f"{'t (h)':>8}  {'S_bayes':>8}  {'IC 94%':>14}  {'S_empir':>8}")
for j, t in enumerate(ts):
    print(f"{t:8.0f}  {S_med[j]:8.3f}  [{S_lo[j]:.3f}, {S_hi[j]:.3f}]  {S_emp[j]:8.3f}")

print(f"\nReponse a l'ingenieur : P(T > 1500 h) ≈ {S_med[list(ts).index(1500)]:.3f}",
      f"(IC 94% [{S_lo[list(ts).index(1500)]:.3f}, {S_hi[list(ts).index(1500)]:.3f}])")




=== Fonction de survie S(t) : bayesienne vs empirique ===
   t (h)   S_bayes          IC 94%   S_empir
       0     1.000  [1.000, 1.000]     1.000
     250     0.739  [0.683, 0.791]     0.767
     500     0.547  [0.467, 0.625]     0.483
    1000     0.299  [0.218, 0.391]     0.367
    1500     0.163  [0.102, 0.245]     0.117
    2000     0.089  [0.047, 0.153]     0.083
    3000     0.027  [0.010, 0.060]     0.033

Reponse a l'ingenieur : P(T > 1500 h) ≈ 0.163 (IC 94% [0.102, 0.245])








Lecture


    

  • Accord bayesien / empirique : les deux courbes se suivent. La bayesienne est parametrique et lisse (elle projette la forme exponentielle), l’empirique est en escalier (fraction brute). Leur proximite valide que l’exponentielle decrit bien ce jeu.
    • 

  • Incertitude quantifiee : l’intervalle de credibilite 94 % sur \(S(t)\) elargit la queue, la ou peu de données constrainent – c’est précisément la region fiabilite critique, et le bayesien y donne une reponse honnete (pas un chiffre unique trompeur).
    • 









4. Value-add PyMC : inferer DIRECTEMENT la forme Weibull k


L’exponentielle suppose un taux constant. La loi de Weibull generalise : \(T \sim \text{Weibull}(k, \eta)\), \(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 des 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 5 (cellule de LOO) retrouve par comparaison de modeles.




Infer-19 (Infer.NET) evite d’inferer la forme k : la vraisemblance de Weibull n’est conjuguee a aucun prior usuel, et l’EP sur la forme est notoirement instable. L’astuce d’Infer-19 est de fixer k, transformer \(U = T^k \sim \text{Exp}\) (conjugue), puis choisir k par balayage sur une grille (7 compilations de modèle).


PyMC n’a pas cette contrainte : NUTS echantillonne la forme k directement via pm.Weibull(k, eta) avec des priors sur k et eta. C’est l’apport d’un echantillonneur generique : pas de conjugaison requise. On inferere simultanement k et eta, et le postérieur sur k nous dit si le regime est a usure (\(k > 1\)), constant (\(k = 1\)) ou mortalite infantile (\(k < 1\)).




# --- Value-add : Weibull avec forme k INFEREE directement (NUTS) ---
with pm.Model() as modele_wei:
    # Priors : k positif (HalfNormal), eta positif large.
    k = pm.HalfNormal("k", sigma=5.0)
    eta = pm.HalfNormal("eta", sigma=2000.0)
    pm.Weibull("T", alpha=k, beta=eta, observed=lifetimes_wei)
    trace_wei = pm.sample(1000, tune=1000, chains=2, target_accept=0.9,
                          random_seed=42, progressbar=False)

# Log-vraisemblance pour az.compare (LOO), cf. cellule modele exponentiel.
pm.compute_log_likelihood(trace_wei, model=modele_wei)

post_k = trace_wei.posterior["k"]
post_eta = trace_wei.posterior["eta"]
print("=== Posterieurs Weibull (k ET eta inferes directement par NUTS) ===")
print(f"  k   : mediane = {float(post_k.median()):.3f}  (vrai k = {k_true})")
print(f"        IC 94% : [{float(post_k.quantile(0.03)):.3f}, {float(post_k.quantile(0.97)):.3f}]")
print(f"  eta : mediane = {float(post_eta.median()):.1f} h  (vrai eta = {eta_true})")
print(f"        IC 94% : [{float(post_eta.quantile(0.03)):.1f}, {float(post_eta.quantile(0.97)):.1f}]")
print("\nContraste : Infer-19 fixait k puis balayait 7 valeurs (EP instable sur la forme).")
print("PyMC inferere k directement -- le postérieur recouvre la vraie valeur 1.8 (usure).")















=== Posterieurs Weibull (k ET eta inferes directement par NUTS) ===
  k   : mediane = 1.964  (vrai k = 1.8)
        IC 94% : [1.623, 2.346]
  eta : mediane = 1043.9 h  (vrai eta = 1000.0)
        IC 94% : [915.0, 1189.7]

Contraste : Infer-19 fixait k puis balayait 7 valeurs (EP instable sur la forme).
PyMC inferere k directement -- le postérieur recouvre la vraie valeur 1.8 (usure).








Lecture : le MCMS generalise ce que l’EP conjuguée ne peut pas


Le postérieur sur k centrer sur la vraie valeur 1.8 : le modèle detecte le regime d’usure (k > 1) sans aucun balayage. La où Infer-19 devait fixer k puis transformer pour ramener le problème au cas conjugue (fragilite de l’EP sur la forme), PyMC echantillonne k et eta jointement – c’est précisément le type de vraisemblance non-conjuguee que le MCMC gere nativerment et que le message-passing EP fuir.


Compromis documente (G.2 honnete) : cette generalite se paie en cout (NUTS echantillonne des milliers de fois, ~secondes ici, la ou l’EP conjugue est quasi-instantanee). Sur un modèle conjugue simple (exponentiel), l’EP d’Infer.NET est plus rapide ; des qu’on veut inferer la forme ou comparer des modèles non-conjugues, NUTS devient l’outil naturel.








5. Value-add PyMC : sélection de modèle par LOO cross-validation


Infer-19 choisit k par balayage manuel (log-vraisemblance predictive sur une grille de 7 valeurs). PyMC + arviZ offrent la leave-one-out cross-validation approximée (PSIS-LOO) : on calcule le score predictive de chaque modèle (exponentiel k=1 vs Weibull k libre) sur tout le jeu, et arviZ donne le poids de chaque modèle (model weighting). C’est systématique et automatique, pas un balayage ad hoc.




# --- Selection Exp (k=1) vs Weibull (k libre) par LOO-PSIS (arviZ) ---
# Sur le jeu Weibull (vrai k=1.8) : Weibull doit gagner. Sur le jeu Exp : ex aequo.
df_loo = az.compare({"Weibull": trace_wei, "Exponentiel_sur_Weibull": trace_exp},
                    method="stacking")  # arviz>=1.0 : LOO est le critere par defaut (kwarg ic retire)
print("=== Comparaison de modeles par LOO-PSIS (arviZ) sur le jeu Weibull ===")
print(df_loo.to_string())
print("\nLe poids (weight) indique la confiance accordee a chaque modele. Weibull (k libre)")
print("doit dominer sur le jeu Weibull (vrai k=1.8) -- la selection automatique retrouve le verdict")
print("qu'Infer-19 obtenait par balayage manuel.")




=== Comparaison de modeles par LOO-PSIS (arviZ) sur le jeu Weibull ===
                         rank   elpd    p  elpd_diff  weight   se  dse  warning
Weibull                     0 -450.0  2.0        0.0    0.65  5.1  0.0    False
Exponentiel_sur_Weibull     1 -460.0  1.0      -10.0    0.35  7.6  9.6    False

Le poids (weight) indique la confiance accordee a chaque modele. Weibull (k libre)
doit dominer sur le jeu Weibull (vrai k=1.8) -- la selection automatique retrouve le verdict
qu'Infer-19 obtenait par balayage manuel.










6. 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^\star\), 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) = e^{-\lambda c_i}\).


Nous executons la comparaison complete sur le jeu exponentiel de la section 2, censuré a \(c^\star = 800\) h — environ 45 % de censures, l’ordre de grandeur courant en essais de duree de vie :


    

  • Modele naif : les durees enregistrees (les censures figees a \(c^\star\)) sont traitees comme des evenements exacts. Chaque censure est comptee comme une mort precoce : le taux estime monte, la survie estimee s’effondre — l’ingenieur conclut a tort que ses composants sont bien moins fiables qu’en realite.
    • 

  • Modele censure : \(f(t_i)\) pour les evenements, \(S(c_i)\) pour les censures. En conjugaison exponentielle, chaque censure ne change pas la forme du posterieur mais ajoute \(c_i\) au denominateur de taux : la correction dit exactement « le temps de survie non observe compte comme du temps a risque ».
    • 

  • Kaplan–Meier : reference non parametrique (Kaplan & Meier, 1958) qui ne suppose aucune loi — l’escalier empirique adapte a la censure, la ou l’empirique naif de la section 3 est fausse.
    • 



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 7) l’etendent.




# --- Censure administrative a c* = 800 h sur le jeu exponentiel (meme protocole qu'Infer-19) ---
c_star = 800.0
observed_mask = lifetimes_exp <= c_star                          # evenement observe avant la fin du test
t_recorded = np.where(observed_mask, lifetimes_exp, c_star)      # duree enregistree (censure figee a c*)
n_obs = int(observed_mask.sum()); n_cens = N - n_obs

print(f"Censure administrative a c* = {c_star:.0f} h sur N = {N} composants")
print(f"  Evenements observes : {n_obs}   (durees exactes)")
print(f"  Censures a droite   : {n_cens}   (T_i > {c_star:.0f}, seule l'information 'survit a c*' est connue)")
print(f"  Fraction censuree   : {n_cens / N:.1%}   (theorie : exp(-lambda c*) = {np.exp(-lambda_true * c_star):.1%})")




Censure administrative a c* = 800 h sur N = 60 composants
  Evenements observes : 36   (durees exactes)
  Censures a droite   : 24   (T_i > 800, seule l'information 'survit a c*' est connue)
  Fraction censuree   : 40.0%   (theorie : exp(-lambda c*) = 44.9%)








# --- Modele NAIF : toutes les durees enregistrees traitees comme des evenements exacts ---
with pm.Model() as modele_naif:
    lam_naif = pm.Gamma("lam", alpha=0.001, beta=0.001)
    # ERREUR du naif : les censures, figees a c*, comptent comme des morts a c*.
    pm.Exponential("T", lam=lam_naif, observed=t_recorded)
    trace_naif = pm.sample(1000, tune=1000, chains=2, target_accept=0.9,
                           random_seed=42, progressbar=False)

lam_naif_mean = float(trace_naif.posterior["lam"].mean())
print("=== Modele NAIF (censures traitees comme des evenements a c*) ===")
print(f"  lambda_naif   = {lam_naif_mean:.6f}   (vrai lambda = {lambda_true:.6f}, biais de {lam_naif_mean / lambda_true:.2f}x)")
print(f"  S_naif(1500)  = {np.exp(-lam_naif_mean * 1500):.3f}   (vraie survie = {np.exp(-lambda_true * 1500):.3f})")




=== Modele NAIF (censures traitees comme des evenements a c*) ===
  lambda_naif   = 0.001988   (vrai lambda = 0.001000, biais de 1.99x)
  S_naif(1500)  = 0.051   (vraie survie = 0.223)








# --- Modele CENSURE : f(t_i) pour les evenements, S(c_i) pour les censures ---
t_events = lifetimes_exp[observed_mask]     # durees exactes observees
c_cens = np.full(n_cens, c_star)            # instants de censure (tous a c*)

with pm.Model() as modele_cens:
    lam_cens = pm.Gamma("lam", alpha=0.001, beta=0.001)
    # (a) evenements observes : contribution de densite f(t_i) = lambda * exp(-lambda * t_i)
    pm.Exponential("T_obs", lam=lam_cens, observed=t_events)
    # (b) censures : contribution de survie S(c_i) = exp(-lambda * c_i), i.e. log S(c_i) = -lambda * c_i,
    #     ecrit a la main via pm.Potential (l'API empaquetee pm.Censored(lower/upper=...) fait la meme correction).
    pm.Potential("censure_logsurvie", -lam_cens * c_cens.sum())
    trace_cens = pm.sample(1000, tune=1000, chains=2, target_accept=0.9,
                           random_seed=42, progressbar=False)

lam_cens_mean = float(trace_cens.posterior["lam"].mean())
# Verification conjuguee : posterieur exact = Gamma(a + n_obs, b + somme(t_obs) + somme(c_cens))
a_exact, b_exact = 0.001 + n_obs, 0.001 + t_events.sum() + c_cens.sum()
print("=== Modele CENSURE (f(t_i) pour les evenements, S(c_i) pour les censures) ===")
print(f"  lambda_censure = {lam_cens_mean:.6f}   (vrai lambda = {lambda_true:.6f})")
print(f"  Forme fermee conjuguee = {a_exact / b_exact:.6f}  -> NUTS retrouve la conjugaison")
print(f"  S_censure(1500) = {np.exp(-lam_cens_mean * 1500):.3f}   (vraie survie = {np.exp(-lambda_true * 1500):.3f})")
print(f"\n  Le terme S(c_i) ajoute {c_cens.sum():.0f} h de temps-a-risque sans ajouter d'evenement :")
print("  chaque composant censure a 'travaille' au-dela de c* sans qu'on puisse le compter comme mort.")




=== Modele CENSURE (f(t_i) pour les evenements, S(c_i) pour les censures) ===
  lambda_censure = 0.001187   (vrai lambda = 0.001000)
  Forme fermee conjuguee = 0.001189  -> NUTS retrouve la conjugaison
  S_censure(1500) = 0.169   (vraie survie = 0.223)

  Le terme S(c_i) ajoute 19200 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) + tableau de confrontation ---
def kaplan_meier(t_rec, event):
    """Estimateur en escalier de S(t) sous censure a droite (Kaplan & Meier, 1958)."""
    order = np.argsort(t_rec)
    t_o, e_o = np.asarray(t_rec)[order], np.asarray(event)[order]
    S, out_t, out_S = 1.0, [0.0], [1.0]
    for tj, ej in zip(t_o, e_o):
        at_risk = int((t_o >= tj).sum())
        if ej:                      # seuls les evenements observes reduisent la survie
            S *= 1.0 - 1.0 / at_risk
        out_t.append(tj); out_S.append(S)
    return np.array(out_t), np.array(out_S)

ts_surv = np.array([0, 250, 500, 1000, 1500, 2000, 3000.0])
km_t, km_S = kaplan_meier(t_recorded, observed_mask.astype(int))
S_km = np.array([km_S[km_t <= t][-1] for t in ts_surv])

print("=== S(t) : verite simulee vs naive vs censuree vs Kaplan-Meier ===")
print(f"{'t (h)':>8}  {'verite':>7}  {'naif':>7}  {'censure':>8}  {'K-M':>6}")
for j, t in enumerate(ts_surv):
    print(f"{t:8.0f}  {np.exp(-lambda_true * t):7.3f}  {np.exp(-lam_naif_mean * t):7.3f}"
          f"  {np.exp(-lam_cens_mean * t):8.3f}  {S_km[j]:6.3f}")
print(f"\nEmpirique NAIF (comme en section 3) : S(1500) = {(t_recorded > 1500).mean():.3f}"
      f"  <- fausse aussi : les censures a {c_star:.0f} h y sont comptees mortes.")




=== 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.608     0.743   0.767
     500    0.607    0.370     0.552   0.483
    1000    0.368    0.137     0.305   0.400
    1500    0.223    0.051     0.169   0.400
    2000    0.135    0.019     0.093   0.400
    3000    0.050    0.003     0.028   0.400

Empirique NAIF (comme en section 3) : S(1500) = 0.000  <- fausse aussi : les censures a 800 h y sont comptees mortes.








Lecture : pourquoi le naif s’effondre et la vraisemblance restaure


    

  • Le biais naif est un biais de comptage : les 24 censures sont enregistrees comme des morts a \(c^\star = 800\) h alors que leurs vraies durees depassent 800 h. Le taux naif double presque la verite (1,99x) alors que l’estimateur censure ne la depasse que de l’alea d’echantillonnage (1,19x) ; le facteur de gonflement naif/censure est exactement \(N / n_{obs} = 60/36\) — deux fois moins de fiabilite estimee sans qu’aucun composant supplementaire soit mort.
    • 

  • Le role exact du terme \(S(c_i)\) : en conjugaison exponentielle, \(\log S(c_i) = -\lambda c_i\) est lineaire en \(\lambda\) — chaque censure ajoute \(c_i\) heures de temps a risque au denominateur du posterieur sans ajouter d’evenement au numerateur (ici +19 200 h). C’est precisement la comptabilite « a risque » que le naif oublie, et NUTS retrouve la forme fermee \(\mathrm{Gamma}(a + n_{obs},\ b + \sum t_{obs} + \sum c_i)\) chiffres en main.
    • 

  • Kaplan–Meier s’arrete ou finit l’observation : l’escalier suit la verite tant que des evenements restent observables, puis devient plat a 0,400 au-dela de \(c^\star = 800\) h — ce n’est pas une estimation de \(S(3000)\) : avec une censure administrative commune, plus aucune duree observee ne depasse \(c^\star\), 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) — c’est la complementarite non-parametrique/parametrique.
    • 

  • Pont moteur : Infer-19 mene la meme comparaison cote EP — la ou PyMC ecrit S(c_i) a la main via pm.Potential (le MCMC accepte des facteurs arbitraires), Infer.NET n’a pas d’operateur pour la forme contrainte T_i > c^\star et branche la meme vraisemblance par sa forme exacte en statistiques suffisantes. Meme protocole, meme forme fermee, deux paradigmes.
    • 









7. Exercices


Les exercices etendent le modele de duree au-dela de l’exemple execute de la section 6. 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 : \(\lambda_i = \lambda_0 \cdot \exp(\beta x_i)\) ou \(x_i\) est la covariable centree. - Indice : declarez lambda0 ~ Gamma et beta ~ Normal(0, grand), avec lam_i = lambda0 * pm.math.exp(beta * x_i) et pm.Exponential("T", lam=lam_i, observed=...). Le postérieur de beta donne l’effet de la contrainte (signe + amplitude).




# Exercice 1 a completer
# TODO etudiant : lambda_i = lambda0 * exp(beta * x_i), inferer beta (effet covariable).
print("Exercice 1 a completer : regression de temps accelere (AFT).")




Exercice 1 a completer : regression de temps accelere (AFT).










Exercice 2 – Exponentiel vs Weibull : a partir de quel N ?


Reprenez le jeu Weibull (\(k = 1.8\)) en faisant varier \(N \in \{20, 50, 100, 200\}\). Pour chaque \(N\), calculez le poids LOO du modèle Weibull. A partir de quel \(N\) le Weibull est-il prefere de facon decisive (poids \(> 0.95\)) ? - Indice : bouclez sur les valeurs de N, re-echantillonnez les deux modèles, appelez az.compare. C’est la generalisation de la question 3 d’Infer-19 (qui utilisait un facteur de Bayes manuel).




# Exercice 2 a completer
# TODO etudiant : boucler sur N, calculer le poids LOO du Weibull, trouver le seuil N.
print("Exercice 2 a completer : a partir de quel N le Weibull est-il decisiivement prefere ?")




Exercice 2 a completer : a partir de quel N le Weibull est-il decisiivement prefere ?










Exercice 3 — Sensibilite au taux de censure


L’exemple de la section 6 fixe \(c^\star = 800\) h (~45 % de censures). Refaites tourner le trio naif / censure / Kaplan–Meier pour \(c^\star \in \{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 c_star change — encapsulez la construction du jeu et les deux modeles dans une fonction comparer_censure(c_star) et bouclez. A \(c^\star = 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^\star\).
    • 

  • Etape 2 : rapporter fraction censurée et biais relatif pour chaque \(c^\star\).
    • 





# Exercice 3 a completer
# TODO etudiant : boucler sur c_star 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.
print("Exercice 3 a completer : sensibilite du biais naif au taux de censure.")




Exercice 3 a completer : sensibilite du biais naif au taux de censure.












Ponts avec la serie






Direction
Lien
Relation






<-> Infer-19
Infer-19 – Analyse de survie
Jumeau .NET (Infer.NET, EP). Ce notebook ajoute : inference directe de la forme k (NUTS) + sélection LOO automatique (arviZ).




<-> PyMC-18
PyMC-18 – Change-Point
Le change-Point est une rupture du taux ; la survie modelise un delai unique. Combiner les deux = survie a rupture (taux qui change a un instant latent).




<-> PyMC-14
PyMC-14 – Sequences (HMM)
Complete la famille temporelle : HMM (discret recurrent), Kalman (continu recurrent), survie (continu, duree unique).










Conclusion


L’analyse de survie complete la famille des modèles temporels : la quantite centrale est la fonction de survie \(S(t) = P(T > t)\), dans la queue, pas au centre.


Dualite des moteurs (jumeaux .NET / Python) : - Infer.NET (EP) : exponentiel conjugue exact (forme fermee \((B/(B+t))^A\)), mais evite d’inferer la forme Weibull (fragilite de l’EP sur la forme → transformee + balayage de grille). - PyMC (NUTS) : inferere k directement (vraie valeur recouvre), sélection de modèle automatique par LOO. Plus general, au prix d’un cout MCMC.


La lecon : sur un modèle conjugue simple (exponentiel), l’EP d’Infer.NET est inegalable en rapidite et exactitude fermee. Des qu’on veut inferer une forme non-conjuguee (Weibull k) ou comparer des modèles, le MCMC generique de PyMC devient l’outil naturel – c’est la complementarite des deux moteurs sur la même famille de problemes.


Pour aller plus loin. La censure a gauche/par intervalle, les modèles a risques proportionnels (Cox, 1972, JRSS Series B 34(2):187-220), les modèles a risques concurrents (competing risks) et les modèles de fragilite partagee (frailty, analogue hiérarchique de PyMC-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 (section 4).
    • 

  • 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-220 — modèle à risques proportionnels, l’extension de référence pour la régression sur la survie (au-delà du présent 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 / censure.
    • 

  • Salvatier J., Wiecki T. V., Fonnesbeck C. (2016), « Probabilistic programming in Python using PyMC3 », PeerJ Computer Science 2:e55 — moteur NUTS utilisé ici.
    • 

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









 Retour au sommet