1. Outils statistiques

1.1. Incertitudes

1.1.1. Généralités

L’incertitude sur une mesure \(x\) est une estimation du « manque de précision » de cette mesure.

• Ce « manque de précision » peuvent être d’origine intrinsèque (p.ex. bruit de Poisson) ou extrinsèque (liés aux performances de l’instrument de mesure) à la quantité mesurée.

• Une mesure sans incertitude ne permet pas d’estimer la précision de la mesure, et est donc sans intérêt quantitatif. Malheureusement, une estimation valable de l’incertitude est toujours plus délicate que la mesure elle-même.

• Si de multiples mesures \(\{x_{i}\}\) sont disponibles, il est possible de les combiner pour construire un estimateur \(\hat{x}\) le plus « précis » possible.

1.1.2. Interprétations fréquentiste et bayésienne

Les mesures et leurs incertitudes doivent permettre d’évaluer la distribution de probabilité de la valeur vraie. Il en existe deux grandes interprétations:

Interprétation fréquentiste:

elle décrit la distribution statistique que suivent/suivraient les mesures si elles sont/étaient répétées un grand nombre de fois. Exemple: le lancer de dé.

Interprétation bayésienne:

elle décrit la plausibilité de la valeur vraie, compte-tenu de toutes les informations et mesures disponibles. Exemple: probabilité d’un accident nucléaire.

Le plus souvent, l’incertitude citée représente la déviation standard (erreur à \(1\sigma\)) des mesures (réelles ou imaginées) autour du meilleur estimateur disponible. Dans l’hypothèse d’une distribution normale des erreurs, \(x \pm \sigma_{x}\) signifie que l’on estime avoir 68% de chance de trouver la « vraie » valeur \(X\) – dont \(x = \hat{X}\) est un estimateur – dans l’intervalle \([x \pm \sigma_{x}]\).

1.1.3. Precision & Accuracy: précision et exactitude

Il faut distinguer deux types de précision:

Precision:

une évaluation de la variabilité statistique de la mesure, i.e. l’erreur statistique;

Accuracy:

une évaluation du biais statistique de la mesure, i.e. l’erreur systématique.

Les deux notions sont indépendantes:

• une mesure fiable est en moyenne juste, mais peut s’écarter ponctuellement de la valeur vraie;

• une mesure précise retourne toujours sensiblement la même valeur, qui peut être significativement différente de la valeur vraie.

Pour améliorer la mesure, il convient donc:

• d’estimer les deux types d’erreur, statistique et systématique;

• de diminuer l’erreur dominante:

  • l’erreur statistique diminue avec le nombre de mesures indépendantes, p.ex. avec la taille de l’échantillon;

  • l’erreur systématique décroît avec une meilleure compréhension du processus de mesure, p.ex. un meilleur modèle ou une correction des biais.

1.2. Statistiques fréquentistes

1.2.1. Propagation des erreurs

1.2.1.1. Introduction

Soit \(f(x, y)\) une quantité dérivée de deux mesures \(\hat{x}\) et \(\hat{y}\) (avec des incertitudes \(\sigma_{x}\) et \(\sigma_{y}\)). Quelle est l’incertitude \(\sigma_{f}\) sur \(\hat{f}\)?

Développons \(f\) à proximité de la valeur \(\hat{f} = f(\hat{x}, \hat{y})\):

\[\begin{split}\begin{aligned} \d f &= \partial_{x}f\,\d x + \partial_{y}f\,\d y \\ (\d f)^{2} &= (\partial_{x}f)^{2}(\d x)^{2} + (\partial_{y}f)^{2}(\d y)^{2} + 2 \partial_{x}f\,\partial_{y}f\,\d x\,\d y \end{aligned}\end{split}\]

\(\d x\) représentant les fluctuations de \(x\) autour de la valeur mesurée \(\hat{x}\), donc \(\moy{\d x} = 0\) (pour une distribution d’erreur symétrique) et \(\moy{(\d x)^{2}} = \sigma_{x}^{2}\), soit:

\[\sigma_{f}^{2} = (\partial_{x}f)^{2}\sigma_{x}^{2} + (\partial_{y}f)^{2}\sigma_{y}^{2} + 2 \partial_{x}f\,\partial_{y}f\,\Cov(x, y)\]

