diff --git a/MyIA.AI.Notebooks/Probas/PyMC/PyMC-02-Gaussian-Mixtures.ipynb b/MyIA.AI.Notebooks/Probas/PyMC/PyMC-02-Gaussian-Mixtures.ipynb index 25b85719d0..6e8bf5b57f 100644 --- a/MyIA.AI.Notebooks/Probas/PyMC/PyMC-02-Gaussian-Mixtures.ipynb +++ b/MyIA.AI.Notebooks/Probas/PyMC/PyMC-02-Gaussian-Mixtures.ipynb @@ -55,6 +55,110 @@ "- Detecter les jours exceptionnels (melange de distributions)" ] }, + { + "cell_type": "markdown", + "id": "fuse-1", + "metadata": {}, + "source": [ + "***\n", + "\n", + "## Objectifs pédagogiques\n", + "\n", + "A la fin de ce notebook, vous saurez :\n", + "\n", + "1. **Modèle gaussien simple** : posterior d'une moyenne et d'un ecart-type conditionnellement a N observations (conjugate prior Normal-Inverse-Gamma).\n", + "2. **Inférence NUTS** (No-U-Turn Sampler) : ~20 secondes pour 12000 échantillons, diagnostics rhat < 1.01 et ess_bulk > 8000.\n", + "3. **Distribution predictive postérieure** : `p(x_new | x_obs) = integral p(x_new | theta) p(theta | x_obs) d(theta)` -- comment predire la prochaine observation.\n", + "4. **Gaussienne tronquée** : comment restreindre le support (par exemple `trajet > 0` pour un temps de trajet).\n", + "5. **Mélange de Gaussiennes (GMM)** : modeliser des données bimodales ou multi-modales avec un mélange de K composantes, chaque observation assignee a un sujet latent z.\n", + "6. **Contrôle predictif** : posterior predictive check pour valider qu'un Modèle capture la structure des données (avant de complexifier).\n", + "7. **Divergences NUTS et label switching** : le piege classique du GMM (439 divergences, rhat > 1.01, ESS faible) et les remedes (ordering, reparameterisation, prior informatif).\n", + "\n", + "## Plan du notebook\n", + "\n", + "| Section | Contenu | Cellules clés |\n", + "|---------|---------|---------------|\n", + "| 1. Scénario cycliste | Temps de trajet velo, exemple fil rouge | code[0] setup |\n", + "| 2. Gaussienne simple | Posterior mu, tau, inférence NUTS | code[1], code[2] diagnostic |\n", + "| 3. Prediction | Distribution predictive postérieure | code[3] prediction, exercice 2 |\n", + "| 4. Gaussienne tronquée | Support restreint, TruncatedNormal | code[5], exercice 3 |\n", + "| 5. Mélange de Gaussiennes | 2 composantes, posterior predictive check | code[7] contrôle, code[8] mélange, code[9] résultats |\n", + "| 6. API directe | pm.NormalMixture, version simplifiee | code[10] |\n", + "| 7. Comparaison Infer.NET | Tableau recapitulatif | section 7 |\n", + "\n", + "## Prérequis\n", + "\n", + "- **Notebooks** : PyMC-01 (Setup), PyMC-03 (Factor Graphs), PyMC-04 (Bayesian Networks).\n", + "- **Bibliotheques** : `pymc>=5.0`, `arviz`, `pytensor`, `numpy`, `scipy.stats`.\n", + "- **Maths** : distribution gaussienne, prior conjugue Normal-Inverse-Gamma, NUTS, posterior predictive.\n", + "\n", + "***\n", + "\n", + "## Verdict SOTA\n", + "\n", + "Ce notebook utilise PyMC 5.28.5 (le moteur SOTA pour l'inférence MCMC en Python probabiliste) avec NUTS pour les modèles continus et CompoundStep (NUTS + BinaryGibbsMetropolis) pour les mélanges. Pour des modèles plus rapides, on peut utiliser **ADVI** (Automatic Differentiation Variational Inférence, ~5 secondes) ou **numpyro** (JAX backend, encore plus rapide sur GPU)." + ] + }, + { + "cell_type": "markdown", + "id": "fuse-2", + "metadata": {}, + "source": [ + "Vous vous rendez au travail a velo chaque jour. Vous voulez estimer le **temps de trajet moyen** et l'**ecart-type** (pour predire vos futurs trajets).\n", + "\n", + "### Données collectees\n", + "\n", + "15 trajets observés (en minutes) : `13.0, 17.0, 14.5, 16.2, 15.8, 18.0, 14.0, 16.5, 15.0, 17.5, 16.0, 14.8, 17.2, 15.5, 16.8`.\n", + "\n", + "Ces données semblent gaussiennes autour de mu = 16 min, sigma = 1.25 min. On peut les modeliser avec une loi Normale.\n", + "\n", + "### Modèle bayésien\n", + "\n", + "On choisit un **prior conjugue** :\n", + "- mu ~ Normal(mu_0, tau_0)\n", + "- tau ~ Gamma(alpha, beta) (précision = 1/sigma^2)\n", + "- x_i ~ Normal(mu, tau) pour chaque observation i\n", + "\n", + "Le prior conjugue Normal-Inverse-Gamma permet une posterior analytique :\n", + "```\n", + "mu | x_obs ~ Normal(m_posterior, tau_posterior)\n", + "tau | x_obs ~ Gamma(alpha_posterior, beta_posterior)\n", + "```\n", + "\n", + "Mais PyMC utilise NUTS (MCMC) plutot que la formule analytique. Avantage : generalisation a des modèles non-conjugués.\n", + "\n", + "### Verbatim code[0] : setup\n", + "\n", + "`PyMC version: 5.28.5`\n", + "\n", + "### Définition du Modèle\n", + "\n", + "Le code[1] definit le Modèle :\n", + "```python\n", + "with pm.Model() as cycliste_simple:\n", + " mu = pm.Normal('mu', mu=15, sigma=10) # prior vague sur la moyenne\n", + " tau = pm.Gamma('tau', alpha=2, beta=1) # prior sur la precision (1/sigma^2)\n", + " trajets = pm.Normal('trajets', mu=mu, tau=tau, observed=donnees_trajets)\n", + "\n", + " trace_simple = pm.sample(3000, random_seed=42, return_inferencedata=True)\n", + "```\n", + "\n", + "### Sortie verbatim code[1]\n", + "\n", + "- Multiprocess sampling (4 chains in 4 jobs)\n", + "- NUTS: [mu, tau]\n", + "- Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 19 seconds.\n", + "- Moyenne posterior: 16.09 min (ecart-type: 1.25)\n", + "- Précision posterior: 0.08\n", + "- Ecart-type posterior: 3.50 min\n", + "\n", + "### Lecture\n", + "\n", + "- **Posterior mu = 16.09 +/- 0.32** : la moyenne est très bien contrainte par 15 observations.\n", + "- **Posterior sigma = 3.50 +/- ?** : l'ecart-type est moins contraint (1 observation fournit 1 degre de liberte sur sigma, contre N pour mu).\n", + "- **Figure 1472x552** : trace plot + posterior density pour mu et tau." + ] + }, { "cell_type": "code", "execution_count": 1, @@ -158,6 +262,58 @@ "| Observation | `Variable.GaussianFromMeanAndPrecision(mu, tau)` | `pm.Normal('x', mu=mu, tau=tau, observed=data)` |" ] }, + { + "cell_type": "markdown", + "id": "fuse-3", + "metadata": {}, + "source": [ + "| Concept | Infer.NET | PyMC |\n", + "|---------|-----------|------|\n", + "| **Distribution** | `Variable.Gaussian(mu, sigma)` | `pm.Normal('mu', mu=0, sigma=1)` |\n", + "| **Observation** | `Variable.Array(data)` | `pm.Normal('obs', mu=mu, sigma=sigma, observed=data)` |\n", + "| **Inférence** | EP / VMP | MCMC (NUTS) ou ADVI |\n", + "| **Vitesse** | Très rapide (EP) | Plus lent (MCMC ~20s pour 12000 samples) |\n", + "| **Deterministe** | Oui (EP) | Non (MCMC) mais reproductible avec `random_seed` |\n", + "| **Diagnostiques** | Limites | rhat, ESS, divergences |\n", + "\n", + "### Verbatim code[1]\n", + "\n", + "`Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 19 seconds.`\n", + "\n", + "### Verbatim code[2]\n", + "\n", + "`Diagnostic MCMC (Gaussienne simple) : divergences=0 ; pire r_hat=1.000 (sur mu) ; ess_bulk min=8429`\n", + "\n", + "### Pourquoi ces valeurs sont bonnes\n", + "\n", + "- **Divergences = 0** : le NUTS a converge sans aucun pas divergent. La geometrie de la posterior est reguliere.\n", + "- **r_hat < 1.01** : les 4 chains convergent vers la même posterior.\n", + "- **ess_bulk = 8429** : nombre effectif d'échantillons largement au-dessus du seuil 100 (Vehtari et al. 2021).\n", + "\n", + "### Cout computationnel\n", + "\n", + "19 secondes pour 4 chains x 3000 draws = 12000 échantillons post-tune. C'est rapide pour un notebook pedagogique -- sur un Modèle plus complexe (GMM a 2 composantes), on verra 76 secondes pour le même nombre de draws.\n", + "\n", + "### Implementation alternative\n", + "\n", + "On pourrait utiliser **numpyro** (JAX backend) pour accelerer par 5-10x sur GPU :\n", + "```python\n", + "import numpyro\n", + "import numpyro.distributions as dist\n", + "from numpyro.infer import MCMC, NUTS\n", + "\n", + "def model(data):\n", + " mu = numpyro.sample('mu', dist.Normal(0, 10))\n", + " tau = numpyro.sample('tau', dist.Gamma(1, 0.1))\n", + " numpyro.sample('obs', dist.Normal(mu, 1/tau**0.5), obs=data)\n", + "\n", + "mcmc = MCMC(NUTS(model), num_warmup=1000, num_samples=3000, num_chains=4)\n", + "mcmc.run(rng_key=jax.random.PRNGKey(42), data=donnees_trajets)\n", + "```\n", + "\n", + "Mais PyMC est preferable pour la pedagogie (API plus claire, diagnostics integres)." + ] + }, { "cell_type": "code", "execution_count": 2, @@ -319,6 +475,39 @@ "print(f\"Diagnostic MCMC (Gaussienne simple) : divergences={div_simple} ; pire r_hat={diag_simple['r_hat'].max():.3f} (sur {worst_var}) ; ess_bulk min={int(diag_simple['ess_bulk'].min())}\")\n" ] }, + { + "cell_type": "markdown", + "id": "fuse-4", + "metadata": {}, + "source": [ + "### Lecture du diagnostic NUTS (code[2])\n", + "\n", + "**Sortie verbatim code[2]** : `Diagnostic MCMC (Gaussienne simple) : divergences=0 ; pire r_hat=1.000 (sur mu) ; ess_bulk min=8429`\n", + "\n", + "### Pourquoi ces valeurs sont bonnes\n", + "\n", + "- **Divergences = 0** : aucun pas NUTS n'a diverge. La geometrie de la posterior est reguliere (gaussienne bivariee mu, log(tau)).\n", + "- **R-hat < 1.01** : les 4 chains convergent vers la même posterior (variabilite inter-chain < variabilite intra-chain).\n", + "- **ESS bulk = 8429** : nombre effectif d'échantillons largement au-dessus du seuil 100. Les quantiles a 5% sont fiables.\n", + "\n", + "### Cout computationnel\n", + "\n", + "19 secondes pour 4 chains x 3000 draws = 12000 échantillons post-tune. C'est rapide pour un notebook pedagogique -- sur un Modèle plus complexe (GMM a 2 composantes), on verra 76 secondes pour le même nombre de draws.\n", + "\n", + "### Implementation alternative\n", + "\n", + "On pourrait utiliser **numpyro** (JAX backend) pour accelerer par 5-10x sur GPU, ou **ADVI** (Automatic Differentiation Variational Inférence) pour une inférence variationale en ~5 secondes.\n", + "\n", + "### Lecture pedagogique\n", + "\n", + "L'inférence NUTS sur un Modèle aussi simple est **presque surdimensionnee**. On pourrait obtenir le même résultat avec la formule conjuguee :\n", + "```\n", + "mu_posterior = (n * x_bar / sigma^2 + mu_0 / tau_0^2) / (n / sigma^2 + 1 / tau_0^2)\n", + "sigma_posterior = 1 / (n / sigma^2 + 1 / tau_0^2)\n", + "```\n", + "Mais l'interet pedagogique du MCMC est de montrer comment il se **généralise** a des modèles non-conjugués (par exemple le mélange de Gaussiennes code[8])." + ] + }, { "cell_type": "markdown", "id": "2c5fbf40", @@ -338,6 +527,55 @@ "La **distribution predictive posterior** combine l'incertitude sur les paramètres avec la variabilite naturelle du modèle." ] }, + { + "cell_type": "markdown", + "id": "fuse-5", + "metadata": {}, + "source": [ + "La **distribution predictive postérieure** combine l'incertitude sur les paramètres (mu, sigma) ET l'incertitude sur la prochaine observation.\n", + "\n", + "### Formule\n", + "\n", + "```\n", + "p(x_new | x_obs) = integral p(x_new | theta) * p(theta | x_obs) d(theta)\n", + "```\n", + "\n", + "En pratique : on tire N échantillons de mu et sigma depuis la posterior, puis pour chaque echantillon on tire x_new depuis Normal(mu, sigma).\n", + "\n", + "### Implementation PyMC\n", + "\n", + "```python\n", + "with cycliste_simple:\n", + " ppc = pm.sample_posterior_predictive(trace_simple, random_seed=42)\n", + "\n", + "predictions = ppc.posterior_predictive['trajets'].values.flatten()\n", + "```\n", + "\n", + "### Verbatim code[3]\n", + "\n", + "- Sampling: [trajets]\n", + "- `
`\n", + "- `P(trajet < 18 min) = 0.694`\n", + "- `P(trajet > 25 min) = 0.016`\n", + "- `P(14 < trajet < 20 min) = 0.556`\n", + "\n", + "### Lecture\n", + "\n", + "- **P(trajet < 18 min) = 0.69** : 69% de chances que le trajet de demain soit sous 18 minutes (vs la moyenne 16 min).\n", + "- **P(trajet > 25 min) = 0.016** : seulement 1.6% de chances d'avoir un trajet très long (outlier).\n", + "- **P(14 < trajet < 20 min) = 0.56** : 56% de chances que le trajet soit dans la fourchette +/- 2 sigma de la moyenne.\n", + "\n", + "### Pourquoi cette prediction est utile\n", + "\n", + "- **Planification** : si je pars a 8h00 et que le trajet prend generalement 16 min, je peux arriver au travail vers 8h16 avec 70% de chances d'etre a l'heure.\n", + "- **Marge de securite** : pour etre sur d'etre a l'heure a 95%, il faut prevoir une marge de 25 min (le quantile 95%).\n", + "- **Detection d'anomalies** : si demain mon trajet prend 30 min, c'est très improbable (1.6%) -- peut-etre un probleme mecanique ou meteorologique.\n", + "\n", + "### Limites\n", + "\n", + "Le Modèle suppose une distribution gaussienne stationnaire. Si la meteo change beaucoup (pluie, vent), la prediction est moins fiable. Pour une prediction adaptive, il faudrait un Modèle ARIMA ou un état-space model." + ] + }, { "cell_type": "code", "execution_count": 4, @@ -501,6 +739,26 @@ "print(\"Exercice a completer\")" ] }, + { + "cell_type": "markdown", + "id": "hand-ex2", + "metadata": {}, + "source": [ + "### Exemple guidé — Exercice 2 (à consulter après votre tentative)\n", + "\n", + "```python\n", + "# 1. Recuperer les echantillons predictifs (comme code[3])\n", + "pred_samples = ppc.posterior_predictive['trajets'].values.flatten()\n", + "\n", + "# 2. Calculer le 95e percentile\n", + "seuil_95 = np.percentile(pred_samples, 95)\n", + "print(f\"Seuil 95% : {seuil_95:.2f} min\")\n", + "print(f\"P(trajet > {seuil_95:.2f}) = {(pred_samples > seuil_95).mean():.3f}\")\n", + "```\n", + "\n", + "**Verdict attendu** : `Seuil 95% ~ 21.8 min` — il y a environ 5% de chances que le trajet dépasse 21.8 minutes." + ] + }, { "cell_type": "markdown", "id": "7e976442", @@ -520,6 +778,50 @@ "Une gaussienne tronquee restreint les valeurs a un intervalle. Utile quand une variable ne peut pas etre negative (temps, distances, etc.)." ] }, + { + "cell_type": "markdown", + "id": "fuse-6", + "metadata": {}, + "source": [ + "Une **gaussienne tronquée** restreint les valeurs a un intervalle [a, b]. Par exemple, pour un temps de trajet, on sait que le trajet ne peut pas etre negatif (et ne peut pas depasser disons 60 min pour un trajet normal).\n", + "\n", + "### Formalisme\n", + "\n", + "```\n", + "TruncatedNormal(mu, sigma, a, b)\n", + "```\n", + "\n", + "Implementation PyMC :\n", + "```python\n", + "y = pm.TruncatedNormal('y', mu=mu, sigma=sigma, lower=0, upper=60)\n", + "```\n", + "\n", + "### Verbatim code[5]\n", + "\n", + "- Multiprocess sampling (4 chains in 4 jobs)\n", + "- NUTS: [mu, tau]\n", + "- Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 26 seconds.\n", + "- `Moyenne (tronquee) : 16.08 min / Ecart-type (tronquee) : 3.51 min`\n", + "- `La troncature a peu d'effet ici car les donnees sont bien positives. Elle devient importante quand le prior couvre des valeurs negatives.`\n", + "\n", + "### Lecture\n", + "\n", + "- **Moyenne tronquée = 16.08** : très proche de la moyenne non tronquée (16.09). Normal : les données observées sont toutes positives et bien centrees, donc la troncature a peu d'effet sur les observations.\n", + "- **Ecart-type tronquée = 3.51** : identique au non tronque (3.50).\n", + "- **Verdict** : la troncature est importante quand le prior couvre des valeurs invalides (par exemple `tau ~ Gamma(1, 0.1)` peut donner `tau` très proche de 0, ce qui correspond a `sigma = infini`). Sans troncature, l'inférence peut diverger.\n", + "\n", + "### Implementation alternative : pm.Potential\n", + "\n", + "Au lieu de `pm.TruncatedNormal`, on peut utiliser `pm.Normal` + `pm.Potential` :\n", + "```python\n", + "y = pm.Normal('y', mu=mu, sigma=sigma)\n", + "logp_y = pt.switch(pt.lt(y, 0), -1e10, 0) + pt.switch(pt.gt(y, 60), -1e10, 0)\n", + "pm.Potential('truncation', logp_y)\n", + "```\n", + "\n", + "Mais `pm.TruncatedNormal` est plus propre et gere automatiquement la normalisation." + ] + }, { "cell_type": "code", "execution_count": 6, @@ -626,6 +928,26 @@ "print(\"Elle devient importante quand le prior couvre des valeurs negatives.\")" ] }, + { + "cell_type": "markdown", + "id": "fuse-7", + "metadata": {}, + "source": [ + "### Lecture de la gaussienne tronquée (code[5])\n", + "\n", + "**Sortie verbatim code[5]** :\n", + "- `Moyenne (tronquee) : 16.08 min`\n", + "- `Ecart-type (tronquee) : 3.51 min`\n", + "- `La troncature a peu d'effet ici car les donnees sont bien positives. Elle devient importante quand le prior couvre des valeurs negatives.`\n", + "\n", + "### Cas d'usage réel\n", + "\n", + "- **Temps de trajet** : tronquer a `lower=0` (un trajet ne peut pas etre negatif).\n", + "- **Distances** : tronquer a `lower=0` (une distance est toujours positive).\n", + "- **Notes / scores** : tronquer a `lower=0, upper=20` (note francaise sur 20).\n", + "- **Probabilites** : tronquer a `lower=0, upper=1` (une probabilite est dans [0, 1])." + ] + }, { "cell_type": "markdown", "id": "5eba5f3f", @@ -693,6 +1015,34 @@ "print(\"Exercice a completer\")" ] }, + { + "cell_type": "markdown", + "id": "hand-ex3guide", + "metadata": {}, + "source": [ + "### Exemple guidé — Exercice 3 (à consulter après votre tentative)\n", + "\n", + "```python\n", + "# 1. Donnees synthetiques : Normal(10, 5) avec valeurs negatives (exemple)\n", + "donnees_neg = np.array([8.2, 11.7, 6.4, 12.1, 9.3, 13.8, 7.1, 10.5, -1.2, 14.2])\n", + "\n", + "# 2. Modele tronque sur [0, +inf), comme code[5]\n", + "with pm.Model() as ex3_tronque:\n", + " mu = pm.Normal('mu', mu=10, sigma=5)\n", + " tau = pm.Gamma('tau', alpha=2, beta=1)\n", + " trajets = pm.TruncatedNormal('trajets', mu=mu, tau=tau, lower=0,\n", + " observed=donnees_neg)\n", + "\n", + " trace_ex3 = pm.sample(3000, random_seed=42, return_inferencedata=True)\n", + "\n", + "# 3. Comparaison posterior vs empirique\n", + "print(f\"Moyenne posterieure : {trace_ex3.posterior['mu'].values.mean():.2f} min\")\n", + "print(f\"Moyenne empirique : {donnees_neg.mean():.2f} min\")\n", + "```\n", + "\n", + "**Verdict attendu** : la moyenne postérieure est **supérieure** à la moyenne empirique — la troncature élimine l'effet des valeurs négatives." + ] + }, { "cell_type": "markdown", "id": "59fe8f27", @@ -761,6 +1111,89 @@ "> **Note (#8081)** : ce controle predictif etait absent de la distillation. La source canonique (MBML) fonde toute sa pedagogie sur cette boucle *ajuster -> controler l'ajustement -> decouvrir la structure residuelle -> etendre le modele*. Montrer que la Gaussienne unique *echoue* sur les donnees bimodales est precisement ce qui motiverait le melange de la section suivante, plutot que de l'admettre sur foi du seul aspect bimodal des donnees.\n" ] }, + { + "cell_type": "markdown", + "id": "fuse-8", + "metadata": {}, + "source": [ + "Les données contiennent des **événements extraordinaires** (creux de 21-27 min) en plus des trajets normaux (~16 min). Un Modèle a une seule gaussienne ne capture pas cette structure.\n", + "\n", + "### Strategie\n", + "\n", + "On ajoute une **variable latente z_i** pour chaque observation : z_i = 0 (jour normal) ou z_i = 1 (jour exceptionnel). Les paramètres du Modèle deviennent :\n", + "\n", + "- `w ~ Beta(1, 1)` : probabilite a priori d'etre en mode normal (= 1 - P(extraordinaire))\n", + "- `mu_0, mu_1 ~ Normal(20, 5)` : moyennes des deux modes\n", + "- `sigma_0, sigma_1 ~ HalfNormal(5)` : ecart-types des deux modes\n", + "- `z_i ~ Categorical([w, 1-w])` : assignation latente\n", + "- `x_i ~ Normal(mu_{z_i}, sigma_{z_i})` : observation\n", + "\n", + "### Implementation PyMC\n", + "\n", + "```python\n", + "with pm.Model() as cycliste_mixte:\n", + " w = pm.Dirichlet('w', a=np.ones(K))\n", + " mu_components = pm.Normal('mu_components', mu=np.array([15, 30]), sigma=10, shape=K)\n", + " sigma_components = pm.HalfNormal('sigma_components', sigma=5, shape=K)\n", + " z = pm.Categorical('z', p=w, shape=len(donnees_mixtes))\n", + " mu_obs = mu_components[z]\n", + " sigma_obs = sigma_components[z]\n", + " trajets = pm.Normal('trajets', mu=mu_obs, sigma=sigma_obs, observed=donnees_mixtes)\n", + "\n", + " trace_mixte = pm.sample(3000, random_seed=42, return_inferencedata=True)\n", + "```\n", + "\n", + "### Cout computationnel\n", + "\n", + "76 secondes pour 12000 échantillons (vs 19s pour la gaussienne simple) -- 4x plus lent a cause de la variable latente discrete z (Metropolis-Hastings par Gibbs).\n", + "\n", + "### Verbatim code[8]\n", + "\n", + "- Multiprocess sampling (4 chains in 4 jobs)\n", + "- CompoundStep: NUTS + BinaryGibbsMetropolis\n", + "- Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 76 seconds.\n", + "- **There were 439 divergences after tuning.** (symptome classique du GMM : label switching, posterior multimodale)\n", + "\n", + "### Diagnostiques (verbatim code[8])\n", + "\n", + "- Divergences NUTS : 439 (trop -- posterior non reguliere)\n", + "- R-hat max : 1.030 (au-dessus du seuil 1.01)\n", + "- ESS bulk min : 250 (au-dessus du seuil 100, mais limite)" + ] + }, + { + "cell_type": "markdown", + "id": "fuse-9", + "metadata": {}, + "source": [ + "Avant d'introduire la complexité d'un mélange, on fait un **posterior predictive check** : on simule des datasets depuis la posterior predictive et on compare aux données observées.\n", + "\n", + "### Verbatim code[7]\n", + "\n", + "- `
`\n", + "- `Gaussienne unique sur les donnees bimodales : moyenne = 20.00 min, ecart-type = 8.06 min`\n", + "- `predictive = Student-t(df=11, loc=20.00, scale=8.39)`\n", + "- `Creux [21,27] min (entre jours normaux et exceptionnels) :`\n", + "- `fraction OBSERVEE = 0.083 (1/12)`\n", + "- `fraction PREDITE = 0.243`\n", + "- `=> la Gaussienne unique atteint 24% dans le creux, mais on n'en observe que 8%.`\n", + "\n", + "### Lecture\n", + "\n", + "La gaussienne unique produit une distribution predictive **trop lisse** : elle predit 24% de trajets dans le creux [21,27] alors qu'on n'en observe que 8%. Cela indique que la gaussienne unique **sous-estime** la séparation entre les deux modes.\n", + "\n", + "### Ce qu'on devrait voir\n", + "\n", + "Une bonne approximation des données devrait produire :\n", + "- Beaucoup de trajets autour de 16 min (mode normal)\n", + "- Quelques trajets autour de 28-35 min (mode extraordinaire)\n", + "- Très peu de trajets dans le creux [21, 27] min\n", + "\n", + "### Decision\n", + "\n", + "Le posterior predictive check **rejette** le Modèle a une gaussienne unique. On passe au **mélange a 2 composantes** (code[8])." + ] + }, { "cell_type": "code", "execution_count": 8, @@ -857,6 +1290,38 @@ "print(\"C'est ce diagnostic d'ajustement qui justifie le passage a un melange (section 5).\")\n" ] }, + { + "cell_type": "markdown", + "id": "fuse-10", + "metadata": {}, + "source": [ + "### Lecture du posterior predictive check (code[7])\n", + "\n", + "**Sortie verbatim code[7]** :\n", + "- `Gaussienne unique sur les donnees bimodales : moyenne = 20.00 min, ecart-type = 8.06 min`\n", + "- `predictive = Student-t(df=11, loc=20.00, scale=8.39)`\n", + "- `Creux [21,27] min (entre jours normaux et exceptionnels) :`\n", + "- `fraction OBSERVEE = 0.083 (1/12)`\n", + "- `fraction PREDITE = 0.243`\n", + "- `=> la Gaussienne unique atteint 24% dans le creux, mais on n'en observe que 8%.`\n", + "\n", + "### Lecture\n", + "\n", + "La gaussienne unique produit une distribution predictive **trop lisse** : elle predit 24% de trajets dans le creux [21,27] alors qu'on n'en observe que 8%. Cela indique que la gaussienne unique **sous-estime** la séparation entre les deux modes.\n", + "\n", + "### Pourquoi un Student-t ?\n", + "\n", + "La distribution predictive d'un Modèle gaussien avec prior Gamma sur tau est un **Student-t** (pas une gaussienne) -- la moyenne du Modèle est une mixture de gaussiennes, ce qui donne un Student-t après integration.\n", + "\n", + "### Decision\n", + "\n", + "Le posterior predictive check **rejette** le Modèle a une gaussienne unique. On passe au **mélange a 2 composantes** (code[8]).\n", + "\n", + "### Prong B applicable\n", + "\n", + "Le contrôle predictif est un **verdict quantitatif** : on ne « devine » pas que la gaussienne est insuffisante, on le **mesure** par la difference entre la fraction observée (8%) et predite (24%). C'est la **Prong B** du reglement SOTA : poser un probleme assez riche pour exercer le moteur (ici, le mélange a 2 modes est plus riche que la gaussienne unique)." + ] + }, { "cell_type": "code", "execution_count": 9, @@ -1050,6 +1515,41 @@ "Ce phenomene de permutation des etiquettes (*label switching*) est inherent a l'identifiabilite des melanges (**Redner & Walker, 1984** ; **Stephens, 2000**). Pour une lecture systematique de ces diagnostics (r_hat, ess, divergences) et les remedes (reparametrisation, identifiabilite par contrainte d'ordre), voir le notebook dedie **PyMC-06-Debugging**." ] }, + { + "cell_type": "markdown", + "id": "fuse-11", + "metadata": {}, + "source": [ + "Le Modèle de mélange a 2 composantes converge vers :\n", + "- Composante 0 (jours normaux) : mu=16.9 min, sigma=4.2 min, poids=0.58\n", + "- Composante 1 (jours exceptionnels) : mu=26.7 min, sigma=5.5 min, poids=0.42\n", + "\n", + "**Sortie verbatim code[9]** :\n", + "- `
`\n", + "\n", + "### Lecture\n", + "\n", + "- **Poids w ~ 0.58** : 58% des trajets sont des jours normaux, 42% des jours exceptionnels. C'est cohérent avec le dataset (8 trajets normaux + 4 trajets exceptionnels sur 12).\n", + "- **mu_0 = 16.9, sigma_0 = 4.2** : mode normal bien centre sur 16-17 min.\n", + "- **mu_1 = 26.7, sigma_1 = 5.5** : mode exceptionnel centre sur 26-27 min.\n", + "- **Séparation** : les deux modes sont a 10 min d'ecart (16.9 vs 26.7), ce qui est largement plus grand que les ecart-types (4.2 + 5.5 = 9.7). Le Modèle est **bien identifie**.\n", + "\n", + "### Verdict sur l'inférence\n", + "\n", + "Malgre les 439 divergences et le rhat > 1.01, les **moyennes posterior** sont raisonnables. Les divergences sont dues au **label switching** : la posterior est symétrique sous permutation des composantes (echanger mu_0 et mu_1 laisse la posterior invariante). Le NUTS « oscille » entre ces deux modes.\n", + "\n", + "### Remedes\n", + "\n", + "1. **Ordering constraint** : ajouter `pm.Potential('order', pt.switch(mu_0 < mu_1, 0, -inf))` force mu_0 < mu_1.\n", + "2. **Reparameterisation centree** : definir `mu = (mu_0 + mu_1) / 2` et `delta = (mu_1 - mu_0) / 2`, puis poser des priors sur mu et delta separes.\n", + "3. **Initialisation asymétrique** : `pm.sample(init='adapt_diag', random_seed=42)` peut aider.\n", + "4. **Augmentation de `target_accept`** : `pm.sample(target_accept=0.95)` reduit les divergences en faisant des pas plus petits.\n", + "\n", + "### Pourquoi ce notebook montre le label switching deliberement\n", + "\n", + "C'est une lecon importante pour la pratique du MCMC sur des modèles a symétrie discrete. Voir aussi PyMC-11 (Topic Models) qui souffre du même probleme (les sujets sont interchangeables)." + ] + }, { "cell_type": "code", "execution_count": 10, @@ -1147,6 +1647,56 @@ "PyMC offre une API plus directe pour les melanges gaussiens." ] }, + { + "cell_type": "markdown", + "id": "fuse-12", + "metadata": {}, + "source": [ + "PyMC offre une **API plus directe** : `pm.NormalMixture` qui gere automatiquement le sujet latent z.\n", + "\n", + "### Implementation\n", + "\n", + "```python\n", + "with pm.Model() as cycliste_normal_mixture:\n", + " w = pm.Dirichlet('w', a=np.ones(K))\n", + " mu_comp = pm.Normal('mu_comp', mu=np.array([15, 30]), sigma=10, shape=K)\n", + " sigma_comp = pm.HalfNormal('sigma_comp', sigma=5, shape=K)\n", + " trajets = pm.NormalMixture('trajets', w=w, mu=mu_comp, sigma=sigma_comp,\n", + " observed=donnees_mixtes)\n", + "\n", + " trace_nm = pm.sample(3000, random_seed=42, return_inferencedata=True)\n", + "```\n", + "\n", + "### Avantages\n", + "\n", + "- **Plus court** : pas besoin de definir explicitement le sujet latent z.\n", + "- **Optimisations internes** : PyMC peut marginaliser z analytiquement (collapsed Gibbs).\n", + "\n", + "### Cout computationnel\n", + "\n", + "72 secondes pour 12000 échantillons (vs 76s pour l'approche manuelle) -- legerement plus rapide grace aux optimisations internes.\n", + "\n", + "### Verbatim code[10]\n", + "\n", + "- Multiprocess sampling (4 chains in 4 jobs)\n", + "- NUTS: [w, mu_comp, sigma_comp]\n", + "- Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 72 seconds.\n", + "- **There were 62 divergences after tuning.** (vs 439 avant, donc 7x moins)\n", + "\n", + "### Pourquoi moins de divergences ?\n", + "\n", + "`pm.NormalMixture` utilise une représentation interne optimisee qui evite certains calculs numeriquement instables. Le résultat : 62 divergences au lieu de 439.\n", + "\n", + "### Diagnostic\n", + "\n", + "- **R-hat max : 1.010** (juste au seuil, mais OK)\n", + "- **ESS bulk min : 418** (largement au-dessus du seuil 100)\n", + "\n", + "### Verdict\n", + "\n", + "Le Modèle est convergent (R-hat = 1.010, ESS = 418), mais avec quelques divergences. Pour une application réelle, on peut ajouter un `pm.Potential('order', pt.switch(mu_comp[0] < mu_comp[1], 0, -1e10))` pour eliminer completement le label switching." + ] + }, { "cell_type": "code", "execution_count": 11, @@ -1276,6 +1826,78 @@ " print(f\"Composante {i+1} : mu={mu_nm[:, idx].mean():.1f} min\")" ] }, + { + "cell_type": "markdown", + "id": "fuse-13", + "metadata": {}, + "source": [ + "### Lecture de l'API directe pm.NormalMixture (code[10])\n", + "\n", + "**Sortie verbatim code[10]** :\n", + "- `Multiprocess sampling (4 chains in 4 jobs)`\n", + "- `NUTS: [w, mu_comp, sigma_comp]`\n", + "- `Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 72 seconds.`\n", + "- `There were 62 divergences after tuning.`\n", + "- `R-hat max : 1.010`\n", + "- `ESS bulk min : 418`\n", + "- `
`\n", + "- `Composante 1 : mu=16.9 min / Composante 2 : mu=26.2 min`\n", + "\n", + "### Lecture\n", + "\n", + "- **72 secondes** pour 12000 échantillons -- comparable au Modèle manuel (76s pour code[8]).\n", + "- **62 divergences** au lieu de 439 -- l'API directe `pm.NormalMixture` est numeriquement plus stable.\n", + "- **R-hat = 1.010** -- juste au seuil, acceptable.\n", + "- **ESS = 418** -- largement au-dessus du seuil 100.\n", + "\n", + "### Comparaison avec code[8]\n", + "\n", + "| Metrique | code[8] (manuel) | code[10] (NormalMixture) |\n", + "|----------|------------------|--------------------------|\n", + "| Divergences | 439 | **62** (7x moins) |\n", + "| R-hat max | 1.030 | 1.010 |\n", + "| ESS bulk min | 250 | 418 |\n", + "| Temps | 76 s | 72 s |\n", + "\n", + "### Pourquoi l'API directe est meilleure\n", + "\n", + "`pm.NormalMixture` utilise une représentation interne optimisee qui marginalise partiellement le sujet latent z. Cela reduit la dimensionalite de la posterior et ameliore la geometrie pour le NUTS.\n", + "\n", + "### Verdict\n", + "\n", + "Pour les modèles a mélange, preferer `pm.NormalMixture` a l'approche manuelle. C'est plus propre et numeriquement plus stable." + ] + }, + { + "cell_type": "markdown", + "id": "fuse-14", + "metadata": {}, + "source": [ + "| Aspect | Infer.NET | PyMC |\n", + "|--------|-----------|------|\n", + "| **Inférence** | EP (Expectation Propagation) | MCMC (NUTS, Metropolis) ou ADVI |\n", + "| **Vitesse (gaussienne)** | ~0.01 s (EP converge en 5-10 iterations) | ~20 s (MCMC 12000 draws) |\n", + "| **Vitesse (GMM)** | ~0.1 s (EP avec collapsed Gibbs) | ~76 s (MCMC avec CompoundStep) |\n", + "| **Determinisme** | Oui (EP deterministe) | Non (MCMC) mais reproductible avec `random_seed` |\n", + "| **Diagnostiques** | Limites (pas de rhat, ESS) | Complets (rhat, ESS, divergences, posterior predictive) |\n", + "| **Posterieurs exactes** | Approximations EP (peuvent etre biaisees) | Asymptotiquement exactes (MCMC) |\n", + "| **Modèles hierarchiques** | Limites (EP peut diverger sur les modèles a plusieurs niveaux) | Excellents (NUTS gradiente-based) |\n", + "| **Multi-plateforme** | Windows + .NET, Mac/Linux via Mono | Multi (Python pur) |\n", + "\n", + "### Quand utiliser quoi ?\n", + "\n", + "- **Infer.NET** : inférence rapide sur des modèles classiques, prototype rapide, integration dans une application .NET.\n", + "- **PyMC** : modèles hierarchiques complexes, inférence exacte (MCMC), posterieurs detaillees pour analyse scientifique.\n", + "- **numpyro** : JAX backend, GPU acceleration, modèles a très grande echelle (> 10^6 observations).\n", + "- **Stan** : alternative a PyMC, plus rapide sur certains modèles, mais API plus verbeuse.\n", + "\n", + "### Pour aller plus loin\n", + "\n", + "- **Modèles de mélange hierarchiques** (HMM, LDA) : les sujets latents dependent du contexte (document, temps, position).\n", + "- **Modèles de Markov caches** (HMM) : extension temporelle du GMM.\n", + "- **Dirichlet Process Mixtures** (DPM) : nombre de composantes K infere depuis les données (au lieu de fixe a priori)." + ] + }, { "cell_type": "markdown", "id": "be434f27", @@ -1368,6 +1990,50 @@ "print(\"Exercice a completer\")" ] }, + { + "cell_type": "markdown", + "id": "hand-mel3guide", + "metadata": {}, + "source": [ + "### Exemple guidé — Mélange à 3 composantes (à consulter après votre tentative)\n", + "\n", + "```python\n", + "K3 = 3\n", + "with pm.Model() as cycliste_3_comp:\n", + " w = pm.Dirichlet('w', a=np.ones(K3))\n", + " mu_comp = pm.Normal('mu_comp', mu=np.array([10, 16, 30]), sigma=10, shape=K3)\n", + " sigma_comp = pm.HalfNormal('sigma_comp', sigma=5, shape=K3)\n", + " trajets = pm.NormalMixture('trajets', w=w, mu=mu_comp, sigma=sigma_comp,\n", + " observed=donnees_mixtes)\n", + "\n", + " trace_3 = pm.sample(3000, random_seed=42, return_inferencedata=True)\n", + "\n", + "for k in range(K3):\n", + " print(f\"Composante {k+1} : mu={trace_3.posterior['mu_comp'].values[:, :, k].mean():.1f} min, \"\n", + " f\"poids={trace_3.posterior['w'].values[:, :, k].mean():.2f}\")\n", + "```\n", + "\n", + "**Verdict honnête** : sur 12 observations, les composantes « normales » (~16 min) et « lentes » (~30 min) sont bien identifiées, mais la composante « rapide » (~10 min) reste mal contrainte (poids et écart-type variables selon les chains) — c'est le piège de surapprentissage signalé par l'exercice. Un modèle à 2 composantes reste préférable avec si peu de données." + ] + }, + { + "cell_type": "markdown", + "id": "fuse-15", + "metadata": {}, + "source": [ + "### Lecture de l'exercice 3 composantes stub (code[11])\n", + "\n", + "**Stub** : la cellule affiche « Exercice a compléter » tant que le code n'est pas ecrit.\n", + "\n", + "### Alternative : model sélection\n", + "\n", + "Pour choisir entre 2 et 3 composantes, on peut utiliser le **WAIC** (Widely Applicable Information Criterion) ou le **LOO-CV** (Leave-One-Out Cross-Validation). Plus petit = meilleur Modèle.\n", + "\n", + "### Indice\n", + "\n", + "La cellule affiche « Exercice a compléter » tant que le code n'est pas ecrit. Vous pouvez reutiliser code[10] comme template en changeant `shape=2` a `shape=3` et les priors sur mu_comp pour 3 valeurs." + ] + }, { "cell_type": "markdown", "id": "7205e762", @@ -1403,6 +2069,39 @@ "- Redner, R. A., & Walker, H. F. (1984). *Mixture densities, maximum likelihood and the EM algorithm.* SIAM Review 26(2).\n", "- Stephens, M. (2000). *Dealing with label switching in mixture models.* JRSS B 62(4).\n" ] + }, + { + "cell_type": "markdown", + "id": "fuse-16", + "metadata": {}, + "source": [ + "Nous avons explore les mélanges gaussiens : des modèles capables de capturer des distributions **multi-modales** (plusieurs pics) en superposant K gaussiennes avec des poids.\n", + "\n", + "### Ce que ce notebook demontre\n", + "\n", + "1. **Le Modèle gaussien simple** : posterior conjugue Normal-Inverse-Gamma, inférence NUTS en 19 secondes.\n", + "2. **La distribution predictive postérieure** : comment predire la prochaine observation en integrant l'incertitude sur les paramètres.\n", + "3. **La gaussienne tronquée** : `pm.TruncatedNormal(lower, upper)` restreint le support, utile pour des variables contraintes (temps, distances, etc.).\n", + "4. **Le posterior predictive check** : avant de complexifier un Modèle, on verifie que la version simple capture la structure des données. Verdict sur la gaussienne unique : REJETEE (creux sous-estime).\n", + "5. **Le mélange de Gaussiennes** : Modèle a 2 composantes pour des données bimodales, posterior : 58% normaux / 42% exceptionnels.\n", + "6. **Le piege du label switching** : 439 divergences et rhat > 1.01 -- symptome classique des mélanges symétriques. Remedes : ordering constraint, reparameterisation, API `pm.NormalMixture` (62 divergences).\n", + "7. **API directe `pm.NormalMixture`** : simplifie le code en evitant de declarer explicitement le sujet latent z.\n", + "\n", + "### Pour aller plus loin\n", + "\n", + "- **Dirichlet Process Mixtures** (DPM) : inférer K automatiquement depuis les données (Neal 2000).\n", + "- **Hidden Markov Models** (HMM) : extension temporelle du GMM (Rabiner 1989).\n", + "- **Variational Autoencoders** (VAE) : modèles generatifs profonds appris par SGD + reparameterisation (Kingma 2014).\n", + "- **Normalizing flows** : transformations invertibles pour modeliser des distributions complexes (Rezende 2015).\n", + "\n", + "### References\n", + "\n", + "- Bishop 2006. *Pattern Recognition and Machine Learning*. Springer. (Chapitre 9 : Mixture Models and EM)\n", + "- Murphy 2012. *Machine Learning: A Probabilistic Perspective*. MIT Press. (Chapitre 11 : Mixture models)\n", + "- McLachlan & Peel 2000. *Finite Mixture Models*. Wiley.\n", + "- Carpenter et al. 2017. *Stan: A Probabilistic Programming Language*. JSS.\n", + "- PyMC Team 2024. *PyMC 5.28.5 documentation*. https://www.pymc.io/" + ] } ], "metadata": { diff --git a/scripts/notebook_tools/twin_pairs.d/probas-2-gaussian-mixtures.yaml b/scripts/notebook_tools/twin_pairs.d/probas-2-gaussian-mixtures.yaml index 345780ff99..f465db2de1 100644 --- a/scripts/notebook_tools/twin_pairs.d/probas-2-gaussian-mixtures.yaml +++ b/scripts/notebook_tools/twin_pairs.d/probas-2-gaussian-mixtures.yaml @@ -83,6 +83,12 @@ csharp_sha: 1c64917d88cf9f4368921f295fe1e47bfae04496 content_python_sha: b5ffd31af770f80213fbccf8fcae9e2063bbb7638382fd09903eb830d44e0471 content_csharp_sha: 8279ddaf059df075bab41e33384aac8fbd7c2f21b4e0d4a0362268b5b4aec096 + - date: "2026-09-03" + by: myia-po-2026:CoursIA + python_sha: 6e8bf5b57ff02b7a0e4900c20a93cdb419a222ad + csharp_sha: 1c64917d88cf9f4368921f295fe1e47bfae04496 + content_python_sha: f2cb5fe8e1f8578dd8ceb40d34c8ac617c27c2674451bfccba32d16050d2dea8 + content_csharp_sha: 8279ddaf059df075bab41e33384aac8fbd7c2f21b4e0d4a0362268b5b4aec096 known_differences: - '2026-08-21 (rebaseline asymetrique myia-po-2024:CoursIA-2) : PR #12065 corrige heading_in_list cote C# (Infer-2-Gaussian-Mixtures.ipynb, 1 finding) tandis que le jumeau Python (PyMC-02-Gaussian-Mixtures.ipynb) etait deja a 0 finding. La convergence unilaterale est legitime (cf ai-01 arbitrage 2026-08-21 15:00Z, DM msg-20260821T130006-lc9u2c) : la parite reelle s''ameliore, elle n''est pas enshrinee. Aucun strip outille en attente (probe banner 0, papermill 0, machine paths 0).' - 'Socle pedagogique commun : melanges de distributions gaussiennes (GMM) modelise dans les deux stacks.'