Théorèmes asymptotiques

HAX603X: Modélisation stochastique

Joseph Salmon

joseph.salmon@inria.fr

Inria

Université de Montpellier (Imag)

CNRS

Jean-Baptiste Fermanian

jean-baptiste.fermanian@inria.fr

Université de Montpellier (Imag)

CNRS

Inria

Table des matières

• Loi des grands nombres
• Théorème central limite

Loi des grands nombres

Loi forte des grands nombres

• Résultat fondamental: 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.

Exemple 1 (Cas Bernouilli: pile ou face) La probabilité d’apparition du côté pile (noté \(p\)) peut-être estimée 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 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 =\mathbb{E}(X_1) \]

• Membre de gauche : la fréquence empirique de piles
• Membre de droite : la fréquence théorique de piles

Visualisation de l’exemple du pile ou face

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| echo: false
#| 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.h2("Loi des grands nombres"),
    ui.row(
        ui.column(4,
                ui.input_action_button(
                    "seed",
                    "Nouveau tirage",
                    class_="btn-primary"
                ),
            ),
        ui.column(4,
            ui.input_slider( "p", "Espérance: p", 0.01, 0.99, value=0.35, step=0.01, ticks=False,
                ),
            ),
        ui.column(4,
            ui.input_slider( "n_samples", "Iterations: n", 2, 1000, value=15, step=1, ticks=False)
                ),
            ),
    ui.row(
        ui.card(
            output_widget("my_widget"),
            )
        )
)


def server(input, output, session):
    seed = reactive.Value(41)

    @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=1.18,
            x=0.85,
            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)

Théorème central limite

Au delà de la loi des grands nombres

• 1er ordre d’approximation de la convergence de \(\bar{X}_n\): loi des grands nombres
• 2ème ordre d’approximation: théorème central limite

Enjeu: quantifier les variations de \(\bar X_n - \mu\)

Réponse: théorème central limite (TCL), avec la convergence en loi d’une transformation affine de la moyenne empirique

Théorème central limite

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 = {\rm 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)\).

Interprétation: la moyenne empirique de v.a. i.i.d de variance \(\sigma^2\) se comporte asymptotiquement comme une loi normale \(\mathcal{N}(\mu, \tfrac{\sigma^2}{n})\): \(\quad \bar X_n \approx \mathcal{N}(\mu, \frac{\sigma^2}{n})\).

Note

Hypothèses du théorème plutôt faibles: variance finie uniquement

Formulation de la convergence

Convergence en loi \(\iff\) convergence des fonctions de répartition (aux pts de continuité)

Notations:

• \(\varphi\) : la densité d’une loi normale centrée réduite \(\varphi(x) = \tfrac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi}}\)
• \(\Phi\) : la fonction de répartition d’une loi normale centrée réduite \(\Phi(x) = \displaystyle \int_{-\infty}^{x}\varphi(u) du\)

Ré-écriture du TCL: pour tout \(a < b\) on a alors

\[ \begin{align} \mathbb{P} \left(\bar X_n \in [ \mu + \tfrac{a \sigma}{\sqrt{n}}, \mu + \tfrac{ b \sigma}{\sqrt{n}}] \right) & = \mathbb{P} \left(\tfrac{\bar X_n - \mu}{\sigma} \in [ \tfrac{a}{\sqrt{n}},\tfrac{b}{\sqrt{n}}]\right) \nonumber\\ & \class{fragment}{{} = \mathbb{P} \bigg( a \leq \sqrt n \left(\tfrac{\bar X_n - \mu}{\sigma} \right) \leq b\bigg) \nonumber} \\ & \class{fragment}{{} \underset{n \to \infty}{\longrightarrow} \int_a^b \varphi(x) \, dx = \Phi(b) - \Phi(a) \nonumber}\\ \end{align} \]

Lien intervalle de confiance et TCL

