from scipy.stats import norm
q = norm.ppf(1-0.05/2)
print(f"Gaussienne centrée réduite,\nQuantile de niveau (1-α/2):\nq = {q:.2f}")Gaussienne centrée réduite,
Quantile de niveau (1-α/2):
q = 1.96Théorèmes asymptotiques
{{
}}
Loi des grands nombres
Le premier résultat fondamental en probabilités concerne le comportement asymptotique de la moyenne empirique: \bar X_n = \dfrac{X_1 + \cdots + X_n}{n} \enspace. quand on observe n variables aléatoires i.i.d X_1,\dots,X_n, ayant une espérance finie.
Théorème 1 (Loi forte des grands nombres) Soit (X_n)_{n \geq 1} une suite de variables aléatoires indépendantes et identiquement distribuées (i.i.d.) dans L^1(\Omega, \mathcal{F}, \mathbb{P}). Notons \mu = \mathbb{E}[X_1]. Alors \bar X_n converge vers \mu presque sûrement : \mathbb{P}\bigg( \dfrac{X_1 + \cdots + X_n}{n} \underset{n \to \infty}{\longrightarrow} \mu \bigg) = 1\,.
Interprétation: Intuitivement, la probabilité d’un événement A correspond à la fréquence d’apparition de A quand on répète une expérience qui fait intervenir cet événement. Par exemple, si on dispose une pièce truquée, on estimera la probabilité d’apparition du côté pile en lançant la pièce un grand nombre de fois et en comptant le nombre de pile obtenu. La loi des grands nombres justifie a posteriori cette intuition : si X_1, \ldots, X_n sont i.i.d. de loi de Bernoulli de paramètre p, alors \dfrac{X_1 + \cdots + X_n}{n} \xrightarrow[n \to \infty]{p.s.} p \enspace. Le membre de gauche correspond au nombre empirique de pile obtenu, celui de droite à la valeur théorique.
Remarque: Bien qu’assez intuitif, ce théorème est difficile à démontrer, cf.(Ouvrard 2008; Barbe et Ledoux 2006) ou encore (Williams 1991) pour une version de preuve avec des martingales.
Williams, D. 1991. Probability with martingales. Cambridge Mathematical Textbooks. Cambridge: Cambridge University Press.
#| standalone: true
#| viewerHeight: 550
import numpy as np
from shiny import ui, render, App, reactive
from shinywidgets import output_widget, register_widget
from plotly.subplots import make_subplots
import plotly.graph_objs as go
n_init = 42
app_ui = ui.page_fluid(
ui.tags.head(
ui.tags.style("""
.bslib-sidebar-layout > .sidebar > .sidebar-content {
display: flex;
flex-direction: column;
padding: 0.5rem;
}
.irs--shiny .irs-line {
top: 27px;
height: 4px;
}
.irs.irs--shiny .irs-bar {
top: 27px;
height: 4px;
}
.irs--shiny .irs-handle {
top: 23px;
width: 22px;
height: 10px;
""")
),
ui.panel_title("Loi des grands nombres"),
ui.layout_sidebar(
ui.panel_sidebar(
ui.input_action_button("seed", "Nouveau tirage", class_="btn-primary"),
ui.input_slider( "p", "Espérance: p", 0.01, 0.99, value=0.35, step=0.01, ticks=False,
),
ui.input_slider( "n_samples", "Échantillons: n", 2, 1000, value=15, step=1, ticks=False),
width=3),
ui.panel_main(output_widget("my_widget"), width = 9)
)
)
def server(input, output, session):
seed = reactive.Value(42)
@reactive.Effect
@reactive.event(input.seed)
def _():
seed.set(np.random.randint(0, 1000))
subplots = make_subplots(
rows=2,
cols=1,
vertical_spacing=0.45,
horizontal_spacing=0.04,
row_heights=[5, 1],
subplot_titles=(
f"Moyenne empirique: loi de Bernoulli",
"Tirages aléatoires <span style='color:rgb(66, 139, 202)'>bleu: 0</span>, <span style='color:rgb(255, 0, 0)'>rouge: 1</span> (seed="
+ f"{n_init:03}"
+ ")",
),
)
fig = go.FigureWidget(subplots)
fig.add_trace(
go.Scatter(
mode="lines",
line=dict(color="black", width=3),
x=[],
y=[],
name=r"Moyenne <br> empirique",
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
mode="lines",
line=dict(dash="dash", color="black", width=1),
marker={},
x=[],
y=[],
name=r"p",
),
row=1,
col=1,
)
fig.add_trace(
go.Heatmap(
x=[],
z=[],
colorscale=[[0, "rgb(66, 139, 202)"], [1, "rgb(255,0,0)"]],
showscale=False,
),
row=2,
col=1,
)
fig.update_yaxes(range=[0, 1.1], row=1, col=1)
fig.update_xaxes(matches="x1", row=2, col=1)
fig.update_yaxes(visible=False, row=2, col=1)
fig.update_xaxes(visible=False, row=2, col=1)
fig.update_layout(
template="simple_white",
showlegend=True,
xaxis_title="Échantillons: n",
)
fig.update_layout(autosize=True)
fig.update_layout(
legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.65,
bgcolor="rgba(0,0,0,0)",
)
)
fig.update_layout(
margin=dict(l=0, r=0, b=10, t=70),
)
register_widget("my_widget", fig)
@reactive.Effect
def _():
p = input.p()
n_samples = input.n_samples()
rng = np.random.default_rng(seed())
iterations = np.arange(1, n_samples + 1)
samples = rng.binomial(1, p, size=n_samples)
means_samples = np.cumsum(samples) / np.arange(1, n_samples + 1)
# Update data in fig:
fig.data[0].x = iterations
fig.data[0].y = means_samples
fig.data[1].x = iterations
fig.data[1].y = np.full((n_samples), p)
fig.data[2].x = iterations
fig.data[2].z = [samples]
fig.update_xaxes(range=[1, n_samples + 1])
# Update the subplot titles:
fig.layout.annotations[1].update(
text=f"Tirages aléatoires (seed="
+ f"{seed():03}"
+ ") <br> <span style='color:rgb(66, 139, 202)'>bleu: 0</span>, <span style='color:rgb(255, 0, 0)'>rouge: 1</span> "
)
app = App(app_ui, server)
Pour aller plus loin:
Quant p varie, à n fixé…les signaux générés sont très très proches, ce qui ne devrait pas être le cas sans structuration particulière de la génération. L’aléa est imparfait!
Théorème central limite (TCL)
Une fois la loi des grands nombres établie, on peut se demander quel est l’ordre suivant dans le développement asymptotique de \bar X_n - \mu, ou de manière équivalente de S_n - n \mu, où S_n = X_1 + \cdots + X_n. Le théorème suivant répond à cette question, en donnant une convergence en loi d’une transformation affine de la moyenne empirique:
Théorème 2 (Théorème central limite) Soit X_1, \ldots, X_n une suite de variables aléatoires i.i.d de variance \sigma^2 = {\mathrm var}(X_1) \in ]0, \infty[. On note \mu = \mathbb{E}[X_1] leur espérance. Alors \sqrt n \left(\tfrac{\bar X_n - \mu}{\sigma} \right) \xrightarrow[n \to +\infty]{\mathcal{L}} N\enspace, où N suit une loi normale centrée réduite : N \sim\mathcal{N}(0,1).
Preuve: cf.(Ouvrard 2008; Barbe et Ledoux 2006).
Ouvrard, J.-Y. 2008. Probabilités : Tome 1, Licence - CAPES. 2ᵉ éd. Enseignement des mathématiques. Cassini.
Barbe, Philippe, et Michel Ledoux. 2006. Probabilités.
On peut interpréter ce théorème grossièrement de la façon suivante: la moyenne empirique de variables aléatoires i.i.d de variance \sigma^2 se comporte asymptotiquement comme une loi normale \mathcal{N}(\mu, \tfrac{\sigma^2}{n}), ce que l’on écrit avec un abus de notation:
\bar X_n \approx \mathcal{N}(\mu, \frac{\sigma^2}{n}) \enspace.
En termes de somme cumulée empirique, la convergence se réécrit
\tfrac{S_n - n \mu}{\sqrt n \sigma} \xrightarrow[n \to +\infty]{\mathcal{L}} N \enspace.
Les hypothèses de ce théorème sont plutôt faibles (il suffit de supposer une variance finie). Pourtant, le résultat est universel : la loi de départ peut être aussi farfelue que l’on veut, elle se rapprochera toujours asymptotiquement d’une loi normale.
On rappelle que la convergence en loi est équivalente à la convergence des fonctions de répartition en tout point de continuité de la limite. Ainsi, le théorème central limite se réécrit de la manière suivante : pour tout a < b, notons \alpha_n=\mathbb{P} \left(\bar X_n \notin [ \mu + \tfrac{a \sigma}{\sqrt{n}}, \mu + \tfrac{ b \sigma}{\sqrt{n}}] \right). Ainsi, \begin{align} 1-\alpha_n& = \mathbb{P} \left(\bar X_n \in [ \mu + \tfrac{a \sigma}{\sqrt{n}}, \mu + \tfrac{ b \sigma}{\sqrt{n}}] \right)\nonumber\\ & = \mathbb{P} \left(\tfrac{\bar X_n - \mu}{\sigma} \in [ \tfrac{a}{\sqrt{n}},\tfrac{b}{\sqrt{n}}]\right) \nonumber\\ & = \mathbb{P} \bigg( a \leq \sqrt n \left(\tfrac{\bar X_n - \mu}{\sigma} \right) \leq b\bigg) \nonumber\\ & \underset{n \to \infty}{\longrightarrow} \int_a^b \varphi(x) \, dx\,. \nonumber\\ \end{align} où l’on note \varphi (resp. \Phi) la densité (resp. la fonction de répartition) d’une loi normale centrée réduite, définie pour tout x\in\mathbb{R} par \varphi(x)=\tfrac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi}} (resp. \Phi(x)= \int_{-\infty}^{x}\varphi(u) du).
Dans le cas classique d’un intervalle de confiance à 95%, c’est-à-dire quand \alpha_n=0.05, et en prenant un intervalle de confiance symétrique (alors a=-t et b=q) on obtient 1-\alpha_n= \int_{-q}^q \varphi(x) \, dx=\Phi(q)-\Phi(-q)=2 \Phi(q)-1 \implies \boxed{q=\Phi^{-1}(1-\tfrac{\alpha_n}{2})} et q est donc le quantile de niveau 1-\tfrac{\alpha_n}{2} de la loi normale centrée réduite. Numériquement on peut facilement évaluer q et vérifier que q\approx 1.96 avec scipy:
Exemple 1 (Loi de Bernoulli) On considère des variables aléatoires X_1, \ldots, X_n i.i.d. suivant une loi de Bernoulli de paramètre p \in ]0,1[, dont l’espérance et la variance sont respectivemenbt p et p(1-p). Le théorème central limite donne alors \sqrt n \left(\frac{\bar X_n - p}{p (1-p)} \right) \xrightarrow[n \to +\infty]{\mathcal{L}} N\,, avec N \sim \mathcal{N}(0,1). Cette convergence est illustrée dans le widget ci-dessous. Le contexte est le suivant. On répète t fois le processus, qui consiste à afficher (\bar{X}_k)_{k \in [n]}, où les n variables aléatoires sont i.i.d. et suivent une loi de Bernoulli de paramètre p.
#| standalone: true
#| viewerHeight: 600
import numpy as np
from scipy.stats import norm
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from shiny import ui, render, App, reactive
from shinywidgets import output_widget, render_widget
app_ui = ui.page_fluid(
ui.tags.head(
ui.tags.style("""
.bslib-sidebar-layout > .sidebar > .sidebar-content {
display: flex;
flex-direction: column;
padding: 0.1rem;
}
.irs--shiny .irs-line {
top: 27px;
height: 4px;
}
.irs.irs.irs--shiny .irs-bar {
top: 27px;
height: 4px;
}
.irs--shiny .irs-handle {
top: 23px;
width: 22px;
height: 10px;
""")
),
ui.panel_title("TCL"),
ui.layout_sidebar(
ui.panel_sidebar(
ui.input_action_button("seed", "Nouveau tirage",class_="btn-primary"),
ui.input_slider("p", "Espérance: p", 0.01, 0.99, value=0.5, step=0.01, ticks=False),
ui.input_slider("n_samples", "Échantillons: n", 1, 200, value=100, step=1, ticks=False),
ui.input_slider("n_repetitions", "Répétitions: t", 1, 300, value=200, step=1, ticks=False),
width=3
),
ui.panel_main(output_widget("my_widget"), width = 9)
)
)
def server(input, output, session):
seed = reactive.Value(42)
@reactive.Effect
@reactive.event(input.seed)
def _():
seed.set(np.random.randint(0, 1000))
@output
@render_widget
def my_widget():
rng = np.random.default_rng(seed())
p = input.p()
n_repetitions = input.n_repetitions()
n_samples = input.n_samples()
iterations = np.arange(1, n_samples + 1)
samples = rng.binomial(1, p, size=(n_repetitions, n_samples))
means_samples = np.cumsum(samples, axis=1) / np.arange(1, n_samples + 1)
x_hist = np.linspace(0, 1, num=300)
# Create figure
fig = make_subplots(
rows=1,
cols=3,
# vertical_spacing=0.5,
horizontal_spacing=0.02,
column_widths=[20, 2, 3],
subplot_titles=("t = " + str(n_repetitions) + " répétitions","","")
)
for i in range(n_repetitions):
fig.add_trace(
go.Scatter(
mode='lines',
line=dict(color="rgba(0,0,0,0.05)", width=1),
x=iterations,
y=means_samples[i,:],
),
row=1, col=1,
)
fig.add_trace(
go.Scatter(
mode='lines',
line=dict(dash="dash", color="blue", width=1),
marker={},
x=iterations,
y=np.full((n_samples), p),
name='p'),
row=1, col=1,
)
fig.update_layout(
template="simple_white",
showlegend=False,
)
fig.add_trace(
go.Scatter(
mode='markers',
marker=dict(color="rgba(0,0,0,0.05)", size=4),
x=np.zeros(n_samples),
y=means_samples[:,-1],
),
row=1, col=2,
)
y_hist, bins = np.histogram(means_samples[:,-1], bins=int(np.sqrt(n_repetitions)), density=True)
fig.add_trace(
go.Bar(x=y_hist, y=bins[:-1] + np.diff(bins)/2,
opacity=0.75,
marker_color = 'black',
orientation='h',
width=np.diff(bins),
name="Tirages de moyennes empiriques",
),
row=1, col=3,
)
fig.add_trace(
go.Scatter(x=norm.pdf(x_hist, p, np.sqrt(p*(1-p) / n_samples)),
y=x_hist,
mode='lines',
line=dict(color="red"),
legendgroup='1',
name="TCL"
),
row=1, col=3,
)
fig.update_xaxes(range=[1, n_samples + 1])
fig.update_yaxes(range=[-.05, 1.05], row=1, col=1)
fig.update_yaxes(matches="y1",visible = False, row=1, col=2)
fig.update_xaxes(range=[-0.2, 0.2], visible = False, row=1, col=2)
fig.update_yaxes(matches="y1", row=1, col=3, visible=False)
fig.update_xaxes(range=[0, 1.1 / np.sqrt(2*np.pi* p*(1-p) / n_samples)], row=1, col=3)
fig.update_xaxes(visible=False, row=1, col=3)
for trace in fig['data']:
print(trace)
if (trace['name'] != 'TCL') and (trace['name'] != 'p'):
trace['showlegend'] = False
fig.update_layout(
margin=dict(l=0, r=0, b=10, t=100),
)
fig.update_layout(
title=dict(text="Distribution de la moyenne empirique<br> (cas loi de Bernoulli)", yanchor="top", y=0.95),
title_x=0.5,
showlegend=True,
)
fig.update_layout(
legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.89,
bgcolor="rgba(0,0,0,0)",
)
)
fig['layout']['xaxis']['title']='Échantillons: n'
fig.update_layout(autosize=True)
return fig
app = App(app_ui, server)
Une autre illustration possible de la convergence donnée par le TCL est celle qui correspond au point de vue donnée par l’analyse. Pour cela supposons que l’on ait une suite de variables aléatoires réelles X_1, \dots, X_n, i.i.d. dont la fonction de densité commune est notée par f.
On rappelle quelques éléments de probabilités concernant les densités. Pour cela on rappelle la définition de la convolution deux fonctions. Pour cela prenons deux fonctions f et g définies sur \mathbb{R} et qui sont intégrables au sens de Lebesgue. La convolution de f par g est alors la fonction f*g suivante: \begin{align*} \mathbb{R} &\mapsto \mathbb{R} \nonumber\\ x &\to \int_{-\infty}^{+\infty} f(x-y)g(y) dy \enspace.\nonumber \end{align*}
Note
On peut aussi obtenir f*g(x) en calculant \int_{\mathbb{R}^2} f(u)g(v) {1\hspace{-3.8pt} 1}_{u+v=x} du dv.
Théorème 3 (Loi de la somme et convolutions) Soient X et Y des v.a. indépendantes de densités f et g respectivement, la loi de X+Y est donnée par la convolution f*g.
Rappel: pour un scalaire \alpha\neq 0, la densité de \alpha X est donnée par la fonction x \mapsto \frac{1}{|\alpha|} \cdot f(\frac{x}{\alpha}).
Corollaire 1 (Loi de la moyenne) Soient X_1,\dots,X_n des v.a. i.i.d. de densité f, la densité de \bar{X}_n est donnée par la fonction x \mapsto n \cdot [f*\dots*f](n \cdot x).
Dessous, pour X_1, \dots, X_n, i.i.d., de densité f, on affiche la densité de la loi de \bar{X}_n.
#| standalone: true
#| echo: false
#| viewerHeight: 550
import numpy as np
from shiny import ui, render, App, reactive
from shinywidgets import output_widget, register_widget
from plotly.subplots import make_subplots
import plotly.graph_objs as go
from scipy import signal
app_ui = ui.page_fluid(
ui.tags.head(
ui.tags.style("""
.bslib-sidebar-layout > .sidebar > .sidebar-content {
display: flex;
flex-direction: column;
padding: 0.5rem;
}
.irs--shiny .irs-line {
top: 27px;
height: 2px;
}
.irs.irs--shiny .irs-bar {
top: 27px;
height: 2px;
}
.irs--shiny .irs-handle {
top: 23px;
width: 22px;
height: 10px;
""")
),
ui.panel_title("TCL et convolutions"),
ui.layout_sidebar(
ui.panel_sidebar(
ui.input_select(
"loi",
"Loi sous-jacente",
{'uniforme': 'Uniforme', 'laplace' : 'Laplace'},
),
ui.input_slider(
"n_iter",
"Échantillons: n",
1,
10,
value=1,
step=1,
ticks=False,
), width = 3
),
ui.panel_main(
output_widget("my_widget"), width = 9
)
)
)
nnzeros = 10001
x_min = -20
x_max = 20
x = np.linspace(x_min, x_max, nnzeros)
y = np.zeros(nnzeros)
mask = np.where(np.abs(x) <= 0.5, 1, 0)
y[mask == 1] = 1
delta = (x_max - x_min) / nnzeros
var = np.sum(y * x**2 *(delta)) - (np.sum(y * x * delta))**2
def convolve(signal_in, n_convolutions, delta):
output = np.zeros(len(signal_in))
if n_convolutions == 0:
return output
elif n_convolutions == 1:
return signal_in
else:
output = signal_in.copy()
for i in range(n_convolutions - 1):
output = signal.fftconvolve(
output * delta, signal_in, mode="same"
)
return output
def server(input, output, session):
fig = go.FigureWidget()
fig.add_trace(
go.Scatter(
mode="lines",
line=dict(color="black", width=3),
x=[],
y=[],
name="loi de de la moyenne empirique<br>(variance adéquate)",
)
)
fig.add_trace(
go.Scatter(
x=x,
y=np.exp(-(x**2) / (2 * var)) / np.sqrt(2 * var * np.pi),
mode="lines",
line=dict(dash="dash", color="red", width=2),
name=f"Loi normale<br>(variance adéquate)",
)
)
fig.update_xaxes(range=[-3, 3], position=0.)
fig.update_yaxes(range=[0, 2 * np.max(y)], position=0.5, showticklabels=False)
fig.update_layout(
template="simple_white",
showlegend=True,
autosize=True,
title=dict(text="Densité : <br> moyenne de n variables aléatoires", yanchor="top", y=0.95),
title_x=0.5,
)
fig.update_layout(legend=dict(
yanchor="top",
y=0.95,
xanchor="left",
x=0.8,
font=dict(size= 18)
))
fig.update_layout(
height=250,
margin=dict(l=0, r=0, b=10, t=100),
)
register_widget("my_widget", fig)
@reactive.Effect
def _():
if str(input.loi()) == 'uniforme':
y = np.zeros(nnzeros)
mask = np.where(np.abs(x) <= 0.5, 1, 0)
y[mask == 1] = 1
else:
y=np.exp(-np.abs(x)) / 2
y = y / (np.sum(y) * delta)
var = np.sum(y * x**2 * delta) - (np.sum(y * x * delta))**2
y_display = convolve(y, input.n_iter(), delta)
# Update data in fig:
fig.data[0].x = x / np.sqrt(input.n_iter())
fig.data[0].y = y_display * np.sqrt(input.n_iter())
fig.data[1].y = np.exp(-(x**2) / (2 * var)) / np.sqrt(2 * var * np.pi)
fig.update_yaxes(range=[0, 2 * np.max(y)], position=0.5, showticklabels=False)
app = App(app_ui, server)Pour aller plus loin sur les convolutions, voir la vidéo de 3Blue1Brown à ce sujet: Convolutions | Why X+Y in probability is a beautiful mess