Travaux pratiques avec Python

Cours complet · mathématiques approfondies (ECG 2e année), chapitre 13 · prépa ECG, 2e année

Travailler ce chapitre sur Adloun Exercices corrigés de ce chapitre

Les douze chapitres précédents ont construit des objets et démontré des théorèmes. Ce dernier chapitre les met en machine. Non pour vérifier ce qui est déjà démontré — une simulation ne démontre rien — mais pour voir ce que les énoncés disent, et pour traiter les cas où le calcul exact n'aboutit pas.

Le programme découpe ces travaux en quatre thèmes, que ce chapitre suit dans l'ordre : statistiques descriptives bivariées, fonctions de plusieurs variables, simulation de lois, estimation.

AttentionCe qu'un TP n'est pas

Le programme est explicite : « L'objectif de ces travaux pratiques n'est pas l'écriture de longs programmes mais l'assimilation de savoir-faire ». On n'attend donc de vous ni une bibliothèque, ni une interface : quinze lignes qui répondent à une question précise, et une phrase qui interprète le résultat.

Autre point de règlement, souvent mal connu : « On rappellera dans les sujets toutes les syntaxes des commandes non exigibles. » Les seules commandes à savoir par cœur sont celles de première année. Tout le reste — et ce chapitre en emploie — vous sera donné le jour de l'épreuve.

13.1 Thème 1 : statistiques descriptives bivariées

13.1.1 Le nuage et son point moyen

Définition 13.1Série double, nuage, point moyen

Une série statistique double est une liste de couples $(x_{1}, y_{1}), \dots, (x_{n}, y_{n})$ relevée sur $n$ individus. Le nuage de points est l'ensemble de ces couples dans le plan.

Le point moyen du nuage est $G = (x, y)$ , où $x = \frac{1}{n} \sum x_{i}$ et $y = \frac{1}{n} \sum y_{i}$ .

Définition 13.2Covariance et corrélation empiriques

La covariance empirique de la série est

$s_{x y} = \frac{1}{n} i = 1 \sum n (x_{i} - x) (y_{i} - y) = \frac{1}{n} i = 1 \sum n x_{i} y_{i} - x y,$

et le coefficient de corrélation empirique est $r = \frac{s _{x y}}{s _{x} s _{y}}$ , où $s_{x}$ et $s_{y}$ sont les écarts types des deux séries simples.

ImportantCe sont les formules du chapitre 8, appliquées à des données

La covariance empirique est la covariance du couple $(X, Y)$ pour la loi uniforme sur les $n$ observations. La formule de Huygens et l'encadrement $∣ r ∣ \leq 1$ ne sont donc pas des résultats nouveaux : ce sont ceux du cours, lus sur cette loi-là. En particulier $∣ r ∣ = 1$ si et seulement si les points sont exactement alignés, ce qui est le cas d'égalité de Cauchy-Schwarz du chapitre 4.

Prenons huit entreprises d'un même secteur, dont on relève le budget publicitaire $x$ et le chiffre d'affaires $y$ , tous deux en milliers d'euros.

$x$ (budget)	12	18	24	30	36	42	48	54
$y$ (chiffre)	46	52	61	64	73	77	86	89

Python : Le nuage, la covariance, la droite


import numpy as np
import matplotlib.pyplot as plt

x = np.array([12, 18, 24, 30, 36, 42, 48, 54])
y = np.array([46, 52, 61, 64, 73, 77, 86, 89])

xb, yb = x.mean(), y.mean()
cov = ((x - xb)*(y - yb)).mean()
r = cov/(x.std()*y.std())
a = cov/x.var()
b = yb - a*xb
print(round(cov, 3), round(r, 4), round(a, 4), round(b, 4))
# 198.0  0.9956  1.0476  33.9286

plt.scatter(x, y)
plt.plot(x, a*x + b)
plt.plot(xb, yb, "o")        # le point moyen
plt.show()

