DecPyMC-1-Utility-Foundations : Axiomes et Fondements

Serie : Programmation Probabiliste avec PyMC (1/7)
Duree estimee : 50 minutes
Prerequis : Avoir explore les notebooks 1-8 (modèles probabilistes de base)

---

Objectifs

• Comprendre les loteries comme representation des choix en environnement stochastique
• Maitriser les axiomes de Von Neumann-Morgenstern
• Deriver la fonction d’utilite par calibration
• Comprendre la maximisation de l’utilite esperee comme fondement de l’agent rationnel

---

Navigation

Précédent	Suivant
PyMC-02b-Debugging-Python	DecPyMC-2-Utility-Money

---

1. Introduction : Pourquoi l’Utilite ?

Le problème de la decision sous incertitude

Jusqu’ici, nous avons utilise PyMC pour calculer des distributions de probabilite. Mais comment agir face a l’incertitude ?

Considerez ce dilemme :

Option	Description
A	Gagner 100 EUR avec certitude
B	50% de chance de gagner 200 EUR, 50% de ne rien gagner

Valeur esperee : E[A] = 100, E[B] = 100. Identiques !

Pourtant, la plupart des gens preferent A. Pourquoi ?

La valeur esperee ne suffit pas

Le Paradoxe de Saint-Petersbourg illustre cette limite. Formule par Nicolas Bernoulli en 1713 (lettre a Pierre Raymond de Montmort), il doit sa resolution a son cousin Daniel Bernoulli (1738) dans Specimen theoriae novae de mensura sortis (ou il introduit l’utilite logarithmique) :

Une piece est lancee jusqu’a obtenir Face. Si Face apparait au tour n, vous gagnez \(2^n\) euros. Combien paieriez-vous pour jouer ?

\[E[\text{gain}] = \sum_{n=1}^{\infty} \frac{1}{2^n} \times 2^n = \sum_{n=1}^{\infty} 1 = \infty\]

Valeur esperee infinie, mais personne ne paierait un million pour jouer. Il faut une autre approche.

2. Loteries : Representation Formelle des Choix

Definition

Une loterie est une distribution de probabilite sur des outcomes (résultats).

Notation : L = [p1:o1, p2:o2, …, pn:on]

Ou : - pi est la probabilite de l’outcome oi - la somme des pi = 1

Exemples

Loterie	Description
[1.0 : 100EUR]	100 EUR avec certitude (outcome degenere)
[0.5 : 200EUR, 0.5 : 0EUR]	Pile ou face pour 200 EUR
[0.6 : Succes, 0.4 : Echec]	Examen avec 60% de reussite

Loteries composees

Une loterie composee a d’autres loteries comme outcomes :

L = [0.5 : L1, 0.5 : L2]

Elle peut etre reduite en loterie simple par le calcul des probabilites totales.

Preparation de l’environnement

import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import List, Tuple, Generic, TypeVar
import warnings

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)

RANDOM_SEED = 42
rng = np.random.default_rng(RANDOM_SEED)

print(f"PyMC version : {pm.__version__}")
print(f"ArviZ version : {az.__version__}")
print("PyMC pret pour la theorie de la decision !")

PyMC version : 5.28.5
ArviZ version : 0.23.4
PyMC pret pour la theorie de la decision !

Packages charges

Package	Rôle
pymc	Modèles probabilistes, inference bayesienne
arviz	Diagnostics et visualisation des posteriors
numpy	Calculs numériques et generation de données
dataclasses	Structure de données pour les loteries

Les packages PyMC (Salvatier et al., 2016) et ArviZ (Kumar et al., 2019) fournissent l’infrastructure d’inference bayesienne. Ce notebook introduit les fondements de la théorie de la decision. Nous utiliserons PyMC pour l’inference bayesienne et des constructions Python pures pour les calculs d’utilite.

Implementation de la classe Loterie

Nous definissons une classe Loterie qui encapsule une distribution discrete sur des outcomes numériques. Cette implementation : - Valide que les probabilites somment a 1 - Fournit un affichage lisible de la forme [p1:o1, p2:o2, ...] - Permet le calcul de l’utilite esperee

# Representation d'une loterie simple
@dataclass
class Loterie:
    """Loterie : distribution de probabilite sur des outcomes numeriques."""
    outcomes: List[Tuple[float, float]]  # Liste de (probabilite, outcome)
    
    def __post_init__(self):
        total = sum(p for p, _ in self.outcomes)
        if abs(total - 1.0) > 1e-6:
            raise ValueError(f"Probabilites doivent sommer a 1, got {total}")
    
    def __str__(self):
        parts = [f"{p:.0%}:{o}" for p, o in self.outcomes]
        return f"[{', '.join(parts)}]"
    
    def expected_value(self):
        """Valeur esperee de la loterie."""
        return sum(p * o for p, o in self.outcomes)
    
    def expected_utility(self, utility_func):
        """Utilite esperee : E[U(X)]."""
        return sum(p * utility_func(o) for p, o in self.outcomes)


# Exemples
certainty = Loterie([(1.0, 100)])
coin_flip = Loterie([(0.5, 200), (0.5, 0)])

print(f"Certitude : {certainty}")
print(f"Pile ou face : {coin_flip}")
print(f"E[certitude] = {certainty.expected_value()}")
print(f"E[coin_flip] = {coin_flip.expected_value()}")

Certitude : [100%:100]
Pile ou face : [50%:200, 50%:0]
E[certitude] = 100.0
E[coin_flip] = 100.0

3. Axiomes de Préférences Rationnelles

Von Neumann et Morgenstern (1944) ont etabli les axiomes qui definissent un agent rationnel. Ces quatre axiomes ne sont pas arbitraires : ils capturent les conditions minimales pour qu’un système de préférences soit coherent et non-exploitable. Un agent qui les viole peut etre “money-pumped” (arnaque par cycle de préférences).

