Serie : Programmation Probabiliste avec PyMC (4/7)
Duree estimee : 55 minutes
Prerequis : Notebooks 3 (Graphes probabilistes), 14-16 (Utilite)
| Précédent | Suivant |
|---|---|
| DecPyMC-3-Multi-Attribute | DecPyMC-5-Value-Information |
Rappel : Reseaux Bayesiens
Un reseau bayesien represente des relations probabilistes entre variables :
[Maladie] ----> [Symptome]
\--------> [Test]Un reseau de decision ajoute :
[Maladie] ----> [Test Result]
| |
v v
<Traitement> ----> ((Utilite))Nous chargeons les bibliotheques necessaires : PyMC pour la modelisation probabiliste, ArviZ pour l’analyse des résultats, et numpy/scipy pour les calculs numériques.
PyMC version: 5.28.5
ArviZ version: 0.23.4
Environnement pret !Les reseaux de decision utilisent une notation graphique standardisee (Howard & Matheson, 1984) qui distingue clairement les différents types de noeuds. Cette convention visuelle permet de lire immediatement la structure du problème decisionnel.
| Type | Forme | Symbole | Description | Exemple |
|---|---|---|---|---|
| Chance | Ovale/Cercle | [X] | Variable aleatoire non controlee | Maladie, Marche, Meteo |
| Decision | Rectangle | Choix contrôle par l’agent | Traiter, Investir, Acheter | |
| Utilite | Losange/Hexagone | ((U)) | Fonction de valeur/payoff | Profit, Sante, Satisfaction |
Chaque type d’arc encode une relation spécifique dans le problème :
| Type d’arc | Signification | Interpretation |
|---|---|---|
| Chance -> Chance | Dépendance probabiliste | P(B|A) - A influence B |
| Chance -> Decision | Arc informationnel | L’agent observe A avant de choisir D |
| Chance -> Utilite | La variable affecte l’utilite | U depend de la valeur de A |
| Decision -> Chance | La decision influence le résultat | L’action D modifie P(B) |
| Decision -> Utilite | Cout/benefice direct de la decision | U inclut le cout de D |
Dans un reseau de decision bien forme : 1. Les noeuds de decision sont totalement ordonnes (D1 avant D2 avant D3…) 2. Un arc informationnel vers Di signifie que la variable est connue avant Di 3. Pas de cycles impliquant des noeuds de decision
La cellule suivante définit une structure de données simple pour representer les noeuds d’un reseau de decision et illustre la notation graphique sur un exemple medical.
class NodeType(Enum):
CHANCE = "chance"
DECISION = "decision"
UTILITY = "utility"
@dataclass
class DecisionNode:
name: str
node_type: NodeType
parents: list[str] = field(default_factory=list)
@property
def symbol(self) -> str:
if self.node_type == NodeType.CHANCE:
return f"[{self.name}]"
elif self.node_type == NodeType.DECISION:
return f"<{self.name}>"
else:
return f"(({self.name}))"
# Exemple : Decision medicale
network = [
DecisionNode("Maladie", NodeType.CHANCE),
DecisionNode("Test", NodeType.CHANCE, parents=["Maladie"]),
DecisionNode("Traiter", NodeType.DECISION, parents=["Test"]),
DecisionNode("Sante", NodeType.CHANCE, parents=["Maladie", "Traiter"]),
DecisionNode("Utilite", NodeType.UTILITY, parents=["Sante", "Traiter"]),
]
print("Reseau de decision : Diagnostic Medical\n")
for node in network:
parents_str = f" <- {', '.join(node.parents)}" if node.parents else ""
print(f" {node.symbol}{parents_str}")Reseau de decision : Diagnostic Medical
[Maladie]
[Test] <- Maladie
<Traiter> <- Test
[Sante] <- Maladie, Traiter
((Utilite)) <- Sante, TraiterLe reseau ci-dessus represente un problème de diagnostic medical classique :
| Noeud | Type | Rôle |
|---|---|---|
[Maladie] |
Chance | Etat du patient (inconnu) |
[Test] |
Chance | Résultat observable, depend de Maladie |
<Traiter> |
Decision | Choix du medecin, base sur le résultat du test |
[Sante] |
Chance | Etat futur, depend de Maladie ET de la decision |
((Utilite)) |
Utilite | Valeur finale, depend de Sante et du cout du traitement |
Point cle : L’arc
Test -> Traiterest un arc informationnel. Il indique que le medecin connait le résultat du test avant de decider. Sans cet arc, le medecin devrait decider “en aveugle”.
Les arcs informationnels sont fondamentaux pour modeliser ce que l’agent sait au moment de prendre une decision. Ils capturent l’asymetrie d’information qui rend les problemes de decision interessants.
Un arc informationnel vers un noeud de decision indique ce que l’agent sait au moment de prendre cette decision. C’est la différence entre agir en aveugle et agir en connaissance de cause.
| Situation | Politique | Complexite |
|---|---|---|
| Aucune information | Action unique | 1 decision |
| Test observe | Action conditionnelle au test | 2^k decisions (k = bits d’information) |
| Information parfaite | Action conditionnelle a l’etat | Action optimale par etat |
L’ajout d’un arc informationnel peut augmenter l’utilite esperee. La différence est la valeur de l’information que nous etudierons dans le notebook suivant (DecPyMC-5).
# Impact des arcs informationnels
print("Scenarios d'information :\n")
# Prior sur la maladie
p_maladie = 0.1
sensibilite = 0.95 # P(test+|malade)
specificite = 0.90 # P(test-|sain)
# Utilites
U_traiter_malade = 90 # Guerison
U_traiter_sain = 70 # Effets secondaires inutiles
U_pas_traiter_malade = 10 # Maladie non traitee
U_pas_traiter_sain = 100 # Parfait
# Scenario 1 : Pas d'information (decision a priori)
EU_traiter = p_maladie * U_traiter_malade + (1 - p_maladie) * U_traiter_sain
EU_pas_traiter = p_maladie * U_pas_traiter_malade + (1 - p_maladie) * U_pas_traiter_sain
print("Scenario 1 : Decision SANS information (pas de test)")
print(f" P(malade) = {p_maladie:.0%}")
print(f" E[U(traiter)] = {EU_traiter:.1f}")
print(f" E[U(pas traiter)] = {EU_pas_traiter:.1f}")
print(f" => Decision : {'TRAITER' if EU_traiter > EU_pas_traiter else 'PAS TRAITER'}\n")
# Scenario 2 : Decision AVEC information (apres test)
# Calculer P(malade|test+) et P(malade|test-) par Bayes
p_test_pos = p_maladie * sensibilite + (1 - p_maladie) * (1 - specificite)
p_malade_test_pos = (p_maladie * sensibilite) / p_test_pos
p_malade_test_neg = (p_maladie * (1 - sensibilite)) / (1 - p_test_pos)
print("Scenario 2 : Decision AVEC information (apres test)")
print(f" P(test+) = {p_test_pos:.1%}")
print(f" P(malade|test+) = {p_malade_test_pos:.1%}")
print(f" P(malade|test-) = {p_malade_test_neg:.1%}")
# Decision si test positif
EU_traiter_tp = p_malade_test_pos * U_traiter_malade + (1 - p_malade_test_pos) * U_traiter_sain
EU_pas_traiter_tp = p_malade_test_pos * U_pas_traiter_malade + (1 - p_malade_test_pos) * U_pas_traiter_sain
print(f" Si test+ : traiter={EU_traiter_tp:.1f}, pas traiter={EU_pas_traiter_tp:.1f} "
f"=> {'TRAITER' if EU_traiter_tp > EU_pas_traiter_tp else 'PAS TRAITER'}")
# Decision si test negatif
EU_traiter_tn = p_malade_test_neg * U_traiter_malade + (1 - p_malade_test_neg) * U_traiter_sain
EU_pas_traiter_tn = p_malade_test_neg * U_pas_traiter_malade + (1 - p_malade_test_neg) * U_pas_traiter_sain
print(f" Si test- : traiter={EU_traiter_tn:.1f}, pas traiter={EU_pas_traiter_tn:.1f} "
f"=> {'TRAITER' if EU_traiter_tn > EU_pas_traiter_tn else 'PAS TRAITER'}")Scenarios d'information :
Scenario 1 : Decision SANS information (pas de test)
P(malade) = 10%
E[U(traiter)] = 72.0
E[U(pas traiter)] = 91.0
=> Decision : PAS TRAITER
Scenario 2 : Decision AVEC information (apres test)
P(test+) = 18.5%
P(malade|test+) = 51.4%
P(malade|test-) = 0.6%
Si test+ : traiter=80.3, pas traiter=53.8 => TRAITER
Si test- : traiter=70.1, pas traiter=99.4 => PAS TRAITERVous etes responsable des recrutements. Un candidat se presente, et vous devez decider de l’embaucher ou non. Avant de decider, vous pouvez demander une lettre de reference qui vous renseignera (imparfaitement) sur la qualite du candidat.
Scénario : - Variable aleatoire : bon_candidat (Oui/Non), avec prior p_bon - Test : lettre de reference, avec sensibilite (recommande un bon candidat) et specificite (deconseille un mauvais candidat) - Decision : embaucher ou refuser - Utilites : productivite d’un bon employe, cout d’un mauvais recrutement, neutralite si pas d’embauche
Objectif : Comparer la decision avec et sans lettre de reference. L’information supplementaire change-t-elle la decision optimale ?
Indices : - Étape 1 : Définir les probabilites a priori et la qualite de la lettre de reference - Étape 2 : Calculer E[U] sans lettre de reference (decision a priori) - Étape 3 : Appliquer le theoreme de Bayes pour obtenir les posteriors P(bon|lettre+) et P(bon|lettre-) - Étape 4 : Calculer E[U] avec lettre de reference et comparer les deux scénarios
# Exercice : Decision de recrutement avec lettre de reference
# Etape 1 : Parametres du probleme
p_bon = None # TODO etudiant : P(bon candidat), essayez 0.3
sens_lettre = None # TODO etudiant : sensibilite (lettre favorable si bon), essayez 0.80
spec_lettre = None # TODO etudiant : specificite (lettre defavorable si mauvais), essayez 0.75
# Utilites
U_embaucher_bon = None # TODO etudiant : gain si bon recrutement, essayez 100
U_embaucher_mauvais = None # TODO etudiant : cout si mauvais recrutement, essayez -60
U_pas_embaucher = 0 # Neutre
# Etape 2 : Decision sans lettre de reference (a priori)
# TODO etudiant : calculer E[U(embaucher)] et E[U(pas embaucher)]
eu_embaucher_sans = None # TODO etudiant
eu_pas_embaucher_sans = None # TODO etudiant
# Etape 3 : Posterior apres lecture de la lettre (theoreme de Bayes)
# TODO etudiant : calculer P(lettre+), P(bon|lettre+), P(bon|lettre-)
p_lettre_pos = None # TODO etudiant
p_bon_lpos = None # TODO etudiant : P(bon candidat|lettre favorable)
p_bon_lneg = None # TODO etudiant : P(bon candidat|lettre defavorable)
# Etape 4 : Decision avec lettre de reference
# TODO etudiant : pour chaque resultat de lettre, calculer E[U(embaucher)] et E[U(pas embaucher)]
# puis calculer E[U] totale avec information
eu_avec_lettre = None # TODO etudiant
# Resultat : comparer les deux scenarios
result = None # TODO etudiant : eu_avec_lettre - max(eu_embaucher_sans, eu_pas_embaucher_sans)
print("Exercice a completer : decision de recrutement avec/sans lettre de reference")Exercice a completer : decision de recrutement avec/sans lettre de referenceLes résultats ci-dessus illustrent le theoreme fondamental de la valeur de l’information :
| Scénario | Decision | E[U] | Commentaire |
|---|---|---|---|
| Sans test | Pas traiter | 91.0 | Prior faible (10%), risque acceptable |
| Avec test+ | Traiter | 80.3 | Posterior eleve (51.4%), traitement justifie |
| Avec test- | Pas traiter | 99.4 | Posterior très faible (0.6%), rassurant |
A priori : P(malade) = 10% est faible. Le cout des effets secondaires (U=70 vs 100) domine.
Test positif : Le theoreme de Bayes concentre la probabilite : \[P(\text{malade}|\text{test}+) = \frac{0.10 \times 0.95}{0.185} = 51.4\%\] Le risque devient majoritaire, justifiant le traitement.
Test negatif : La probabilite residuelle (0.6%) est negligeable.
Lecon : L’information ne change pas toujours la decision (si le prior etait de 50%, on traiterait dans tous les cas). La valeur de l’information depend de la distance entre le prior et les seuils de decision.
L’algorithme standard pour resoudre un reseau de decision est la backward induction (ou programmation dynamique de Bellman, 1957), formalisee pour les diagrammes d’influence par Shachter (1986). L’idee est de resoudre le problème de la fin vers le debut.
Étape 1 : Dernière decision (Dn) - Pour chaque configuration possible des parents informationnels de Dn : - Calculer E[U | parents, action] pour chaque action possible - Sélectionner l’action qui maximise E[U]
Étape 2 : Remonter (Dn-1, Dn-2, …) - Utiliser les valeurs optimales des decisions suivantes
Étape 3 : Première decision (D1) - La dernière itération donne la decision optimale initiale
| Nombre de decisions | Complexite | Commentaire |
|---|---|---|
| 1 | O(|A| x |S|) | A = actions, S = etats |
| n séquentielles | O(|A|^n x |S|^n) | Exponentielle en n |
| Avec info parfaite | O(|S| x |A|) | Lineaire car decision par etat |
@dataclass
class SimpleDecisionNetwork:
"""Solveur de reseau de decision par backward induction."""
# Noeud chance : Maladie (M)
p_m: float = 0.1 # P(malade)
# Noeud chance : Test (T) | M
p_tpos_given_m: float = 0.95 # Sensibilite
p_tpos_given_notm: float = 0.10 # 1 - Specificite
# Noeud utilite : U(M, D)
u_traiter_malade: float = 90
u_traiter_sain: float = 70
u_pas_traiter_malade: float = 10
u_pas_traiter_sain: float = 100
def get_utility(self, malade: bool, traiter: bool) -> float:
if traiter and malade:
return self.u_traiter_malade
if traiter and not malade:
return self.u_traiter_sain
if not traiter and malade:
return self.u_pas_traiter_malade
return self.u_pas_traiter_sain
def solve(self) -> tuple[dict[bool, bool], float]:
"""Resout le reseau par backward induction.
Returns:
policy: {test_result: should_treat}
expected_utility: utilite esperee globale
"""
policy = {}
# P(T+)
p_tpos = self.p_m * self.p_tpos_given_m + (1 - self.p_m) * self.p_tpos_given_notm
p_tneg = 1 - p_tpos
# Pour T+ : decision optimale
p_m_tpos = (self.p_m * self.p_tpos_given_m) / p_tpos
eu_traiter_tp = p_m_tpos * self.u_traiter_malade + (1 - p_m_tpos) * self.u_traiter_sain
eu_pas_traiter_tp = p_m_tpos * self.u_pas_traiter_malade + (1 - p_m_tpos) * self.u_pas_traiter_sain
policy[True] = eu_traiter_tp > eu_pas_traiter_tp
best_u_tp = max(eu_traiter_tp, eu_pas_traiter_tp)
# Pour T- : decision optimale
p_m_tneg = (self.p_m * (1 - self.p_tpos_given_m)) / p_tneg
eu_traiter_tn = p_m_tneg * self.u_traiter_malade + (1 - p_m_tneg) * self.u_traiter_sain
eu_pas_traiter_tn = p_m_tneg * self.u_pas_traiter_malade + (1 - p_m_tneg) * self.u_pas_traiter_sain
policy[False] = eu_traiter_tn > eu_pas_traiter_tn
best_u_tn = max(eu_traiter_tn, eu_pas_traiter_tn)
# Utilite esperee totale
total_eu = p_tpos * best_u_tp + p_tneg * best_u_tn
return policy, total_eu
dn = SimpleDecisionNetwork()
policy, eu = dn.solve()
print("=== Resolution du Reseau de Decision ===")
print()
print("Politique optimale :")
print(f" Si test POSITIF : {'TRAITER' if policy[True] else 'PAS TRAITER'}")
print(f" Si test NEGATIF : {'TRAITER' if policy[False] else 'PAS TRAITER'}")
print()
print(f"Utilite esperee : {eu:.2f}")=== Resolution du Reseau de Decision ===
Politique optimale :
Si test POSITIF : TRAITER
Si test NEGATIF : PAS TRAITER
Utilite esperee : 95.90La classe SimpleDecisionNetwork implemente l’algorithme de backward induction de maniere explicite. Le résultat est une politique conditionnelle :
| Observation | P(malade|obs) | Action optimale | E[U|obs, action*] |
|---|---|---|---|
| Test positif | 51.4% | TRAITER | 80.3 |
| Test negatif | 0.6% | PAS TRAITER | 99.4 |
L’utilite esperee globale de 95.90 est calculee comme :
\[E[U] = P(\text{test}+) \times E[U|\text{test}+, \text{action}^*] + P(\text{test}-) \times E[U|\text{test}-, \text{action}^*]\]
Remarque : Cette utilite (95.90) est superieure a celle sans test (91.0). La différence (4.90) represente la valeur de l’information apportee par le test.
Un investisseur considere un projet. Le marche peut etre favorable ou non. Il peut commander une étude de marche avant de decider.
[Marche] -----> [Étude] -----> <Investir>
| |
+----------------------------> ((Profit))# Reseau de decision : Investissement
# Marche (M) : Favorable (F) ou Defavorable (D)
p_favorable = 0.4
# Qualite de l'etude
p_epos_given_f = 0.8 # Detecte correctement marche favorable
p_epos_given_d = 0.3 # Faux positif
# Profits
profit_investir_favorable = 500_000
profit_investir_defavorable = -200_000
profit_pas_investir = 0
cout_etude = 20_000
print("=== Decision d'Investissement ===")
print(f"P(marche favorable) = {p_favorable:.0%}")
print(f"Cout de l'etude = {cout_etude:,.0f} EUR\n")
# Option 1 : Pas d'etude, decision directe
EU_investir_direct = (p_favorable * profit_investir_favorable
+ (1 - p_favorable) * profit_investir_defavorable)
EU_pas_investir_direct = profit_pas_investir
best_eu_sans_etude = max(EU_investir_direct, EU_pas_investir_direct)
print("Option 1 : Decision directe (sans etude)")
print(f" E[Profit(investir)] = {EU_investir_direct:,.0f} EUR")
print(f" E[Profit(pas investir)] = {EU_pas_investir_direct:,.0f} EUR")
print(f" => Decision : {'INVESTIR' if EU_investir_direct > EU_pas_investir_direct else 'NE PAS INVESTIR'}")
print(f" => E[Profit] = {best_eu_sans_etude:,.0f} EUR\n")
# Option 2 : Avec etude
p_epos = p_favorable * p_epos_given_f + (1 - p_favorable) * p_epos_given_d
p_eneg = 1 - p_epos
# Si etude positive
p_f_epos = (p_favorable * p_epos_given_f) / p_epos
EU_investir_epos = p_f_epos * profit_investir_favorable + (1 - p_f_epos) * profit_investir_defavorable
investir_si_epos = EU_investir_epos > profit_pas_investir
best_eu_epos = max(EU_investir_epos, profit_pas_investir)
# Si etude negative
p_f_eneg = (p_favorable * (1 - p_epos_given_f)) / p_eneg
EU_investir_eneg = p_f_eneg * profit_investir_favorable + (1 - p_f_eneg) * profit_investir_defavorable
investir_si_eneg = EU_investir_eneg > profit_pas_investir
best_eu_eneg = max(EU_investir_eneg, profit_pas_investir)
EU_avec_etude = p_epos * best_eu_epos + p_eneg * best_eu_eneg - cout_etude
print("Option 2 : Avec etude de marche")
print(f" P(etude positive) = {p_epos:.1%}")
print(f" Si E+ : P(favorable|E+) = {p_f_epos:.1%} => {'INVESTIR' if investir_si_epos else 'NE PAS INVESTIR'}")
print(f" Si E- : P(favorable|E-) = {p_f_eneg:.1%} => {'INVESTIR' if investir_si_eneg else 'NE PAS INVESTIR'}")
print(f" => E[Profit avec etude] = {EU_avec_etude:,.0f} EUR\n")
print("=== Comparaison ===")
print(f"Sans etude : {best_eu_sans_etude:,.0f} EUR")
print(f"Avec etude : {EU_avec_etude:,.0f} EUR")
print(f"=> Valeur de l'etude : {EU_avec_etude - best_eu_sans_etude:,.0f} EUR")
print(f"=> Decision : {'COMMANDER L\'ETUDE' if EU_avec_etude > best_eu_sans_etude else 'DECIDER DIRECTEMENT'}")=== Decision d'Investissement ===
P(marche favorable) = 40%
Cout de l'etude = 20,000 EUR
Option 1 : Decision directe (sans etude)
E[Profit(investir)] = 80,000 EUR
E[Profit(pas investir)] = 0 EUR
=> Decision : INVESTIR
=> E[Profit] = 80,000 EUR
Option 2 : Avec etude de marche
P(etude positive) = 50.0%
Si E+ : P(favorable|E+) = 64.0% => INVESTIR
Si E- : P(favorable|E-) = 16.0% => NE PAS INVESTIR
=> E[Profit avec etude] = 104,000 EUR
=== Comparaison ===
Sans etude : 80,000 EUR
Avec etude : 104,000 EUR
=> Valeur de l'etude : 24,000 EUR
=> Decision : COMMANDER L'ETUDEVous envisagez d’acheter une maison. Avant de signer, vous pouvez commander une expertise structurelle pour detecter d’eventuels problemes (fissures, fondations, humidite). L’expertise a un cout, mais peut vous eviter un achat desastreux.
Scénario : - Variable aleatoire : problemes_structurels (Oui/Non), prior p_problemes - Test : expertise structurelle avec sensibilite sens_expert et specificite spec_expert - Decision : acheter ou ne pas acheter - Utilites : gain/perte selon presence de problemes
Objectif : Calculer l’esperance d’utilite avec et sans expertise, en deduire la valeur de l’information (EVPI/EVSI).
Indices : - Étape 1 : Définir les probabilites et utilites du problème - Étape 2 : Calculer E[U] sans expertise (decision a priori) - Étape 3 : Calculer les posteriors P(problemes|expertise+) et P(problemes|expertise-) via Bayes - Étape 4 : Calculer E[U] avec expertise et la valeur nette de l’information
# Exercice : Achat immobilier et expertise structurelle
# Etape 1 : Parametres du probleme
p_problemes = None # TODO etudiant : P(problemes structurels), essayez 0.15
sens_expert = None # TODO etudiant : sensibilite de l'expertise, essayez 0.85
spec_expert = None # TODO etudiant : specificite de l'expertise, essayez 0.90
cout_expertise = None # TODO etudiant : cout de l'expertise en kEUR, essayez 3
# Utilites (en kEUR)
U_acheter_sain = None # TODO etudiant : gain si achat sans problemes, essayez 50
U_acheter_problemes = None # TODO etudiant : perte si achat avec problemes, essayez -80
U_pas_acheter = 0 # Ne pas acheter = pas de gain ni perte
# Etape 2 : Decision sans expertise (a priori)
# TODO etudiant : calculer E[U(acheter)] et E[U(pas acheter)] sans information
eu_acheter_sans = None # TODO etudiant
eu_pas_acheter_sans = None # TODO etudiant
# Etape 3 : Posterior apres expertise (theoreme de Bayes)
# TODO etudiant : calculer P(expertise+), P(problemes|expertise+), P(problemes|expertise-)
p_expert_pos = None # TODO etudiant
p_prob_epos = None # TODO etudiant : P(problemes|expertise+)
p_prob_eneg = None # TODO etudiant : P(problemes|expertise-)
# Etape 4 : Decision avec expertise et valeur de l'information
# TODO etudiant : calculer E[U] optimal si expertise+, si expertise-
# puis E[U] totale avec expertise (moins le cout)
eu_avec_expertise = None # TODO etudiant
# Resultat : valeur nette de l'information
result = None # TODO etudiant : eu_avec_expertise - max(eu_acheter_sans, eu_pas_acheter_sans)
print("Exercice a completer : valeur de l'information pour l'expertise structurelle")Exercice a completer : valeur de l'information pour l'expertise structurelle| Stratégie | Decision | E[Profit] | Commentaire |
|---|---|---|---|
| Sans étude | Investir directement | 80 000 EUR | Prior P(favorable)=40% suffit |
| Avec étude | Conditionnel au résultat | 104 000 EUR | Evite les grosses pertes |
La valeur de l’étude (24 000 EUR) se decompose ainsi :
Point cle : L’étude est rentable car elle permet de changer de decision dans certains cas. Si P(favorable) etait de 80%, on investirait dans tous les cas (E+ ou E-), et l’étude n’aurait aucune valeur decisionnelle.
Quand il y a plusieurs noeuds de decision, leur ordre temporel doit etre specifie.
<Faire Test?> ----> [Résultat] ----> <Traiter?> ----> ((Utilite))
^
|
[Maladie] ---------------------------------------------+La maladie reste non observée avant les deux décisions : l’arc Maladie -> Utilite indique son effet sur le gain, pas une information disponible pour Traiter. Le résultat dépend de la maladie ; il n’est observé que si le test est fait. Le schéma omet cette dépendance probabiliste et l’arc de coût Faire Test? -> Utilite.
Deux decisions séquentielles : 1. D’abord : Faire le test ou non ? 2. Ensuite : Traiter ou non ? (en connaissant le résultat si test fait ; sinon, selon la probabilité a priori)
# Decision sequentielle : Faire test puis traiter
# Parametres
p_maladie_seq = 0.2
sensibilite_seq = 0.90
specificite_seq = 0.85
cout_test = 50
# Utilites (sans cout test)
U_traiter_malade = 100
U_traiter_sain = 60
U_pas_traiter_malade = 20
U_pas_traiter_sain = 100
# Resolution par backward induction
print("=== Decision Sequentielle : Test puis Traitement ===\n")
# Politique pour D2 (traiter) sachant test fait
p_tpos_seq = p_maladie_seq * sensibilite_seq + (1 - p_maladie_seq) * (1 - specificite_seq)
p_m_tpos = (p_maladie_seq * sensibilite_seq) / p_tpos_seq
eu_traiter_tp = p_m_tpos * U_traiter_malade + (1 - p_m_tpos) * U_traiter_sain
eu_pas_traiter_tp = p_m_tpos * U_pas_traiter_malade + (1 - p_m_tpos) * U_pas_traiter_sain
best_u_tp = max(eu_traiter_tp, eu_pas_traiter_tp) - cout_test
p_m_tneg = (p_maladie_seq * (1 - sensibilite_seq)) / (1 - p_tpos_seq)
eu_traiter_tn = p_m_tneg * U_traiter_malade + (1 - p_m_tneg) * U_traiter_sain
eu_pas_traiter_tn = p_m_tneg * U_pas_traiter_malade + (1 - p_m_tneg) * U_pas_traiter_sain
best_u_tn = max(eu_traiter_tn, eu_pas_traiter_tn) - cout_test
eu_faire_test = p_tpos_seq * best_u_tp + (1 - p_tpos_seq) * best_u_tn
# Si pas de test : decision a priori
eu_traiter_ap = p_maladie_seq * U_traiter_malade + (1 - p_maladie_seq) * U_traiter_sain
eu_pas_traiter_ap = p_maladie_seq * U_pas_traiter_malade + (1 - p_maladie_seq) * U_pas_traiter_sain
eu_pas_test = max(eu_traiter_ap, eu_pas_traiter_ap)
print("Etape 2 : Politique de traitement si test fait")
print(f" Si T+ : P(M|T+)={p_m_tpos:.1%} => {'Traiter' if eu_traiter_tp > eu_pas_traiter_tp else 'Pas traiter'}")
print(f" Si T- : P(M|T-)={p_m_tneg:.1%} => {'Traiter' if eu_traiter_tn > eu_pas_traiter_tn else 'Pas traiter'}\n")
print("Etape 1 : Decision de faire le test")
print(f" E[U | faire test] = {eu_faire_test:.1f}")
print(f" E[U | pas test] = {eu_pas_test:.1f} ({'traiter' if eu_traiter_ap > eu_pas_traiter_ap else 'pas traiter'})\n")
print("=== Politique Optimale Globale ===")
print(f"D1 : {'FAIRE LE TEST' if eu_faire_test > eu_pas_test else 'NE PAS FAIRE LE TEST'}")
if eu_faire_test > eu_pas_test:
print(f"D2 : Si T+ => {'TRAITER' if eu_traiter_tp > eu_pas_traiter_tp else 'PAS TRAITER'}")
print(f" Si T- => {'TRAITER' if eu_traiter_tn > eu_pas_traiter_tn else 'PAS TRAITER'}")
else:
print(f"D2 : {'TRAITER' if eu_traiter_ap > eu_pas_traiter_ap else 'PAS TRAITER'} (sans test)")=== Decision Sequentielle : Test puis Traitement ===
Etape 2 : Politique de traitement si test fait
Si T+ : P(M|T+)=60.0% => Traiter
Si T- : P(M|T-)=2.9% => Pas traiter
Etape 1 : Decision de faire le test
E[U | faire test] = 43.6
E[U | pas test] = 84.0 (pas traiter)
=== Politique Optimale Globale ===
D1 : NE PAS FAIRE LE TEST
D2 : PAS TRAITER (sans test)Ce résultat est contre-intuitif : malgre un test informatif (sensibilite 90%, specificite 85%), la politique optimale est de ne pas faire le test.
Le cout du test (50 points) est applique dans les deux branches (T+ et T-), ce qui reduit significativement l’utilite esperee.
\[E[U|\text{faire test}] = P(T+) \times (E[U^*|T+] - 50) + P(T-) \times (E[U^*|T-] - 50)\]
Lecon importante : Dans les decisions séquentielles, le timing des couts est crucial. Un cout fixe en amont peut rendre toute la branche sous-optimale, même si l’information obtenue est de haute qualite.
PyMC (Salvatier et al., 2016) n’a pas de support natif pour les noeuds de decision, mais nous pouvons :
pm.BernoulliL’avantage de PyMC est de pouvoir effectuer l’inference exacte du posterior P(maladie|test) et de propager l’incertitude de maniere rigoureuse.
# Reseau de decision avec PyMC
# Probleme : Diagnostic et traitement
# Parametres
p_maladie_prior = 0.15
sensibilite = 0.92 # P(test+|malade)
specificite_pymc = 0.88 # P(test-|sain)
# Utilites
u_traiter_malade = 85
u_traiter_sain = 65
u_pas_traiter_malade = 15
u_pas_traiter_sain = 100
print("=== Reseau de Decision avec PyMC ===\n")
for obs_test in [True, False]:
with pm.Model() as diag_model:
# Variable de chance : Maladie
maladie = pm.Bernoulli("maladie", p=p_maladie_prior)
# Variable de chance : Test | Maladie
p_test = pm.math.switch(
maladie,
sensibilite, # P(test+|malade)
1 - specificite_pymc # P(test+|sain)
)
test_positif = pm.Bernoulli("test_positif", p=p_test, observed=int(obs_test))
# Inference
trace = pm.sample(5000, return_inferencedata=True, progressbar=False, random_seed=42, chains=4)
# Extraire P(maladie|test)
p_m = float(trace.posterior["maladie"].mean())
# Calculer E[U] pour chaque decision
eu_traiter = p_m * u_traiter_malade + (1 - p_m) * u_traiter_sain
eu_pas_traiter = p_m * u_pas_traiter_malade + (1 - p_m) * u_pas_traiter_sain
decision = "TRAITER" if eu_traiter > eu_pas_traiter else "PAS TRAITER"
print(f"Observation : test = {'POSITIF' if obs_test else 'NEGATIF'}")
print(f" P(maladie|test) = {p_m:.1%}")
print(f" E[U(traiter)] = {eu_traiter:.1f}")
print(f" E[U(pas traiter)] = {eu_pas_traiter:.1f}")
print(f" => Decision optimale : {decision}\n")=== Reseau de Decision avec PyMC ===
Observation : test = POSITIF
P(maladie|test) = 58.1%
E[U(traiter)] = 76.6
E[U(pas traiter)] = 50.6
=> Decision optimale : TRAITER
Observation : test = NEGATIF
P(maladie|test) = 1.7%
E[U(traiter)] = 65.3
E[U(pas traiter)] = 98.6
=> Decision optimale : PAS TRAITER
L’implementation PyMC utilise l’echantillonnage MCMC pour calculer le posterior P(maladie|test). Contrairement au calcul analytique, cette approche :
| Test observe | P(maladie|test) | E[U(traiter)] | E[U(pas traiter)] | Decision |
|---|---|---|---|---|
| POSITIF | ~57% | ~76 | ~51 | TRAITER |
| NEGATIF | ~2% | ~65 | ~99 | PAS TRAITER |
| Avantage | Description |
|---|---|
| Modèles complexes | Variables continues, dependances multiples, hiérarchies |
| Diagnostics | ArviZ fournit R-hat, ESS, trace plots |
| Ecosysteme | Integration avec scipy, matplotlib, xarray |
Note : Pour ce problème simple (2 variables binaires), le calcul analytique de Bayes est exact et plus rapide. L’approche PyMC prend tout son sens pour des modèles plus complexes.
La cellule précédente modelisait une seule variable binaire (maladie d’un patient) et laissait PyMC l’echantillonner par BinaryGibbsMetropolis – c’est-a-dire compter des 0 et des 1. Pour un problème a 2 variables binaires, le calcul analytique de Bayes donne la même reponse plus rapidement (cf. note ci-dessus) : le moteur MCMC n’etait pas reellement mis a contribution.
Pour faire valoir l’inference bayesienne par MCMC, considerons un problème non conjugue ou aucune forme close n’existe. Scénario : un reseau de cliniques applique le même test diagnostique (imparfait : sensibilite 0.92, specificite 0.88). On observe, par clinique, le nombre de patients testes \(n_i\) et de tests positifs \(k_i\) – mais jamais le vrai etat maladie (latent). Objectif : inferer la vraie prevalence de chaque clinique.
Piege du taux apparent : le taux de test positif naif \(k_i/n_i\) surestime la vraie prevalence, car il melange les vrais positifs (malades bien detectes) et les faux positifs (patients sains malgre tout testes positifs, \(1-\text{specificite}=0{,}12\)). Le taux apparent verifie \(p_{\text{app}} = p\cdot\text{se} + (1-p)\cdot(1-\text{sp})\).
Modèle hiérarchique non-centre (recette #3801) :
\[\text{logit}(p_i) = \mu_{pop} + \tau \cdot z_i, \quad z_i \sim \mathcal{N}(0,1), \quad k_i \sim \text{Binomial}(n_i,\; p_i\cdot\text{se} + (1-p_i)(1-\text{sp}))\]
Le prior logit-Normal hiérarchique + le lien non-lineaire entre vraie prevalence et taux apparent cassent toute conjugaison : MCMC est requis, et aucun calcul analytique par clinique (Beta-Binomial indépendant) ne peut produire cette correction du test imparfait ni le shrinkage entre cliniques.
# Prevalence hierarchique avec test imparfait (cas NON-CONJUGUE : MCMC requis)
# K cliniques, vraie prevalence inconnue, test imparfait connu (se, sp).
# On n'observe que les comptes de tests positifs k_i : l'etat maladie est LATENT.
K = 6
n_patients = np.array([60, 90, 40, 200, 70, 130])
prevalences_vraies = np.array([0.12, 0.25, 0.08, 0.40, 0.15, 0.50]) # inconnues de l'inference
se, sp = 0.92, 0.88 # sensibilite / specificite du test (specs constructeur, cf cellule precedente)
# Taux apparent reel (simule) = vraie prevalence corrigee par le test imparfait
p_apparent_vrai = prevalences_vraies * se + (1 - prevalences_vraies) * (1 - sp)
rng_obs = np.random.default_rng(7)
k_positifs = rng_obs.binomial(n_patients, p_apparent_vrai)
taux_apparent = k_positifs / n_patients
print("Taux de test positif OBSERVES (apparent, naif k/n) par clinique :", taux_apparent.round(3))
print("(ces taux apparents SURESTIMENT la vraie prevalence a cause des faux positifs)\n")
with pm.Model() as model_prevalence:
# Niveau population (parametrisation NON-CENTREE : evite le funnel, recette #3801)
mu_pop = pm.Normal("mu_pop", mu=-1.5, sigma=1.5) # logit-prevalence moyenne
tau = pm.HalfNormal("tau", sigma=1.0) # dispersion entre cliniques
z = pm.Normal("z", 0.0, 1.0, shape=K)
p_vraie = pm.Deterministic("p_vraie", pm.math.sigmoid(mu_pop + tau * z))
# Lien NON-LINEAIRE vraie prevalence -> taux apparent (casse la conjugaison)
p_apparent = pm.Deterministic(
"p_apparent", p_vraie * se + (1 - p_vraie) * (1 - sp)
)
# Vraisemblance binomiale sur les comptes de tests positifs observes
k = pm.Binomial("k", n=n_patients, p=p_apparent, observed=k_positifs)
trace_prev = pm.sample(
draws=2000, tune=1500, chains=4, target_accept=0.99,
random_seed=42, progressbar=False, return_inferencedata=True
)
# Diagnostics MCMC (honnetete des sorties)
div = int(trace_prev.sample_stats["diverging"].sum())
resume = az.summary(trace_prev, var_names=["mu_pop", "tau", "p_vraie"])
print("Inference hierarchique de la vraie prevalence (MCMC, test imparfait non-conjugue)")
print("=" * 64)
print(f"Divergences NUTS : {div} (cible = 0)")
print(f"R-hat max : {resume['r_hat'].max():.3f} (cible < 1.01)")
print(f"ESS bulk min : {int(resume['ess_bulk'].min())}")
print()
print("Comparaison : taux apparent (naif k/n) vs vraie prevalence inferiee (corrigee)")
print(f"{'clinique':>9} {'n':>5} {'vrai':>6} {'apparent':>9} {'vraie inferiee':>16} correction")
print("-" * 64)
p_post = resume.loc[[f"p_vraie[{i}]" for i in range(K)], "mean"].values
for i in range(K):
corr = abs(taux_apparent[i] - p_post[i])
marqueur = " <-- corrige (faux positifs retires)" if corr > 0.02 else ""
print(f"{i:>9} {n_patients[i]:>5} {prevalences_vraies[i]:>6.2f} {taux_apparent[i]:>9.3f} {p_post[i]:>16.3f}{marqueur}")
print()
print("Le moteur MCMC a inferer la VRAIE prevalence (etat maladie latent) en inversant")
print("le test imparfait, avec shrinkage entre cliniques -- ce qu'aucun calcul analytique")
print("(Beta-Binomial independant par clinique) ne peut produire : il donnerait seulement")
print("la distribution du taux apparent, sans correction ni emprunt d'information.)")Taux de test positif OBSERVES (apparent, naif k/n) par clinique : [0.233 0.378 0.225 0.45 0.3 0.592]
(ces taux apparents SURESTIMENT la vraie prevalence a cause des faux positifs)
Inference hierarchique de la vraie prevalence (MCMC, test imparfait non-conjugue)
================================================================
Divergences NUTS : 0 (cible = 0)
R-hat max : 1.000 (cible < 1.01)
ESS bulk min : 1791
Comparaison : taux apparent (naif k/n) vs vraie prevalence inferiee (corrigee)
clinique n vrai apparent vraie inferiee correction
----------------------------------------------------------------
0 60 0.12 0.233 0.169 <-- corrige (faux positifs retires)
1 90 0.25 0.378 0.317 <-- corrige (faux positifs retires)
2 40 0.08 0.225 0.174 <-- corrige (faux positifs retires)
3 200 0.40 0.450 0.405 <-- corrige (faux positifs retires)
4 70 0.15 0.300 0.230 <-- corrige (faux positifs retires)
5 130 0.50 0.592 0.568 <-- corrige (faux positifs retires)
Le moteur MCMC a inferer la VRAIE prevalence (etat maladie latent) en inversant
le test imparfait, avec shrinkage entre cliniques -- ce qu'aucun calcul analytique
(Beta-Binomial independant par clinique) ne peut produire : il donnerait seulement
la distribution du taux apparent, sans correction ni emprunt d'information.)Contrairement a Infer.NET qui genere des factor graphs, PyMC permet de visualiser le DAG (Directed Acyclic Graph) du modèle. Cette representation montre les dependances entre variables.
Le diagramme ci-dessous represente le reseau de decision medical (exemple fil rouge) avec la notation standardisee de Howard & Matheson (1984) rappelee en section 2 : les cercles sont des noeuds de chance (variables aleatoires), le rectangle est le noeud de decision, et l’hexagone est le noeud d’utilite. L’arc en pointilles Test -> Traiter est l’arc informationnel : le medecin observe le résultat du test avant de choisir le traitement.
flowchart TD
Mal(("Maladie"))
Test(("Test"))
Treat["Traiter"]
Sante(("Sante"))
Util{{"Utilite"}}
Mal --> Test
Mal --> Sante
Test -.-> Treat
Treat --> Sante
Sante --> Util
Treat -.-> Util
classDef chance fill:#e8f0fe,stroke:#1a73e8,stroke-width:2px,color:#174ea6
classDef decision fill:#fef7e0,stroke:#f9ab00,stroke-width:2px,color:#b06000
classDef utility fill:#e6f4ea,stroke:#188038,stroke-width:2px,color:#0d652d
class Mal,Test,Sante chance
class Treat decision
class Util utilityLecture.
Maladie(incertaine) influenceTest(observable) etSante(future) ;Traiterest decide en connaissantTest(arc en pointilles) ;Utilitecombine laSanteobtenue et le cout duTraiter. La cellule de code qui suit montre comment PyMC n’encode que la partie probabiliste (chance) de ce DAG – les decisions et l’utilite restent specifiees hors modèle, comme enumere dans la classeDecisionNodeplus haut.
# Visualisation du DAG du modele de diagnostic
import matplotlib.pyplot as plt
# Creer le modele pour la visualisation
with pm.Model() as diag_viz:
maladie = pm.Bernoulli("maladie", p=0.15)
p_test = pm.math.switch(maladie, 0.92, 0.12)
test_positif = pm.Bernoulli("test_positif", p=p_test)
# Afficher la structure du modele
print("Structure du DAG du modele de diagnostic :")
print(diag_viz)
print()
print("Variables :")
for rv in diag_viz.free_RVs:
print(f" {rv.name} ~ {rv.owner.op.name}")
print()
print("Structure du reseau de decision complet :")
print(" [Maladie] -----> [TestResult]")
print(" | |")
# v v
print(" <Traiter> ------> ((Utilite))")
print()
print("Dans le DAG PyMC :")
print(" - Noeuds = Variables aleatoires (maladie, test_positif)")
print(" - Arcs = Dependances probabilistes (maladie -> test_positif)")
print(" - Decisions = Enumerees manuellement hors du modele")Structure du DAG du modele de diagnostic :
<pymc.model.core.Model object at 0x000002902AE5EA50>
Variables :
maladie ~ bernoulli
test_positif ~ bernoulli
Structure du reseau de decision complet :
[Maladie] -----> [TestResult]
| |
<Traiter> ------> ((Utilite))
Dans le DAG PyMC :
- Noeuds = Variables aleatoires (maladie, test_positif)
- Arcs = Dependances probabilistes (maladie -> test_positif)
- Decisions = Enumerees manuellement hors du modeleModelisez le reseau de decision suivant :
Un etudiant doit decider s’il revise ou pas pour un examen. - Variable aleatoire : Difficulte de l’examen (Facile/Difficile) - Decision : Reviser (cout en temps) ou pas - Utilite : Note obtenue moins cout de revision
Rappel de la formule :
\[E[U(a)] = \sum_{s} P(s) \times U(a, s)\]
ou \(s \in \{\text{facile}, \text{difficile}\}\) et \(a \in \{\text{reviser}, \text{pas reviser}\}\).
# Exercice : Reseau de decision Etudiant
# Parametres
p_difficile = 0.4 # P(examen difficile)
# Notes esperees
note_revise_facile = 18
note_revise_difficile = 14
note_pas_revise_facile = 14
note_pas_revise_difficile = 8
# Cout de revision (en points d'utilite)
cout_revision = 2
# A completer 1 : Calculer l'utilite esperee de "reviser"
# Formule : E[U(reviser)] = P(facile) * (note_revise_facile - coutRevision) + P(difficile) * (note_revise_difficile - coutRevision)
# A completer 2 : Calculer l'utilite esperee de "ne pas reviser"
# Formule : E[U(pas reviser)] = P(facile) * note_pas_revise_facile + P(difficile) * note_pas_revise_difficile
# A completer 3 : Comparer et afficher la decision optimale
# Indice : utilisez un if/else sur les deux utilites esperees
# A completer 4 : Calculer le cout seuil (cout de revision au-dela duquel il ne faut plus reviser)
# Indice : resolvez E[U(reviser)] = E[U(pas reviser)] pour coutRevision
print("Exercice a completer")Exercice a completerAprès avoir implemente le code, verifiez :
Un comite doit choisir l’emplacement d’un nouvel aeroport parmi 3 sites candidats. Les critères d’evaluation sont incertains car bases sur des études preliminaires :
| Critere | Description | Incertitude |
|---|---|---|
| Couts | Construction, infrastructure | Depend du terrain, geologie |
| Securite | Trafic aerien, conditions meteo | Études locales necessaires |
| Nuisances | Bruit pour riverains, pollution | Impact environnemental |
Cet exemple combine Decision Networks et MAUT (Keeney & Raiffa, 1976) avec des attributs incertains.
# Reseau de Decision Multi-Attribut : Site d'Aeroport
# 3 sites candidats, 3 attributs incertains modelises avec scipy
n_sites = 3
sites = ["Site A (cotier)", "Site B (rural)", "Site C (periurbain)"]
# Priors sur les attributs bases sur etudes preliminaires
# Couts (millions EUR) - Gaussian
cout_mean = [500, 300, 400]
cout_std = [100, 71, 89] # sqrt(variance)
# Securite (score 0-1) - Beta
sec_alpha = [8, 5, 6]
sec_beta = [2, 5, 4]
# Nuisances (score 1-10) - Gaussian
nuis_mean = [3.0, 1.0, 5.0]
nuis_std = [1.0, 0.7, 1.4]
# Poids MAUT
w_cout = 0.35
w_securite = 0.40
w_nuisances = 0.25
print("=== Decision Multi-Attribut : Site d'Aeroport ===")
print(f"Poids : Cout={w_cout:.0%}, Securite={w_securite:.0%}, Nuisances={w_nuisances:.0%}\n")
# Calculer les esperances des distributions
print("Distributions des attributs par site :\n")
print(f"{'Site':<22} | {'Cout (MEUR)':>14} | {'Securite':>10} | {'Nuisances':>10}")
print("-" * 70)
eu = np.zeros(n_sites)
np.random.seed(42)
for s in range(n_sites):
# Esperances des distributions
e_cout = cout_mean[s]
e_sec = sec_alpha[s] / (sec_alpha[s] + sec_beta[s]) # Esperance Beta
e_nuis = nuis_mean[s]
print(f"{sites[s]:<22} | {e_cout:>7.0f} +/- {cout_std[s]:>4.0f} | {e_sec:>8.1%} | {e_nuis:>6.1f} +/- {nuis_std[s]:.1f}")
# Normalisation (0-1)
norm_cout = 1 - (e_cout - 300) / 200 # 300-500 MEUR -> 1-0
norm_sec = e_sec # Deja 0-1
norm_nuis = 1 - (e_nuis - 1) / 9 # 1-10 -> 1-0
eu[s] = w_cout * norm_cout + w_securite * norm_sec + w_nuisances * norm_nuis
print("\n=== Calcul de l'Utilite Esperee par Site ===\n")
for s in range(n_sites):
norm_cout = 1 - (cout_mean[s] - 300) / 200
norm_sec = sec_alpha[s] / (sec_alpha[s] + sec_beta[s])
norm_nuis = 1 - (nuis_mean[s] - 1) / 9
print(f"{sites[s]:<22} :")
print(f" v_cout={norm_cout:.3f}, v_sec={norm_sec:.3f}, v_nuis={norm_nuis:.3f}")
print(f" => EU = {eu[s]:.3f}")
# Decision optimale
best_idx = np.argmax(eu)
print(f"\n=== Decision Optimale : {sites[best_idx]} (EU = {eu[best_idx]:.3f}) ===")=== Decision Multi-Attribut : Site d'Aeroport ===
Poids : Cout=35%, Securite=40%, Nuisances=25%
Distributions des attributs par site :
Site | Cout (MEUR) | Securite | Nuisances
----------------------------------------------------------------------
Site A (cotier) | 500 +/- 100 | 80.0% | 3.0 +/- 1.0
Site B (rural) | 300 +/- 71 | 50.0% | 1.0 +/- 0.7
Site C (periurbain) | 400 +/- 89 | 60.0% | 5.0 +/- 1.4
=== Calcul de l'Utilite Esperee par Site ===
Site A (cotier) :
v_cout=0.000, v_sec=0.800, v_nuis=0.778
=> EU = 0.514
Site B (rural) :
v_cout=1.000, v_sec=0.500, v_nuis=1.000
=> EU = 0.800
Site C (periurbain) :
v_cout=0.500, v_sec=0.600, v_nuis=0.556
=> EU = 0.554
=== Decision Optimale : Site B (rural) (EU = 0.800) ===Plutot que d’utiliser uniquement les esperances, nous pouvons propager l’incertitude complete via Monte Carlo pour obtenir la distribution de l’utilite de chaque site.
# Simulation Monte Carlo pour propager l'incertitude
n_samples = 10_000
np.random.seed(42)
print("=== Simulation Monte Carlo (10 000 echantillons) ===\n")
eu_samples = np.zeros((n_sites, n_samples))
for s in range(n_sites):
# Echantillonner depuis les distributions
cout_samples = np.random.normal(cout_mean[s], cout_std[s], n_samples)
sec_samples = np.random.beta(sec_alpha[s], sec_beta[s], n_samples)
nuis_samples = np.random.normal(nuis_mean[s], nuis_std[s], n_samples)
# Normaliser
norm_cout = np.clip(1 - (cout_samples - 300) / 200, 0, 1)
norm_sec = np.clip(sec_samples, 0, 1)
norm_nuis = np.clip(1 - (nuis_samples - 1) / 9, 0, 1)
eu_samples[s] = w_cout * norm_cout + w_securite * norm_sec + w_nuisances * norm_nuis
# Probabilite que chaque site soit le meilleur
best_per_sample = np.argmax(eu_samples, axis=0)
p_best = [np.mean(best_per_sample == s) for s in range(n_sites)]
print("Probabilite d'etre le meilleur site :")
for s in range(n_sites):
mean_eu = eu_samples[s].mean()
std_eu = eu_samples[s].std()
print(f" {sites[s]:<22} : P(meilleur) = {p_best[s]:.1%}, EU = {mean_eu:.3f} +/- {std_eu:.3f}")
print(f"\nSite dominant : {sites[np.argmax(p_best)]} (P = {max(p_best):.1%})")=== Simulation Monte Carlo (10 000 echantillons) ===
Probabilite d'etre le meilleur site :
Site A (cotier) : P(meilleur) = 12.5%, EU = 0.584 +/- 0.112
Site B (rural) : P(meilleur) = 76.0%, EU = 0.742 +/- 0.095
Site C (periurbain) : P(meilleur) = 11.5%, EU = 0.554 +/- 0.138
Site dominant : Site B (rural) (P = 76.0%)| Site | Cout (normalise) | Securite | Nuisances | EU |
|---|---|---|---|---|
| A (cotier) | 0.0 (500M = max) | 0.80 | 0.78 | ~0.55 |
| B (rural) | 1.0 (300M = min) | 0.50 | 1.00 | ~0.78 |
| C (periurbain) | 0.5 | 0.60 | 0.56 | ~0.55 |
La simulation Monte Carlo revelle que Site B a environ 76% de chances d’etre le meilleur choix, confirmant la robustesse de la decision face a l’incertitude des attributs.
| Extension | Description |
|---|---|
| Analyse de sensibilite | Varier les poids pour voir si la decision change |
| Dominance stochastique | Comparer les distributions completes |
| Information supplementaire | Modeliser une étude de terrain comme source d’information |
| Distribution | Usage | Equivalent PyMC |
|---|---|---|
Bernoulli |
Variables binaires (maladie, test) | pm.Bernoulli |
Gaussian |
Attributs continus incertains (couts) | np.random.normal / pm.Normal |
Beta |
Scores bornes [0,1] (securite) | np.random.beta / pm.Beta |
| Pattern | Description | Section |
|---|---|---|
| Modèle conditionnel | pm.math.switch pour CPT |
7 |
| Observation | observed= pour conditionner |
7 |
| Enumeration | Boucle sur les observations possibles | 7 |
| Monte Carlo | Propagation d’incertitude | 8bis |
| Concept | Description | Section |
|---|---|---|
| Arcs informationnels | Ce que l’agent sait avant de decider | 3 |
| Backward induction | Resolution de la fin vers le debut | 4 |
| Valeur de l’information | Gain d’utilite grace a l’observation | 5 |
| Decisions séquentielles | Plusieurs decisions dans le temps | 6 |
| MAUT probabiliste | Multi-attributs avec incertitude | 8bis |
| Concept | Description |
|---|---|
| Reseau de decision | Extension des reseaux bayesiens avec decisions et utilites |
| Noeuds de decision | Choix contrôles par l’agent |
| Noeuds d’utilite | Consequences des choix |
| Arcs informationnels | Ce que l’agent sait au moment de decider |
| Politique | Fonction observations -> decisions |
| Backward induction | Resolution du dernier au premier noeud de decision |
| Si vous voulez… | Consultez… |
|---|---|
| Calculer la valeur de l’information | DecPyMC-5-Value-Information |
| Decisions robustes | DecPyMC-6-Expert-Systems |
| Decisions séquentielles (MDPs) | DecPyMC-7-Sequential |
Dans DecPyMC-5-Value-Information, nous verrons :
Conclusion : Ce notebook a presente les reseaux de decision (diagrammes d’influence), les arcs informationnels, la resolution par backward induction, et les decisions séquentielles avec PyMC.
Retour au sommaire : Index Probas
Navigation : << DecPyMC-3 | DecPyMC-5 >>