2.2.1 Loi binomiale

\[{\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)

2.2.2 Loi de Poisson

\[{\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)

2.2.3 Loi Exponentielle

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="")

2.2.4 Loi Gaussienne

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})}\).

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:sm':
## 
##     muscle

fitdistr(Davis$height,"normal")

##       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)

2.2.5 Loi du Chi-Deux

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="")

2.2.6 Loi de Student

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)\).

2.2.7 Loi de Fisher

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\).

2.2.8 Vecteur Gaussien

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})\):

• on simule \(d\) réalisations indépendantes d’une \(\mathcal{N}(0,1)\) que l’on stocke dans \(\boldsymbol{y}\), i.e. y=rnorm(d)
• on diagonalise la matrice des vecteurs propres \(\boldsymbol{P}\) et la matrice diagonale des racines carrées des valeurs propres \(\boldsymbol{D}^{1/2}\), puis on calcule \(\boldsymbol{P} \boldsymbol{D}^{1/2} \boldsymbol{y}\).
• la réalisation désirée est: \(\boldsymbol{x} = \boldsymbol{\mu} + \boldsymbol{P} \boldsymbol{D}^{1/2} \boldsymbol{y}\).

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

A0.5%*% matrix(rnorm(2*3),nr=2)

##            [,1]      [,2]     [,3]
## [1,] -0.8129428 0.6786791 2.119525
## [2,] -1.4852915 1.5903142 1.664039

m = 10^4
y = matrix(rnorm(2*m),nr=2)
X0=y
apply(X0,1,var)

## [1] 1.008306 1.032374

cov(X0[1,],X0[2,])

## [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")

par(mfrow = c(1,1))

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,

• si \(\boldsymbol{A}\) est une matrice idempotente, \(\boldsymbol{Y}^\top \boldsymbol{A} \boldsymbol{Y} \sim \chi^2(r)\) où \(r=\operatorname{rang}(\boldsymbol{A})\).
• si \(\boldsymbol{Y}^\top \boldsymbol{A} \boldsymbol{Y} \sim \chi^2(s)\) alors \(\boldsymbol{A}\) est une matrice idempotente et \(s=\operatorname{rang}(\boldsymbol{A})\).

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,

• si \(\boldsymbol{A} \boldsymbol{\Sigma} \boldsymbol{A}=\boldsymbol{A}\), \(\boldsymbol{Y}^\top \boldsymbol{A} \boldsymbol{Y} \sim \chi^2(r,\lambda)\) où \(r=\operatorname{rang}(\boldsymbol{A})\) et \(\lambda=\boldsymbol{\mu}^\top \boldsymbol{A} \boldsymbol{\mu}\).
• si \(\boldsymbol{Y}^\top \boldsymbol{\Sigma} \boldsymbol{Y} \sim \chi^2(s,\lambda)\) alors \(\boldsymbol{A} \boldsymbol{\Sigma} \boldsymbol{A}=\boldsymbol{A}\), \(s=\operatorname(\boldsymbol{A})\) et \(\lambda=\boldsymbol{\mu}^\top \boldsymbol{A} \boldsymbol{\mu}\).

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

• Les vecteurs \(\boldsymbol{P}_{F_1} \boldsymbol{X}, \dots,\boldsymbol{P}_{F_m} \boldsymbol{X}\) sont deux à deux indépendants et de loi \(\mathcal{N}(\boldsymbol{P}_{F_i}\boldsymbol{\mu}, \sigma^2 \boldsymbol{P}_{F_i})\) pour \(1\leq i \leq m\).
• Les variables aléatoires \(\|\boldsymbol{P}_{F_i} (\boldsymbol{X}-\boldsymbol{\mu})\|^2/\sigma^2\) pour \(1\leq i\leq m\) sont 2 à 2 indépendantes et de loi \(\chi^2(d_i)\).

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