Notation

• A > B : l’agent prefere strictement A a B
• A < B : l’agent prefere strictement B a A
  
• A ~ B : l’agent est indifferent entre A et B
• A >= B : l’agent prefere A ou est indifferent (préférence faible)

Les 4 Axiomes

Axiome 1 : Completude (Orderability)

Pour toute paire de loteries A et B, exactement une relation est vraie : A > B, ou B > A, ou A ~ B

Interpretation : L’agent peut toujours comparer deux options. Sans cet axiome, l’agent pourrait etre paralyse face a certains choix.

Axiome 2 : Transitivite

Si A > B et B > C, alors A > C

Interpretation : Les préférences sont coherentes, pas de cycles. Si A > B > C > A, on peut vous faire payer pour echanger C->B, B->A, A->C et vous revenez au point de depart avec moins d’argent (“money pump”).

Axiome 3 : Continuite

Si A > B > C, alors il existe p dans (0,1) tel que : B ~ [p:A, (1-p):C]

Interpretation : Même le pire outcome peut etre compense par une probabilite suffisante du meilleur.

Axiome 4 : Indépendance

Si A > B, alors pour toute loterie C et toute probabilite p : [p:A, (1-p):C] > [p:B, (1-p):C]

Interpretation : Les préférences ne dependent pas des alternatives non pertinentes. C’est l’axiome le plus controverse : les expériences de Maurice Allais (1953) (Le comportement de l’homme rationnel devant le risque, Econometrica) ont montre que des agents rationnels le violent systematiquement, et Daniel Ellsberg (1961) a exhibit une aversion a l’ambiguite similaire. La Prospect Theory de Kahneman & Tversky (1979) proposera une refonte descriptive de ces comportements.

Demonstration : Verification de la transitivite

Illustrons l’axiome de transitivite avec un agent simule utilisant une fonction d’utilite concave U(x) = sqrt(x). Nous comparons trois loteries et verifions que les préférences sont transitives.

# Verification de la transitivite des preferences
# Agent avec fonction d'utilite concave : aversion au risque

def utilite_sqrt(x):
    """Fonction d'utilite concave U(x) = sqrt(x)."""
    return np.sqrt(x)

# Trois loteries
A = Loterie([(1.0, 100)])            # 100 certain
B = Loterie([(0.5, 196), (0.5, 0)])  # E[B] = 98
C = Loterie([(1.0, 49)])             # 49 certain

UA = A.expected_utility(utilite_sqrt)
UB = B.expected_utility(utilite_sqrt)
UC = C.expected_utility(utilite_sqrt)

print(f"U(A) = {UA:.3f} (100 certain)")
print(f"U(B) = {UB:.3f} (50% de 196)")
print(f"U(C) = {UC:.3f} (49 certain)")
print()

print(f"A > B ? {UA > UB}")
print(f"B ~ C ? {abs(UB - UC) < 1e-10} (U(B) = {UB}, U(C) = {UC})")
print(f"A > C ? {UA > UC} (transitivite verifiee)")

U(A) = 10.000 (100 certain)
U(B) = 7.000 (50% de 196)
U(C) = 7.000 (49 certain)

A > B ? True
B ~ C ? True (U(B) = 7.0, U(C) = 7.0)
A > C ? True (transitivite verifiee)

Exercice : Verifier l’axiome d’indépendance

Objectif : L’axiome d’indépendance dit que si A ≻ B, alors pour toute loterie C et tout α ∈ (0,1), αA + (1-α)C ≻ αB + (1-α)C.

On considere : - Loterie A : 80% chance de 500 EUR, 20% de 0 EUR - Loterie B : 30% chance de 1000 EUR, 70% de 0 EUR - Loterie C : 100% chance de 200 EUR

Travail : 1. Calculez l’utilite esperee de A et B avec u(x) = sqrt(x) (fonction CRRA, rho=0.5) 2. Verifiez que A ≻ B 3. Construisez les loteries composees αA + (1-α)C et αB + (1-α)C pour α = 0.5 4. Calculez l’utilite esperee de chaque loterie composee et verifiez que l’ordre est preserve

Indices : - La loterie composee αA + (1-α)C a pour outcomes : les outcomes de A ponderes par α, plus l’outcome de C pondere par (1-α) - αA + (1-α)C = [0.50.8 : 500, 0.50.2 : 0, 0.5*1.0 : 200]

# Exercice : Verifier l'axiome d'independance
A_outcomes = {500: 0.8, 0: 0.2}
B_outcomes = {1000: 0.3, 0: 0.7}
C_outcomes = {200: 1.0}
alpha = 0.5

# TODO etudiant : definissez u(x) = np.sqrt(x)
# u = lambda x: ...

# TODO etudiant : calculez E[u(A)] et E[u(B)], verifiez A > B
# eu_A = sum(p * u(o) for o, p in A_outcomes.items())
# eu_B = ...

# TODO etudiant : construisez les loteries composees
# composee_A = {o: alpha * p for o, p in A_outcomes.items()}
# composee_A[200] = composee_A.get(200, 0) + (1 - alpha) * 1.0
# composee_B = ...

# TODO etudiant : calculez E[u(composee_A)] et E[u(composee_B)]
# verifiez que l'ordre est preserve : composee_A > composee_B

print("Exercice a completer")

Exercice a completer

Analyse des résultats : Transitivite et Indifference

Loterie	Utilite esperee	Interpretation
A (100 certain)	U(A) = 10.000	sqrt(100) = 10
B (50% de 196)	U(B) = 7.000	0.5 x sqrt(196) + 0.5 x sqrt(0) = 7
C (49 certain)	U(C) = 7.000	sqrt(49) = 7

Observation cle : U(B) = U(C), donc B ~ C (indifference).

