\documentclass{pbsheet}

\TITRE{Feuille de TP n°7 \\ Vecteurs aléatoires et modèle linéaire gaussiens}
\FORMATION{Épreuve de modélisation - Agrégation Externe de Mathématiques}
\ETABLISSEMENT{Université Paul Sabatier Toulouse III}
\ANNEE{2005}
\MEL{chafai@math.ups-tlse.fr}
\AUTEUR{B. Bercu \& D. Chafaï}
\WEB{http://www.lsp.ups-tlse.fr/Chafai/agregation.html}
\DATE{Novembre 2004}

\begin{document}

%%
\section{Vecteurs aléatoires gaussiens}
%%

\begin{defi}
  Soit $X$ un vecteur aléatoire défini sur un espace de probabilité $(\Om,
  \cA, \dP)$, à valeurs dans $\dR^d$ avec $d\geq 1$.  $X$ est dit gaussien si
  toute combinaison linéaire de ses composantes est une v.a. gaussienne.
\end{defi}

\begin{thm} 
  La loi d'un vecteur aléatoire gaussien est entièrement déterminée par son
  vecteur espérance $m=\dE(X)\in\dR^d$ et sa matrice de covariance
  $\Ga=\dE((X-m)(X-m)^\top)\in\dS_d^+$ où $\dS_d^+$ est le cône convexe des
  matrices carrées $d\times d$ symétriques semi-définies positives (pas
  forcément inversibles). Soit $X\sim \cN(m,\Ga)$. Alors, pour tout $u\in
  \dR^d$, on a 
  $$
  \dE\PAR{\exp(i\DP{u,X}} =
  \exp\PAR{i\DP{u,m}-\frac{1}{2}u^\top\Ga u}. 
  $$  
\end{thm}

\begin{thm}
  Soit $X\sim \cN(m,\Ga)$. $X$ admet une densité $f_X$ par rapport à la mesure
  de Lebesgue de $\dR^d$ si et seulement si $\Ga$ est inversible et l'on a 
  $$
  f_X(x)
  =((2\pi)^d\det\Ga)^{-1/2}\exp\PAR{-\frac{1}{2}(x-m)^\top\Ga^{-1}(x-m)}.
  $$
\end{thm}

La loi particulière $\cN(0,\rI_d)$ où $\rI_d$ est la matrice identité de
$\dR^d$ est appelée \emph{gaussienne standard}. Dans tout ce texte,
l'espérence des vecteurs et matrices est prise composante par composante.
Voici quelques propriétés fondamentales des vecteurs gaussiens, cf. par
exemple \cite[IV.4 et VI.4]{barbe-ledoux}

\begin{enumerate}
\item[\textbf{P1.}] Soit $X=(X_1,\ldots,X_d)$ un vecteur gaussien de matrice
  de covariance $\Ga$. Les trois propositions suivantes sont équivalentes :
  (a) $(X_1,\ldots,X_d)$ sont deux à deux indépendantes, (b)
  $(X_1,\ldots,X_d)$ sont indépendantes dans leur ensemble, (c) La matrice
  $\Ga$ est diagonale.
\item[\textbf{P2.}] Soit $m \in \dR^d$, $A\in\cM_d(\dR)$ et
  $X\sim\cN(0,\rI_d)$, alors\footnote{Réciproquement, une matrice de
    covariance $\Ga\in\dS_d^+$ peut toujours s'écrire $AA^\top$ où
    $A\in\cM_d(\dR)$. La matrice $A$ n'est pas unique. La méthode de Choleski
    en fournit une (rapide), la racine carrée matricielle obtenue par
    diagonalisation en base orthonormée de $\Ga$ également (lent). Ce procédé
    ne requière pas l'inversibilité de $\Ga$. } $m+AX\sim\cN(m,AA^\top)$.
\item[\textbf{P3.}] Soit $Z=(X,Y)$ un vecteur gaussien de $\dR^{d+1}$ avec
  $X=(X_1,\ldots,X_d)$ d'espérance $m$ et de matrice de covariance inversible
  $\Ga$. Alors, la loi conditionnelle de $Y$ sachant $X$ est gaussienne
  d'espérance affine en $X$ $\dE(Y|X)=a+b^\top X=\dE(Y)+b^\top(X-m)$ et de
  variance $\var{}{Y|X}=\var{}{Y}-b^\top\Ga b$ avec $a=\dE(Y)-b^\top m$ et
  $b=\Ga^{-1}\cov{}{X}{Y}$. De plus, $\veps=Y-\dE(Y|X)$ est indépendante de
  $X$.
\end{enumerate}

\begin{exo}[Algorithme de Box-Muller]\label{exo:bm}
  Soit $(X,Y)$ un vecteur aléatoire de $\dR^2$. Montrer que $(X,Y)$ suit la
  loi normale $\cN(0,\rI_2)$ si et seulement si $X=r\cos\te$ et $Y=r\sin\te$
  où $r$ et $\te$ sont deux variables aléatoires indépendantes avec $r^2$ de
  loi exponentielle $\cE(1/2)$ et $\te$ de loi uniforme $\cU([0,2\pi])$. En
  déduire un algorithme permettant de générer des réalisations de variables
  aléatoires indépendantes gaussiennes $\cN(m,\Ga)$. Cf. \cite{ycart}.
\end{exo}

\begin{exo}
  Étudier et expliquer le programme suivant
\begin{lstlisting}[language=Matlab]{}
N=input('Entrez la taille de l'échantillon N : '); 
max=round(3*N/2);
m=input('Précisez la valeur de la moyenne m : ');
sigma=input('Précisez la valeur de écart type : ');
X=2*rand(max,1)-ones(max,1); Y=2*rand(max,1)-ones(max,1);
S=X.^2+Y.^2; X=X(find(S<1)); Y=Y(find(S<1));
r=sqrt(X.^2+Y.^2); R=2*sqrt(-log(r))./r; Z=R(1:N).*X(1:N); 
T=m*ones(N,1)+sqrt(sigma^2)*Z;
\end{lstlisting}
% 
\end{exo}

\section{Bref rappel sur les échantillons gaussiens}

\noindent
\textbf{Loi du $\chi^2$.} La loi $\chi^2(n)$ est la loi de
$Z_1^2+\cdots+Z_n^2=\NRM{(Z_1,\ldots,Z_n)}_2^2$ où $Z_1,\ldots,Z_n$ sont des
v.a.r. i.i.d. de loi $\cN(0,1)$. Sa moyenne est $n$, sa variance $2n$, et sa
densité est
$x\mapsto(2^{n/2}\Ga(n/2))^{-1}x^{n/2-1}\exp\PAR{-x/2}\,\rI_{\dR_+}(x)$. C'est
donc une loi Gamma particulière de paramètres $\la=1/2$ et $a=n/2$. Il est
clair que $\chi^2(n)*\chi^2(m)=\chi^2(n+m)$.

\noindent
\textbf{Loi de Student.} Soit $A$ et $B$ des v.a.r. indépendantes de loi
respectives $\cN(0,1)$ et $\chi^2(n)$, alors la loi de Student $t(n)$ est la
loi de $\sqrt{n}A/\sqrt{B}$. La loi $t(n)$ est centrée, de densité
$x\mapsto\frac{\Ga((n+1)/2)}{\Ga(n/2)\sqrt{n\pi}}(1+x^2/n)^{-(n+1)/2}$ sur
$\dR$, et de variance $n/(n-2)$. On a enfin que $t(n)$ converge étroitement
vers $\cN(0,1)$ quand $n$ tend vers $+\infty$.

\noindent
\textbf{Loi de Fisher(-Snedecor).} Soit $A$ et $B$ deux v.a.r. indépendantes
de lois respectives $\chi^2(n)$ et $\chi^2(m)$, alors la loi de
Fisher(-Snedecor) $F(n,m)$ est la loi de $(A/n)/(B/m)$. Sa moyenne et sa
variance valent respectivement $m/(m-2)$ et
$\frac{2m^2(m+n-2)}{n(m-4)(m-2)^2}$ et sa densité est
$x\mapsto\frac{\Ga((n+m)/2)n^{n/2}m^{m/2}}{\Ga(n/2)\Ga(m/2)}
\frac{x^{2/n-1}}{(m+nx)^{(n+m)/2}}\,\rI_{\dR_+}(x)$. Il est immédiat que si
$Z\sim t(n)$, alors $Z^2\sim F(1,n)$. De plus, si $Z\sim F(n,m)$ alors $1/Z
\sim F(m,n)$.

\noindent
\textbf{Moyenne et variance empiriques.}  Soit $Z_1,\ldots,Z_n$ des
v.a.r. i.i.d. de loi $\cN(m,\si^2)$. Alors les estimateurs de la
moyenne empirique $\WH{m}$ et de la variance empirique $\WH{\si^2}$
définis par
\begin{equation*}
\WH{m}:=\frac{Z_1+\cdots+Z_n}{n}
\text{\quad et \quad}
\WH{\si^2}:=\frac{1}{n-1}\sum_{i=1}^n (Z_i-\WH{m})^2
\end{equation*}
sont indépendants et sans biais. De plus, on a
\begin{equation*}
  \frac{1}{\si^2}\sum_{i=1}^n (Z_i-m)^2 \sim \chi^2(n)
  \text{\quad et \quad}
  (n-1)\,\frac{\WH{\si^2}}{\si^2} \sim \chi^2(n-1).
\end{equation*}   
et
\begin{equation*}
  \sqrt{n}\,\frac{\WH{m}-m}{\si}\sim\cN(0,1) 
  \text{\quad et \quad}
  \sqrt{n}\,\frac{\WH{m}-m}{\sqrt{\WH{\si^2}}}\sim t(n-1)
  \text{\quad et \quad}
  n\,\frac{(\WH{m}-m)^2}{\WH{\si^2}} \sim F(1,n-1).
\end{equation*}   
Pour tout ce qui concerne ces lois et les échantillons gaussiens, on pourra
consulter par exemple \cite[Chap. 5]{dacunha-castelle-duflo}.

\begin{thm}[de Cochran]\label{th:cochran} 
  Soit $Y=(Y_1,\ldots,Y_n)$ des v.a.r. i.i.d. de loi $\cN(0,\si^2)$, alors
  l'image de $(Y_1,\ldots,Y_n)$ par toute transformation orthogonale de
  $\dR^n$ est encore un vecteur gaussien de loi $\cN(0,\si^2\rI_n)$. De plus,
  les projections orthogonales $P_{E_1}(Y),\ldots,P_{E_k}(Y)$ de $Y$ sur des
  sous-espaces vectoriels $E_1,\ldots,E_k$ de $\dR^n$ deux à deux orthogonaux
  sont indépendantes, et $\si^{-2}\NRM{P_{E_i}(Y)}_2^2$ suit la loi
  $\chi^2(\dim(E_i))$ pour tout $1\leq i\leq p$.
\end{thm}

\begin{cor}[Du théorème de Cochran]\label{co:cochran}
  Soit $\veps_1,\ldots,\veps_n$ des v.a.r. i.i.d. de loi $\cN(0,1)$ et $V$ un
  sous-espace vectoriel de dimension $p$ de $\dR^n$. On observe $Y=m+\si\veps$
  où $\veps:=(\veps_1,\ldots,\veps_n)$ et où $m\in V$ et $\si^2>0$ sont
  inconnus. Si $P_V\in\mathrm{End}(\dR^n)$ désigne la projection orthogonale
  sur $V$, on a
  \begin{enumerate}
  \item $\WH{m}:=P_V(Y)$ est un estimateur sans biais de $m$;
  \item La variance résiduelle $\WH{\si^2}:=(n-p)^{-1}\NRM{Y-P_V(Y)}_2^2$ est
    un estimateur sans biais de $\si^2$;
  \item Les vecteurs aléatoires $Y-P_V(Y)$ et $P_V(Y)$ sont indépendants, et
    on a $\si^{-2}\NRM{P_V(Y)-m}_2^2\sim\chi^2(p)$ et
    $\si^{-2}\NRM{Y-P_V(Y)}_2^2\sim\chi^2(n-p)$.
  \end{enumerate}
\end{cor}

%%
\section{Modèle linéaire gaussien}
%%

Soit $(X,Y)$ un vecteur aléatoire de $\dR^{p-1}\times\dR$. On suppose que la
loi de $X$ admet une densité de probabilité $\vphi$ par rapport à une mesure
positive $\mu$ sur $\dR^{p-1}$. On pensera par exemple aux cas particuliers où
$\mu$ est la mesure de Lebesgue sur un borélien de $\dR^{p-1}$, ou la mesure
de comptage sur un sous-ensemble dénombrable de $\dR^{p-1}$. On considère une
modélisation gaussienne de la loi conditionnelle $\cL(Y\,\vert\,X)$ sous la
forme
\begin{equation*}
  \cL\PAR{Y\,\vert\,X=x}
  = \cN\PAR{\be_1 f_1(x_1)+\cdots+\be_{p-1} f_{p-1}(x_{p-1})+\be_p\,,\,\si^2}
\end{equation*}
où les fonctions réelles de la variable réelle $f_1,\ldots,f_{p-1}$ sont
connues, et où les paramètres $\be:=(\be_1,\ldots,\be_p)\in\dR^p$ et
$\si^2\in\dR_+^*$ sont inconnus. La moyenne est linéaire en $\be$, d'où le
qualificatif de \emph{modèle linéaire gaussien}. On peut écrire
$\cL\PAR{Y\,\vert\,X=x}=\cN(\mathbf{x}^\top \be,\si^2)$ où $\mathbf{x} :=
(f_1(x_1),\ldots,f_{p-1}(x_{p-1}),1)^\top$. Partant de l'observation d'un
échantillon de taille $n\geq p$ noté $(x_1,y_1),\ldots,(x_n,y_n)$ du vecteur
$(X,Y)$, l'objectif est d'estimer $\be$ et $\si^2$ et d'effectuer des tests
d'hypothèse sur $\be$. On ne s'intéresse pas à la loi de $X$ ici. En
particulier, $\vphi$ restera inconnue et inutile. L'espérence conditionnelle
de $Y$ sachant $X$ est donnée pour tout $x\in\dR^{p-1}$ par
\begin{equation*}
  \moy{}{Y\,\vert\,X=x}
  = \be_1 f_1(x_1) + \cdots + \be_{p-1} f_{p-1}(x_{p-1}) + \be_p.
\end{equation*} 
L'ensemble paramétrique associé est
$\Te:=\{(\be,\si^2)\in\dR^p\times\dR_+^*\}$. Le nombre entier $p$ correspond à
la dimension du paramètre $\be$. On a coutume en statistique appliquée
d'écrire ce modèle conditionnel sous la forme
\begin{equation*}
  Y_i = \be_1 f_1(x_{i,1})+\cdots+\be_{p-1} f_{p-1}(x_{i,p-1})+\be_p+\veps_i,
\end{equation*} 
où les $\veps_i$ sont i.i.d. de loi $\cN(0,\si^2)$ et où les $x_{i,j}$ sont
des «variables déterministes» connues. Cette écriture n'est pas abusive si
l'on considère que les $Y_i$ sont des v.a.r. indépendantes, la loi de $Y_i$
étant précisément $\cN\PAR{\be_1
  f_1(x_1)+\cdots+\be_{p-1}f_{p-1}(x_{p-1})+\be_p\,,\,\si^2}$. Le paramètre
$\be_p$ est appelé \emph{intercept}. On parle de \emph{modèle linéaire simple}
lorsque $p=2$, et on écrit alors $Y_i=\be_1 f(x_i)+\be_2+\veps_i$. On parle de
\emph{modèle linéaire multiple} quand $p\geq 2$. De manière équivalente, on
parle de \emph{régression linéaire simple} et de \emph{régression linéaire
  multiple}.

Pour le modèle linéaire simple ($p=2$), l'estimateur des moindres carrés pour
$\be$ consiste à minimiser $(\be_1,\be_2)\mapsto\sum_{i=1}^n (y_i-\be_1 f(x_i)
-\be_2)^2$. Cela donne en notant $u_i=f(x_i)$, $\OL{u}:=(u_1+\cdots+u_n)/n$ et
$\OL{y}:=(y_1+\cdots+y_n)/n$
\begin{equation*}
  \WH{\be_1} = 
  \frac{\sum_{i=1}^n (u_i-\OL{u})(y_i-\OL{y})}
       {\sum_{i=1}^n \PAR{u_i-\OL{u}}^2}
  \text{\quad et \quad}
  \WH{\be_2} = \OL{y} - \WH{\be_1}\,\OL{u}.
\end{equation*} 
Ces estimateurs sont sans biais. Les \emph{valeurs ajustées}
$\WH{y_i}:=\WH{\be_1}f(x_i)+\WH{\be_2}$ sont les estimées de $\moy{}{Y_i\vert
  X=x_i}$. Les \emph{résidus} correspondants sont les $e_i:=y_i-\WH{y_i}$. La
courbe d'équation $y=\be_1 f(x)+ \be_2$ est appelée \emph{courbe de
  regression}. On observera que la nature gaussienne du problème n'intervient
pas dans la déduction des estimateurs des moindres carrés. En revanche, c'est
elle qui explique qu'il s'agit des estimateurs du maximum de vraisemblance.

\noindent
\textbf{Expression des estimateurs.} Le corollaire \ref{co:cochran} donne des
estimateurs de $\be$ et $\si^2$. En effet, il suffit de considérer le
sous-espace vectoriel $V$ de $\dR^n$ engendré par les colonnes de la matrice
$\bX\in\cM_{n,p}(\dR)$ définie par
\begin{equation*}
  \bX=
  \left(
  \begin{array}{cccc}
    f_1(x_{1,1}) & \cdots & f_{p-1}(x_{1,p-1}) & 1 \\
    \vdots       &        & \vdots       & \vdots \\
    f_1(x_{n,1}) & \cdots & f_{p-1}(x_{n,p-1}) & 1
  \end{array}
  \right).
\end{equation*}
On a alors $V=\{\bX\be,\,\be\in\dR^p\}=\mathrm{Im}(\bX)$. Dans la base
canonique de $\dR^n$, la matrice de la projection orthogonale $P_V$ s'écrit
$\bX(\bX^\top\bX)^{-1}\bX^\top$. L'estimateur $\WH{m}$ de $m=\bX\be$ s'écrit
alors $\WH{m}=P_V(Y)=\bX(\bX^\top\bX)^{-1}\bX^\top Y$. Lorsque la matrice
$\bX$ est de rang plein $p$, et le paramètre $\be$ peut être estimé par
$\WH{\be}:=(\bX^\top\bX)^{-1}\bX^\top Y$. Comme $m=\bX\be$ et comme
$\bX\WH{\be}=P_V(Y)$, on a par linéarité que
$\WH{\be}\sim\cN(\be,(\bX^\top\bX)^{-1}\si^2)$. On a également
$(n-p)\WH{\si^2}/\si^2\sim\chi^2(n-p)$. Les estimateurs $\WH{\be}$ et
$\WH{\si^2}$ sont indépendants et sans biais. On note $\WH{\si}$ la racine
carrée de $\WH{\si^2}$. Lorsque $p=2$, on retrouve bien entendu les formules
du modèle linéaire simple pour $\WH{\be_1}$ et $\WH{\be_2}$.

\noindent
\textbf{Test sur une hypothèse linéaire.} Soit $W$ un sous-espace vectoriel de
$V$, de dimension $p-q<p$. On définit les hypothèses antagonistes
$\text{$\cH_0$:«$m\in W$»}$ et $\text{$\cH_1$:«$m\in V\setminus W$»}$. Soit
$U$ le sous-espace vectoriel de $V$ tel que $V=W\perp U$. Le théorème de
Cochran \ref{th:cochran} appliqué au vecteur gaussien
$Y-m\sim\cN(0,\si^2\,\rI_n)$ et à la décomposition $\dR^n=W\oplus U\oplus
V^\perp$ entraîne alors en particulier que les vecteurs aléatoires
$P_{V^\perp}(Y-m)$ et $P_U(Y-m)$ sont indépendants, et que
$\si^{-2}\NRM{P_{V^\perp}(Y-m)}_2^2$ et $\si^{-2}\NRM{P_U(Y-m)}_2^2$ ont pour
lois respectives $\chi^2(n-p)$ et $\chi^2(q)$. Or on a $P_V-P_W=P_U$ et
$Id-P_V=P_{V^\perp}$. En particulier, $P_{V^\perp}(m)=0$ donc
$P_{V^\perp}(Y-m)=Y-P_V(Y)$, et sous $\cH_0$, $P_U(m)=0$ donc
$P_U(Y-m)=P_V(Y)-P_W(Y)$. On en déduit que sous $\cH_0$
\begin{equation}\label{eq:statQ}
  Q=\frac{(n-p)}{q}
  \frac{\NRM{P_V(Y)-P_W(Y)}_2^2}{\NRM{Y-P_V(Y)}_2^2}\sim F(q,n-p).
\end{equation}
Cela donne un test de niveau $\al$ consistant à rejetter $\cH_0$ dès que
$Q>F_\al(q,n-p)$ où $F_\al(q,n-p)$ est le quantile $\al$ de la loi de Fisher
$F(q,n-p)$. En considérant une définition matricielle de $W$, il est possible
d'exprimer très simplement $Q$ en fonction de $\be$, cf. \cite[Chap.
6.4.3]{paul-toulouse}. La statistique $Q$ est à une transformation croissante
près celle du rapport de maximum de vraisemblance, cf. \cite[Chap.
6.4.2]{paul-toulouse}.

\noindent
\textbf{Test de validité du modèle.} La première question que l'on se pose
dans la pratique est de savoir si les $x_i$ influent sur la moyenne de $Y$.
Cela revient exactement à considérer l'hypothèse linéaire
$\cH_0$:«$\be_1=\cdots=\be_{p-1}=0$» qui correspond à $q=p-1$ et donc à
considérer la statistique de test \eqref{eq:statQ}. Comme dans ce cas
$W=\dR\mathbf{1}$ où $\mathbf{1}:=(1,\ldots,1)^\top$, on a
$p_W(Y)=\frac{1}{n}(Y_1+\cdots+Y_n)\mathbf{1}=:\OL{Y}\mathbf{1}$ et
$p_V(Y)=\bX\WH{\be}=:\WH{Y}$. Cela conduit à la statistique de test
\begin{equation*}
  Q=\frac{(n-p)}{(p-1)}
  \frac{\sum_{i=1}^n (\WH{Y}_i-\OL{Y})^2}{\sum_{i=1}^n (Y_i-\WH{Y}_i)^2}
  =\frac{(n-p)}{(p-1)}
  \frac{\text{variance expliquée}}{\text{variance totale}}.
\end{equation*}
Le théorème de Pythagore correspond donc à la fameuse «formule de
décomposition de la variance»
\begin{equation*}
  \underbrace{
    \frac{1}{n}
    \sum_{i=1}^n 
    (Y_i-\OL{Y})^2}_{\text{variance totale}}
  =    
  \underbrace{
    \frac{1}{n}
    \sum_{i=1}^n 
    (\WH{Y}_i-\OL{Y})^2}_{\text{variance «expliquée»}}
  +  
  \underbrace{
    \frac{1}{n}
    \sum_{i=1}^n 
    (Y_i-\WH{Y}_i)^2}_{\text{variance «résiduelle»}},
\end{equation*}
et en terme de lois : $\chi^2(n-1)=\chi^2(p-1)*\chi^2(n-p)$. Le praticien
préfère parfois disposer d'un indice plutôt que d'un test. Un indice
fréquement utilisé pour quantifier la qualité de la régression est l'indice
$R^2:=1-(\text{variance résiduelle}/\text{variance totale})=\text{variance
  expliquée}/\text{variance totale}$. On a toujours $R^2\in[0,1]$, et $R^2$
est d'autant plus proche de $1$ que la régression est «explicative» au sens de
la variance.

\noindent
\textbf{Intervalle de confiance pour une composante $\be_i$.} Comme $\WH{\be}$
et $(n-p)\WH{\si^2}/\si^2$ sont indépendants et de lois respectives
$\cN(\be,(\bX^\top\bX)^{-1}\si^2)$ et $\chi^2(n-p)$, on en déduit que
\begin{equation*}
  \frac{(\WH{\be}_j-\be_j)}{\sqrt{(\bX^\top\bX)^{-1}_{j,j}\,\WH{\si^2}}}
  \sim t(n-p).
\end{equation*}  
On peut ainsi effectuer des tests d'hypothèse sur une composante donnée de
$\be$. On a par exemple pour $\be_j$ l'intervalle de confiance exact de niveau
$(1-\al)$ suivant
\begin{equation*}
  I_j(\al):= \WH{\be}_j 
  +\PAR{t_{\al/2}(n-p)\sqrt{(\bX^\top\bX)^{-1}_{j,j}\,\WH{\si^2}}}\SBRA{-1,+1},
\end{equation*}  
où $t_{\al/2}(n-p)$ est le quantile $\al/2$ de la loi de Student $t(n-p)$.
Ainsi, on peut construire un test de niveau $\al$ pour $\cH_0:\{\be_j=0\}$ en
acceptant $\cH_0$ dès que $0\in I_j(\al)$. Sous $\cH_0$, $\dP(0\not\in
I_j(\al))=1-\al$. La $p$-value de ce test, i.e. la plus grande valeur de
$\alpha$ pour laquelle $\cH_0$ est acceptée, se calcule très facilement en
distinguant deux cas (selon le signe de $\WH{\be_j}$). On prendra garde au
fait que ces intervalles de confiance $(I_j(\al))_{1\leq j\leq p}$ ne sont pas
indépendants, et il en va de même pour les tests associés, d'où l'intérêt de
la statistique \eqref{eq:statQ}.

\noindent
\textbf{Intervalle de confiance pour $\si^2$.} Comme
$\si^{-2}\NRM{Y-P_V(Y)}_2^2\sim\chi^2(n-p)$, on en déduit qu'avec probabilité
$1-\al$, le paramètre $\si^2$ se trouve dans l'intervalle
\begin{equation*}
  \NRM{Y-P_V(Y)}_2^2\,
  \SBRA{\frac{1}{\chi^2_{1-\al/2}(n-p)}\,,\,\frac{1}{\chi^2_{\al/2}(n-p)}}
\end{equation*}
où $\chi^2_{\tau}(n-p)$ est le quantile $\tau$ de la loi $\chi^2(n-p)$.

\noindent
\textbf{Région de confiance pour $\bH\be$.} Rappelons que $m=\bX\be$. Comme
$\si^{-2}\NRM{P_V(Y)-m}_2^2$ et $\si^{-2}\NRM{P_V(Y)-Y}_2^2$ sont
indépendantes et de lois respectives $\chi^2(p)$ et $\chi^2(n-p)$, on en
déduit que
\begin{equation*}
  \frac{(n-p)}{p}\frac{\NRM{P_V(Y)-m}_2^2}{\NRM{Y-P_V(Y)}_2^2}
  \sim F(p,n-p).
\end{equation*}
Ainsi, avec probabilité $1-\al$, le paramètre $m$ est dans la boule
$B_2(P_V(Y);r_\al)$ de rayon
\begin{equation*}
  r_\al=\NRM{Y-P_V(Y)}_2\sqrt{\frac{p}{n-p}F_\al(p,n-p)}
\end{equation*}
où $F_\al(p,n-p)$ est le quantile $\al$ de la loi $F(p,n-p)$. La traduction de
cette région en terme de $\be$ et $\bX$ est immédiate puisque $m=\bX\be$ et
$P_V(Y)=\bX(\bX^\top\bX)^{-1}\bX^\top Y=\bX\WH{\be}$.
Si $\bH$ est une matrice $q\times p$ de rang $q<p$, la région de confiance
pour $\bH\be$ est donnée par
\begin{equation*}
  \BRA{\be\in\dR^p,\,
    (\WH{\be}-\be)^\top\bH^\top
    (\bH(\bX^\top\bX)^{-1}\bH^\top)^{-1}\bH
    (\WH{\be}-\be)
    \leq q\,\WH{\si^2}\,F_\al(q,n-p)}.
\end{equation*}

\noindent
\textbf{Lien avec le maximum de vraisemblance et formulation conditionnelle.}
La vraisemblance s'écrit
\begin{equation*}  
  L_n(\be,\si^2)
  =(2\pi\si^2)^{-n/2}\exp\PAR{-\frac{1}{2\si^2}\NRM{\bY-\bX\be}_2^2}
  \prod_{i=1}^n\vphi(x_{i,1},\ldots,x_{i,p}),
\end{equation*}
où $\bY:=(y_1,\ldots,y_n)^\top$ et où $\bX:=(f_j(x_{i,j}))_{1\leq i\leq
  n,1\leq j\leq p}$. L'estimateur du maximum de vraisemblance
$(\WH{\be},\tilde{\si^2})\in\Om$ de $(\be,\si^2)$ est défini par
\begin{equation*}
  (\WH{\be},\tilde{\si^2}) 
  := \arg\sup_{(\be,\si^2)\in\Om} L_n(\be,\si^2)
  = \arg\sup_{(\be,\si^2)\in\Om} \log L_n(\be,\si^2).
\end{equation*}
On montre sans difficulté qu'il est bien défini et qu'il ne dépend pas
de $\vphi$, i.e. de la loi de $X$.  Lorsque la matrice $\bX$ est de
rang plein, on a les formules explicites suivantes obtenues en
annulant le gradient de $L_n$
\begin{equation*}  
  \WH{\be} = (\bX^\top\bX)^{-1}\bX\bY 
  \text{\quad et \quad}
  \tilde{\si^2} = n^{-1}\NRM{\bY-\bX\WH{\be}}_2^2.
\end{equation*} 
La matrice $\bX(\bX^\top\bX)^{-1}\bX$ correspond à la projection orthogonale
dans $\dR^n$ sur l'espace vectoriel $\mathrm{Im}(\bX)$ des moyennes
admissibles. Ainsi, $\bX\WH{\be}$ est la projection orthogonale de $\bY$ sur
$\mathrm{Im}(\bX)$. L'estimateur $\WH{\be}$ coïncide avec l'estimateur des
moindres carrés obtenu en minimisant $z\in\dR^p\to\NRM{\bY-\bX z}_2^2$. On a
de plus $\WH{\be}\sim\cN(\be,\si^2(\bX^\top\bX)^{-1})$, et $\WH{\be}$ est donc
sans biais. En revanche, l'estimateur $\tilde{\si^2}$ fait intervenir
$\WH{\be}$ au lieu de $\be$, ce qui explique qu'il soit biaisé. Sa version
débiaisée est donnée par $\WH{\si^2} := (n/(n-p))\tilde{\si^2}$. On a alors
$(n-p)\WH{\si^2}/\si^2\sim\chi^2(n-p)$. L'approche conditionnelle et le lien
avec le maximum de vraisemblance à le mérite de placer le modèle linéaire
gaussien dans le cadre de la statistique asymptotique paramétrique.

\noindent
\textbf{Intervalle de confiance sur l'espérance conditionnelle et intervalle
  de prédiction sur la loi conditionnelle.} On sait que pour tout
$x\in\dR^{p-1}$, $\cL(Y\,\vert\,X=x)=\cN(\mathbf{x}^\top\be,\,\si^2)$ et
$\moy{}{Y\,\vert\,X=x}=\mathbf{x}^\top\be$, où
$\mathbf{x}:=(f_1(x_1),\ldots,f_{p-1}(x_{p-1}),1)^\top$. Pour un niveau $\al$,
on peut donner un intervalle de confiance qui contiendra la moyenne
$\mathbf{x}^\top\be$ avec probabilité $1-\al$. D'un autre côté, on peut
également donner un intervalle de prédiction qui sera de probabilité $1-\al$
pour la loi conditionnelle $\cL(Y\,\vert\,X=x)$. Ces deux intervalles sont
différents\footnote{On peut par exemple retenir que les intervalles de
  confiance concernent des paramètres (moyenne,\ldots) tandisque les
  intervalles de prédiction concernent des lois (ou des v.a.). Ces intervalles
  ne sont pas uniques pour un niveau donné. Idéalement, on cherchera toujours
  pour un niveau donné les intervalles les plus petits.}.

Pour construire un intervalle de confiance sur la moyenne, on écrit que
$\mathbf{x}^\top\WH{\be}\sim\cN(\mathbf{x}^\top\be,\si^2\mathbf{x}^\top(\bX^\top\bX)\mathbf{x})$.
Cette v.a. est indépendante de $\WH{\si^2}$. On obtient donc que
$(x^\top\WH{\be}-x^\top\be)/\WH{\si}\sim t(n-p)$. Ainsi, avec probabilité
$1-\al$, $\moy{}{Y\,\vert\,X=x}$ se trouve dans l'intervalle
\begin{equation*}
  \mathbf{x}^\top\WH{\be}+
  t_{\al/2}(n-p)\sqrt{(\mathbf{x}^\top(\bX^\top\bX)^{-1}\mathbf{x})\WH{\si^2}}
  \,\,[-1,+1],
\end{equation*}
où $t_{\al/2}(n-p)$ est le quantile $\al/2$ de la loi de Student $t(n-p)$.

Pour construire un intervalle de prédiction sur la loi conditionnelle, on
considère une variable aléatoire $Y_x$ de loi $\cL(Y\,\vert\,X=x)$
indépendante de l'échantillon. On a alors
\begin{equation*}
  \mathbf{x}^\top\WH{\be}-Y_x
  \sim\cN(0,\si^2 \mathbf{x}^\top(\bX^\top\bX)^{-1}\mathbf{x}) * \cN(0,\si^2)
  =\cN(0,\si^2(1+\mathbf{x}^\top(\bX^\top\bX)^{-1}\mathbf{x})).
\end{equation*}
Par suite, par indépendance de $\WH{\be}$, $\WH{\si^2}$, et $Y_x$, il vient
\begin{equation*}
  \frac{\mathbf{x}^\top\WH{\be}-Y_x}
  {\sqrt{(1+\mathbf{x}^\top(\bX^\top\bX)^{-1}\mathbf{x})\WH{\si^2}}}
  \sim t(n-p).
\end{equation*}
Ainsi, avec une probabiltié $1-\al$, la variable aléatoire $Y_x$ prend sa
valeur dans l'intervalle
\begin{equation*}
  \mathbf{x}^\top\WH{\be}
  +t_{\al/2}(n-p)\sqrt{(1+\mathbf{x}^\top(\bX^\top\bX)^{-1}\mathbf{x})\WH{\si^2}}
  \,\,[-1,+1].
\end{equation*}

\noindent
\textbf{Résidus et contrôle.} Si le modèle est exact, le vecteur des résidus
$\bY-\bX\WH{\be}$ est gaussien centré de matrice de covariance
$\si^2(\rI_n-\bX^\top(\bX^\top\bX)^{-1}\bX)$. Nous avons supposé que les
composantes de $\bY-\bX\be$ étaient de moyenne nulle, indépendantes, de
variance constante, et de loi gaussienne. Ces hypothèses sont rarement toutes
satisfaites dans la pratique malheureusement. On peut par exemple tracer un
\texttt{qqplot} des quantiles de la la gaussienne versus les quantiles
empiriques des résidus (droite de Henri), cf. par exemple \cite[Chap.
15.4.1.C]{saporta}. On peut alors s'intéresser aux résidus studentisés, cf.
par exemple \cite[Chap. 17.3]{saporta}.

\noindent
\textbf{Références bibliographiques.} Par exemple \cite[Chap. 15 et
16]{saporta}, \cite[Chap. 5]{dacunha-castelle-duflo}, \cite[Chap.
6]{paul-toulouse}, \cite{guyon}, \cite{scheffe}.

\noindent
\textbf{Les données et outils de Stixbox.} Stixbox possède une collection de
jeux de données empiriques. La commande \texttt{getdata} permet d'y accéder,
et en donne également un bref descriptif. Stixbox fournit également trois
fonctions utiles pour le modèle linéaire, qui sont
\begin{itemize}
\item \texttt{identify}. Permet d'itentifier des points sur un nuage de point
  issu d'un échantillon. Un click sur le bouton du mileu fait sortir.
\item \texttt{linreg}. Effectue une \emph{régression linéaire
    simple}\footnote{ \texttt{linreg} permet également de faire de la
    régression linéaire polynomiale.}. Un graphique est automatiquement
  produit. Ce graphique comporte la courbe de régression de $y$ en fonction de
  $x$, les zônes de confiance et de prédiction du niveau désiré, ainsi que le
  nuage de points du jeu de données.
\item \texttt{lsfit}. Effectue une \emph{régression linéaire multiple}. La
  syntaxe est identique à \texttt{linreg}, pour les trois paramètres
  obligatoires. Aucun graphique n'est produit.
\item \texttt{cmpmod}. Quantifie la différence entre un modèle et un
  sous-modèle pour une \emph{régression linéaire multiple}, au moyen d'un test
  de Fisher.
\end{itemize}
Il est très instructif de lire le code de ces fonctions en utilisant la
commande \texttt{type}.

\begin{exo}[Régression linéaire simple]  
  Créer un échantillon de taille $100$ du couple $(X,Y)$ de telle sorte que
  $\cL(X)=\cU([0,1])$ et $\cL(Y\,\vert\,X=x)=\cN(1+x,\,1)$. Étudier la
  régression linéaire simple $Y_i=1+x_i+\veps_i$ au moyen de \texttt{linreg}.
  Les résultats sont-ils compatibles avec les vraies valeurs ? Voici un
  programme à titre d'exemple et de correction. \MFILE{exolinreg}
\end{exo}

\FIG{linreg}{0.4}{htbp}{Graphique de la régression linéaire de l'exercice
  \ref{code:exolinreg}, comprenant la droite de regression, le nuage de
  points, et les régions de confiance et de prédiction (cette dernière étant
  la plus large).}

\FIG{residus}{0.4}{htbp}{Résidus de la régression linéaire de l'exercice
  \ref{code:exolinreg}.}

\begin{exo}[Régression linéaire simple]
  On s'intéresse au jeu de données numéro 13 de Stixbox. Il correspond à un
  échantillon de taille 14 sur deux variables aléatoires : le nombre de
  composants électriques en panne dans un ordinateur et la durée de la
  réparation en minutes. Visualisez les données avec \texttt{identify}.
  Effectuez une régression linéaire simple $\text{durée}_i=\be_2+\be_1
  \text{nombre}_i + \veps_i$ au moyen de \texttt{linreg}. Tester l'hypothèse
  de nullité de $\be_1$ à $5\%$. Idem pour $\be_2$. Conclusion ? Que vous
  inspire le graphe des résidus ? Sont-ils corrélés ? Nous avons supposé que
  la durée de réparation sachant le nombre de composant en panne était
  gaussienne, cela vous semble-t-il plausible ? Peut-on effectuer la
  régression inverse ?
\end{exo}

\begin{exo}[Régression linéaire multiple]
  On souhaite étudier la variation du taux d'hémoglobine dans le sang au cours
  d'une opération chirurgicale en fonction de la durée de l'opération et du
  volume de sang perdu pendant l'opération. On dispose des résultats suivants
  où $y_i$ représente la valeur observée en pourcentage de la variation du
  taux d'hémoglobine, $x_{i,1}$ est la durée de l'opération en heures
  décimales et $x_{i,2}$ est le volume en litres de sang perdu.
\par\medskip
\begin{tabular}{|c|cccccccc|}\hline
  $y_i$ & -1.70 & -4.61 & -5.82 & -1.17 & -4.23 & -3.31 & +0.42 & -2.98\\
  $x_{i,1}$ & 1.75 & 1.33 & 1.43 & 1.86 & 1.81 & 1.66 & 1.60 & 2.00\\
  $x_{i,2}$ & 0.52 & 0.59 & 0.61 & 0.50 & 0.54 & 0.49 & 0.27 & 0.47\\
  \hline
\end{tabular}
\par\medskip
On suppose que $y_i$ est une réalisation d'une variable aléatoire $Y_i$ de loi
$\cN(\be_3+\be_2x_{i,2}+\be_1 x_{i,1},\si^2)$. Étudier cette régression
linéaire multiple grâce à \texttt{lsfit}. Tester l'hypothèse suivant laquelle
la variation du taux d'hémoglobine ne dépend ni de la durée de l'opération ni
du volume de sang perdu ou encore l'hypothèse suivant laquelle la variation du
taux d'hémoglobine ne dépend pas de la durée de l'opération.
\end{exo}

\begin{exo}[Analyse de la variance à un facteur]
  L'analyse de la variance est la comparaison de $k$ échantillons gaussiens
  indépendants de même variance. On se ramène à la formule de décomposition de
  la variance de la regression linéaire, cf. \cite[Chap. 15.6]{saporta} et
  \cite[Chap. 5.1.7]{dacunha-castelle-duflo} par exemple.
\end{exo}

{\tiny
\bibliographystyle{smfplain}
\bibliography{biblio}
}

\end{document}