Soient \(m,n\geq 1\). Une matrice \(\boldsymbol{A}\) de taille \((m,n)\) ou \(m\times n\) (à coefficients réels) est une application de \(\{1,\dots,m\} \times \{1,\dots,n\}\) dans \(\mathbb{R}\).
Plus simplement, il s’agit d’un tableau de nombres rééls ayant \(m\) lignes et \(n\) colonnes. On note \(A_{ij}\) l’élément sur la ligne \(i\) et sur la colonne \(j\) de \(\boldsymbol{A}\).
## [,1] [,2] [,3] [,4]
## [1,] 7 6 5 4
## [2,] 8 3 2 1
A
est une matrice \(2\times4\) et \(A_{1,3}\) (par exemple) vaut 5 - correspondant à A[1,3]
.
## [1] 2 4
## [1] 5
## [1] 7 6 5 4
## [1] 5 2
Soit \(\boldsymbol{A}\) une matrice réelle de taille \((m,n)\). La matrice transposée notée \(\boldsymbol{A}^\top\) de taille \((n,m)\) est définie par \((\boldsymbol{A}^\top)_{ij}=A_{ji}\) pour \(i=1,\dots,n\) et \(j=1,\dots,m\).
(bien sûr \((\boldsymbol{A}^\top)^\top=\boldsymbol{A}\))
## [,1] [,2]
## [1,] 7 8
## [2,] 6 3
## [3,] 5 2
## [4,] 4 1
Une matrice carrée \(\boldsymbol{A}\) de taille \((n,n)\) est dite symétrique si \(\boldsymbol{A}=\boldsymbol{A}^\top\).
Si \(\boldsymbol{x}\) et \(\boldsymbol{y}\) sont deux vecteurs de dimension \(n\) (deux matrices de taille \((n,1)\)), on note \(\boldsymbol{x}^\top y\) le produit scalaire entre \(\boldsymbol{x}\) et \(\boldsymbol{y}\) défini par \(\boldsymbol{x}^\top \boldsymbol{y}=\sum_{i=1}^n x_i y_i\).
## [,1]
## [1,] 88
On note \(\| \boldsymbol{x}\|\) la norme euclidienne du vecteur \(\boldsymbol{x}\), donnée par \(\|\boldsymbol{x}\|^2 = \boldsymbol{x}^\top \boldsymbol{x} = \sum x_i^2\).
## [1] 7 6 5 4
## [1] 126
Pour une matrice carrée de taille \(n\), on appelle diagonale de \(\boldsymbol{A}\), notée \(\operatorname{diag}(\boldsymbol{A})\) la matrice de taille \((n,n)\) définie par \((\operatorname{diag}(\boldsymbol{A}))_{ii}=A_{ii}\) pour \(i=1,\dots,n\) et \(0\) sinon.
Soient \(\boldsymbol{A}\) et \(\boldsymbol{B}\) deux matrices réelles de taille \((m,n)\) et \(\alpha,\beta\in R\) alors la matric \(\boldsymbol{C}= \alpha \boldsymbol{A}+\beta\boldsymbol{B}\) est une matrice réelle de taille \((m,n)\) et donnée par \(C_{ij} = \alpha A_{ij} +\beta B_{ij}\), \(i=1,\dots,m\), \(j=1,\dots,n\).
Soient \(\boldsymbol{A}\) et \(\boldsymbol{B}\) deux matrices réelles de taille \((m,n)\) et \((n,p)\) respectivement alors la matrice \(\boldsymbol{C}\) de taille \((m,p)\) définie par le produit matriciel entre \(\boldsymbol{A}\) et \(\boldsymbol{B}\) est bien définie et on a \(\boldsymbol{C}=\boldsymbol{A} \boldsymbol{B}\), \(C_{ij} = \sum_{k=1}^n A_{ik} B_{kj}\) pour \(i=1,\dots,m\), \(j=1,\dots,p\).
Le produit matriciel n’est pas commutatif pour deux matrices quelconque de même taille: \(\boldsymbol{A} \boldsymbol{B} \neq \boldsymbol{B} \boldsymbol{A}\).
Soit \(\mathbb{I}_n\) la matrice de taille \((n,n)\) composée de 1 sur la diagonale et de 0 ailleurs. Alors, pour \(\boldsymbol{A}\) de taille \((n,n)\), \(\mathbb{I}_n\) est l’élément neutre tel que \(\boldsymbol{A} \mathbb{I}_n = \mathbb{I}_n \boldsymbol{A} = \boldsymbol{A}\).
Soient \(\boldsymbol{A}\), \(\boldsymbol{B}\) et \(\boldsymbol{C}\) trois matrices réelles de dimension concordante, alors
## [,1] [,2]
## [1,] 1 3
## [2,] 2 4
## [,1] [,2]
## [1,] 1 3
## [2,] -2 1
## [,1] [,2]
## [1,] 3 -1
## [2,] -1 2
## [,1] [,2]
## [1,] -5 6
## [2,] -6 10
## [,1] [,2]
## [1,] -5 11
## [2,] -4 16
## [,1] [,2]
## [1,] -5 11
## [2,] -4 16
Une matrice symétrique réelle \(\boldsymbol{A}\) de taille \((n,n)\) est dite définie positive (resp. semi-définie positive) si \(\forall \boldsymbol{x} \in \mathbb{R}^n\), \(\boldsymbol{x}^\top \boldsymbol{A} \boldsymbol{x} >0\) (resp. \(\geq 0\))
Toute matrice \(\boldsymbol{A}\) qui peut s’écrire sous la forme \(\boldsymbol{A} =\boldsymbol{B}^\top \boldsymbol{B}\) est semi-définie positive positive
Soit \(\boldsymbol{A}\) une matrice carrée de taille \((n,n)\). La trace d’une matrice est définie par \[ \operatorname{trace} (\boldsymbol{A}) = \sum_{i=1}^n A_{ii} \] On a
Le déterminant d’une matrice carrée \(\boldsymbol{A}\) de taille \((n,n)\) est noté \(\mathrm{det}(\boldsymbol{A})\) ou \(|\boldsymbol{A}|\) est le réel dont la valeur absolue mesure le volume du parallélépipède engendré par les colonnes de \(\boldsymbol{A}\). La formule générale est \[ |\boldsymbol{A}| = \sum_{\sigma \in S_n} \mathrm{sgn}(\sigma) \prod_{i=1}^n A_{i\sigma_i} \] où \(S_n\) est l’ensemble des permutations de \(\{1,\dots,n\}\) et où \(\mathrm{sgn}(\sigma)\) désigne la signature de \(\sigma \in \{-1,1\}\) (la signature d’une permutation vaut 1 si le nombre de transpositions est pair et vaut -1 sinon).
## [1] 17
## [1] 17
## [1] 17
Soit \(\boldsymbol{A}^{ij}\) la matrice de taille \((n-1,n-1)\) correspondant à la matrice \(\boldsymbol{A}\) à laquelle ont été supprimées les \(i\)ème ligne et \(j\)ème colonne, alors \(|\boldsymbol{A}| = \sum_{i=1}^n (-1)^{i+j} A_{ij} |\boldsymbol{A}^{ij}|\)
## [1] 0
## [1] 17
## [1] 0
## [,1] [,2] [,3]
## [1,] 1 4 0
## [2,] 2 3 8
## [3,] 0 -1 5
## [1] -17
Théorème de Sylvester: soient \(\boldsymbol{A}\) et \(\boldsymbol{B}\) deux matrices de taille \((m,n)\) et \((n,m)\) respectivement alors \[ \det( \mathbb{I}_m + \boldsymbol{A} \boldsymbol{B}) = \det( \mathbb{I}_n + \boldsymbol{B} \boldsymbol{A}). \]
Soit \(\boldsymbol{A}\) une matrice carrée de taille \((n,n)\) dont le déterminant est non nul, alors \(\boldsymbol{A}\) est dite non singulière et il existe une matrice inverse (de même taille) notée \(\boldsymbol{A}^{-1}\) vérifiant \(\boldsymbol{A} \boldsymbol{A}^{-1} = \boldsymbol{A}^{-1}\boldsymbol{A} = \mathbb{I}_n\). Ses coefficients sont donnés par \[ (\boldsymbol{A}^{-1})_{ij} = (-1)^{i+j} \frac{|\boldsymbol{A}^{ij}|}{|\boldsymbol{A}|} \]
En dimension 2, on a la formule simple \[ \left( \begin{array}{cc} A_{11} & A_{12} \\ A_{21} &A_{22}\end{array} \right)^{-1} = (A_{11}A_{22}-A_{12}A_{21})^{-1} \left(\begin{array}{cc} A_{22} & -A_{12} \\ -A_{21} & A_{11} \end{array} \right) \]
## [1] 36
## [,1] [,2] [,3]
## [1,] 1 -1.040834e-17 1.387779e-17
## [2,] 0 1.000000e+00 0.000000e+00
## [3,] 0 -5.551115e-17 1.000000e+00
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
Soient \(\boldsymbol{A}\) et \(\boldsymbol{B}\) deux matrices inversibles de taille \((n,n)\) alors
## [1] 0.02777778
## [1] 0.02777778
Soient \(n_1\) et \(n_2\) tels que \(n=n_1+n_2\). Une matrice \((n,n)\) peut être écrite de la façon suivante \[ \boldsymbol{A} =\left(\begin{array}{cc} \boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{A}_{21} &\boldsymbol{A}_{22}\end{array}\right) \] où \(\boldsymbol{A}_{11}\), \(\boldsymbol{A}_{12}\), \(\boldsymbol{A}_{21}\) et \(\boldsymbol{A}_{22}\) sont de taille \((n_1,n_1)\), \((n_1,n_2)\), \((n_2,n_1)\) et \((n_2,n_2)\).
Les vecteurs \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p\) de \(\mathbb{R}^n\) sont dit linéairement indépendants si \(\forall \mathbf a \in \mathbb{R}^p\), \(\sum_{i=1}^p a_i \boldsymbol{x}_i=0 \Leftrightarrow a_1=\dots=a_p=0\). Si \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p\) sont linéairement indépendants alors \(p\leq n\).
On appelle rang d’une matrice \(\boldsymbol{A}\) et on note \(\operatorname{rang}(\boldsymbol{A})\) le nombre maximal de colonnes qui sont linéairement indépendantes.
Une matrice \(\boldsymbol{A}\) de taille \((n,p)\) est dite de plein rang si \(\operatorname{rang}(\boldsymbol{A})=\min(n,p)\).
Une matrice \(\boldsymbol{A}\) de taille \((p,p)\) est dite inversible ssi \(\operatorname{rang}(\boldsymbol{A})=p\).
Si \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p\) sont des vecteurs linéairement indépendants de \(\mathbb{R}^n\), ils sont la base de l’espace vectoriel \[ \mathcal V(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p) = \{ \boldsymbol{y} \in \mathbb{R}^n: \boldsymbol{y} = \sum_{i=1}^p a_i \boldsymbol{x}_i, \mathbf a\in \mathbb{R}^p\}. \]
Les valeurs propres d’une matrice carrée de taille \((p,p)\), \(\boldsymbol{A}\), sont les solutions de l’équation \(\chi_\boldsymbol{A}(\lambda)=|\boldsymbol{A}-\lambda \mathbb{I}_p|=0\). En fait, \(\chi_\boldsymbol{A}(\lambda)\) s’appelle le polynôme caractéristique de degré \(p\) en \(\lambda\); les racines de ce polynôme peuvent donc être complexes. De plus, certaines valeurs propres peuvent avoir une multiplicité supérieure à 1.
A chaque valeur propre, on peut associer un vecteur propre \(\boldsymbol{u}_i\) vérifiant \(\boldsymbol{A} \boldsymbol{u}_i = \lambda_i \boldsymbol{u}_i\). Le vecteur propre n’est pas défini de manière unique, puisque \(c \boldsymbol{u}_i\) est aussi un vecteur propre de \(\lambda_i\) pour tout réel \(c\) non nul.
Une matrice réelle carrée \(\boldsymbol{A}\) est dite diagonalisable s’il existe une matrice inversible \(\boldsymbol{P}\) et une matrice diagonale \(\boldsymbol{D}\) à coefficients réels satisfaisant la relation~: \(\boldsymbol{A} = \boldsymbol{P} \boldsymbol{D} \boldsymbol{P}^{-1}\). Et dans ce cas chaque vecteur colonne de \(\boldsymbol{P}\) est un vecteur propre de \(\boldsymbol{A}\).
Si une matrice carrée \(\boldsymbol{A}\) diagonalisable admet une valeur propre nulle, alors nécessairement \(\det(\boldsymbol{A})=0\).
Soit \(\boldsymbol{A}\) une matrice carrée de taille \((p,p)\) diagonalisable alors, \(\operatorname{rang}(\boldsymbol{A})= p - m_0\) où \(m_0\leq p\) est la multiplicité de la valeur propre 0 (éventuellement nulle).
Une matrice carrée de taille \((p,p)\), \(\boldsymbol{P}\) est orthogonale si et seulement si : \(\boldsymbol{P}\) est à coefficients réels, est inversible et son inverse est égale à sa transposée : \(\boldsymbol{P}^{-1}= \boldsymbol{P}^\top\), et donc \(\boldsymbol{P} \boldsymbol{P}^\top =\mathbb{I}_n\).
Les vecteurs colonnes d’une matrice orthogonale sont orthogonaux.
Une matrice orthogonale est une transformation rigide; elle préserve les longueurs et les angles: \(\|\boldsymbol{P} \boldsymbol{x}\|= \|\boldsymbol{x}\|\) et \((\boldsymbol{P}\boldsymbol{x})^\top (\boldsymbol{P} \boldsymbol{y})=\boldsymbol{x}^\top \boldsymbol{y}\).
Soit \(\boldsymbol{A}\) une matrice symétrique réelle de taille \((p,p)\), alors il existe une matrice orthogonale \(\boldsymbol{P}\) (c’est-à-dire \(\boldsymbol{P}^{-1}=\boldsymbol{P}^\top\)) et une matrice \(\boldsymbol{D}\) diagonale \(\boldsymbol{D}=\operatorname{diag}(\lambda_i, i=1,\dots,p)\), telle que \[ \boldsymbol{A} = \boldsymbol{P} \boldsymbol{D} \boldsymbol{P}^\top \] Et dans ce cas, les éléments de \(\boldsymbol{D}\) sont les valeurs propres de \(\boldsymbol{A}\) et les vecteurs colonnes de \(\boldsymbol{P}\) les vecteurs propres associés.
Théorème d’encadrement des valeurs propres d’une matrice symétrique: Soit \(\boldsymbol{A}\) une matrice réelle symétrique de taille \((p,p)\) et soient \(\lambda_1\leq \lambda_2\leq \dots \leq \lambda_p\) ses valeurs propres ordonnées. Alors, pour tout \(\boldsymbol{x}\in \mathbb{R}^p\) \[ \lambda_1 \leq \frac{\boldsymbol{x}^\top \boldsymbol{A} \boldsymbol{x}}{\boldsymbol{x}^\top \boldsymbol{x}} \leq \lambda_p. \]
Une matrice \(\boldsymbol{A}\) réelle symétrique est définie positive (resp. semi-définie positive) si et seulement si toutes ses valeurs propres sont toutes positives (resp. positives ou nulles).
## [,1] [,2] [,3]
## [1,] -2 -3 -4
## [2,] -1 -8 -3
## [,1] [,2] [,3]
## [1,] 5 14 11
## [2,] 14 73 36
## [3,] 11 36 25
## eigen() decomposition
## $values
## [1] 9.566163e+01 7.338365e+00 -1.444161e-15
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.1900915 0.4585891 0.86807906
## [2,] -0.8624949 -0.5004045 0.07548514
## [3,] -0.4690073 0.7343646 -0.49065338
## [,1] [,2] [,3]
## [1,] 5 14 11
## [2,] 14 73 36
## [3,] 11 36 25
## [,1] [,2] [,3]
## [1,] -0.1900915 -0.86249486 -0.4690073
## [2,] 0.4585891 -0.50040445 0.7343646
## [3,] 0.8680791 0.07548514 -0.4906534
## [,1] [,2] [,3]
## [1,] -0.1900915 -0.86249486 -0.4690073
## [2,] 0.4585891 -0.50040445 0.7343646
## [3,] 0.8680791 0.07548514 -0.4906534
Soit \(\boldsymbol{y}\) un vecteur de \(\mathbb{R}^n\). Sa projection sur \(\mathcal V(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p)\) est le vecteur \(\hat{\boldsymbol{y}} = \boldsymbol{x}\hat{\mathbf a}\) minimisant \(\|\boldsymbol{y} -\boldsymbol{x} \mathbf a\|\). La solution de ce problème donne \(\hat{\mathbf a} = (\boldsymbol{x}^\top \boldsymbol{x})^{-1} \boldsymbol{x}^\top \boldsymbol{y}\) et le projeté est \(\hat{\boldsymbol{y}} = \boldsymbol{x} \hat{\mathbf a}\).
La matrice \(\mathcal{P} = \boldsymbol{x} ( \boldsymbol{x}^\top \boldsymbol{x})^{-1} \boldsymbol{x}\) est appelée projecteur ou matrice de projection sur \(\{\boldsymbol{x}_1,\dots,\boldsymbol{x}_p\}\) et on a \(\hat{ \boldsymbol{y}} = \mathcal{P} \boldsymbol{y}\).
Déterminons le projeté orthogonal de \(\boldsymbol{y}=(1,-1,0)^\top\) sur les vecteurs \(\boldsymbol{x}_1=(1,1,0)^\top\) et \(\boldsymbol{x}_2=(3,0,1)^\top\) :
f = function(a,y=c(1,-1,0)){ (y[1]-(a[1]+3*a[2]))^2+(y[2]-a[1])^2+(y[3]-a[2])^2 }
optim(par=c(0,0),f)$par
## [1] -0.8181057 0.5454202
## [,1]
## [1,] -0.8181818
## [2,] 0.5454545
Soient \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_p\), \(p\) vecteurs de \(\mathbb{R}^n\) et soit \(Q\subset\{1,\dots,p\}\) un sous-ensemble d’indices. On note \(\hat y_Q\) le projeté orthogonal de \(y\) sur l’espace \(\mathcal V_Q\) engendré par les vecteurs \(\boldsymbol{x}_j\) pour \(j\in Q\) et on note \(\boldsymbol{x}_Q=(\boldsymbol{x}_j, j\in Q)\) la matrice de taille \((n,\#Q)\). Alors on a les propriétés suivantes:
## [,1] [,2]
## [1,] 1 3
## [2,] 0 4
## [3,] 1 2
## [,1] [,2] [,3]
## [1,] 0.5151515 0.1212121 0.4848485
## [2,] 0.1212121 0.9696970 -0.1212121
## [3,] 0.4848485 -0.1212121 0.5151515
## [,1] [,2] [,3]
## [1,] 0.5151515 0.1212121 0.4848485
## [2,] 0.1212121 0.9696970 -0.1212121
## [3,] 0.4848485 -0.1212121 0.5151515
## [,1] [,2] [,3]
## [1,] 0.4848485 -0.12121212 -0.4848485
## [2,] -0.1212121 0.03030303 0.1212121
## [3,] -0.4848485 0.12121212 0.4848485
## [,1] [,2] [,3]
## [1,] -8.326673e-17 -1.387779e-16 5.551115e-17
## [2,] 9.714451e-17 2.081668e-16 2.775558e-17
## [3,] 0.000000e+00 -2.012279e-16 -5.551115e-17
Si \(\boldsymbol{A}\) est une matrice idempotente alors
Théorème de Cochran: soit \(\boldsymbol{A}=\boldsymbol{A}_1+\dots+\boldsymbol{A}_k\). Alors les deux énoncés suivants sont équivalents
Davis = read.table(
"http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)] = Davis[12,c(3,2)]
Un vecteur aléatoire \(\boldsymbol{X}=(X_1,\dots,X_d)\) de dimension \(d\) admet pour fonction de répartition \[ F(\boldsymbol{x})=F(x_1,\dots,x_d)=\mathbb{P}[X_1\leq x_1,\dots,X_d\leq x_d] \] Si les composantes sont absolument continues, la densité associée est \[ f(\boldsymbol{x})=f(x_1,\dots,x_d)=\frac{\partial^n F(x_1,\dots,x_d)}{\partial x_1\dots\partial x_d} \] Soit \(\boldsymbol{X}\) un vecteur aléatoire admettant une densité ou fonction de masse, de dimension \(d\), alors la densité marginale de \(X_i\) (\(i=1,\dots,d\)) est définie par \[ f_{X_i}(x_i ) = \int_{\mathbb{R}^{n-1}} f_{\boldsymbol{X}}(x_1,\dots,x_d) d x_1 \dots d x_{i-1} d x_{i+1} \dots d x_d \] dans le cas absolument continu.
Soit \(I,J \subset \{1,\dots,d\}\) tels que \(I\cap J =\emptyset\) et \(I\cup J=\{1,\dots,d\}\). Pour \(\boldsymbol{x}\in \mathbb{R}^d\) on note \(\boldsymbol{X}_I\) et \(\boldsymbol{x}_J\) les vecteurs associées aux composantes d’indices \(I\) et \(J\) respectivement.
La distribution conditionnelle de \(\boldsymbol{X}_I\) sachant \(\boldsymbol{X}_J=\boldsymbol{x}_J\) est donnée par \[ f_{\boldsymbol{X}_I}(\boldsymbol{x}_I) = \frac{f_\boldsymbol{X}(\boldsymbol{x})}{f_{\boldsymbol{X}_J}(\boldsymbol{x}_J)} \] pour \(\boldsymbol{x}_J\) tel que \(f_{\boldsymbol{X}_J}(\boldsymbol{x}_J)>0\), dans le cas absolument continu.
Soient \(\boldsymbol{X}\) et \(\boldsymbol{Y}\) deux vecteurs aléatoires de dimension \(d\) et \(d^\prime\).
La matrice de variance-covariance est nécessairement une matrice symétrique semi-définie positive.
Soit \(\boldsymbol{Y}\) un vecteur aléatoire de dimension \(d\), de moyenne \(\boldsymbol{\mu}\) et de matrice de covariance \(\boldsymbol{Sigma}\). Soient \(\boldsymbol{A}\) et \(\boldsymbol{B}\) deux matrices réeeles de taille \((d,p)\) et \((d,q)\) et enfin soit \(\mathbf a \in \mathbb{R}^p\) alors
Soit \(\boldsymbol{A}\) une matrice réelle symétrique de taille \((d,d)\) et \(\boldsymbol{Y}\) un vecteur aléatoire de moyenne \(\boldsymbol{\mu}\) et de matrice de covariance \(\boldsymbol{\Sigma}\), alors \[ \mathbb{E} \left( \boldsymbol{Y}^\top \boldsymbol{A} \boldsymbol{Y}\right) = \boldsymbol{\mu}^\top \boldsymbol{A} \boldsymbol{\mu} + \operatorname{trace}( \boldsymbol{A} \boldsymbol{\Sigma}). \]
\[{\displaystyle f(k)=\mathbb {P} (X=k)={n \choose k}\,p^{k}(1-p)^{n-k}} \] de moyenne \(\mathbb{E}[X]=np\) et de variance \(\text{Var}[X]=np(1-p)\)
par(mfrow=c(1,2))
p=dbinom(0:10,size = 10,prob = .3)
names(p)=0:10
barplot(p,col=rgb(0,0,1,.4))
u=seq(-1,11,length=251)
plot(u,pbinom(u,size = 10,prob = .3),type="l",xlab="",ylab="")
lines(u,pnorm(u,10*.3,sqrt(10*.3*(1-.3))),col="red",lty=2)
\[{\displaystyle f(k)= {\frac {\lambda ^{k}}{k!}}\mathrm {e} ^{-\lambda }} \] de moyenne \(\mathbb{E}[X]=\lambda\) et de variance \(\text{Var}[X]=\lambda\)
par(mfrow=c(1,2))
p=dpois(0:10,lambda = 3)
names(p)=0:10
barplot(p,col=rgb(0,0,1,.4))
u=seq(-1,11,length=251)
plot(u,ppois(u,lambda = 3),type="l",xlab="",ylab="")
lines(u,pnorm(u,3,sqrt(3)),col="red",lty=2)
La fonction de répartition est donnée par : \[ {\displaystyle F(x)=\left\{{\begin{matrix}1-e^{-\lambda x}&{\text{si}}\;x\geqslant 0\\0&{\text{si}}\;x<0\end{matrix}}\right.} \] et la densité \({\displaystyle \lambda e^{-\lambda x}\boldsymbol{1}_{\mathbb {R} _{+}}(x)}\). La moyenne est \({\displaystyle \mathbb {E} (X)={\frac {1}{\lambda }}}\) et la variance \({\displaystyle \text{Var}(X)={\dfrac {1}{\lambda ^{2}}}}\)
par(mfrow=c(1,2))
u=seq(-1,11,length=251)
plot(u,dexp(u,1/3),col="white",xlab="",ylab="")
polygon(c(u,rev(u)),c(dexp(u,1/3),rep(0,length(u))),col=rgb(0,0,1,.4), border=NA)
lines(u,dexp(u,1/3))
plot(u,pexp(u,1/3),type="l",xlab="",ylab="")
La densité de la loi normale (Gaussienne) \(\mathcal{N}(\mu,\sigma^2)\) est donnée par : \[ {\displaystyle f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}\exp\big[-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}}\big] \] La moyenne est \({\displaystyle \mathbb {E} (X)=\mu}\) et la variance \({\displaystyle \text{Var}(X)=\sigma^2}\)
par(mfrow=c(1,2))
u=seq(-1,11,length=251)
plot(u,dnorm(u,3,1.5),col="white",xlab="",ylab="")
polygon(c(u,rev(u)),c(dnorm(u,3,1.5),rep(0,length(u))),col=rgb(0,0,1,.4), border=NA)
lines(u,dnorm(u,3,1.5))
plot(u,pnorm(u,3,1.5),type="l",xlab="",ylab="")
Si \({\displaystyle X\sim {\mathcal {N}}(\mu ,\sigma ^{2})}\) alors \({\displaystyle aX+b\sim {\mathcal {N}}(a\mu +b,a^{2}\sigma ^{2})}\).
##
## Attaching package: 'MASS'
## The following object is masked from 'package:sm':
##
## muscle
## mean sd
## 170.5650000 8.9098695
## ( 0.6300229) ( 0.4454935)
hist(Davis$height,probability = TRUE, col=rgb(0,0,1,.4), main="")
x=seq(min(Davis$height)-2,max(Davis$height)+2,length=251)
y=dnorm(x,fitdistr(Davis$height,"normal")$estimate[1],fitdistr(Davis$height,"normal")$estimate[2])
lines(x,y,col="red",lwd=2)
Soient \(\{X_1, \cdots , X_k\}\), \(k\) variables aléatoires indépendantes suivant des lois normales centrées (de moyennes 0) et réduites (de variance 1), alors \[ {\displaystyle X=\sum _{i=1}^{k}X_{i}^{2}}\sim\chi^2(k) \] suit une loi du chi-deux à \(k\) degrés de liberté. La densité de \(X\) est \[ {\displaystyle f(x)={\frac {1}{2^{\frac {k}{2}}\Gamma ({\frac {k}{2}})}}x^{{\frac {k}{2}}-1}e^{-{\frac {x}{2}}}}\boldsymbol{1}_{\mathbb{R}_+}(x) \] On peut noter que \(\mathbb{E}[X]=k\) et \(\operatorname{var}[X]=2k\).
par(mfrow=c(1,2))
u=seq(-1,11,length=251)
plot(u,dchisq(u,3),col="white",xlab="",ylab="")
polygon(c(u,rev(u)),c(dchisq(u,3),rep(0,length(u))),col=rgb(0,0,1,.4), border=NA)
lines(u,dchisq(u,3))
plot(u,pchisq(u,3),type="l",xlab="",ylab="")
Si \(X\) est une variable aléatoire de loi normale centrée et réduite et \(Y\) suit une loi \(\chi^2(n)\), telles que \(X\) et \(Y\) soient indépendantes, alors \({\displaystyle {\frac {X}{\sqrt {Y/n}}}}\) suit une loi de Student à \(n\) degrés de liberté, \(\mathcal{S}td(n)\)
par(mfrow=c(1,2))
u=seq(-4,4,length=251)
plot(u,dt(u,3),col="white",xlab="",ylab="",ylim=c(0,.4))
polygon(c(u,rev(u)),c(dt(u,3),rep(0,length(u))),col=rgb(0,0,1,.4), border=NA)
lines(u,dt(u,3))
lines(u,dnorm(u),col="red",lty=2)
plot(u,pt(u,3),type="l",xlab="",ylab="")
lines(u,pnorm(u),col="red",lty=2)
La densité de \(T\) est donnée par : \[ {\displaystyle f(t)={\frac {1}{\sqrt {k\pi }}}{\frac {\Gamma ({\frac {k+1}{2}})}{\Gamma ({\frac {k}{2}})}}\left(1+{\frac {t^{2}}{k}}\right)^{-{\frac {k+1}{2}}}} \] qui est définie pour \(n\in\mathbb{R}_{+\star}\) et pas juste \(n\in\mathbb{N}_{\star}\). Si \(n>1\), \(\mathbb{E}[T]=1\), et si \(n>2\), \(\operatorname{var}[T]=n/(n-2)\).
Le quotient de deux variables aléatoires indépendantes, \(Q_1\) et \(Q_2\), distribuées chacune selon des lois \(\chi^2(d_1)\) et \(\chi^2(d_2)\) respectivement suit une loi de Fisher \[ {\displaystyle \mathrm {F} = {\frac {U_{1}/d_{1}}{U_{2}/d_{2}}}}\sim\mathcal{F}(d_1,d_2) \] La densité d’une loi de Fisher est donnée par \[ {\displaystyle f(x)={\frac {1}{x\;\mathrm {B} (d_{1}/2,d_{2}/2)}}\left({\frac {d_{1}\,x}{d_{1}\,x+d_{2}}}\right)^{d_{1}/2}\;\left(1-{\frac {d_{1}\,x}{d_{1}\,x+d_{2}}}\right)^{d_{2}/2}} \] pour \(x>0\). L’espérance de cette variable aléatoire est \[ \mathbb{E}[\mathrm{F}] =\frac{d_2}{d_{2}-2} \] dès lors que \(d_2>2\).
La densité d’un vecteur Gaussien \(\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})\) est \[ {\displaystyle f(x_{1},\ldots ,x_{k})={\frac {\exp \left(-{\frac {1}{2}}({\mathbf {x} }-{\boldsymbol {\mu }})^{\top }{\boldsymbol {\Sigma }}^{-1}({\mathbf {x} }-{\boldsymbol {\mu }})\right)}{\sqrt {(2\pi )^{k}|{\boldsymbol {\Sigma }}|}}}} \]
plot(Davis$height,Davis$weight,xlab="Height (cm)",ylab="Weight (kg)")
library(mnormt)
f=function(x,y) dmnorm(cbind(x,y),c(mean(Davis$height),mean(Davis$weight)),var(Davis[,c("height","weight")]))
x=seq(min(Davis$height),max(Davis$height),length=101)
y=seq(min(Davis$weight),max(Davis$weight),length=101)
z=outer(x,y,f)
contour(x,y,z,add=TRUE,col="red")
Let \[ {\displaystyle \mathbf {X} ={\begin{bmatrix}\mathbf {x} _{1}\\\mathbf {x} _{2}\end{bmatrix}}{\text{ with sizes }}{\begin{bmatrix}q\times 1\\(N-q)\times 1\end{bmatrix}}}\]
and accordingly μ and Σ are partitioned as follows \[ {\boldsymbol {\mu }}={\begin{bmatrix}{\boldsymbol {\mu }}_{1}\\{\boldsymbol {\mu }}_{2}\end{bmatrix}}{\text{ with sizes }}{\begin{bmatrix}q\times 1\\(N-q)\times 1\end{bmatrix}}\] and \[{\displaystyle {\boldsymbol {\Sigma }}={\begin{bmatrix}{\boldsymbol {\Sigma }}_{11}&{\boldsymbol {\Sigma }}_{12}\\{\boldsymbol {\Sigma }}_{21}&{\boldsymbol {\Sigma }}_{22}\end{bmatrix}}{\text{ with sizes }}{\begin{bmatrix}q\times q&q\times (N-q)\\(N-q)\times q&(N-q)\times (N-q)\end{bmatrix}}} \]
\[{\displaystyle {\bar {\boldsymbol {\mu }}_{\mathbf {a}}}={\boldsymbol {\mu }}_{1}+{\boldsymbol {\Sigma }}_{12}{\boldsymbol {\Sigma }}_{22}^{-1}\left(\mathbf {a} -{\boldsymbol {\mu }}_{2}\right)}\]
and covariance matrix \[{\displaystyle {\overline {\boldsymbol {\Sigma }}_{\mathbf {a}}={\boldsymbol {\Sigma }}_{11}-{\boldsymbol {\Sigma }}_{12}{\boldsymbol {\Sigma }}_{22}^{-1}{\boldsymbol {\Sigma }}_{21}.}}\] In the bivariate case\[ {\displaystyle {\begin{pmatrix}X_{1}\\X_{2}\end{pmatrix}}\sim {\mathcal {N}}\left({\begin{pmatrix}\mu _{1}\\\mu _{2}\end{pmatrix}},{\begin{pmatrix}\sigma _{1}^{2}&\rho \sigma _{1}\sigma _{2}\\\rho \sigma _{1}\sigma _{2}&\sigma _{2}^{2}\end{pmatrix}}\right)}\] xx \[{\displaystyle X_{1}\mid X_{2}=x_{2}\ \sim \ {\mathcal {N}}\left(\mu _{1}+{\frac {\sigma _{1}}{\sigma _{2}}}\rho (x_{2}-\mu _{2}),\,(1-\rho ^{2})\sigma _{1}^{2}\right).} \] so that \[{\displaystyle \mathbb {E} (X_{1}\mid X_{2}=x_{2})=\mu _{1}+\rho {\frac {\sigma _{1}}{\sigma _{2}}}(x_{2}-\mu _{2})} \] while \[ {\displaystyle \operatorname {var} (X_{1}\mid X_{2}=x_{2})=(1-\rho ^{2})\sigma_1^2<\sigma_1^2= \operatorname {var} (X_{1})} \]
Supposons qu’on veuille simuler \(\boldsymbol{X} \sim \mathcal{N}(\boldsymbol{\mu} ,\boldsymbol{\Sigma})\), un vecteur gaussien de taille \(d\).
\(\boldsymbol{X} \stackrel{\mathcal{L}}{=} \boldsymbol{X}^\prime +\boldsymbol{\mu}\), avec \(\boldsymbol{X}^\prime \sim \mathcal{N}(\boldsymbol{0} ,\boldsymbol{\Sigma})\): il suffit donc de simuler \(\boldsymbol{X}^\prime\).
Puisque \(\boldsymbol{\Sigma}\) est réelle symétrique, \(\exists \boldsymbol{P}\) orthogonale et \(\boldsymbol{D}\) diagonale (dont les éléments sont \(\geq 0\)), telle que \(\boldsymbol{\Sigma} = \boldsymbol{P} \boldsymbol{D} \boldsymbol{P}^\top\), alors \[ \boldsymbol{P} \boldsymbol{D}^{1/2} \boldsymbol{Y} \stackrel{\mathcal{L}}{=} \boldsymbol{X}^\prime \] où \(\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I}_d)\), puisque \(\operatorname{Var}(\boldsymbol{P} \boldsymbol{D}^{1/2} \boldsymbol{Y}) = \boldsymbol{P} \boldsymbol{D}^{1/2} \mathbb{E}(\boldsymbol{Y} \boldsymbol{Y}^\top) \boldsymbol{D}^{1/2} \boldsymbol{P}^\top = \boldsymbol{\Sigma}\)
Aussi, pour simuler une réalisation de \(\boldsymbol{X} \sim \mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})\):
y=rnorm(d)
A = function(rho){
Sigma=cbind(c(1,rho),c(rho,1))
eig=eigen(Sigma)
eig$vector %*% diag(sqrt(eig$val))
}
A0.5 = A(0.5)
Y = rnorm(2,0,1)
A0.5%*%Y
## [,1]
## [1,] 1.417492
## [2,] 1.290876
## [,1] [,2] [,3]
## [1,] -0.8129428 0.6786791 2.119525
## [2,] -1.4852915 1.5903142 1.664039
## [1] 1.008306 1.032374
## [1] 0.01081344
par(mfrow = c(1,2))
library(plot3D)
xc=cut(X0[1,],20);yc=cut(X0[2,],20);z=table(xc,yc)
hist3D(z=z, border="black")
xc=cut(X0[1,],20);yc=cut(X0[2,],20);z=table(xc,yc)
image2D(z=z, border="black")
X0.5=A(0.5)%*%y
par(mfrow = c(1,2))
library(plot3D)
xc=cut(X0.5[1,],20);yc=cut(X0.5[2,],20);z=table(xc,yc)
hist3D(z=z, border="black")
xc=cut(X0.5[1,],20);yc=cut(X0.5[2,],20);z=table(xc,yc)
image2D(z=z, border="black")
X0.9=A(0.9)%*%y
par(mfrow = c(1,2))
library(plot3D)
xc=cut(X0.9[1,],20);yc=cut(X0.9[2,],20);z=table(xc,yc)
hist3D(z=z, border="black")
xc=cut(X0.9[1,],20);yc=cut(X0.9[2,],20);z=table(xc,yc)
image2D(z=z, border="black")
Soit \(\boldsymbol{Y}\sim \mathcal{N}(\boldsymbol{0}, \mathbb{I}_d)\) et soit \(\boldsymbol{A}\) une matrice réelle symétrique de rang \(r\). Alors,
Soit \(\boldsymbol{Y} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbb{I}_d)\), et définissons \(\boldsymbol{X} = \boldsymbol{Y}^\top \boldsymbol{Y} = \|\boldsymbol{Y}\|^2\). On dit que \(\boldsymbol{X}\) suit une loi du \(\chi^2(d)\) et de paramètre de non centralité \(\lambda = \boldsymbol{\mu}^\top \boldsymbol{\mu}\), et on note \(\boldsymbol{X} \sim \chi^2(d,\boldsymbol{\mu}^\top \boldsymbol{\mu})\).
Soit \(\boldsymbol{Y}\sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) et soit \(\boldsymbol{A}\) une matrice réelle symétrique de rang \(r\). Alors,
Théorème de Cochrane: Soit \(\boldsymbol{X}\sim \mathcal{N}(\boldsymbol{\mu},\sigma^2 \mathbb{I}_d)\) où \(\boldsymbol{\mu}\in \mathbb{R}^d\) et \(\sigma>0\). Soient \(F_1,\dots,F_m\) des sous-espaces vectoriels de \(\mathbb{R}^d\) deux à deux orthogonaux, et \(\boldsymbol{P}_{F_i}\) la matrice de projection orthogonale sur \(F_i\) et \(d_i\) la dimension de \(F_i\) pour \(1\leq i \leq m\), alors
Autrement, si \(\boldsymbol{X}\sim \mathcal{N}(\boldsymbol{ 0},\mathbb{I}_d)\) et \(F\) un sous-espace vectoriel de \(\mathbb{R}^d\) de dimension \(p< d\), si \(F^\bot\) désigne son orthogonal, et si \(\boldsymbol{P}_F\) et \(\boldsymbol{P}_{F^\bot}\) les matrices de projection sur \(F\) et \(F^\bot\), alors
Une conséquence est que si \(Y_1,\dots,Y_n\) \(n\) sont des variables i.i.d. de loi \(N(0,1)\) alors \(\overline Y=n^{-1} \sum Y_i\) et \(S^2=(n-1)^{-1}\sum_i (Y_i-\overline Y)\) sont deux variables indépendantes et de plus, \((n-1)S^2\sim\chi^2_{n-1}\).
On dira que la loi de \(Y\) est dans la famille exponentielle si sa masse de probabilité / densité s’écrit \[ f(y|\theta,\varphi) = \exp\left(\frac{y\theta-b(\theta)}{ \varphi }+c(y,\varphi)\right) \] où \(a(\cdot)\), \(b(\cdot)\) et \(c(\cdot)\) sont des fonctions, et où \(\theta\) est appelé . On peut aussi écrire \[ f(y|\theta,\varphi) = A(y)B(\theta)\cdot\exp\left(\frac{y\theta}{ \varphi }\right) \] Plus généralement, on peut considérer \[ f(y|\theta,\varphi) = A(y)B(\theta)\cdot\exp\left(\frac{1}{ \varphi }\sum_{j=1}^r T_j(y)\eta_j(\theta)\right) \] (par la suite, la version simple précédante suffira). \(\boldsymbol{\eta}=(\eta_1(\theta),\dots,\eta_r(\theta))\) est le paramètre canonique.
A sequence \(\{X1, \dots, Xn, \cdots\}\) of random variables is said to converge in distribution to a random variable \(X\) if \[ {\displaystyle \lim _{n\to \infty }F_{n}(x)=F(x)} \] for every $ {x }$ at which \(F\) is continuous.
We will write \(X_n{\xrightarrow {\mathcal{L}}}X\).
Let \(\{X1, \cdots, Xn, \cdots\}\) denote a sample of size \(n\), i.e. independent and identically distributed (i.i.d.) random variables drawn from a distribution of expected value \(\mu\) and finite variance \(\sigma^2\). Let \[ \displaystyle{S_{n}={\frac {X_{1}+\cdots +X_{n}}{n}}} \] as \(n\) goes to infinity, the random variables \(\sqrt{n}(Sn − \mu)\) converge in distribution to a normal \(\mathcal{N}(0,\sigma^2)\) \[ {\displaystyle {\sqrt {n}}\left(S_{n}-\mu \right)\ {\xrightarrow {\mathcal{L}}}\mathcal{N}\left(0,\sigma ^{2}\right).} \]
m=100
mean_samples=function(n=10){
X=matrix(rnorm(n*m),nrow=m,ncol=n)
return(apply(X,1,mean))
}
B=matrix(NA,100,20)
for(i in 1:20){
B[,i] = mean_samples(i*10)
}
colnames(B) = as.character(seq(10,200,by=10))
boxplot(B)
u = seq(0,21,by=.2)
v = sqrt(u*10)
lines(u,1.96/v,col="red")
lines(u,-1.96/v,col="red")
For a univariate parameter \(\theta\) such that \[{\displaystyle {{\sqrt {n}}[X_{n}-\theta ]\,{\xrightarrow {\mathcal{L}}}\,{\mathcal {N}}(0,\sigma ^{2})}}\] then \[ {\displaystyle {{\sqrt {n}}[g(X_{n})-g(\theta )]\,{\xrightarrow {\mathcal L}}\,{\mathcal {N}}(0,\sigma ^{2}\cdot [g'(\theta )]^{2})}} \] for any function \(g\) such that \(g'(\theta)\) exists and is non-zero.
For a multivariate parameter \(\boldsymbol\theta\) such that \[{\displaystyle {{\sqrt {n}}[\boldsymbol{X}_{n}-\boldsymbol{\theta} ]\,{\xrightarrow {\mathcal{L}}}\,{\mathcal {N}}(\boldsymbol0,\boldsymbol\Sigma )}}\] then \[ {\displaystyle {{\sqrt {n}}[g(\boldsymbol{X}_{n})-g(\boldsymbol\theta )]\,{\xrightarrow {\mathcal{L}}}\,{\mathcal {N}}(0,\nabla g(\boldsymbol\theta )^\top\cdot\boldsymbol\Sigma\cdot \nabla g(\boldsymbol\theta ))}} \] for any function \(g\) such that $ g())$ exists and is invertible.
On appelle échantillon le vecteur \(\boldsymbol{x}= (x_1,\dots,x_n)\) vu comme la réalisation de \(n\) variables aléatoires identiquement distribuées, \(\boldsymbol{X}= (X_1,\dots,X_n)\). On supposera que les \(X_i\) suivent une loi \(F_{\theta}\in\mathcal{F}\) où le paramètre \(\theta\) est inconnu, appartenant à un ensemble \(\Theta\).
On parle de modèle identifiable si \(F_{\theta}=F_{\theta'} \Longrightarrow\theta=\theta'\).
Une statistique est une variable aléatoire \(T\) dépendant de l’échantillon \(\boldsymbol{X}\) ne faisant pas intervenir le paramètre \(\theta\). On notera \(T=T(\boldsymbol{X})\).
On note \(\mathcal{L}(\theta; \boldsymbol{x})\) la distribution jointe de \((\boldsymbol{X}_1,\dots,\boldsymbol{X}_n)\), appelée vraisemblance.
On dit que \(T\) est exhaustive pour \(\theta\), si la loi conditionnelle de \(\boldsymbol{X}\) sachant \(T=t\) ne dépend pas de \(\theta\) pour tout \(t\). Par le théorème de factorisation, \(T\) est exhaustive s’il existe deux applications mesurables \(g\) et \(h\) telles que \[ \mathcal{L}(\theta;\boldsymbol{x}) = g(\boldsymbol{x}) \cdot h(T(\boldsymbol{x});\theta) \]
Soit \(\boldsymbol{X}\) un échantillon de \(n\) variables aléatoires i.i.d. de loi de Bernoulli \(\mathcal B(p)\), alors la statistique \(T(\boldsymbol{X})=\sum_j X_j\) est exhaustive pour \(p\)
Soit \(\boldsymbol{X}\) un échantillon de \(n\) variables aléatoires i.i.d. de loi normale \(\mathcal{N}(\mu,\sigma^2)\), alors \(T(\boldsymbol{X})=(\sum_j X_j,\sum_j X_j^2)\) est exhaustive pour le couple \((\mu,\sigma^2)\).
Soit une statistique \(T\). Si pour \(\boldsymbol{x}=(x_1,\dots,x_n)\) et \(\boldsymbol{y}=(y_1,\dots,y_n)\) on a l’équivalence \[ T(\boldsymbol{x}) = T(\boldsymbol{y}) \Leftrightarrow \theta \mapsto \frac{\mathcal{L}(\theta,\boldsymbol{x})}{\mathcal{L}(\theta,\boldsymbol{y})} \text{ ne dépend pas de } \theta. \] alors \(T\) est une statistique exhaustive minimale.
Une statistique exhaustive \(T\)t (à valeurs dans \(\mathbb{R}^d\)) d’un modèle paramétrique est dite complète si pour toute fonction mesurable \(g:\mathbb{R}^d\to\mathbb{R}\) telle que \(g(T)\) soit intégrable on ait \[ \forall \theta \in \Theta, \mathbb{E}_\theta (f(T))=0 \Rightarrow g(T)=0 \;\; P_\theta-p.s. \]
Soit un modèle paramétrique. Si \(T\) est une statistique exhaustive complète
Dans un modèle statistique paramétrique, il existe toujours une statistique exhaustive minimale mais pas toujours de statistique exhaustive complète.
Example d’un pile ou face
## [1] 0 1 0 1 1 0 1 1 1 0 1 0 1 1 0 1 0 0 0 1
p=seq(0,1,by=.01)
logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
LL=sapply(p,logL)
plot(p,LL,type="l",col="red",lwd=2)
p0=.5
points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
abline(v=p0,lty=2)
p_hat = mean(X)
points(p_hat,logL(p_hat),pch=3,cex=1.5,lwd=2,col="blue")
abline(v=p_hat,lty=2,col="blue")
Recall that if \(X\sim F_\theta\) \[ {\displaystyle \mathbb{E}\left(\frac{\partial \log \mathcal{L}(X)}{\partial \theta}\right)}=0 \] while \[ {\displaystyle \operatorname{var}\left(\frac{\partial \log \mathcal{L}(X)}{\partial \theta}\right)}=\mathbb{E}\left(\frac{\partial \log \mathcal{L}(X)}{\partial \theta}\frac{\partial \log \mathcal{L}(X)}{\partial \theta^\top}\right)=I_n(\theta) \] Observe that \[ I_n(\theta)=n\cdot \operatorname{var}\left(\frac{\partial \log f_\theta(X)}{\partial \theta}\right)=\mathbb{E}\left(\frac{\partial \log f_\theta(X)}{\partial \theta}\frac{\partial \log f_\theta(X)}{\partial \theta^\top}\right)=n\cdot I(\theta) \] Given a vector \(\boldsymbol{x}=(x_1,\dots,x_n)\), the score is the vector \[ s_n(\boldsymbol{x},\theta)=\nabla\log\mathcal{L}(\theta;\boldsymbol{x}) \]
To compute derivatives, use
n = 77
S = rep(NA,1e5)
for(s in 1:length(S)){
Xs = sample(0:1,size=n,replace=TRUE)
logL = function(p){sum(log(dbinom(Xs,size=1,prob=p)))}
S[s] = slope_h(logL,p0)
}
mean(S)
## [1] -0.03784
hist(S,probability = TRUE)
u=seq(min(S),max(S),length=251)
v=dnorm(u,0,sd(S))
lines(u,v,lty=2,col="red")
n = 77
S = rep(NA,1e5)
for(s in 1:length(S)){
Xs = sample(0:1,size=n,replace=TRUE)
logL = function(p){sum(log(dbinom(Xs,size=1,prob=p)))}
dlogL = function(p){slope_h(logL,p)}
S[s] = -slope_h(dlogL,p0)
}
mean(S)
## [1] 308.0002
## [1] 308
## [1] 0 1 0 1 1 0 1 1 1 0 1 0 1 1 0 1 0 0 0 1
neglogL = function(p){-sum(log(dbinom(X,size=1,prob=p)))}
(pml = optim(fn=neglogL,par=p0,method="BFGS")$par)
## Warning in dbinom(X, size = 1, prob = p): NaNs produced
## Warning in dbinom(X, size = 1, prob = p): NaNs produced
## [1] 0.5499996
## [1] 0.55
On appelle estimateur (ponctuel) de \(\theta\) toute statistique \(T\) prenant ses valeurs dans \(\Theta\).
Un estimateur, \(T\), de \(\theta\) est dit sans biais, ou non biaisé, si \(\mathbb{E}_\theta(T)=\theta\), où \[ \mathbb{E}_{\color{red}{\theta}}(T)=\int_{\mathbb{R}^n}T(x_1,\dots,x_n)dF_{\color{red}{\theta}}(x_1)\dots dF_{\color{red}{\theta}}(x_n) \] On définit le biais par \(b(\theta)=\mathbb{E}_\theta(T)-\theta\)
Le risque quadratique d’un estimateur \(T\) de \(\theta\) est \(R(T,\theta) = \mathbb{E}_\theta[ (T-\theta)^2]\). On a \[ R(T,\theta) = b(\theta)^2 \, + \, \operatorname{Var}_\theta(T) \] (pour un estimateur sans biais, \(R(T,\theta)=\operatorname{Var}_\theta(T)\))
Soient \(T_1\) et \(T_2\) deux estimateurs de \(\theta\). On dira que \(T_1\) est plus efficace que \(T_2\) si \(R(T_1,\theta)\leq R(T_2,\theta)\).
On dit que la suite d’estimateurs \((T_n)_{n\geq 1}\) d’estimateurs de \(\theta\) est
convergente si \(T_n \stackrel{\mathbb{P}}{\to} \theta\) pour tout \(\theta\in \Theta\).
fortement convergente si \(T_n \stackrel{p.s.}{\to} \theta\) pour tout \(\theta\in \Theta\).
asymptotiquement normale si pour tout \(\theta\in \Theta\), il existe une matrice de covariance \(\boldsymbol{\Sigma}(\theta)\) telle que \(\sqrt{n}(T_n-\theta) \to \mathcal{N}(0,\boldsymbol{\Sigma}(\theta))\)
Sous certaines conditions de régularité (que nous n’énoncerons pas ici), alors l’estimateur du maximum de vraisemblance est fortement convergent et asymptotiquement efficace \[ \sqrt{n} \left( \hat\theta_n - \theta\right) \stackrel{\mathcal{L}}{\to } \mathcal{N}(0,I^{-1}(\theta)). \]
options(warn=-1)
set.seed(123)
n = 200
X = sample(0:1,size=n,replace=TRUE)
neglogL = function(p){-sum(log(dbinom(X,size=1,prob=p)))}
(pml = optim(fn=neglogL,par=p0,method="BFGS")$par)
## [1] 0.485
S = rep(NA,1e3)
for(s in 1:length(S)){
Xs = sample(X,size=n,replace=TRUE)
neglogL = function(p){-sum(log(dbinom(Xs,size=1,prob=p)))}
S[s] = optim(fn=neglogL,par=p0,method="BFGS")$par
}
mean(S)
## [1] 0.4842846
## [1] 0.001293088
## [1] 0.00125
hist(S,probability = TRUE)
u=seq(min(S),max(S),length=251)
v=dnorm(u,mean(S),sd(S))
lines(u,v,lty=2,col="red")
X = Davis$height
par(mfrow=c(1,2))
plot(ecdf(X),main="")
u=seq(min(X)-2,max(X)+2,length=251)
v=ecdf(X)(u)
plot(u,v,type="s")
m=100
inf_sample=function(n=10){
X=matrix(rnorm(n*m),nrow=m,ncol=n)
Xs=t(apply(X,1,sort))
Pe_inf=matrix(rep((0:(n-1))/n,
each=m),nrow=m,ncol=n)
Pe_sup=matrix(rep((0:n)/n,each=m),
nrow=m,ncol=n)
Pt=pnorm(Xs)
D1=abs(Pe_inf-Pt)
D2=abs(Pe_sup-Pt)
Df=(D1+D2)/2+abs(D2-D1)/2
return(apply(Df,1,max))
}
B=matrix(NA,100,20)
for(i in 1:20){
B[,i]=inf_sample(i*10)
}
colnames(B)=as.character(seq(10,200,by=10))
boxplot(B)
u=seq(-3,3,by=.1)
plot(u,u,ylim=c(0,1),col="white")
M=matrix(NA,length(u),1000)
for(m in 1:1000){
n=100
x=rnorm(n)
Femp=Vectorize(function(t) mean(x<=t))
v=Femp(u)
M[,m]=v
lines(u,v,col='light blue',type="s")
}
lines(u,apply(M,1,mean),col="red",type="l")
lines(u,apply(M,1,function(x) quantile(x,.05)),
col="red",type="s")
lines(u,apply(M,1,function(x) quantile(x,.95)),
col="red",type="s")
x0=-1
y=M[which(u==x0),]
hist(y,probability=TRUE,
breaks=seq(.015,0.55,by=.01))
vu=seq(0,1,by=.001)
lines(vu,dnorm(vu,pnorm(x0),
sqrt((pnorm(x0)*(1-pnorm(x0)))/100)),
col="red")
library(MASS)
hist(X,col=rgb(0,0,1,.4), border="white",proba=TRUE,xlab="",main="")
(param <- fitdistr(X,"normal")$estimate)
## mean sd
## 170.56500 8.90987
f1 <- function(x) dnorm(x,param[1],param[2])
x=seq(100,210,by=.2)
lines(x,f1(x),lty=2,col="red")
lines(density(X))
## mixtools package, version 1.1.0, Released 2017-03-10
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.
## One of the variances is going to zero; trying new starting values.
## number of iterations= 237
## [1] 0.1502809 163.7585980 171.7687831 2.5671638 9.0894373
f2 <- function(x){ param12[1]*dnorm(x,param12[2],param12[4])+
(1-param12[1])*dnorm(x,param12[3],param12[5]) }
hist(X,col=rgb(0,0,1,.4), border="white",proba=TRUE,xlab="",main="",ylim=c(0,.05))
(param <- fitdistr(X,"normal")$estimate)
## mean sd
## 170.56500 8.90987
logdf <- function(x,parameter){
p <- parameter[1]
m1 <- parameter[2]
s1 <- parameter[4]
m2 <- parameter[3]
s2 <- parameter[5]
return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2)))
}
logL <- function(parameter) -sum(logdf(X,parameter))
Amat <- matrix(c(1,-1,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,0,0,0,1), 4, 5)
bvec <- c(0,-1,0,0)
constrOptim(c(.5,160,180,10,10), logL, NULL, ui = Amat, ci = bvec)$par
## [1] 0.5996263 165.2690084 178.4991624 5.9447675 6.3564746
## [1] 0.44
## mean sd
## 164.714286 5.633808
## mean sd
## 178.011364 6.404001
hist(X,col=rgb(0,0,1,.4), border="white",proba=TRUE,xlab="",main="")
(param <- fitdistr(X,"normal")$estimate)
## mean sd
## 170.56500 8.90987
f4 <- function(x) pM*dnorm(x,paramM[1],paramM[2])+(1-pM)*dnorm(x,paramF[1],paramF[2])
lines(x,f4(x),lwd=3,col="blue")
## [1] 0 1 0 1 1 0 1 1 1 0 1 0 1 1 0 1 0 0 0 1
neglogL = function(p){-sum(log(dbinom(X,size=1,prob=p)))}
(pml = optim(fn=neglogL,par=p0,method="BFGS")$par)
## [1] 0.5499996
nx=sum(X==1)
f = expression(nx*log(p)+(n-nx)*log(1-p))
Df = D(f, "p")
Df2 = D(Df, "p")
p=p0
score=eval(Df)
(IF=-eval(Df2))
## [1] 80
## [1] 80
On veut tester ici si la probabilité d’avoir pile peut-être égale à 50%, \(H_0:p=1/2\).
Considérons un test \(H_0:\theta=\theta_0\), contre \(H_1:\theta\neq\theta_0\), alors la statistique de Wald est \[ {\displaystyle W={\frac {({\widehat {\theta }}-\theta _{0})^{2}}{\operatorname {var} ({\hat {\theta }})}}}\] qui suit asymptotiquement, sous \(H_0\) une loi \(\chi^2(1)\).
## [1] 0.199997
## [1] TRUE
u=seq(0,qchisq(.99,df=1),length=601)
v=dchisq(u,df=1)
plot(u,v,col="white",xlab="",ylab="")
text(2,2,"ACCEPT",col="green")
text(5.5,2,"REJECT",col="red")
idx=which(u<=qchisq(1-alpha,df=1))
polygon(c(u[idx],rev(u[idx])),c(v[idx],rep(0,length(idx))),col="green",border=NA)
idx=which(u>qchisq(1-alpha,df=1))
polygon(c(u[idx],rev(u[idx])),c(v[idx],rep(0,length(idx))),col="red",border=NA)
abline(v=qchisq(1-alpha,df=1),lty=2)
abline(v=T,col="blue")
or we can compute the p-value
## [1] 0.6547233
u=seq(0,qchisq(.99,df=1),length=601)
v=dchisq(u,df=1)
plot(u,v,col="white",xlab="",ylab="")
text(3,2,"p-VALUE",col="blue")
idx=which(u<=T)
polygon(c(u[idx],rev(u[idx])),c(v[idx],rep(0,length(idx))),col="yellow",border=NA)
idx=which(u>T)
polygon(c(u[idx],rev(u[idx])),c(v[idx],rep(0,length(idx))),col="blue",border=NA)
abline(v=qchisq(1-alpha,df=1),lty=2)
abline(v=T,col="blue")
Si on cherche à tester \({\displaystyle H_{0}\,:\,\theta \in \Theta _{0}}\) contre \({\displaystyle H_{1}\,:\,\theta \notin \Theta _{0}}\), la statistique de test est \[ {\displaystyle LR=-2\log \left[{\frac {\sup _{\theta \in \Theta _{0}}{\mathcal {L}}(\theta )}{\sup _{\theta \in \Theta }{\mathcal {L}}(\theta )}}\right]} \] qui suit asymptotiquement, sous \(H_0\), une loi \(\chi^2(1)\) si \(\theta\) est univarié.
## [1] 0.2003347
## [1] TRUE
Considérons un test \(H_0:\theta=\theta_0\), contre \(H_1:\theta\neq\theta_0\), et notons \({\displaystyle U(\theta )}\) le score, \[ {\displaystyle U(\theta )={\frac {\partial \log L(\theta \mid x)}{\partial \theta }}.} \] et \(I(\theta)\) l’information de Fisher, \[ {\displaystyle I(\theta )=-\mathbb{E} \left[{\frac {\partial ^{2}}{\partial \theta ^{2}}}\log \mathcal{L}(\theta )\right]} \] alors la statistique du score (ou le multiplicateur de Lagrande) est \[ {\displaystyle S(\theta _{0})={\frac {U(\theta _{0})^{2}}{I(\theta _{0})}}} \] qui suit asymptotiquement, sous \(H_0\), une loi \(\chi^2(1)\).
## [1] 0.2
## [1] TRUE
## shape rate
## 1.3562877 0.8207035
## (0.1732663) (0.1263275)
## [1] 1.016686 1.695890
log_lik = function(theta){
a = theta[1]
b = theta[2]
logL = sum(log(dgamma(x,a,b)))
return(-logL)
}
optim(c(1,1),log_lik)
## $par
## [1] 1.3558113 0.8206505
##
## $value
## [1] 147.6097
##
## $counts
## function gradient
## 45 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
options(warn=-1)
prof_log_lik = function(a){
b = (optim(1,function(z) -sum(log(dgamma(x,a,z)))))$par
return(-sum(log(dgamma(x,a,b))))
}
optim(1,prof_log_lik)
## $par
## [1] 1.356445
##
## $value
## [1] 147.6097
##
## $counts
## function gradient
## 26 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
vx = seq(.5,3,length=101)
vl = -Vectorize(prof_log_lik)(vx)
plot(vx,vl,type="l")
abline(v=optim(1,prof_log_lik)$par,lty=2)
abline(h=-optim(1,prof_log_lik)$value)
abline(h=-optim(1,prof_log_lik)$value-qchisq(.95,1)/2)
segments(F$estimate[1]-1.96*F$sd[1],
-170,F$estimate[1]+1.96*F$sd[1],
-170,lwd=3,col="blue")
borne = -optim(1,prof_log_lik)$value-qchisq(.95,1)/2
(b1 = uniroot(function(z) Vectorize(prof_log_lik)(z)+borne,c(.5,1.5))$root)
## [1] 1.046505
## [1] 1.727274
Avec les notations usuelles, on suppose disposer d’observations individuelles \(y_{i,j}\) appartenant à un groupe \(j\in\{1,2,\cdots,k\}\). On suppose que \[ y_{i,j} = \beta_0+\beta_j+\varepsilon_{i,j} \] où (classiquement) les \(\varepsilon_{i,j}\) sont i.i.d. de loi \(\mathcal{N}(0,\sigma^2)\).
Sur la forme, on a un soucis car ce modèle n’est pas identiable. On a trois options possible.
Ces trois modèles sont équivalents.
L’estimateur du maximum de vraisemblance du modèle est ici \[ \widehat{\alpha}_j=\frac{1}{n_j}\sum_{i\in\mathcal{I}_j}y_{i,j}=\overline{y}_{\cdot,j} \] \[ \widehat{\beta}_0=\frac{1}{n_1}\sum_{i\in\mathcal{I}_1}y_{i,j}=\overline{y}_{\cdot,1}~~~~~~~~ \widehat{\beta}_j=\overline{y}_{\cdot,j}-\overline{y}_{\cdot,1} \] ou \[ \widehat{\gamma}_0=\frac{1}{n}\sum_{j=1}^k\sum_{i\in\mathcal{I}_j}y_{i,j}~~~~~~~~ \widehat{\gamma}_j=\overline{y}_{\cdot,j}-\overline{y} \] On va alors chercher à tester \(H_0:\alpha_1=\cdots=\alpha_k\), ou encore \(H_0:\beta_2=\cdots=\beta_k=0\) ou enfin \(H_0:\gamma_1=\cdots=\gamma_k=0\). Techniquement, on a ici un test multiple.
Notons pour commencer que \[
{\displaystyle SCR_{\text{total}}=SCR_{\text{facteur}}+SCR_{\text{résidus}}~}
\] où \[
\underbrace{{\displaystyle \sum _{j=1}^{k}\sum _{i\in\mathcal{I}_j}(y_{i,j}-{\overline {y}})^{2}}}_{SCR_{\text{total}}} =
\underbrace{{\displaystyle \sum _{j=1}^{k}(\overline{y}_{\cdot,j}-{\overline {y}})^{2}}}_{SCR_{\text{facteur}}} +
\underbrace{{\displaystyle \sum _{j=1}^{k}\sum _{i\in\mathcal{I}_j}(y_{i,j}-{\overline {y}_{\cdot,j}})^{2}}}_{SCR_{\text{résidus}}}
\] Sous une hypothèse de normalité des résidus \(\varepsilon_{i,j}\), on peut montrer que si \(H_0\) est vraie \[
{\displaystyle SCR_{\text{facteur}}=\sum _{j=1}^{k}n_{j}({\overline {y_{\cdot,j}}}-{\overline {y}})^{2}\sim \chi ^{2}(DL_{\text{facteur}})\qquad {\text{avec}}~DL_{\text{facteur}}=\sum _{j=1}^{k-1}1=k-1}
\] \[
{\displaystyle SCR_{\text{résidus}}=\sum _{j=1}^{k}\sum _{i\in\mathcal{I}_j}(y_{i,j}-{\overline {y_{\cdot,j}}})^{2}\sim \chi ^{2}(DL_{\text{résidus}})\quad {\text{avec}}~DL_{\text{résidus}}=\sum_{j=1}^{k}(n_{j}-1)=(n_{1}-1)+(n_{2}-1)+\cdots +(n_{k}-1)=n-k}
\] Alors \[
{\displaystyle F={\frac {\dfrac {SCR_{\text{facteur}}}{k-1}}{\dfrac {SCR_{\text{résidus}}}{n-k}}}\sim \mathcal F(k-1,n-k)}
\] Dans les données Davis
, on a deux groupes (correspondant au genre).
## Call:
## aov(formula = weight ~ sex, data = Davis)
##
## Terms:
## sex Residuals
## Sum of Squares 17799.20 17522.79
## Deg. of Freedom 1 198
##
## Residual standard error: 9.407389
## Estimated effects may be unbalanced
ou
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 17799 17799 201.1 <2e-16 ***
## Residuals 198 17523 88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Le genre sex
dans la base Davis
est ici une variable catégorielle, ou factorielle, prenant deux modalités F
et M
(dans cet ordre)
## Factor w/ 2 levels "F","M": 2 1 1 2 1 2 2 2 2 2 ...
On peut changer l’ordre (et la modalité dite de référence) en utilisant
## Factor w/ 2 levels "M","F": 1 2 2 1 2 1 1 1 1 1 ...
mais pour ne pas embrouiller, laissons F
comme modalité de référence par la suite
Pour créer les variables indicatrices associées, on peut utiliser
## sexF sexM
## 1 2 0 1
## 2 1 1 0
## 3 1 1 0
## 4 2 0 1
## 5 1 1 0
## 6 2 0 1
## 7 2 0 1
## 8 2 0 1
## 9 2 0 1
## 10 2 0 1
\[ y_i = \beta_0+\beta_1 x_i+\varepsilon_i \] where
Consider \(n\) pairs of observations \((x_i,y_i)\).
Hypothese 1 les observations \((y_i,x_i)\) sont des realisations de variables aleatoires \((Y_i,X_i)\).
Hypothese 2 les erreurs \(\varepsilon_i\) sont des realisations de variables aleatoires \(\varepsilon_i\), independantes, centrees \(\mathbb{E}[\varepsilon_i]=0\) et de variance constante \(\text{Var}[\varepsilon_i]=\sigma^2\).
Hypothese 2’ les erreurs \(\varepsilon_i\) sont des réalisations de variables aleatoires \(\varepsilon_i\), independantes, de loi \(\mathcal{N}(0,\sigma^2)\).
On suppose ici que \((Y|X=x)\sim\mathcal{N}( \beta_0+\beta_1 x,\sigma^2)\)
##
## Call:
## lm(formula = weight ~ height, data = Davis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.658 -5.381 -0.555 4.807 42.894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -130.91040 11.52792 -11.36 <2e-16 ***
## height 1.15009 0.06749 17.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.505 on 198 degrees of freedom
## Multiple R-squared: 0.5946, Adjusted R-squared: 0.5925
## F-statistic: 290.4 on 1 and 198 DF, p-value: < 2.2e-16
La somme des carres de residus est \[ SCR(\beta_0,\beta_1)=\sum_{i=1}^n \big(y_i - [ \beta_0+\beta_1 x_i]\big)^2 \] On note \[ (\widehat{\beta}_0,\widehat{\beta}_1)=\text{argmin}\left\lbrace SCR(\beta_0,\beta_1)\right\rbrace \] First order conditions are here \[ \left.\frac{\partial SCR(\beta_0,\beta_1)}{\partial\beta_0}\right|_{(\widehat{\beta}_0,\widehat{\beta}_1)}=0 \] while \[ \left.\frac{\partial SCR(\beta_0,\beta_1)}{\partial\beta_1}\right|_{(\widehat{\beta}_0,\widehat{\beta}_1)}=0 \] Then \[ {\hat {\beta}_0}={\bar {y}}-{\hat {\beta }_1}{\bar {x}} \] while \[ {\hat {\beta }_1}={\frac {\sum_{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum_{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}={\frac {s_{xy}}{s_{x}}}= \text{corr}(x,y)\cdot{\frac {s_{y}}{s_{x}}} \] or \[ {\hat {\beta }_1}=\frac{\sum_{i=1}^{n}x_{i}y_{i}-\bar {x}\bar {y}}{\sum_{i=1}^{n}x_{i}^2-n{\bar {x}}} \]
## [1] 1.150092
## [1] -130.9104
par(mfrow=c(1,2))
plot(X,Y)
abline(v=mean(X),col="red",lty=2)
abline(h=mean(Y),col="red",lty=2)
abline(a = Ahat, b=Bhat, lwd=3, col="red")
Y_hat = Ahat + Bhat*X
E_hat = Y-Y_hat
plot(X,Y)
abline(a = Ahat, b=Bhat, lwd=3, col="red")
segments(X,Y,X,Y_hat,col="light blue")
SCR = function(B){sum((Y - (Ahat + B*X))^2)}
x = seq(-10,10)
y = Vectorize(SCR)(x)
plot(x,y,type="l",ylab="Somme des carres des erreurs")
abline(v=Bhat,col="red",lty=2)
The fitted values (or predictions) are \[ \widehat{y}_i = \hat {\beta }_0+\hat {\beta }_1x_i \] and residuals are \[ \widehat{\varepsilon}_i={y}_i-\widehat{y}_i \] Observe (first first order condition) that \[ \sum_{i=1}^n \widehat{\varepsilon}_i=0 \] (provided that there is an intercept term - \(\beta_0\) - in the model) and (second first order condition) \[ \sum_{i=1}^n x_i\widehat{\varepsilon}_i=0 \] Notons que \[ \underbrace{\sum_{i=1}^n (y_i-\bar y)^2}_{\text{total sum of squares}} = \underbrace{\sum_{i=1}^n (y_i-\hat y_i)^2}_{\text{residual sum of squares}} + \underbrace{\sum_{i=1}^n (\hat y_i-\bar y)^2}_{\text{explained sum of squares}} \]
Le coefficient de determination \(R^2\) est \[R^2=1-\frac{x}{x} \]
On peut montrer que \[ R^2=\widehat{\beta}_1^2\frac{s_x}{s_y}=\frac{s_{xy}^2}{s_xs_y}=\text{cor}[x,y]^2 \]
Considerons l’estimateur suivant de \(\sigma^2\), \[ s^2=\frac{1}{n-2}\sum_{i=1}^n\hat\varepsilon_i^2 \]
Parmi les autres interprétations des estimateurs, notons que les paramètres sont des fonctions linéaires des \(y_i\): \[ \widehat{\beta}_1 = \sum_{i=1}^n\color{blue}{\omega_{1,i}}\cdot y_i~~~\color{blue}{\omega_{1,i}}=\frac{x_i-\overline{x}}{s_x^2} \] et \[ \widehat{\beta}_0 = \sum_{i=1}^n\color{red}{\omega_{0,i}}\cdot y_i~~~\color{red}{\omega_{0,i}}=\frac{1}{n}-\overline{x}\color{blue}{\omega_{1,i}} \] Attention: la notation \(\omega\) ne signifie pas vraiment que l’on ait ici des poids : les \(\omega_i\) peuvent être négatifs. Par exemple, on notera que \[ \sum_{i=1}^n \color{blue}{\omega_{1,i}}=0,~~ \sum_{i=1}^n \color{blue}{\omega_{1,i}}\cdot x_i=1,~~ \sum_{i=1}^n \color{blue}{\omega_{1,i}}^2=\frac{1}{s_x^2} \]
\(\widehat{\beta}_0\) et \(\widehat{\beta}_1\) sont des estimateurs sans biais de \(\beta_0\) et \(\beta_1\) respectivement (cf section simulations), \[ \mathbb{E}[\widehat{\beta}_0]={\beta}_0~~,~~ \mathbb{E}[\widehat{\beta}_1]={\beta}_1 \] Les variances sont respectivement \[ \text{Var}[\widehat{\beta}_0]=\sigma^2\left(\frac{1}{n}+\frac{\overline{x}^2}{s_x^2}\right)~~,~~\text{Var}[\widehat{\beta}_1]=\frac{\sigma^2}{s_x^2} \] mais comme \(\sigma\) est ici inconnue, on estime ces variances par \[ \widehat{\text{Var}}[\widehat{\beta}_0]=\widehat{\sigma}^2\left(\frac{1}{n}+\frac{\overline{x}^2}{s_x^2}\right)~~,~~\widehat{\text{Var}}[\widehat{\beta}_1]=\frac{\widehat{\sigma}^2}{s_x^2} \] soit pour le second un ecart-type \[ {\displaystyle s_{\hat {\beta }_1}={\sqrt {\frac {{\frac {1}{n-2}}\sum_{i=1}^{n}{\hat{\varepsilon }}_{i}^{,2}}{\sum_{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}}} \] et pour le premier, une reecriture sous la forme \[ {\displaystyle s_{\hat {\beta }_0}=s_{\hat {\beta }_1}{\sqrt {{\frac {1}{n}}\sum_{i=1}^{n}x_{i}^{2}}}={\sqrt {{\frac{1}{n(n-2)}}\left(\sum_{i=1}^{n}{\hat{\varepsilon}}_{j}^{,2}\right){\frac {\sum_{i=1}^{n}x_{i}^{2}}{\sum_{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}}}} \]
Rappelons que \[ \underbrace{\sum_{i=1}^n (y_i-\bar y)^2}_{\text{TSS}} = \underbrace{\sum_{i=1}^n (y_i-\hat y_i)^2}_{\text{RSS}} + \underbrace{\sum_{i=1}^n (\hat y_i-\bar y)^2}_{\text{ESS}} \] La variance de \(y\) est estimée par \[ s^2_y=\frac{1}{n-1}\sum_{i=1}^n (y_i-\bar y)^2=\frac{\text{TSS}}{n-1} \] La variance de \(\varepsilon\) est estimée par \[ \widehat{\sigma}^2=\frac{1}{n-2}\sum_{i=1}^n (y_i-\widehat{y}_i)^2=\frac{1}{n-2}\sum_{i=1}^n \widehat{\varepsilon}_i^2=\frac{\text{RSS}}{n-2} \]
Si \(Y|X\) suit une loi normale, avec \(\varepsilon\) suit une loi normale, et en supposant l’indépendance des résidus, on obtient la normalité de \(\widehat{\beta}_0\) et \(\widehat{\beta}_1\). Aussi, une statistique naturelle de test de \(H_0:\beta_1=\color{blue}{0}\) contre \(H_1:\beta_0\neq\color{blue}0\) est basé sur \[ T=\frac{\widehat{\beta}_1-\color{blue}{0}}{\sqrt{\widehat{\text{Var}}[\widehat{\beta}_1]}} \] qui suit, sous \(H_0\)une loi de Student \(\mathcal{S}td(n-2)\). Comme c’est un test bi-latéral, on utilise alors une des deux méthodes - région critique : si \(t^{-1}_{n-2}(1-\alpha)\) désigne le quantile de niveau \(1-\alpha\) de la loi \(\mathcal{S}td(n-2)\), on rejette \(H_0\) si \(|T|>t^{-1}_{1,n-2}(1-\alpha/2)\), - p-value : on rejette \(H_0\) si \(p=\mathbb{P}[|Y|>|T||Y\sim\mathcal{S}td(n-2)]<\alpha\).
Considerons le test \(H_0:\beta_1=0\) contre \(H_1:\beta_0\neq0\). La statistique de Fisher (analyse de la variance) vise à comparer les résidus de deux modèles : (0) \(y_i=\beta_0+\varepsilon_i\) et (0) \(y_i=\beta_0+\beta_1x_i+\varepsilon_i\), i.e. \[ F=\frac{(TSS - RSS)/1}{RSS/(n-2)}=\frac{ESS}{RSS/(n-2)} \] qui suit, sous \(H_0\) une loi de Fisher \(\mathcal{F_{1,n-2}}\). On utilise alors une des deux méthodes - région critique : si \(F^{-1}_{1,n-2}(1-\alpha)\) désigne le quantile de niveau \(1-\alpha\) de la loi \(\mathcal{F_{1,n-2}}\), on rejette \(H_0\) si \(F>F^{-1}_{1,n-2}(1-\alpha)\), - p-value : on rejette \(H_0\) si \(p=\mathbb{P}[Y>F|Y\sim\mathcal{F_{1,n-2}}]<\alpha\).
On peut noter que \[ F=(n-2)\frac{R^2}{1-R^2} \]
## [1] 0.5945555
## [1] 290.353
## value numdf dendf
## 290.353 1.000 198.000
Dans le cas de la significativité de la pente, i.e. \(H_0:\beta_1=0\), notons que \[ T^2=\frac{\widehat{\beta}_1^2}{\widehat{\sigma}^2/s_x^2}=W= \frac{\widehat{\beta}_1^2s_x^2}{\widehat{\sigma}^2}= \frac{ESS}{RSS/(n-2)}=F. \]
Comme \(\varepsilon\perp x\), \[ \widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1x_i \] et si on considere une nouvelle observation \(\color{red}{x}\), on aurait \[ \widehat{y}_{\color{red}{x}}=\widehat{\beta}_0+\widehat{\beta}_1\color{red}{x} \] et on aura comme observation \[y_{\color{red}{x}}=\beta_0+\beta_1{\color{red}{x}}+\varepsilon\] ou \(\varepsilon\) est un bruit imprevisible.
Notons que \[ y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}=\big(\beta_0+\beta_1{\color{red}{x}}+\varepsilon\big)-\big(\widehat{\beta}_0+\widehat{\beta}_1\color{red}{x}\big) \] soit \[ y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}=\varepsilon+\big(\beta_0-\widehat{\beta}_0\big)+\big(\beta_1-\widehat{\beta}_1\big){\color{red}{x}} \] de telle sorte que \[ \mathbb{E}[y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}]=0 \] donc \(\widehat{y}_{\color{red}{x}}\) est un estimateur dans biais de \(y_{\color{red}{x}}\).
Si on continue, \[ \text{Var}[y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}] =\text{Var}[\varepsilon]+\text{Var}\big(\widehat{\beta}_0+\widehat{\beta}_1\color{red}{x}\big) \] qui se réécrit \[ \text{Var}[y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}] = \sigma^2 + \left(\frac{\sigma^2}{n}+\frac{\sigma^2(\color{red}{x}-\overline{x})^2}{s_x^2}\right) \] Frees (Regression Modeling with Actuarial and Financial Applications, 2014) appele “standard error of prediction” la racine carrée de l’estimateur de cette grandeur \[ \sqrt{\widehat{\text{Var}}[y_{\color{red}{x}}-\widehat{y}_{\color{red}{x}}] }=\widehat{\sigma}\cdot\sqrt{1+\frac{1}{n}+\frac{(\color{red}{x}-\overline{x})^2}{s_x^2}} \] En notant que si les residus sont supposes Gaussien, on peut alors construire un intervalle de confiance de prédiction pour \(y_{\color{red}{x}}\) : \[ \underbrace{\widehat{\beta}_0+\widehat{\beta}_1\color{red}{x}}_{\widehat{y}_{\color{red}{x}}} \pm t_{n-2,1-\alpha/2}\cdot\widehat{\sigma}\cdot\sqrt{1+\frac{1}{n}+\frac{(\color{red}{x}-\overline{x})^2}{s_x^2}} \]
reg = lm(weight~height, data=Davis)
par(mfrow=c(1,1))
plot(Davis$height,Davis$weight)
u=seq(min(Davis$height),max(Davis$height),lenght=251)
v=predict(reg,newdata = data.frame(height=u),interval = "prediction")
lines(u,v[,1],col="blue",lwd=2)
lines(u,v[,2],col="red",lty=2)
lines(u,v[,3],col="red",lty=2)
On retrouve dans l’expression précédante deux sources d’erreur - l’erreur de modèle, venant du fait que la réalisation \(y\) est \(\beta_0+\beta_1x\) auquel s’ajoute un bruit \(\varepsilon\) - l’erreur d’estimation, venant du fait que \(\beta_0+\beta_1x\) est incertain
Pour l’erreur d’estimation, notons que \[ \text{Var}[\widehat{y}_{\color{red}{x}}] = \left(\frac{\sigma^2}{n}+\frac{\sigma^2(\color{red}{x}-\overline{x})^2}{s_x^2}\right) \] d’où un intervalle de confiance pour \(\widehat{y}_{\color{red}{x}}\) de la forme \[ \underbrace{\widehat{\beta}_0+\widehat{\beta}_1\color{red}{x}}_{\widehat{y}_{\color{red}{x}}} \pm t_{n-2,1-\alpha/2}\cdot\widehat{\sigma}\cdot\sqrt{\frac{1}{n}+\frac{(\color{red}{x}-\overline{x})^2}{s_x^2}} \]
reg = lm(weight~height, data=Davis)
par(mfrow=c(1,1))
plot(Davis$height,Davis$weight)
u=seq(min(Davis$height),max(Davis$height),lenght=251)
v=predict(reg,newdata = data.frame(height=u),interval = "confidence")
lines(u,v[,1],col="blue",lwd=2)
lines(u,v[,2],col="red",lty=2)
lines(u,v[,3],col="red",lty=2)
Deux méthodes de boostrap sont envisageables.
La première consiste a adapter directement la méthode dans le cas univarié en tirant avec remise de paires \((x_i,y_i)\)
library(scales)
BETA = matrix(NA,1000,2)
for(s in 1:1000){
set.seed(s)
idx = sample(1:nrow(Davis),nrow(Davis),replace=TRUE)
reg_sim = lm(weight~height, data=Davis[idx,])
BETA[s,] = reg_sim$coefficients
}
pic_ani = function(k){
plot(Davis$height,Davis$weight,ylab="Weight (kg)",xlab="Height (cm)",col="white")
for(s in 1:100) abline(BETA[s,1],BETA[s,2],col=alpha("light blue",.3))
set.seed(k)
idx = sample(1:nrow(Davis),nrow(Davis),replace=TRUE)
points(Davis$height,Davis$weight,col="grey")
points(Davis$height[idx],Davis$weight[idx],pch=19)
abline(BETA[k,1],BETA[k,2],col="blue",lwd=2)
}
for (k in 1:50) {pic_ani(k)}