La transitivite est preservee car : - A > B (10 > 7) - B ~ C (7 = 7, indifference) - Donc A > C (10 > 7)

Note technique : Avec une utilite concave (sqrt), la loterie B a le même attrait que le montant certain C = 49, bien que E[B] = 98 >> 49. C’est exactement l’aversion au risque (formalisee par Pratt, 1964 et Arrow, 1965) : l’agent prefere 49 certain a une esperance de 98 avec risque.

4. Theoreme de Representation

Le theoreme fondamental

Theoreme (Von Neumann-Morgenstern, 1944)

Si les préférences d’un agent satisfont les 4 axiomes, alors il existe une fonction U (unique a transformation affine pres) telle que :

A > B ssi E[U(A)] > E[U(B)]

Ce theoreme est remarquable : des axiomes qualitatifs sur les préférences impliquent l’existence d’une fonction quantitative (utilite). C’est le pont entre psychologie (préférences) et mathematiques (optimisation).

Consequences

1. L’utilite esperee suffit : Pour prendre des decisions rationnelles, il suffit de maximiser E[U(outcome)]

2. Unicite relative : U est unique a transformation affine pres (U’ = aU + b, a > 0). Seules les différences d’utilite comptent.

3. Fondement de l’agent rationnel : Un agent qui viole ces axiomes peut etre “money-pumped”

L’agent rationnel fonde sur l’utilite

Observation -> [Inference Bayesienne] -> [Calcul E[U]] -> [argmax] -> Action
  (capteurs)     (PyMC)                 (pour chaque      (choix)
                                         action possible)

L’agent : 1. Observe le monde 2. Infere l’etat via inference bayesienne (PyMC) 3. Evalue l’utilite esperee de chaque action possible 4. Choisit l’action qui maximise E[U] 5. Agit dans le monde

Modelisation d’un agent rationnel avec PyMC

Scénario : Un patient a potentiellement une maladie. L’agent (medecin) doit decider : traiter ou ne pas traiter.

Le modèle bayesien calcule d’abord P(malade|test+), puis l’agent utilise ce posterior pour calculer l’utilite esperee de chaque action.

# Modele de diagnostic medical avec PyMC

prior_malade = 0.3
sensibilite = 0.9   # P(test+|malade)
specificite = 0.8   # P(test-|sain)

with pm.Model() as model_diagnostic:
    # Prior sur la maladie
    malade = pm.Bernoulli("malade", p=prior_malade)
    
    # Probabilite du test selon l'etat
    p_test = pm.math.switch(malade, sensibilite, 1 - specificite)
    
    # Observation : test positif
    test = pm.Bernoulli("test", p=p_test, observed=1)
    
    # Inference
    trace_diag = pm.sample(
        draws=5000, chains=4, random_seed=RANDOM_SEED,
        progressbar=False, return_inferencedata=True
    )

# Extraire le posterior
p_malade_post = float(trace_diag.posterior["malade"].mean().values)
print(f"Posterior P(malade|test+) = {p_malade_post:.1%}")

# Verification analytique
p_test_total = prior_malade * sensibilite + (1 - prior_malade) * (1 - specificite)
p_malade_analytique = (prior_malade * sensibilite) / p_test_total
print(f"Analytique P(malade|test+) = {p_malade_analytique:.1%}")
print(f"(Theoreme de Bayes : 0.30 x 0.90 / {p_test_total:.2f} = {p_malade_analytique:.3f})")

Posterior P(malade|test+) = 66.3%
Analytique P(malade|test+) = 65.9%
(Theoreme de Bayes : 0.30 x 0.90 / 0.41 = 0.659)

Interpretation : Le theoreme de Bayes en action

Le résultat P(malade|test+) ~ 65.9% illustre le theoreme de Bayes :

\[P(\text{malade}|\text{test}+) = \frac{P(\text{test}+|\text{malade}) \times P(\text{malade})}{P(\text{test}+)}\]

Terme	Valeur	Signification
P(malade)	0.30	Prior (prevalence)
P(test+ | malade)	0.90	Sensibilite
P(test+ | sain)	0.20	1 - Specificite
P(test+)	0.41	Probabilite totale

\[P(\text{malade}|\text{test}+) = \frac{0.90 \times 0.30}{0.41} = \frac{0.27}{0.41} \approx 0.659\]

Malgre un test positif avec 90% de sensibilite, la probabilite n’est “que” de 66%. C’est parce que la maladie est relativement rare (prior 30%) et que les faux positifs contribuent significativement.

Au-dela du conjugue : diagnostic hiérarchique multi-sites

Le diagnostic a 2 etats précédent est conjugue : le posterior a une forme close, et MCMC redonne exactement le calcul analytique de Bayes (les deux valeurs coincident a 66%). Pour montrer la vraie puissance de l’inference bayesienne par MCMC, considerons un problème non conjugue ou aucune forme close n’existe.

Scénario : un reseau de 8 hopitaux partage le même test, mais chaque hopital a sa propre prevalence de la maladie (liee a sa population). On observe le nombre de patients testes \(n_i\) et de cas positifs \(k_i\) par site. Objectif : inferer la prevalence de chaque site en empruntant l’information de tous les sites (partial pooling).

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)\]

Le lien logit et la structure hiérarchique cassent la conjugaison : MCMC est requis. Le benefice pedagogique concret est le shrinkage (retrecissement) : les sites peu observes sont tires vers la moyenne de population \(\mu_{pop}\), ce qu’une mise a jour Beta-Binomial independante par site ne produirait jamais.

Ce modèle est une instance du cadre hiérarchique bayesien classique (Gelman & Hill, Data Analysis Using Regression and Multilevel/Hierarchical Models, 2006). La parametrisation non-centree \(\mu_{pop} + \tau \cdot z_i\) (plutot que la centree) evite le funnel de Neal qui apparait quand \(\tau\) est petit.