où \(\Cov(x, y) \defeq \moy{(x - \hat{x})(y - \hat{y})}\) désigne la covariance entre \(x\) et \(y\).

• Le développement au 1er ordre de \(\d f\) n’est valable que pour de « petites » erreurs \(\sigma_{x,y}\).

• Même si les distributions initiales en \(x\) et \(y\) sont normales, celle sur \(f\) ne l’est pas nécessairement. Comme les calculs ne sont généralement pas analytiques, il est nécessaire d’estimer le « posterior » de façon numérique (Markov Chain Monte Carlo).

1.2.1.2. Généralisation

Soient

• \(N\) mesures \(\v{y} = (y_1, \ldots y_N)\) avec une matrice \(N\times N\) de covariance \(V_{ij} = \Cov(y_i, y_j)\),

• \(M\) fonctions \(\v{f} = (f_1, \ldots, f_M)\) s’appliquant sur \(\v{y}\).

Au 1er ordre:

\[\Cov(\v{f}(\v{y})) \simeq \m{J}\cdot\m{V}\cdot{\m{J}}^T\]

où la matrice jacobienne \(M\times N\) des dérivées est définie par

\[J_{ij} = \left.\frac{\partial f_i}{\partial y_j}\right|_{\v{y}}\]

Au 1er ordre, cette expression n’est exacte que si \(\v{f}\) est linéaire: \(\v{f}(\v{y}) \equiv \m{J}\cdot\v{y}\). Sinon, il faut s’assurer que les fonctions sont raisonnablement linéaires sur le domaine \(\v{y} \pm \v{\sigma}_{\v{y}}\), ou faire un développement d’ordre plus élevé, ou réaliser la propagation des erreurs numériquement (MCMC).

1.2.1.3. Cas particuliers

Au 1er ordre, sous une forme non-matricielle:

\[\Cov(f_i, f_j) \simeq \sum_{k,l = 1}^{N} \frac{\partial f_i}{\partial y_k}\frac{\partial f_j}{\partial y_l} V_{kl}.\]

Cas décorrélé

Dans le cas particulier où les mesures d’entrée \(\v{y}\) sont totalement décorrélées:

\[\begin{split}\begin{aligned} V_{ij} &= \sigma_i^2 \delta_{ij} \\ \sigma_{f_i}^2 &= \sum_{k=1}^{N} \left(\frac{\partial f_i}{\partial y_k}\right)^2 \sigma_k^2. \end{aligned}\end{split}\]

1.2.2. Maximum de vraisemblance

1.2.2.1. Définitions

Soient \(\v{X}\) des variables aléatoires suivant une densité de probabilité \(P_{\v{\theta}}(\v{x}) = P(\v{x}|\v{\theta})\) décrite par les paramètres \(\v{\theta}\).

Fonction de vraisemblance (likelihood)

\[\mathcal{L}(\v{\theta} | \v{x}) \defeq P(\v{x}|\v{\theta}).\]

L’estimateur de maximum de vraisemblance \(\hat{\v{\theta}}\) des paramètres \(\v{\theta}\) est celui maximisant la fonction de vraisemblance, compte-tenu des valeurs observées \(\v{x}\) (et éventuellement de leurs erreurs \(\v{\sigma}\)):

\[\hat{\v{\theta}} = \argmax \mathcal{L}(\v{\theta} | \v{x}, \v{\sigma}) = \argmin \tilde{\mathcal{L}}(\v{\theta} | \v{x}, \v{\sigma}) \quad\text{avec}\quad \tilde{\mathcal{L}} \defeq -2\ln\mathcal{L}\]

1.2.2.2. Moyenne « optimale »

Soient \(N\) observations indépendantes \(\{x_{i} \pm \sigma_{i}\}\), chacune supposée suivre une loi normale \(\mathcal{N}(\mu, \sigma_{i}^{2})\) de moyenne constante \(\mu\).

Quel est le Maximum Likelihood Estimate \(\hat{\mu}\) de la moyenne \(\mu\), et son erreur \(\sigma_{\hat{\mu}}\)?