• \(\boldsymbol{P}_F \boldsymbol{X} \sim \mathcal{N}(\boldsymbol{ 0}, \boldsymbol{P}_F)\) et \(\boldsymbol{P}_{F^\bot} \boldsymbol{X} \sim \mathcal{N}(\boldsymbol{ 0}, \boldsymbol{P}_{F^\bot})\) et \((\boldsymbol{P}_F \boldsymbol{X}) \bot (\boldsymbol{P}_{F^\bot} \boldsymbol{X})\)
• \(\|\boldsymbol{P}_F \boldsymbol{X} \|^2 \sim \chi^2(p)\), \(\|\boldsymbol{P}_{F^\bot} \boldsymbol{X} \|^2 \sim \chi^2(d-p)\) et ces deux variables sont indépendantes.

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}\).

2.2.9 La famille exponentielle

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.

2.3.1 central limit theorem

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")

2.3.2 Delta Method

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.

2.4.1 Statistique

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

• alors pour toute fonction mesurable \(\varphi\) bijective, \(\varphi(T)\) est une statistique exhaustive complète
• alors \(T\) est une statistique exhaustive minimale
• alors \(T\) est indépendante de toute statistique libre du modèle (théorème de Basu)

Dans un modèle statistique paramétrique, il existe toujours une statistique exhaustive minimale mais pas toujours de statistique exhaustive complète.

2.4.2 Vraisemblance, score et information de Fisher

Example d’un pile ou face

set.seed(123)
n = 20
(X = sample(0:1,size=n,replace=TRUE))

##  [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

slope_h = function(f,x,h=1e-6) (f(x+h)-f(x-h))/(2*h)

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

n/(p0*(1-p0))

## [1] 308

2.4.3 Maximum de Vraisemblance

set.seed(123)
n = 20
(X = sample(0:1,size=n,replace=TRUE))

##  [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

mean(X)

## [1] 0.55

2.4.4 Estimateur (ponctuel)

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))\)

2.4.5 Maximum de Vraisemblance : propriétés asymptotiques

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

var(S)

## [1] 0.001293088

(p0*(1-p0))/n

## [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")

2.4.6 Methode des Moments

2.4.7 Statistique Nonparametrique

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")

par(mfrow=c(1,2))
hist(X, main="")
plot(density(X))

 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")

2.4.8 Melange

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))

library(mixtools) 

## 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.

mix <- normalmixEM(X)

## One of the variances is going to zero;  trying new starting values.
## number of iterations= 237

(param12 <- c(mix$lambda[1],mix$mu,mix$sigma)) 

## [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

lines(x,f2(x),lwd=2, col="red") 
lines(density(X))

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

(pM <- mean(Davis$sex=="M"))

## [1] 0.44

(paramF <- fitdistr(X[Davis$sex=="F"],"normal")$estimate)

##       mean         sd 
## 164.714286   5.633808

(paramM <- fitdistr(X[Davis$sex=="M"],"normal")$estimate)

##       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")

2.5.1 Decision

set.seed(123)
n = 20
(X = sample(0:1,size=n,replace=TRUE))

##  [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/(p0*(1-p0)/n)

## [1] 80

On veut tester ici si la probabilité d’avoir pile peut-être égale à 50%, \(H_0:p=1/2\).

2.5.2 Test de Wald

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)\).

alpha=0.05
(T=(pml-p0)^2*IF)

## [1] 0.199997

T<qchisq(1-alpha,df=1)

## [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-pchisq(T,df=1)

## [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")

2.5.3 Test du Rapport de Vraisemblance

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é.

logL = function(p){sum(log(dbinom(X,size=1,prob=p)))}
(T=2*(logL(pml)-logL(p0)))

## [1] 0.2003347

T<qchisq(1-alpha,df=1)

## [1] TRUE

2.5.4 Test du Score

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)\).

(T=score^2/IF)

## [1] 0.2

T<qchisq(1-alpha,df=1)

## [1] TRUE

4.5.1 Test de Student

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\).

4.5.2 Test de Fisher

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} \]

reg = lm(weight~height, data=Davis)
summary(reg)$r.squared

## [1] 0.5945555

(nrow(Davis)-2)*summary(reg)$r.squared/(1-summary(reg)$r.squared)

## [1] 290.353

summary(reg)$fstatistic

##   value   numdf   dendf 
## 290.353   1.000 198.000

4.5.3 Test de Wald

4.5.4 Lien entre les trois tests

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. \]

4.6.1 Prévision

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}}\).

4.6.2 Incertitude

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)

4.6.3 Simulations

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)}