# Diagnostic hierarchique multi-sites (cas NON-CONJUGUE : MCMC requis)
# 8 hopitaux, prevalences variables. On infer chaque prevalence avec partial pooling.
n_patients = np.array([20, 30, 15, 100, 25, 40, 12, 60])
k_positifs = np.array([3, 7, 2, 43, 8, 14, 2, 10])
empirique = k_positifs / n_patients
print("Prevalences empiriques par site :", empirique.round(3))

with pm.Model() as model_hierarchique:
    # 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 hopitaux
    z = pm.Normal("z", 0.0, 1.0, shape=len(n_patients))
    logit_p = mu_pop + tau * z
    p = pm.Deterministic("p", pm.math.sigmoid(logit_p))
    # Vraisemblance binomiale observee
    k = pm.Binomial("k", n=n_patients, p=p, observed=k_positifs)
    trace_hier = pm.sample(
        draws=2000, tune=1500, chains=4, target_accept=0.95,
        random_seed=RANDOM_SEED, progressbar=False, return_inferencedata=True
    )

# Diagnostics MCMC (honetete des sorties)
div = int(trace_hier.sample_stats["diverging"].sum())
resume = az.summary(trace_hier, var_names=["mu_pop", "tau", "p"])
print("Inference hierarchique des prevalences (MCMC, logit-Binomial non-centre)")
print("=" * 60)
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 : prevalence empirique vs posterior (shrinkage vers la moyenne)")
p_post = resume.loc[[f"p[{i}]" for i in range(len(n_patients))], "mean"].values
for i in range(len(n_patients)):
    diff = abs(p_post[i] - empirique[i])
    marqueur = "  <-- shrinkage (site petit tire vers mu_pop)" if (diff > 0.02 and n_patients[i] < 30) else ""
    print(f"  site {i} (n={n_patients[i]:3d}): empirique={empirique[i]:.3f} -> posterior={p_post[i]:.3f}{marqueur}")
print()
print("Partial pooling : les sites peu observes sont tires vers mu_pop,")
print("ce qu'aucune mise a jour conjuguee independante (Beta-Binomial par site) ne produit.")

Prevalences empiriques par site : [0.15  0.233 0.133 0.43  0.32  0.35  0.167 0.167]

Inference hierarchique des prevalences (MCMC, logit-Binomial non-centre)
============================================================
Divergences NUTS : 0   (cible = 0)
R-hat max         : 1.000   (cible < 1.01)
ESS bulk min      : 2697

Comparaison : prevalence empirique vs posterior (shrinkage vers la moyenne)
  site 0 (n= 20): empirique=0.150 -> posterior=0.203  <-- shrinkage (site petit tire vers mu_pop)
  site 1 (n= 30): empirique=0.233 -> posterior=0.243
  site 2 (n= 15): empirique=0.133 -> posterior=0.206  <-- shrinkage (site petit tire vers mu_pop)
  site 3 (n=100): empirique=0.430 -> posterior=0.400
  site 4 (n= 25): empirique=0.320 -> posterior=0.294  <-- shrinkage (site petit tire vers mu_pop)
  site 5 (n= 40): empirique=0.350 -> posterior=0.319
  site 6 (n= 12): empirique=0.167 -> posterior=0.223  <-- shrinkage (site petit tire vers mu_pop)
  site 7 (n= 60): empirique=0.167 -> posterior=0.193

Partial pooling : les sites peu observes sont tires vers mu_pop,
ce qu'aucune mise a jour conjuguee independante (Beta-Binomial par site) ne produit.

De l’inference a la decision

Maintenant, l’agent utilise le posterior pour prendre une decision. Nous definissons une matrice d’utilite qui specifie les consequences de chaque action dans chaque etat.

# Matrice d'utilite et decision optimale

# Matrice d'utilite (outcomes)
U_traiter_malade = 80      # Guerison (moins effets secondaires)
U_traiter_sain = 60        # Effets secondaires inutiles
U_pas_traiter_malade = 0   # Catastrophe
U_pas_traiter_sain = 100   # Parfait

p_malade = p_malade_post
p_sain = 1 - p_malade

# Utilite esperee de chaque action
EU_traiter = p_malade * U_traiter_malade + p_sain * U_traiter_sain
EU_pas_traiter = p_malade * U_pas_traiter_malade + p_sain * U_pas_traiter_sain

print(f"P(malade) = {p_malade:.1%}, P(sain) = {p_sain:.1%}")
print()
print(f"E[U(traiter)] = {p_malade:.2f} x {U_traiter_malade} + {p_sain:.2f} x {U_traiter_sain} = {EU_traiter:.1f}")
print(f"E[U(pas traiter)] = {p_malade:.2f} x {U_pas_traiter_malade} + {p_sain:.2f} x {U_pas_traiter_sain} = {EU_pas_traiter:.1f}")
print()

decision = "TRAITER" if EU_traiter > EU_pas_traiter else "NE PAS TRAITER"
print(f"Decision optimale : {decision}")

P(malade) = 66.3%, P(sain) = 33.7%

E[U(traiter)] = 0.66 x 80 + 0.34 x 60 = 73.3
E[U(pas traiter)] = 0.66 x 0 + 0.34 x 100 = 33.7

Decision optimale : TRAITER

Synthese : La matrice de decision medicale

Action	Etat : Malade (~66%)	Etat : Sain (~34%)	E[U(action)]
Traiter	U = 80 (guerison)	U = 60 (effets secondaires)	~73
Ne pas traiter	U = 0 (catastrophe)	U = 100 (parfait)	~34

\[E[U(\text{traiter})] = 0.66 \times 80 + 0.34 \times 60 \approx 73.2\] \[E[U(\text{ne pas traiter})] = 0.66 \times 0 + 0.34 \times 100 \approx 34.1\]