\[\begin{split}\begin{aligned} \text{Vraisemblance} && \mathcal{L}(\mu) &\defeq P(\v{x}|\mu, \v{\sigma}) = \prod_{i=1}^{N} \frac{1}{\sigma_{i}\sqrt{2\pi}} \exp\left(-\frac{(x_{i} - \mu)^{2}}{2\sigma_{i}^{2}}\right) \\ \text{Log-vraisemblance} && \tilde{\mathcal{L}}(\mu) &\defeq -2\ln\mathcal{L}(\mu) = \sum_{i=1}^{N} \frac{(x_{i} - \mu)^{2}}{\sigma_{i}^{2}} + \sum_{i} \ln\sigma_{i}^{2} + N\ln 2\pi \\ \text{Extremum} && \left.\dpart{\tilde{\mathcal{L}}}{\mu}\right|_{\hat{\mu}} &= 0 = -2 \sum_{i} \frac{x_{i} - \hat{\mu}}{\sigma_{i}^{2}} \\ \text{Estimateur moyenne} && \hat{\mu} &= \sum_{i} \tilde{w}_{i} x_{i} \quad\text{avec}\quad w_{i} = 1/\sigma_{i}^{2}, \quad\tilde{w}_{i} = w_{i} / \sum_{j} w_{j}\\ \text{Erreur sur la moyenne} && \sigma_{\hat{\mu}}^{2} &= \sum_{i} \tilde{w}_{i}^{2} \sigma_{i}^{2} = 1 \Big/ \sum_{i} 1/\sigma_{i}^{2} \\ \text{Cas particulier} \quad \sigma_{i} = \sigma = \text{cte} && \hat{\mu} &= \moy{x}, \quad \sigma_{\hat{\mu}} = \sigma/\sqrt{N} \end{aligned}\end{split}\]

1.2.3. Minimisation du \(\chi^{2}\)

Soient \(N\) observations indépendantes \(\{x_{i} \pm \sigma_{i}\}\), chacune supposée suivre une loi normale \(\mathcal{N}(\mu_{i}(\v{\theta}), \sigma_{i}^{2})\) dont la moyenne est modélisée par les fonctions \(\mu_{i}(\v{\theta})\) des paramètres \(\v{\theta}\).

\[\begin{split}\begin{aligned} \text{Vraisemblance} && \mathcal{L}(\v{\theta}) &\defeq P(\v{x}|\v{\theta}, \v{\sigma}) = \prod_{i=1}^{N} \frac{1}{\sigma_{i}\sqrt{2\pi}} \exp\left(-\frac{(x_{i} - \mu_{i}(\v{\theta}))^{2}}{2\sigma_{i}^{2}}\right) \\ \text{Log-vraisemblance} && -2\ln\mathcal{L}(\v{\theta}) &\equiv {\color{vert} \sum_{i=1}^{N} \left(\frac{x_{i} - \mu_{i}(\v{\theta})}{\sigma_{i}}\right)^{2} \defeq \chi^{2}(\v{\theta})} \\ \text{Estimateur MLE} && \hat{\v{\theta}} &= \argmin \chi^{2}(\v{\theta}) \end{aligned}\end{split}\]

1.2.3.1. Limites

Attention

L’estimateur de maximum de vraisemblance minimise le \(\chi^{2}\), dans l’hypothèse où les erreurs \(\v{\sigma}\) sont normales et connues a priori et le modèle \(\v{\mu}(\v{\theta})\) est correct.

Si les erreurs ne sont pas gaussiennes (p.ex. bruit de photons), ne sont pas correctement estimées (p.ex. erreurs systématiques ou présence de points aberrants), dépendent des paramètres \(\v{\theta}\), ou bien que le modèle ne reflète pas la réalité (paramètres supplémentaires ou forme fonctionnelle inadéquate), \(\argmin \chi^{2}\) n’est pas un MLE!

1.2.3.2. Estimation des erreurs sur les paramètres

Développons le \(\chi^{2}\) autour de sa valeur minimale \(\chi^{2}_{\min} = \chi^{2}(\hat{\v{\theta}})\):

\[\begin{split}\begin{aligned} \chi^{2} - \chi^{2}_{\min} &= \Delta\chi^{2} \\ &\simeq \sum_{i} \cancel{\left.\frac{\partial\chi^{2}}{\partial\theta_{i}}\right|_{\min}} \Delta\theta_{i} + \sum_{i,j} \frac{1}{2}{\color{bleu} \left. \frac{\partial^{2}\chi^{2}}{\partial\theta_{i}\partial\theta_{j}} \right|_{\min}} \Delta\theta_{i}\Delta\theta_{j} \end{aligned}\end{split}\]