Trois lignes de calcul, une ligne de tracé. Le .std() et le .var() de numpy divisent par $n$ , ce qui est bien la convention de la définition ci-dessus.

13.1.2 Il y a deux droites de régression

◆Théorème 13.3Les deux droites, et ce que leur écart mesure

La droite de régression de $y$ en $x$ est $y = a x + b$ avec $a = \frac{s _{x y}}{s _{x}^{2}}$ et $b = y - a x$ ; celle de $x$ en $y$ est $x = a^{'} y + b^{'}$ avec $a^{'} = \frac{s _{x y}}{s _{y}^{2}}$ et $b^{'} = x - a^{'} y$ .

Ces deux droites passent par le point moyen $G$ , et

$a a^{'} = \frac{s _{x y}^{2}}{s _{x}^{2} s _{y}^{2}} = r^{2} .$

Démonstration

Les deux droites passent par $G$ : c'est immédiat sur $y = a (x - x) + y$ , et sur $x = a^{'} (y - y) + x$ . Le produit $a a^{'}$ est celui des deux quotients, soit $s_{x y}^{2} / (s_{x}^{2} s_{y}^{2}) = r^{2}$ .

ImportantChoisir laquelle, c'est choisir ce qu'on explique

Les deux droites ne coïncident que si $∣ r ∣ = 1$ . Ajuster $y$ en $x$ , c'est minimiser les écarts verticaux : on suppose que $x$ est connu et que $y$ est à expliquer. Ajuster $x$ en $y$ minimise les écarts horizontaux et inverse les rôles.

Si l'on veut prédire le chiffre d'affaires à partir du budget, c'est la première qu'il faut. Prendre l'autre parce qu'elle est déjà tracée est une faute de modèle, pas une approximation.

13.1.3 Redresser avant d'ajuster

Une droite ne résume bien qu'un nuage droit. Quand la croissance est multiplicative, on pré-transforme les données pour se ramener au cas linéaire. Voici le nombre d'abonnés d'un service, en milliers, sur huit trimestres.

$t$	0	1	2	3	4	5	6	7
$y$	12	19	31	47	78	121	195	310

ImportantLe passage au logarithme change le modèle, pas les données

Si $y = α e^{β t}$ , alors $ln y = ln α + β t$ : la série $(t, ln y)$ est exactement alignée. Ajuster une droite sur les logarithmes revient donc à ajuster une exponentielle sur les données.

Le coefficient $r$ calculé après transformation ne mesure plus la même chose : il juge la qualité du modèle exponentiel, pas celle d'un lien linéaire. Ici il passe de $0, 915$ à $0, 9999$ — et c'est le second chiffre qui dit que le modèle est le bon.

Python : Redresser un nuage par le logarithme


import numpy as np

t = np.array([0, 1, 2, 3, 4, 5, 6, 7])
y = np.array([12, 19, 31, 47, 78, 121, 195, 310])
L = np.log(y)

beta = ((t - t.mean())*(L - L.mean())).mean()/t.var()
alpha = np.exp(L.mean() - beta*t.mean())
print(round(alpha, 3), round(beta, 4))       # 11.996  0.4642
print(round(np.log(2)/beta, 3))  # 1.493 : doublement en 1,5 trimestre
print(round(alpha*np.exp(8*beta)))   # 492 : prevision au 9e trimestre

Le modèle s'écrit $y ≃ 12, 0 e^{0, 464 t}$ . Le temps de doublement se lit directement : $ln 2/ β ≃ 1, 49$ trimestre.

AttentionUne corrélation forte ne désigne aucune cause

$r = 0, 9956$ dit que les deux séries varient ensemble. Il ne dit pas que la publicité fait le chiffre d'affaires : les entreprises qui vendent le plus sont aussi celles qui peuvent dépenser le plus. Un troisième facteur — la taille de l'entreprise — explique peut-être les deux.

Le programme demande explicitement un regard critique sur ces méthodes. C'est ici qu'il s'exerce : la statistique descriptive décrit, elle n'explique pas.

13.2 Thème 2 : fonctions de plusieurs variables

13.2.1 Graphe, lignes de niveau, gradient

Python : Les trois représentations d'une même fonction


import numpy as np
import matplotlib.pyplot as plt

f = lambda x, y: x**2 + 2*y**2
X, Y = np.meshgrid(np.linspace(-3, 3, 200), np.linspace(-3, 3, 200))
Z = f(X, Y)

fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection="3d")
ax.plot_surface(X, Y, Z)                    # le graphe
ax2 = fig.add_subplot(1, 2, 2)
ax2.contour(X, Y, Z, levels=[2, 6, 12])     # les lignes de niveau
ax2.quiver(2, 1, 4, 4)                      # le gradient en (2, 1)
ax2.set_aspect("equal")
plt.show()

meshgrid fabrique la grille des couples $(x, y)$ ; contour trace les lignes de niveau demandées ; quiver pose une flèche. Ces trois syntaxes ne sont pas exigibles et vous seraient rappelées.

◆Théorème 13.4Le gradient est orthogonal à la ligne de niveau

Soit $f$ de classe $C^{1}$ sur un ouvert de $R^{2}$ et $A = (x_{0}, y_{0})$ un point où $\nabla f (A) \neq = 0$ . Si $t \mapsto γ (t)$ est un arc de classe $C^{1}$ tracé sur la ligne de niveau de $f$ passant par $A$ , avec $γ (0) = A$ , alors

$⟨ \nabla f (A), γ^{'} (0)⟩ = 0.$

Démonstration

La fonction $t \mapsto f (γ (t))$ est constante sur la ligne de niveau. Sa dérivée en $0$ est donc nulle ; or la règle de dérivation composée du chapitre 5 donne $⟨ \nabla f (A), γ^{'} (0)⟩$ .

13.2.2 Le graphe et son plan tangent

Le développement limité d'ordre $2$ du chapitre 10 dit que, près d'un point critique $A$ ,

$f (A + H) - f (A) = \frac{1}{2}^{t} H \nabla^{2} f (A) H + o (∥ H ∥^{2}),$

donc que la position du graphe par rapport à son plan tangent est gouvernée par le signe de la forme quadratique hessienne — c'est-à-dire par les signes des valeurs propres de la matrice hessienne, qui est symétrique réelle et donc diagonalisable en base orthonormée (chapitre 9).

Python : Classer un point critique par les valeurs propres


import numpy as np

def hessienne(x, y):
    # pour f = x^3 - 3xy + y^3
    return np.array([[6*x, -3.0], [-3.0, 6*y]])

for (x, y) in [(0, 0), (1, 1)]:
    vp = np.linalg.eigvalsh(hessienne(x, y))      # matrice symetrique
    print((x, y), vp.round(4), "selle" if vp[0]*vp[-1] < 0
          else ("minimum" if vp[0] > 0 else "maximum"))
# (0, 0) [-3.  3.] selle
# (1, 1) [ 3.  9.] minimum

eigvalsh est réservé aux matrices symétriques et rend les valeurs propres triées. Le test de signe suffit alors : deux signes opposés donnent un point selle, deux valeurs propres strictement positives un minimum local.

ImportantLe programme dit « lien avec les valeurs propres », et rien d'autre

On classe un point critique par les signes des valeurs propres de la hessienne. Aucune autre quantité associée à la matrice n'intervient, et il n'est pas question d'en calculer une pour obtenir ces signes : ici on les lit, sur machine, avec eigvalsh ; à la main, on cherche un polynôme annulateur comme au chapitre 2.

Le point $(0, 0)$ mérite qu'on regarde ce qui s'y passe. La fonction y vaut $0$ et son gradient y est nul : le plan tangent est le plan $z = 0$ . Deux coupes suffisent alors à conclure.

13.2.3 Extrema sous contrainte linéaire

ImportantUne contrainte linéaire se substitue

Sous une contrainte affine, on exprime une variable en fonction des autres et l'on reporte : le problème contraint à $n$ variables devient un problème libre à $n - 1$ variables, que l'on traite avec les outils du chapitre 10. C'est la seule méthode au programme, et elle suffit.

Cherchons le minimum de $f (x, y) = x^{2} + 2 y^{2}$ sous la contrainte $x + y = 3$ . En posant $y = 3 - x$ , on étudie $g (x) = x^{2} + 2 (3 - x)^{2} = 3 x^{2} - 12 x + 18$ , dont la dérivée $6 x - 12$ s'annule en $x = 2$ . Comme $g^{''} = 6 > 0$ , c'est un minimum, et il vaut $g (2) = 6$ , atteint en $(2, 1)$ .

Python : Vérifier le minimum contraint par balayage


import numpy as np

f = lambda x, y: x**2 + 2*y**2
xs = np.linspace(-2, 6, 100001)
vals = f(xs, 3 - xs)                 # la contrainte est deja injectee
i = vals.argmin()
print(round(xs[i], 5), round(3 - xs[i], 5), round(vals[i], 6))
# 2.0  1.0  6.0

Le balayage corrobore le calcul, il ne le remplace pas : il ne visite qu'une grille, et rien ne garantit qu'il n'existe pas ailleurs un point plus bas. Ici c'est le calcul de dérivée, lui, qui le garantit.

13.3 Thème 3 : simulation de lois

13.3.1 La méthode d'inversion

◆Théorème 13.5Méthode d'inversion

Soit $F$ la fonction de répartition d'une variable à densité, continue et strictement croissante sur un intervalle $I$ où elle prend ses valeurs dans $] 0, 1 [$ . Si $U$ suit la loi uniforme sur $] 0, 1 [$ , alors $X = F^{- 1} (U)$ admet $F$ pour fonction de répartition.

Démonstration

Pour $x \in I$ , l'événement $[F^{- 1} (U) \leq x]$ est égal à $[U \leq F (x)]$ par croissance stricte de $F$ et de sa réciproque. Or $F (x) \in] 0, 1 [$ , donc $P (U \leq F (x)) = F (x)$ . La fonction de répartition de $X$ est donc $F$ .

Pour la loi exponentielle de paramètre $λ$ , $F (x) = 1 - e^{- λ x}$ sur $R_{+}$ , d'où $F^{- 1} (u) = - \frac{ln ( 1 - u )}{λ}$ .

Python : Simuler une loi exponentielle sans numpy


import random, math

random.seed(4)
lam, N = 1.5, 200000
ech = [-math.log(1 - random.random())/lam for _ in range(N)]
print(round(sum(ech)/N, 4))                          # 0.6671    1/lam
# 0.7763    1 - e^{-1.5}
print(round(sum(1 for v in ech if v <= 1)/N, 4))

La moyenne empirique approche $1/ λ = 0, 6667$ et la fréquence de $[X \leq 1]$ approche $1 - e^{- 1, 5} = 0, 7769$ . Deux contrôles valent mieux qu'un : le premier vérifie l'échelle, le second la forme.

13.3.2 Une loi sans espérance : Cauchy

La loi de Cauchy a pour densité $f (x) = \frac{1}{π ( 1 + x ^{2} )}$ sur $R$ , donc pour fonction de répartition $F (x) = \frac{1}{2} + \frac{arctan x}{π}$ et pour réciproque $F^{- 1} (u) = tan (π (u - \frac{1}{2}))$ .

ImportantPourquoi elle n'a pas d'espérance

L'intégrale $\int_{0}^{+ \infty} \frac{x}{π ( 1 + x ^{2} )} d x$ vaut $[\frac{l n ( 1 + x ^{2} )}{2 π}]_{0}^{+ \infty}$ , qui diverge. La condition d'existence de l'espérance n'est pas remplie ; il n'y a donc aucun nombre vers lequel la moyenne empirique aurait à converger. La loi faible des grands nombres ne s'applique pas : son hypothèse manque.

AttentionCe que cette figure prouve, et ce qu'elle ne prouve pas

Elle ne démontre pas que l'espérance n'existe pas — une simulation ne démontre rien. C'est le calcul de l'intégrale divergente qui le démontre.

Ce qu'elle montre, en revanche, est précieux : le symptôme de cette absence, qu'on saura reconnaître ailleurs. Une moyenne empirique qui refuse de se stabiliser quand $n$ est déjà grand doit faire soupçonner une loi à queue lourde, et non un manque de tirages.

13.3.3 Trois façons de simuler une loi géométrique

Python : Bernoulli et while, exponentielle et floor, numpy


import numpy as np
import random, math

p, N = 0.3, 200000
rng = np.random.default_rng(5)

def par_bernoulli(p):               # methode 1 : on compte les essais
    k = 1
    while random.random() >= p:
        k += 1
    return k

ech1 = [par_bernoulli(p) for _ in range(N)]
# methode 2
ech2 = np.floor(rng.exponential(1/(-math.log(1 - p)), N)) + 1
ech3 = rng.geometric(p, N)                                 # methode 3

for e in (ech1, ech2, ech3):
    e = np.array(e)
    print(round(float(e.mean()), 3), round(float(np.mean(e == 1)), 3))
#  3.333 et  0.300 pour les trois : moyenne 1/p, et P(X = 1) = p

La méthode 2 s'explique : si $E$ suit $E (μ)$ avec $μ = - ln (1 - p)$ , alors $P (⌊ E ⌋ \geq k) = e^{- μ k} = (1 - p)^{k}$ , qui est exactement $P (X > k)$ pour $X$ géométrique. C'est la même loi.

ImportantTrois méthodes, trois coûts

La première ne demande qu'une loi de Bernoulli, mais sa durée est aléatoire : en moyenne $1/ p$ tours de boucle, et beaucoup plus quand $p$ est petit. La deuxième et la troisième rendent un tirage en temps constant. Comparer des méthodes de simulation — pas seulement les écrire — fait partie des compétences attendues.

13.3.4 Une loi normale par le théorème limite central

Si $U_{1}, \dots, U_{12}$ sont indépendantes de loi uniforme sur $[0, 1]$ , chacune a pour espérance $1/2$ et pour variance $1/12$ . La somme $S = \sum U_{i}$ vérifie donc $E (S) = 6$ et $V (S) = 1$ : la variable $S - 6$ est centrée réduite, et le théorème limite central du chapitre 11 suggère qu'elle est presque normale.

AttentionLa méthode des douze uniformes ment sur les extrêmes

$P (∣ S - 6∣ > 6) = 0$ exactement, alors qu'une vraie normale centrée réduite dépasse $6$ avec probabilité $2 \times 1 0^{- 9}$ . Sur $300000$ tirages, le plus extrême observé vaut $- 4, 27$ .

Ce n'est négligeable que si l'on s'intéresse au centre. Pour un calcul de risque, qui vit précisément dans la queue, la méthode est fausse — et numpy.random.standard_normal, qui n'a pas cette borne, est alors le bon outil. Choisir sa méthode de simulation en fonction de la question posée est, là encore, une compétence du programme.

13.4 Thème 4 : Monte-Carlo, estimation, intervalles

13.4.1 Le principe et sa garantie

Définition 13.6Méthode de Monte-Carlo

Estimer une quantité $θ = E (φ (X))$ en simulant un $n$ -échantillon $(X_{1}, \dots, X_{n})$ de la loi de $X$ et en retenant

$θ_{n} = \frac{1}{n} i = 1 \sum n φ (X_{i}),$

c'est appliquer la méthode de Monte-Carlo. Le cas où $φ$ est la fonction indicatrice d'un événement donne l'estimation d'une probabilité par une fréquence.

◆Théorème 13.7Garantie d'approximation

Si $φ (X)$ admet une variance $σ^{2}$ , alors $θ_{n}$ est un estimateur sans biais de $θ$ , il est convergent (loi faible des grands nombres), et le théorème limite central donne, pour $n$ grand,

$P (θ_{n} - θ \leq \frac{1 , 96 σ}{n}) ≃ 0, 95.$

ImportantL'erreur décroît en

1/ n

, et cela coûte cher

Le $n$ du dénominateur est la loi d'airain de la méthode : pour diviser l'erreur par $10$ , il faut multiplier le nombre de tirages par $100$ . Aucune astuce de programmation n'y change rien — seule une réduction de $σ$ le peut.

C'est aussi ce qui rend la méthode utilisable en grande dimension, là où les méthodes de quadrature s'effondrent : le $1/ n$ ne dépend pas du nombre de variables.

Python : Une intégrale, et la taille d'échantillon qu'elle exige


import numpy as np

rng = np.random.default_rng(11)
ech = np.exp(-rng.random(1000000)**2)
m, s = ech.mean(), ech.std(ddof=1)
print(round(m, 6), round(1.96*s/np.sqrt(len(ech)), 6))
# 0.746516  0.000394        (valeur exacte : 0.746824)

sigma = 0.20098
print(int(np.ceil((1.96*sigma/0.005)**2)))     # 6207

La dernière ligne répond à la question utile : pour que la demi-largeur tombe sous $0, 005$ — deux décimales sûres — il faut environ $6200$ tirages. On dimensionne l'échantillon avant de lancer le calcul, à partir de $σ$ .

Python : Une probabilité qu'on saurait calculer, et une qu'on ne saurait pas


import numpy as np

rng = np.random.default_rng(13)
X = rng.exponential(1, 1000000)
Y = rng.exponential(1, 1000000)
print(round(float(np.mean(X + Y <= 1)), 6))    # 0.264411
print(round(1 - 2/np.e, 6))                 # 0.264241  (calcul exact)

# pas de forme simple a comparer
print(round(float(np.mean(X*Y <= 1)), 6))

La première ligne se vérifie : $X + Y$ suit une loi $γ$ de paramètre $2$ (chapitre 8), dont on sait intégrer la densité. La seconde ne se vérifie pas, et c'est justement pour ces cas-là que la méthode existe. Toujours l'essayer d'abord sur un cas où l'on connaît la réponse : c'est ainsi qu'on teste le code.

13.4.2 Comparer deux estimateurs

Soit $X$ de loi uniforme sur $[0, θ]$ , $θ > 0$ inconnu. Deux estimateurs se présentent : $T_{n} = 2 X_{n}$ , et $T_{n}^{'} = \frac{n + 1}{n} M_{n}$ où $M_{n} = max (X_{1}, \dots, X_{n})$ .

Exemple 13.8Les deux sont sans biais, l'un est bien meilleur

Comme $E (X) = θ /2$ , on a $E (T_{n}) = θ$ et $V (T_{n}) = \frac{4}{n} \cdot \frac{θ ^{2}}{12} = \frac{θ ^{2}}{3 n}$ .

Pour $M_{n}$ , l'indépendance donne $P (M_{n} \leq t) = (t / θ)^{n}$ sur $[0, θ]$ , donc une densité $n t^{n - 1} / θ^{n}$ , puis $E (M_{n}) = \frac{n θ}{n + 1}$ et $E (M_{n}^{2}) = \frac{n θ ^{2}}{n + 2}$ . On en tire $E (T_{n}^{'}) = θ$ et $V (T_{n}^{'}) = \frac{θ ^{2}}{n ( n + 2 )}$ .

Le rapport des variances vaut $\frac{n + 2}{3}$ : pour $n = 10$ , l'estimateur par le maximum est quatre fois moins dispersé.

Python : Comparer par histogrammes


import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(17)
theta, n, N = 3.0, 10, 100000
ech = rng.random((N, n))*theta
T1 = 2*ech.mean(axis=1)
T2 = (n + 1)/n*ech.max(axis=1)

for T, th in ((T1, theta**2/(3*n)), (T2, theta**2/(n*(n + 2)))):
    print(round(T.mean(), 4), round(T.var(), 5), round(th, 5))
# 2.9961  0.29974  0.3      -> T1
# 2.9994  0.07548  0.075    -> T2

plt.hist(T1, bins=60, alpha=0.5)
plt.hist(T2, bins=60, alpha=0.5)
plt.show()

Les deux histogrammes sont centrés au même endroit — les deux estimateurs sont sans biais — mais le second est nettement plus étroit. C'est la variance, et elle seule, qui les départage : le programme exclut explicitement toute comparaison par une quantité qui mélangerait biais et variance.

13.4.3 Comparer des intervalles de confiance

Exemple 13.9L'espérance d'une loi normale d'écart type connu

Soit $(X_{1}, \dots, X_{n})$ un $n$ -échantillon de loi $N (m, σ^{2})$ avec $σ$ connu. La moyenne $X_{n}$ suit $N (m, σ^{2} / n)$ , ce qui donne un intervalle de confiance exact au niveau $95 %$ :

$[X_{n} - \frac{1 , 96 σ}{n}, X_{n} + \frac{1 , 96 σ}{n}] .$

L'inégalité de Bienaymé-Tchebychev en donne un autre, valable sans hypothèse de loi : $[X_{n} \pm \frac{σ}{0 , 05 n}]$ .

Pour $σ = 4$ et $n = 25$ , les demi-largeurs valent $1, 568$ et $3, 578$ .

Python : Demi-largeur moyenne et taux de couverture


import numpy as np

rng = np.random.default_rng(23)
m, sigma, n, N = 20.0, 4.0, 25, 200000
Xb = m + sigma/np.sqrt(n)*rng.standard_normal(N)

d_exact = 1.96*sigma/np.sqrt(n)
d_tcheb = sigma/np.sqrt(0.05*n)
for d in (d_exact, d_tcheb):
    print(round(d, 4), round(float(np.mean(np.abs(Xb - m) <= d)), 4))
# 1.568   0.9495    -> exact : tient sa promesse au plus juste
# 3.5777  1.0000  -> Tchebychev : couvre toujours, 2x plus large

Les deux intervalles sont honnêtes — aucun ne couvre moins de $95 %$ — mais le second est $2, 28$ fois plus large. Tenir sa promesse exactement vaut mieux que la dépasser : un intervalle trop large ne dit plus rien.

ImportantÀ partir de quel rang l'intervalle asymptotique est-il honnête ?

Pour une loi de Bernoulli, l'intervalle $[F_{n} \pm 1, 96 F_{n} (1 - F_{n}) / n]$ n'est justifié qu'à la limite. En simulant, on mesure sa couverture réelle :

	$n = 10$	$n = 30$	$n = 100$	$n = 300$	$n = 1000$
$p = 0, 5$	$0, 89$	$0, 96$	$0, 94$	$0, 94$	$0, 95$
$p = 0, 1$	$0, 65$	$0, 81$	$0, 93$	$0, 94$	$0, 95$
$p = 0, 02$	$0, 18$	$0, 46$	$0, 87$	$0, 93$	$0, 95$

Le rang à partir duquel l'approximation vaut dépend de $p$ : $n = 30$ suffit pour $p = 0, 5$ , il en faut environ $300$ pour $p = 0, 02$ . Un intervalle qui annonce $95 %$ et n'en couvre que $18$ n'est pas approximatif, il est faux — et seule la simulation le révèle.

Continuer sur Adloun : animation, QCM, fiches, exercices