L’ecart de ~39 unites en faveur du traitement est substantiel. Le cout asymetrique des erreurs (0 vs 60) rend le traitement rationnel même si le patient pourrait etre sain.

5. Calibration par Mise a l’Indifference

Le problème pratique

Comment determiner la fonction d’utilite U d’un agent ? On ne peut pas observer U directement. La calibration permet de reconstruire U a partir des choix.

Méthode de l’indifference (Axiome de continuite)

Étape 1 : Fixer les references - o_best avec U(o_best) = 1 - o_worst avec U(o_worst) = 0

Étape 2 : Pour chaque outcome intermediaire - Option A : outcome o avec certitude - Option B : loterie [p : o_best, (1-p) : o_worst] - Faire varier p jusqu’a indifference

Étape 3 : En deduire l’utilite - U(o) = p* (la probabilite d’indifference)

Exemple : Evaluer une assurance

Vous avez une voiture de 10 000 EUR. Probabilite de vol = 5%.

Outcome	Valeur
Pas de vol	10 000 EUR
Vol	0 EUR

# Calibration de l'utilite pour la decision d'assurance

# Utilite logarithmique (introduite par Daniel Bernoulli, 1738 ; aversion au risque classique)
# U(x) = log(x + 1) normalise sur [0, 10000]
def utilite_log(x):
    """Utilite logarithmique normalisee."""
    if x <= 0:
        return 0
    return np.log(x + 1) / np.log(10001)  # Normalise: U(10000) ~ 1

# Sans assurance : loterie [0.95 : 10000, 0.05 : 0]
p_vol = 0.05
loterie_sans = Loterie([(1 - p_vol, 10000), (p_vol, 0)])
EU_sans = loterie_sans.expected_utility(utilite_log)

# Avec assurance : certitude de (10000 - 600) = 9400
prime_assurance = 600
EU_avec = utilite_log(10000 - prime_assurance)

print(f"E[U(sans assurance)] = 0.95 x U(10000) + 0.05 x U(0)")
print(f"                     = 0.95 x {utilite_log(10000):.4f} + 0.05 x {utilite_log(0):.4f}")
print(f"                     = {EU_sans:.4f}")
print()
print(f"E[U(avec assurance)] = U(9400) = {EU_avec:.4f}")
print()

if EU_avec > EU_sans:
    print(f"Decision : PRENDRE l'assurance (gain utilite = {EU_avec - EU_sans:.4f})")
else:
    print(f"Decision : NE PAS prendre l'assurance")

# Prime maximale acceptable
# U(10000 - p) = EU_sans => 10000 - p = exp(EU_sans * log(10001)) - 1
prime_max = 10000 - (np.exp(EU_sans * np.log(10001)) - 1)
print(f"\nPrime maximale acceptable : {prime_max:.0f} EUR")

E[U(sans assurance)] = 0.95 x U(10000) + 0.05 x U(0)
                     = 0.95 x 1.0000 + 0.05 x 0.0000
                     = 0.9500

E[U(avec assurance)] = U(9400) = 0.9933

Decision : PRENDRE l'assurance (gain utilite = 0.0433)

Prime maximale acceptable : 3691 EUR

Exercice : Modèle PyMC pour une loterie d’assurance

Objectif : L’exemple ci-dessus utilise un calcul analytique pour la decision d’assurance. On peut aussi modeliser la loterie directement dans PyMC.

On considere une maison de 200 000 EUR, avec P(inondation) = 2%. Si inondation, perte totale. L’assurance couvre la perte mais coute 5 000 EUR/an.

Travail : 1. Construisez un modèle PyMC avec une variable inondation ~ Bernoulli(0.02) et perte conditionnelle 2. Echantillonnez la distribution du gain net (sans assurance) : 0 si pas d’inondation, -200 000 si inondation 3. Calculez l’utilite esperee avec U(x) = log(x + 1) pour les deux options (assurance vs pas d’assurance) 4. Comparez et decidez

Indices : - Utilisez pm.math.switch pour le gain conditionnel - Richesse initiale = 200 000 EUR - Sans assurance : loterie [0.98 : 0, 0.02 : -200 000] - Avec assurance : certitude de -5 000 EUR (prime)

# Exercice : Modele PyMC pour une loterie d'assurance
richesse_maison = 200_000
p_inondation = 0.02
prime_assurance_maison = 5_000

# TODO etudiant : construisez un modele PyMC avec la variable inondation ~ Bernoulli(0.02)
# with pm.Model() as model_assurance:
#     inondation = pm.Bernoulli("inondation", p=p_inondation)
#     perte = pm.math.switch(inondation, -richesse_maison, 0.0)
#     gain_net = pm.Deterministic("gain_net", perte)
#     trace_assurance = pm.sample(..., chains=4)

# TODO etudiant : calculez E[U(sans assurance)] et E[U(avec assurance)]
# eu_sans = ...
# eu_avec = ...
# print(f"Decision : {'prendre assurance' if eu_avec > eu_sans else 'ne pas prendre assurance'}")

print("Exercice a completer : modele PyMC pour la loterie d'assurance")

Exercice a completer : modele PyMC pour la loterie d'assurance

Analyse du résultat : Pourquoi une prime si elevee ?

Metrique	Sans assurance	Avec assurance (600 EUR)
Valeur esperee	9 500 EUR	9 400 EUR
Utilite esperee	0.9500	0.9933

L’utilite logarithmique U(x) = log(x+1) est fortement concave : perdre les 10 000 EUR de la voiture serait catastrophique en utilite, tandis que payer 600 EUR de prime a un effet marginal. Jusqu’a combien l’agent est-il alors pret a payer pour eliminer le risque ? La reponse distingue deux grandeurs differentes, souvent confondues : la prime de risque de Pratt et la prime maximale d’assurance.