• Question: comment choisir \(a\) et \(b\) pour obtenir un intervalle de confiance à 95% pour \(\mu\)?
• Notation: \(\alpha_n=\mathbb{P} \left(\bar X_n \notin [ \mu + \tfrac{a \sigma}{\sqrt{n}}, \mu + \tfrac{ b \sigma}{\sqrt{n}}] \right)\), on cherche \(a, b\) tels que \(\alpha_n \approx 0.05\).
• Simplification: choix d’un intervalle symétrique autour de \(\mu \implies q=a=-b\), on cherche alors: \(\mathbb{P} \left( |\bar X_n - \mu| > q \tfrac{\sigma}{\sqrt{n}} \right) \approx 0.05\) \[ \begin{align} & 1-\alpha_n \approx \int_{-q}^q \varphi(x) \, dx=\Phi(q)-\Phi(-q)=2 \Phi(q)-1 \\ \implies & \boxed{q\approx\Phi^{-1}(1-\tfrac{\alpha_n}{2})} \end{align} \]
• Interpretation: \(q\) est (approx.) 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

from scipy.stats import norm  # import du module "norm" de scipy.stats
q = norm.ppf((1-0.05/2))      # Quantile 1-0.05/2 (en: Percent point function)
print(f"{q:.2f}")             # Affichage à 2 décimales

1.96

Visualization du TCL

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| 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.h2("TCL"),
    ui.row(
        ui.column(3,
                ui.input_action_button(
                    "seed",
                    "Nouveau tirage",
                    class_="btn-primary"
                ),
            ),
        ui.column(3,
                ui.input_slider(
                    "p",
                    "Espérance: p",
                    min=0.01,
                    max=0.99,
                    value=0.5,
                    step=0.01
                ),
                ),
            ui.column(3,
                ui.input_slider(
                    "n_samples",
                    "Échantillons: n",
                    min=1,
                    max=200,
                    value=100,
                    step=1
                ),
            ),
            ui.column(3,
                ui.input_slider(
                    "n_repetitions",
                    "Répétitions: t",
                    min=1,
                    max=300,
                    value=200,
                    step=1
                ),
            ),
        ),
    ui.row(
        ui.card(
            output_widget("my_widget"),
        )
    )
)


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=1.24,
                xanchor="left",
                x=0.85,
                bgcolor="rgba(0,0,0,0)",
            )
        )
        fig['layout']['xaxis']['title']='Échantillons: n'

        fig.update_layout(autosize=True)

        return fig


app = App(app_ui, server)

Pour aller plus loin: vision des convolutions

Notation: Soient \(f\) et \(g\) définies sur \(\mathbb{R}\) (intégrables au sens de Lebesgue).

Définition 1 (Convolution) La convolution de \(f\) par \(g\) est la fonction \(f*g\) suivante: \[ \begin{align} f*g: \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\).

Somme et convolutions

Théorème 3 (Loi de la somme et convolutions) Soient \(X\) et \(Y\) des v.a. indépendantes de densités respectives \(f\) et \(g\), alors la densité 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)\).

Convolution et TCL

Pour \(X_1, \dots, X_n\), i.i.d., de densité \(f\), on affiche la loi de \(\bar{X}_n\) (à une constante près)

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 460
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_fillable(
    ui.h2("TCL et convolutions"),
    ui.row(
        ui.column(4,
            ui.input_select(
                "loi",
                "Loi sous-jacente",
                {'uniforme': 'Uniforme', 'laplace' : 'Laplace'},
            ),
            ),
        ui.column(4,
            ui.input_slider(
                "n_iter",
                "Échantillons: n",
                1,
                10,
                value=1,
                step=1,
                ticks=False,
                )
            ),
    ),
    ui.row(
        ui.card(
            output_widget("my_widget"),
            )
        )
)

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 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é : moyenne de n variables aléatoires", yanchor="top", y=0.95),
        title_x=0.5,
    )
    fig.update_layout(legend=dict(
        yanchor="top",
        # y=1.2,
        xanchor="left",
        # x=0.1,
        font=dict(size= 18)
    ))
    fig.update_layout(
    height=200,
    margin=dict(l=0, r=0, b=0, t=20),
    )

    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

Pour plus d’info sur les convolutions, voir la vidéo de 3Blue1Brown : Convolutions | Why X+Y in probability is a beautiful mess, 🇬🇧