soit, dans l'approximation quadratique

\[\Cov(\v{\theta}) = (2\Delta\chi^{2}\,{\color{bleu}H}^{-1})\]

où la matrice hessienne \(H\) du \(\chi^{2}\) est évaluée en \(\hat{\v{\theta}}\):

\[\begin{split}\begin{aligned} H_{ij} &\defeq \left.\frac{\partial^2\chi^2}{% \partial\theta_i\partial\theta_j}\right|_{\v{\theta}=\hat{\v{\theta}}} \\ &= 2 \left.\frac{\partial\v{\mu}}{\partial\theta_i}\right|_{\hat{\v{\theta}}}^T \cdot \m{V}^{-1} \cdot \left.\frac{\partial\v{\mu}}{\partial\theta_j}\right|_{\hat{\v{\theta}}}. \end{aligned}\end{split}\]

\(\Delta\chi^{2}\) dépend du niveau de confiance \(\alpha\) et du nombre \(M' \leq M\) de paramètres pour lesquels l’erreur doit être calculée en marginalisant sur les \(M - M'\) autres paramètres de nuisance:

• pour l’erreur à \(n\,\sigma\) associée à un seul paramètre (\(M' = 1\)): \(\Delta\chi^{2} = n^{2}\)

• pour la matrice de covariance \(M'\times M'\) à \(1\sigma\) (soit \(\alpha = \int_{-1}^{+1}\mathcal{N}(x)\d x = 0.683\)): \(\Delta\chi^{2} = X_{M'}^{-1}(\alpha)\) où \(X_{M'}^{-1}\) est la fonction quantile de la distribution du \(\chi^{2}\) à \(M'\) degrés de liberté.

Cas décorrélé

Dans le cas particulier où la matrice de covariance \(\m{V}\) des données est purement diagonale:

\[\begin{split}\begin{aligned} V_{ij} &= \sigma_i^2 \delta_{ij} \\ \chi^2(\v{\theta}) &= \sum_{i=1}^N \frac{(y_i - \mu_i(\v{\theta}))^2}{\sigma_i^2}, \\ H_{ij} &= 2 \sum_{k=1}^N \frac{1}{\sigma_k^2} \frac{\partial \mu_k}{\partial\theta_i} \frac{\partial \mu_k}{\partial\theta_j}. \end{aligned}\end{split}\]

1.2.3.3. Qualité de l’ajustement

\(\displaystyle \chi^{2}(\hat{\v{\theta}}) = \min_{\v{\theta}}\sum_{i} \left(\frac{x_{i} - \mu_{i}(\v{\theta})}{\sigma_{i}}\right)^{2}\) est par définition une somme de \(N\) variables aléatoires normales \(\mathcal{N}(0,1)\) soumises à \(M\) contraintes \(\partial\chi^{2}/\partial \theta_{j = 1\ldots M} = 0\): c’est donc une variable aléatoire suivant une loi de distribution du \(\chi^{2}\) à \(k = N - M\) degrés de liberté:

\[\chi^{2}_{k}(x) = \frac{1}{2^{k/2}\Gamma(k/2)} x^{k/2-1}e^{-x/2}\]

Qualité de l’ajustement (goodness of fit) Puisque \(E(\chi^{2}_{k}) = k\), on devrait avoir \(\chi^{2}_{\min} \approx k = N - M\). Dans le cas contraire (\(\chi^{2}_{\min} \gg k\)), le résultat de l’ajustement n’est pas fiable: les erreurs de tout ou partie des observations ne sont pas crédibles, ou le modèle est simplement faux.

Si un \(\chi^{2}\) aussi élevé est trop improbable (valeur \(p\) trop faible), le modèle doit être rejeté (dans l’hypothèse où les erreurs sont correctement estimées et normalement distribuées).

1.2.3.4. Cas particulier: modèle linéaire

On suppose que les modèles \(\mu_{i}\) sont linéaires en \(\v{\theta}\) via la matrice jacobienne \(\m{J}\):

\[\v{\mu} = \m{J} \v{\theta}\]

On a alors

\[\chi^{2}(\v{\theta}) = (\v{x} - \m{J}\v{\theta})^{T} \m{W} (\v{x} - \m{J}\v{\theta})\]

où \(\m{W} = [1/\sigma_{i}^{2}] = \m{C}^{-1}\) est la matrice des poids, et

\[\hat{\v{\theta}} = (\m{J}^{T} \m{W} \m{J})^{-1} \m{J}^{T} \m{W} \v{x}.\]

1.2.3.5. Cas général non-linéaire

Dans le cas général, il faut déterminer \(\hat{\v{\theta}}\) par minimisation numérique de \(\chi^{2}(\v{\theta})\), en s’assurant que:

• la méthode de minimisation a convergé vers un minimum local,

• la méthode a bien convergé vers le minimum global,

• le résultat de l’ajustement est crédible.

1.3. Statistiques bayésiennes

Deux paradigmes statistiques:

fréquentiste:

les probabilités sont des fréquences d’apparition;

bayésien:

les probabilités sont des « degrés de plausibilité ».

• Dans la limite d’un grand nombre d’observations, les deux approches sont équivalentes.

• Les probabilités fréquentistes ne s’appliquent pas dans le cas d’un petit nombre d’observations (voire une seule réalisation): seules les probabilités bayésiennes sont alors correctement définies.

• L’inférence bayésienne permet d’estimer la probabilité d’un évènement à partir de celles d’autres évènements déjà évalués et de connaissances a priori.

1.3.1. Estimateurs bayésiens

1.3.1.1. Théorème de Bayes

Soient deux évenements \(A\) (« proposition ») et \(B\) (« évidence »):

\[\begin{split}\begin{aligned} P(A, B) &= P(A | B) P(B) \\ = P(B, A) &= P(B | A) P(A) \\ \text{donc}\qquad P(A | B) &= \frac{P(B | A) P(A)}{P(B)} \quad\text{Théorème de Bayes} \end{aligned}\end{split}\]

où

• \(P(A | B)\) est la probabilité a posteriori (postérieure ou conditionnelle) de \(A\) sachant \(B\), i.e. le degré de plausibilité de \(A\) en tenant compte de \(B\).

• \(P(B | A)\) est la fonction de vraisemblance de \(B\) pour un \(A\) connu;

• \(P(A)\) est la probabilité a priori (antérieure ou marginale) de \(A\) indépendamment de \(B\): elle représente le degré de plausibilité initiale sur \(A\).

• \(P(B)\) est la probabilité inconditionnelle de \(B\), que l’on considère souvent comme une simple constante de normalisation, mais qui peut avoir une interprétation comme évidence de \(B\).

1.3.1.2. Estimation de paramètres

Dans le cadre d’une comparaison entre les données \(d\) et le modèle \(m\) dépendant des paramètres \(\theta\) que l’on cherche à estimer:

\[P(\theta | d, m) = \frac{P(d | m, \theta) P(\theta | m)}{P(d | m)}\]

où

• \(P(\theta) \equiv P(\theta | d, m)\) est la distribution a posteriori des paramètres \(\theta\) une fois acquis les observations \(d\) sous l’hypothèse du modèle \(m\);

• \(P(d | m, \theta) \defeq \mathcal{L}(\theta | d, m)\) est la fonction de vraisemblance, i.e. la probabilité d’observer les données \(d\) compte-tenu du modèle \(m(\theta)\);

• \(P(\theta | m)\) est le prior sur les paramètres \(\theta\) caractérisant une « plausibilité » a priori (p.ex. \(0 \leq \theta \leq 1\));

• \(P(d | m) = \int_{\theta} P(d | m, \theta) P(\theta | m)\d\theta\) est l’évidence bayésienne, ou vraisemblance marginalisée. Si l’on ne considère qu’un modèle \(m\), c’est une constante de normalisation.

1.3.1.3. Facteur de Bayes

Deux modèles \(m_{1}\) et \(m_{2}\) peuvent être comparés et éventuellement discriminés à l’aide du facteur de Bayes, i.e. le rapport des évidences bayésiennes de chacun des modèles compte-tenu des données \(d\):

\[K \defeq \frac{P(d | m_{1})}{P(d | m_{2})}\]

Si \(K > 1\), les données favorisent le modèle \(m_{1}\), mais pas forcément de façon significative:

\(K\)	Évidence
\(< 1\)	négative (\(m_{2}\) est favorisée)
\(1 - 3\)	négligeable
\(3 - 20\)	positive
\(20 - 150\)	forte
\(> 150\)	très forte