Cadre de calcul : richesse initiale \(w = 10\,000\) EUR ; sans assurance, le gain net vaut \(X = 0\) avec probabilite \(0{,}95\) et \(X = -10\,000\) avec probabilite \(0{,}05\), donc \(E[X] = -500\) EUR : la perte esperee vaut 500 EUR. L’equivalent certain de richesse finale, note \(CE\), verifie \(U(CE) = E[U(w+X)] = 0{,}95\), soit \(CE = 10001^{0{,}95} - 1 \approx 6\,309{,}17\) EUR.

1. Prime de risque de Pratt : le prix du risque seul

\[U(w + E[X] - \pi) = E[U(w + X)]\]

L’agent est indifferent entre garder la loterie et recevoir avec certitude \(w + E[X] - \pi\). La prime de risque mesure le cout du risque au-dela de la perte esperee. Comme \(U\) est croissante, \(w + E[X] - \pi = CE\), d’ou :

\[\pi = w + E[X] - CE = 10\,000 - 500 - 6\,309{,}17 = 3\,190{,}83 \text{ EUR}\]

C’est la forme classique de Pratt (1964) et Arrow (1965), \(\pi = E[X] - CE_{gain}\), avec l’equivalent certain exprime en gains : \(CE_{gain} = CE - w = -3\,690{,}83\) EUR.

2. Prime maximale d’assurance : le prix total du contrat

\[U(w - p_{max}) = E[U(w + X)]\]

Cette fois, l’agent compare la richesse certaine \(w - p_{max}\) — il paie la prime et l’assureur couvre le vol — a la loterie sans assurance. On resout le paiement total maximal accepte :

\[p_{max} = w - CE = 10\,000 - 6\,309{,}17 = 3\,690{,}83 \text{ EUR}\]

C’est exactement la quantite affichee par la cellule de calcul de l’assurance (“Prime maximale acceptable : 3691 EUR”, arrondie a l’unite).

Les deux grandeurs s’additionnent

Grandeur	Equation de definition	Valeur
Perte esperee	\(-E[X] = 0{,}05 \times 10\,000\)	500 EUR
Prime de risque (Pratt)	\(U(w + E[X] - \pi) = E[U(w+X)]\)	3 190,83 EUR
Prime maximale d’assurance	\(U(w - p_{max}) = E[U(w+X)]\)	3 690,83 EUR

\[p_{max} = \underbrace{500}_{\text{perte esperee}} + \underbrace{3\,190{,}83}_{\text{prime de risque}} = 3\,690{,}83 \text{ EUR}\]

Note technique : \(U(w - p_{max}) = E[U(w+X)]\) definit la prime d’assurance (paiement total, perte esperee comprise) ; \(U(w + E[X] - \pi) = E[U(w+X)]\) definit la prime de risque (cout du risque au-dela de la perte esperee). Un agent neutre au risque (\(\pi = 0\)) n’accepterait de payer que 500 EUR.

Cette decomposition explique pourquoi l’assurance peut creer un surplus pour un agent averse au risque : son prix de reserve, 3 690,83 EUR, depasse la perte esperee de 500 EUR de la prime de risque de 3 190,83 EUR. Le prix effectivement propose reste determine par le contrat et le marche.

6. Modelisation avec PyMC

Loteries comme modèles probabilistes

Une loterie peut etre naturellement modelisee comme une variable aleatoire dans PyMC. L’avantage est triple :

1. Inference : PyMC peut calculer les posteriors sur des variables latentes
2. Composition : Les loteries composees sont modelisees par des modèles graphiques
3. Integration : Le calcul de E[U] s’appuie sur les distributions inferees

# Loterie modelisee avec PyMC
# 50% de gagner 1000, 50% de perdre 500

with pm.Model() as model_loterie:
    # Variable latente : issue de la loterie
    gagne = pm.Bernoulli("gagne", p=0.5)
    
    # Gain conditionnel
    gain = pm.math.switch(gagne, 1000.0, -500.0)
    
    # Deterministic pour tracer le gain
    gain_det = pm.Deterministic("gain", gain)
    
    # Echantillonnage (pas d'observation = prior predictive)
    trace_loterie = pm.sample(
        draws=5000, chains=4, random_seed=RANDOM_SEED,
        progressbar=False, return_inferencedata=True
    )

# Analyser la distribution du gain
gains = trace_loterie.posterior["gain"].values.flatten()
print(f"Distribution du gain :")
print(f"  E[gain] = {gains.mean():.0f} EUR")
print(f"  P(gain = +1000) = {(gains == 1000).mean():.1%}")
print(f"  P(gain = -500) = {(gains == -500).mean():.1%}")
print()

# Calcul de l'utilite esperee avec U(x) = sqrt(x + 1000)
def U(x):
    return np.sqrt(x + 1000)

EU = 0.5 * U(1000) + 0.5 * U(-500)
U_certain = U(gains.mean())

print(f"E[U(loterie)] = 0.5 x U(1000) + 0.5 x U(-500)")
print(f"             = 0.5 x {U(1000):.2f} + 0.5 x {U(-500):.2f}")
print(f"             = {EU:.2f}")
print()
print(f"U(valeur esperee) = U(250) = {U_certain:.2f}")
print()
print(f"EU < U(EV) ? {EU < U_certain} => agent averse au risque prefere la certitude")

Distribution du gain :
  E[gain] = 254 EUR
  P(gain = +1000) = 50.3%
  P(gain = -500) = 49.7%

E[U(loterie)] = 0.5 x U(1000) + 0.5 x U(-500)
             = 0.5 x 44.72 + 0.5 x 22.36
             = 33.54

U(valeur esperee) = U(250) = 35.42

EU < U(EV) ? True => agent averse au risque prefere la certitude

Interpretation : Equivalent certain et aversion au risque

Concept	Valeur	Interpretation
E[gain]	250 EUR	Valeur esperee de la loterie
U(E[gain])	35.36	Utilite de la valeur esperee
E[U(gain)]	33.54	Utilite esperee de la loterie

Inegalite de Jensen (Jensen, 1906 ; appliquee a l’utilite par Pratt, 1964, et Arrow, 1965) : Pour une fonction d’utilite concave U,

\[E[U(X)] \leq U(E[X])\]

C’est la signature mathematique de l’aversion au risque.

Equivalent certain : Le montant CE tel que U(CE) = E[U(loterie)] :

\[U(CE) = 33.54 \Rightarrow \sqrt{CE + 1000} = 33.54 \Rightarrow CE \approx 125 \text{ EUR}\]

L’agent averse au risque est indifferent entre : - Recevoir 125 EUR avec certitude - Jouer la loterie [50% : +1000, 50% : -500] dont E = 250 EUR

7. Exemple guide : Calibrer Votre Fonction d’Utilite

Imaginez que vous devez choisir entre :

• Option A : Recevoir X EUR avec certitude
• Option B : 50% de chance de recevoir 1000 EUR, 50% de ne rien recevoir

Question : Pour quelle valeur de X etes-vous indifferent entre A et B ?

Votre X*	Profil
X* = 500	Neutre au risque
X* < 500	Averse au risque
X* > 500	Amateur de risque

# Exploration : calibrer le coefficient d'aversion au risque (CRRA)

# Fonction CRRA avec protection contre les valeurs extremes
EPSILON = 1.0  # Richesse minimale pour eviter U(0) = -inf

def crra(x, rho):
    """Utilite CRRA (Arrow 1965 ; Pratt 1964) : U(x) = x^(1-rho) / (1-rho), ou ln(x) si rho ~ 1.
    Modele canonical d'aversion relative constante au risque."""
    x = max(x, EPSILON)
    if abs(rho - 1.0) < 0.01:
        return np.log(x)
    return x ** (1 - rho) / (1 - rho)

def find_rho(ce, high=1000, low=EPSILON, p=0.5):
    """Recherche de rho par dichotomie pour satisfaire U(CE) = E[U(loterie)]."""
    rho_min, rho_max = 0.01, 10.0
    for _ in range(100):
        rho = (rho_min + rho_max) / 2
        U_ce = crra(ce, rho)
        EU_loterie = p * crra(high, rho) + (1 - p) * crra(low, rho)
        if U_ce > EU_loterie:
            rho_max = rho
        else:
            rho_min = rho
        if rho_max - rho_min < 0.0001:
            break
    return (rho_min + rho_max) / 2

# Equivalent certain : modifiez cette valeur selon votre preference
equivalent_certain = 400

rho = find_rho(equivalent_certain)

print(f"Equivalent certain : {equivalent_certain} EUR")
print(f"Loterie : 50% de 1000 EUR, 50% de {EPSILON} EUR")
print(f"Coefficient d'aversion au risque (rho) estime : {rho:.2f}")
print()

if rho < 0.5:
    print("Profil : Faible aversion au risque (proche de la neutralite)")
elif rho < 1.5:
    print("Profil : Aversion au risque moderee (profil equilibre)")
elif rho < 3.0:
    print("Profil : Aversion au risque significative (profil prudent)")
else:
    print("Profil : Tres forte aversion au risque (profil conservateur)")

# Verification
U_ce = crra(equivalent_certain, rho)
EU_loterie = 0.5 * crra(1000, rho) + 0.5 * crra(EPSILON, rho)
print(f"\nVerification de l'indifference :")
print(f"  U({equivalent_certain}) = {U_ce:.4f}")
print(f"  E[U(loterie)] = {EU_loterie:.4f}")
print(f"  Difference : {abs(U_ce - EU_loterie):.6f} (proche de 0)")

# Table de calibration
print("\n--- Exemples de calibration ---")
print(f"{'CE (EUR)':>10} | {'rho':>5} | Interpretation")
print("-" * 10 + "-+-" + "-" * 5 + "-+-" + "-" * 20)
for ce in [200, 300, 400, 450, 480]:
    r = find_rho(ce)
    interp = "Quasi-neutre" if r < 0.5 else "Modere" if r < 1.5 else "Prudent" if r < 3 else "Conservateur"
    print(f"{ce:>10} | {r:>5.2f} | {interp}")

Equivalent certain : 400 EUR
Loterie : 50% de 1000 EUR, 50% de 1.0 EUR
Coefficient d'aversion au risque (rho) estime : 0.25

Profil : Faible aversion au risque (proche de la neutralite)

Verification de l'indifference :
  U(400) = 119.4424
  E[U(loterie)] = 119.4390
  Difference : 0.003425 (proche de 0)

--- Exemples de calibration ---
  CE (EUR) |   rho | Interpretation
-----------+-------+---------------------
       200 |  0.61 | Modere
       300 |  0.44 | Quasi-neutre
       400 |  0.25 | Quasi-neutre
       450 |  0.14 | Quasi-neutre
       480 |  0.06 | Quasi-neutre

Analyse du résultat : Calibration CRRA

Equivalent certain	Coefficient rho	Interpretation
480 EUR	~0.06	Quasi-neutralite
400 EUR	~0.25	Faible aversion
300 EUR	~0.44	Aversion moderee
200 EUR	~0.61	Aversion significative

Experimentez : Modifiez equivalent_certain dans la cellule précédente : - CE = 200 : profil prudent (rho ~ 0.6) - CE = 480 : profil quasi-neutre (rho ~ 0.06) - CE = 500 : neutralite parfaite (rho = 0)

Valeurs typiques en economie empirique :

Population	rho estime
Investisseurs institutionnels	0.5 - 1.5
Menages americains	1.0 - 3.0
Decisions de sante	2.0 - 5.0

8. Exercice : Comparer CARA et CRRA

Objectif

Implementer et comparer deux familles de fonctions d’utilite pour une même loterie :

• CARA (Constant Absolute Risk Aversion, Arrow 1965 ; Pratt 1964) : U(x) = -exp(-alpha * x)
• CRRA (Constant Relative Risk Aversion) : U(x) = x^(1-rho) / (1-rho)

Loterie

50% de chance de gagner 2000 EUR, 50% de perdre 500 EUR (richesse initiale = 10 000 EUR).

Travail a realiser

1. Implementer U_CARA(x, alpha) = -exp(-alpha * x) et calculer CE pour alpha = 0.001
2. Implementer U_CRRA(x, rho) = x^(1-rho)/(1-rho) et calculer CE pour rho = 2.0
3. Calculer les primes de risque : Prime = E[gain] - CE
4. Comparer et interpreter
5. Bonus : Tester avec richesse = 1000 EUR et observer que la prime CARA ne change pas mais CRRA change

# Exercice : Comparer les primes de risque CARA vs CRRA
richesse_initiale = 10000.0
gain_haut = 2000.0
perte_basse = -500.0
proba = 0.5

valeur_esperee = proba * gain_haut + (1 - proba) * perte_basse
print(f"Loterie : {proba:.0%} de +{gain_haut:.0f} EUR, {(1-proba):.0%} de {perte_basse:.0f} EUR")
print(f"Valeur esperee du gain : {valeur_esperee:.0f} EUR")
print(f"Richesse initiale : {richesse_initiale:.0f} EUR")
print()

# A completer 1 : Implementer la fonction d'utilite CARA
# Hint : def u_cara(x, alpha): return -np.exp(-alpha * x)

# A completer 2 : Implementer la fonction d'utilite CRRA
# Hint : utiliser la fonction crra() definie plus haut

# A completer 3 : Calculer l'equivalent certain CE pour chaque fonction
# L'equivalent certain CE verifie : U(richesse + CE) = E[U(richesse + gain)]
# Methode : recherche par dichotomie dans l'intervalle [perte_basse - 100, gain_haut + 100]

# A completer 4 : Calculer et afficher les primes de risque
# Prime de risque = valeur_esperee - CE

# A completer 5 (BONUS) : Repeter avec richesse_initiale = 1000 EUR
# Observer : la prime CARA change-t-elle ? La prime CRRA change-t-elle ?

print("Exercice a completer : fonctions CARA et CRRA, calcul des equivalents certains et primes de risque.")
print("Indices : CARA = aversion absolue constante, CRRA = aversion relative constante.")

Loterie : 50% de +2000 EUR, 50% de -500 EUR
Valeur esperee du gain : 750 EUR
Richesse initiale : 10000 EUR

Exercice a completer : fonctions CARA et CRRA, calcul des equivalents certains et primes de risque.
Indices : CARA = aversion absolue constante, CRRA = aversion relative constante.

9. Resume

Concept	Description
Loterie	Distribution de probabilite sur des outcomes
Axiomes VNM	Completude, Transitivite, Continuite, Indépendance
Theoreme de representation	Existence d’une fonction U telle que préférences = max E[U]
Agent rationnel	Maximise l’utilite esperee de ses actions
Calibration	Méthode de l’indifference pour determiner U
Aversion au risque	E[U(X)] < U(E[X]) pour U concave (Jensen)

Distributions utilisees dans ce notebook

Distribution	Paramètres	Usage typique
Bernoulli	(p,)	Variables binaires (malade/sain, gain/perte)
Deterministic	(expr,)	Calculs derivés (gain conditionnel)

---

Pour aller plus loin

Si vous voulez…	Consultez…
Approfondir l’aversion au risque	DecPyMC-2-Utility-Money
Decisions multi-critères	DecPyMC-3-Multi-Attribute
Visualiser les reseaux de decision	DecPyMC-4-Decision-Networks

References

Théorie de la decision et utilite esperee - Bernoulli, D. (1738) : Specimen Theoriae Novae de Mensura Sortis (paradoxe de Saint-Petersbourg, utilite logarithmique) - Bernoulli, N. (1713) : lettre a Montmort (formulation initiale du paradoxe) - Von Neumann, J. & Morgenstern, O. (1944) : Theory of Games and Economic Behavior (axiomes VNM, theoreme de representation) - Savage, L. J. (1954) : The Foundations of Statistics (axiomatique subjective, sure-thing principle) - Allais, M. (1953) : Le comportement de l’homme rationnel devant le risque, Econometrica (paradoxe d’Allais) - Ellsberg, D. (1961) : Risk, Ambiguity, and the Savage Axioms, QJE (paradoxe d’Ellsberg) - Pratt, J. W. (1964) : Measuring Risk Aversion, Annals of Mathematical Statistics (prime de risque, aversion absolue/relative) - Arrow, K. J. (1965) : Aspects of the Theory of Risk-Bearing (CRRA, CARA) - Kahneman, D. & Tversky, A. (1979) : Prospect Theory, Econometrica (aversion aux pertes)

Implementation - Salvatier, J., Wiecki, T. V. & Fonnesbeck, C. (2016) : Probabilistic programming in Python using PyMC3, PeerJ Computer Science - Kumar, R. et al. (2019) : ArviZ: a unified library for exploratory analysis of Bayesian models in Python, JOSS - Russell, S. & Norvig, P. : Artificial Intelligence: A Modern Approach, chapitre 16 (Making Simple Decisions)

---

Conclusion : Ce notebook a etabli les fondements de la théorie de la decision : loteries, axiomes VNM, theoreme de representation, et calibration de la fonction d’utilite par la méthode de l’indifference.

Retour au sommaire : Index Probas

Navigation : << PyMC-02b | DecPyMC-2 >>