Tag Archives: convergence

Convergence and Asymptotic Results

Last week, in our mathematical statistics course, we’ve seen the law of large numbers (that was proven in the probability course), claiming that

given a collection  of i.i.d. random variables, with

To visualize that convergence, we can use

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

It is possible to visualize also the  bounds (used in the central limit theorem to get a limiting non degenerated distribution)

> u=seq(0,21,by=.2)
> v=sqrt(u*10)
> lines(u,1.96/v,col="red")
> lines(u,-1.96/v,col="red")

Yesterday, we’ve been discussing properties of the empirical cumulative distribution function,

We’ve seen Glivenko-Cantelli theorem, which states that (under mild assumptions)

To visualize that convergence use the following code. Here I use the trick

to get the maximum (componentwise) between two matrices

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

We have also discussed the pointwise asymptotic normality of the empirical cumulative distribution function

Here again, it is possible to visualize it. The first step is to compute several trajectories for empirical cumulative distribution function

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

Note that we can compute (pointwise) confidence bands

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

Now, if we focus on one specific point, we can visualize the asmptotic normality (i.e. the almost normality when we have a sample of size 100)

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

Defining Properly MA(∞) Time Series

In order to properly define https://latex.codecogs.com/gif.latex?MA(\infty) series, we need to get back on some properties of infinite sequences, as briefly mentioned yesterday in the MAT8181 course. Consider some sequence https://latex.codecogs.com/gif.latex?(a_i)_{i\in\mathbb{N}}. The sequence is said to be summable if

https://latex.codecogs.com/gif.latex?S_n=\sum_{i=0}^n%20a_i

is convergent, i.e. if the limit of https://latex.codecogs.com/gif.latex?S_n exists when https://latex.codecogs.com/gif.latex?n\rightarrow\infty.

From Cauchy criterionhttps://latex.codecogs.com/gif.latex?\sum%20a_i converges if and only if for each https://latex.codecogs.com/gif.latex?\eta%3E0, there is https://latex.codecogs.com/gif.latex?m\in\mathbb{N} for which

https://latex.codecogs.com/gif.latex?\vert%20a_i+a_{i+1}+\cdots+a_{j-1}+a_j\vert%3C\eta

when https://latex.codecogs.com/gif.latex?%20i,j%3Em. The sequence https://latex.codecogs.com/gif.latex?(a_i)_{i\in\mathbb{N}} is said to be absolutely summable if

https://latex.codecogs.com/gif.latex?%20\sum_{i=0}^\infty%20\vert%20a_i\vert%20%3C\infty

and square-summable if

https://latex.codecogs.com/gif.latex?%20\sum_{i=0}^\infty%20a_i^2%20%3C\infty

Observe that absolute summability will imply square summability (since for https://latex.codecogs.com/gif.latex?j‘s large enough https://latex.codecogs.com/gif.latex?%20\vert%20a_j\vert%20\leq1, and then https://latex.codecogs.com/gif.latex?%20a_j^2\leq\vert%20a_j\vert)

Consider now some https://latex.codecogs.com/gif.latex?MA(\infty) time series

https://latex.codecogs.com/gif.latex?%20X_t=\sum_{h=0}^\infty%20\theta_h%20\varepsilon_{t-h}

If the sequence of coefficients https://latex.codecogs.com/gif.latex?%20(\theta_i) is square-summable, then

https://latex.codecogs.com/gif.latex?%20S_T%20=%20\sum_{h=0}^T%20\theta_h%20\varepsilon_{t-h}

converges in https://latex.codecogs.com/gif.latex?%20L_2  to some random varible as https://latex.codecogs.com/gif.latex?%20T\rightarrow\infty. This can be proved easily using Cauchy criteria, in the sense that for any https://latex.codecogs.com/gif.latex?\eta%3E0, there is a https://latex.codecogs.com/gif.latex?%20T large enough such that, for any https://latex.codecogs.com/gif.latex?%20h,

In that case, if the sequence of coefficients https://latex.codecogs.com/gif.latex?%20(\theta_i) is square-summable, then https://latex.codecogs.com/gif.latex?%20(X_t) is stationary (in the https://latex.codecogs.com/gif.latex?%20L_2 sense) since the process is centered, and

https://latex.codecogs.com/gif.latex?\gamma(h)=\sigma^2%20\cdot%20\sum_{i=0}^\infty%20\theta_i%20\theta_{i+h}

for all https://latex.codecogs.com/gif.latex?h\in\mathbb{N}.

Further, ergodicity of the time series, define as the absolute summability of the autocovariance sequence, is obtained when the sequence of coefficients https://latex.codecogs.com/gif.latex?%20(\theta_i) is absolutely summable.

Simple and heuristic optimization

This week, at the Rmetrics conference, there has been an interesting discussion about heuristic optimization. The starting point was simple: in complex optimization problems (here we mean with a lot of local maxima, for instance), we do not necessarily need extremely advanced algorithms that do converge extremly fast, if we cannot ensure that they reach the optimum. Converging extremely fast, with a great numerical precision to some point (that is not the point we’re looking for) is useless. And some algorithms might be much slower, but at least, it is much more likely to converge to the optimum. Wherever we start from.
We have experienced that with Mathieu, while we were looking for maximum likelihood of our MINAR process: genetic algorithm have performed extremely well. The idea is extremly simple, and natural. Let us consider as a starting point the following algorithm,

  1. Start from some 
  2. At step , draw a point  in a neighborhood of 
  • either  then 
  • or  then 

This is simple (if you do not enter into details about what such a neighborhood should be). But using that kind of algorithm, you might get trapped and attracted to some local optima if the neighborhood is not large enough. An alternative to this technique is the following: it might be interesting to change a bit more, and instead of changing when we have a maximum, we change if we have almost a maximum. Namely at step ,

  • either then 
  • or  then 

for some . To illustrate the idea, consider the following function

> f=function(x,y) { r <- sqrt(x^2+y^2);
+ 1.1^(x+y)*10 * sin(r)/r }
(on some bounded support). Here, by picking noise and  values arbitrary, we have obtained the following scenarios
> x0=15
> MX=matrix(NA,501,2)
> MX[1,]=runif(2,-x0,x0)
> k=.5
> for(s in 2:501){
+  bruit=rnorm(2)
+  X=MX[s-1,]+bruit*3
+  if(X[1]>x0){X[1]=x0}
+  if(X[1]<(-x0)){X[1]=-x0}
+  if(X[2]>x0){X[2]=x0}
+  if(X[2]<(-x0)){X[2]=-x0}
+  if(f(X[1],X[2])+k>f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=X}
+  if(f(X[1],X[2])+k<=f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=MX[s-1,]}
+}

It does not always converge towards the optimum,

and sometimes, we just missed it after being extremely unlucky

Note that if we run 10,000 scenarios (with different random noises and starting point), in 50% scenarios, we reach the maxima. Or at least, we are next to it, on top.

What if we compare with a standard optimization routine, like Nelder-Mead, or quasi gradient ?Since we look for the maxima on a restricted domain, we can use the following function,

> g=function(x) f(x[1],x[2])
> optim(X0, g,method="L-BFGS-B",
+ lower=-c(x0,x0),upper=c(x0,x0))$par

In that case, if we run the algorithm with 10,000 random starting point, this is where we end, below on the right (while the heuristic technique is on the left),

In only 15% of the scenarios, we have been able to reach the region where the maximum is.

So here, it looks like an heuristic method works extremelly well, if do not need to reach the maxima with a great precision. Which is usually the case actually.

L’espérance est un opérateur linéaire, so what ?

Oui, “l’espérance est un opérateur linéaire”. On n’arrête pas d’insister sur cette propriété dans la plupart des cours de probabilité. En fait, je pense que cette propriété a deux implications importantes. La première est que l’on interprète souvent cette phrase en disant que si http://freakonometrics.blog.free.fr/public/perso5/espl01.gif intégrable, alors pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gifhttp://freakonometrics.blog.free.fr/public/perso5/espl7.gif. La conséquence de cette propriété (qui est juste) est que http://freakonometrics.blog.free.fr/public/perso5/espl09.gif dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).

Mais cette linéarité dit davantage. Elle dit que pour toutes variables aléatoires http://freakonometrics.blog.free.fr/public/perso5/espl01.gifet http://freakonometrics.blog.free.fr/public/perso5/espl2.gif (définies sur le même espace), http://freakonometrics.blog.free.fr/public/perso5/espl3.gif, pour tout http://freakonometrics.blog.free.fr/public/perso5/espl4.gif et http://freakonometrics.blog.free.fr/public/perso5/espl6.gif. Et ça c’est fort. En particulier, on ne dit pas qu’il faut que http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre http://freakonometrics.blog.free.fr/public/perso5/espl01.gif et http://freakonometrics.blog.free.fr/public/perso5/espl2.gif.

Pourtant @Arnaud (qui passe souvent à l’occasion sur le blog) m’a laissé un couple de questions par courriel: “la dépendance modifie-t-elle la moyenne ?[…] sur des exemples simples en R, la moyenne est toujours perturbée”. L’idée était que les espérances étaient approchées numériquement par simulation. Or les méthodes de Monte Carlo reposent sur la loi des grands nombres et sur le théorème central limite. En particulier, si http://freakonometrics.blog.free.fr/public/perso5/espl34.gif, et que l’on simule les couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif, alors

http://freakonometrics.blog.free.fr/public/perso5/espl31.gif

autrement dit

http://freakonometrics.blog.free.fr/public/perso5/espl32.gif

i.e. la vitesse de convergence est affectée par la dépendance, pas la valeur vers laquelle on va tendre (car la moyenne tend vers l’espérance, qui est linéaire). En particulier, la convergence sera d’autant plus lente que http://freakonometrics.blog.free.fr/public/perso5/espl33.gif sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.

Par exemple, si on s’amuse à simuler 100 couples http://freakonometrics.blog.free.fr/public/perso5/espl10.gif Gaussiens, avec

http://freakonometrics.blog.free.fr/public/perso5/espl11.gif

on a les valeurs suivantes pour la moyenne de ce couple (obtenu en générant 10,000 échantillons de 100 paires), avec en abscisse la corrélation entre les deuxvariables

le trait rouge est la moyenne empirique des moyennes obtenues, et les régions sont délimitées par les quantiles à 5% (et 95%), 10% (et 90%) et 25% (et 75%). En moyenne, les simulations donnent la même chose, peu importe la dépendance. Mais plus la dépendance est forte, plus la variance sera grande, et plus la convergence sera lente (et les simulations auront d’autant plus tendance à s’éloigner – ou se disperser autour – de la vraie valeur). Mais cette réponse est partielle, car on a supposé que la variance était finie. Or cette hypothèse n’est pas nécessaire pour assurer la convergence de la moyenne vers l’espérance (et donc pour utiliser des méthodes de simulations). Que se passerait-il avec des lois de variance infinie ?

Pour autre chose que les vecteurs Gaussiens (par exemple des variables de variance infinie, comme des lois de Pareto) on peut utiliser le code suivant (on continue pour l’instant à utiliser une copule Gaussienne),

library(mnormt)
library(copula)
library(evir)
set.seed(1)
ns=100
R=seq(-.95,.95,by=.05)
M=BINF1=BSUP1=rep(NA,length(R))
for(i in 1:length(R)){
r=R[i]
norm.cop = normalCopula(r, dim = 2)
S=rep(NA,20000)
for(j in 1:20000){
U=rcopula(norm.cop,ns)
X=cbind(qgpd(U[,1], .7 ,1),qgpd(U[,2], .5 ,1))
S[j]=mean(X[,1]+X[,2])}
BINF1[i]=quantile(S,.05)
BSUP1[i]=quantile(S,.95)
M[i]=mean(S)
}
plot(R,M,col="white",ylim=range(c(BINF1,BSUP1)))
polygon(c(R,rev(R)),c(BINF1,rev(BSUP1)),
col="light yellow",border=NA)
lines(R,BINF1,lty=2)
lines(R,BSUP1,lty=2)
lines(R,M,lwd=2,col="red")
esperance=1+1/(1-.7)+1+1/(1-.5)
segments(min(R),esperance,max(R),esperance,col="blue")

Sur cette somme de deux lois de Pareto (comme dans le code ci-dessus), on obtient numériquement

Autrement dit, ici, l’effet dépendance semble avoir un impact plus faible sur l’approximation de la moyenne. On peut aussi essayer de changer la copule, par exemple mettre une copule de type Clayton (avec en abscisse cette fois le tau de Kendall),

(qui présente de la dépendance dans la queue inférieure) ou la copule duale (qui présente de la dépendance dans la queue supérieure)

Bref, la structure de la dépendance ne semble pas trop changer l’estimation de la moyenne (ou plutôt la vitesse de la convergence). On peut aussi se demander s’il y a un effet dimension. Par exemple, si au lieu de simuler des couples, on simule des vecteurs de plus grande taille. Par exemple, un vecteur échangeable en dimension 10, avec des marges normales centrées réduites,

norm.cop =normalCopula(r, dim = 10,dispstr="ex")

Moralité ? effectivement, la dépendance a un impact sur la vitesse de convergence quand on travaille sur des vecteurs aléatoires. Mais de manière assez surprenante, c’est surtout le cas pour la loi normale. Et étrangement, cet impact ne semble pas lié à de la dépendance extrême, ou à la dimension…

Convergences de variables aléatoires

Le schéma – ou diagramme – des relations entre les modes convergences ressemble à ce qui est présenté sur la droite. C’est grosso modo ce que j’ai dit en cours (sans le représenter explicitement).

Je ne vais pas démontrer ces implications car ce sont des théorèmes (qui peuvent donc se trouver dans n’importe quel textbook de probabilité).

En revanche, un étudiant m’avait demandé si les réciproques étaient vraies. Comme toujours en maths, si les réciproques étaient toujours vraies, on le mettrait. Donc ici les réciproques sont généralement fausses. On peut toutefois se demander quand ces réciproques peuvent être vraies…et quand elles sont fausses.

  • de la convergence en proba à la convergence en moyenne

Commençons par le plus simple (à mon avis). Il existe un cas particulier où cette réciproque est vraie: le cas de variables aléatoires bornées. Autrement dit, si http://freakonometrics.free.fr/blog/cvvvv20.png et que l’on a http://freakonometrics.free.fr/blog/cvvvvv22.png pour tout http://freakonometrics.free.fr/blog/lgn10.png, alors http://freakonometrics.free.fr/blog/cvvvvvv23.png pour tout http://freakonometrics.free.fr/blog/cvvvvv24.png.
Bon, en fait la condition est un peu trop forte, et il suffit que la suite http://freakonometrics.free.fr/blog/cvvv1.png soit uniformément intégrable.
Maintenant, il existe aussi des cas où cette réciproque sera fausse (et on se doute qu’en construisant une suite non bornée, on devrait y arriver). Par exemple, supposons que pour tout http://freakonometrics.free.fr/blog/lgn10.png,
http://freakonometrics.free.fr/blog/cvvv5.png
On peut montrer que http://freakonometrics.free.fr/blog/cvvvvv7.png, car en effet
http://freakonometrics.free.fr/blog/cvvvv6.png
Toutefois, certes la suite converge dans L1, mais vers une limite différente, étant donné que http://freakonometrics.free.fr/blog/cvvvvv9.png pour tout http://freakonometrics.free.fr/blog/lgn10.png. De plus, http://freakonometrics.free.fr/blog/cvvvv10.png et donc la suite ne converge pas dans L2.

  • de la convergence en loi à la convergence en proba

Un exemple simple ici: si http://freakonometrics.free.fr/blog/cvvv1.png converge en loi vers une constante, alors on aura convergence en probabilité (vers la même constante).
Maintenant la réciproque est relativement simple, l’idée étant que http://freakonometrics.blog.free.fr/public/perso/cvvv555.png n’implique en rien http://freakonometrics.blog.free.fr/public/perso/cvvvvv5666.png, au sens où http://freakonometrics.blog.free.fr/public/perso/cvvvvvv56666.png, pour tout http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png. Supposons maintenant que l’on ait à la fois http://freakonometrics.blog.free.fr/public/perso/cvvvv4646.png et http://freakonometrics.free.fr/blog/cvvvv20.png, où on va supposer que http://freakonometrics.free.fr/wesp30.png est une loi symétrique. La symétrie signifie que http://freakonometrics.blog.free.fr/public/perso/cvvv5656655.png. Notons que l’on a alors forcément http://freakonometrics.blog.free.fr/public/perso/cvvvcv_56.png. Supposons que la réciproque soit toujours vraie, alors on devrait avoir http://freakonometrics.blog.free.fr/public/perso/cv65464565645vv.png, et donc on aurait que http://freakonometrics.blog.free.fr/public/perso/cvwgfrwegr656.png, ce qui n’arrive que si http://freakonometrics.free.fr/wesp30.png est la constante http://freakonometrics.free.fr/blog/cvvvv4.png.

  • de la convergence en proba vers la convergence presque sûre

Je ne vais pas m’étendre sur la construction d’exemples pour lesquels la réciproque est fausse (formellement, il suffit de changer d’espace de probabilité sous-jacent, je peux renvoyer vers les notes de cours de Philippe Briand (ici) et Hélène m’a suggéré également les exercices en ligne ). Par contre, il existe des cas importants où la réciproque est vraie.
Considérons une suite de variables indépendantes et identiquement distribuées http://freakonometrics.free.fr/blog/cvvv1.png. Posons alors http://freakonometrics.free.fr/blog/cvvvv30.png On peut alors montrer que http://freakonometrics.free.fr/blog/cvvvvvv26.png si et seulement si http://freakonometrics.free.fr/blog/cvvvvv35.png (ça doit être du Borel-Cantelli).

  • de la convergence en loi à la convergence presque sûre

Je ne vais pas non plus m’étendre davantage, mais on a un résultat qui ressemble à ce que je disais auparavant. En fait, il existe un théorème (dit de Skorohod) qui permet de dire que si on a une convergence en loi, alors on peut construite une suite de variables aléatoires qui
convergera presque sûrement. Formellement on a que si http://freakonometrics.blog.free.fr/public/perso/cvvvv4646.png, alors il existe une suite http://freakonometrics.free.fr/blog/cvvv11.png et une variable http://freakonometrics.free.fr/blog/cvvv13.png, telles que http://freakonometrics.free.fr/blog/cvvv16.png et http://freakonometrics.free.fr/blog/cvvvvv15.png telles que http://freakonometrics.free.fr/blog/cvvvv17.png. Mais de manière générale, si la suite n’est pas construite comme il faut (http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png par http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png en quelques sortes) alors il est assez rare d’avoir la réciproque.

Échantillonnage, TCL et LGN

Mercredi matin, nous allons aborder l’échantillonnage. Nous allons partir d’un échantillon http://freakonometrics.free.fr/blog/lgn01.png, supposé indépendant et identiquement distribué. Le modèle sous-jacent est la collection de variables aléatoires http://freakonometrics.free.fr/blog/lgn02.pngindépendantes et identiquement distribuées. Dans un premier temps, on s’intéressera à la moyenne de l’échantillon, correspondant à une réalisation de la variable aléatoire http://freakonometrics.free.fr/blog/lgn04.png définie par

http://freakonometrics.free.fr/blog/lgn05.png

Afin d’avoir des résultats sur le comportement de la moyenne, on rappellera deux théorèmes importants: la loi des grands nombres et le théorème central limite. Afin d’illustrer ces théorèmes, considérons le cas de tirages de PILE et FACE. On suppose que http://freakonometrics.free.fr/blog/lgn06.png désigne la variable aléatoire correspondant au fait d’avoir obtenu PILE au http://freakonometrics.free.fr/blog/lgn07.pngème tirage. On remet des pièces à une centaine d’élèves (chaque élève correspond à une couleur différente), et on regarde la moyenne du nombre de PILE obtenus sur les http://freakonometrics.free.fr/blog/lgn10.png premiers lancers. Et on regarde ce qui se passe quand n augmente. La loi de grands nombres permet d’affirmer que

http://freakonometrics.free.fr/blog/lgn11.png

ici, si l’on suppose que la pièce n’est pas biaisée, soit

http://freakonometrics.free.fr/blog/lgn12.png

pour la version faible, ou pour la version forte

http://freakonometrics.free.fr/blog/lgn13.png

soit

http://freakonometrics.free.fr/blog/lgn16.png

C’est ce que l’on observe ci-dessous, chaque trait correspondant à la trajectoire http://freakonometrics.free.fr/blog/lgn21.png obtenue par un élève. Un histogramme est représenté à droite.

http://freakonometrics.free.fr/animation-cv-binomiale.gif

On peut aussi centrer (ici retrancher http://freakonometrics.free.fr/blog/lgn22.png) et réduire (ici mulplier par http://freakonometrics.free.fr/blog/lgn23.png) la moyenne empirique pour étudier une éventuelle convergence. On pose alors

http://freakonometrics.free.fr/blog/lgn24.png

Le théorème central limite garantie que

http://freakonometrics.free.fr/blog/lgn27.png

On a toujours l’histogramme des valeurs obtenues, au bout de http://freakonometrics.free.fr/blog/lgn10.png tirages, par les élèves, avec en rouge la densité de la loi normale centrée réduire. Et là encore, l’expérience semble valider la théorie…

http://freakonometrics.free.fr/animation-cv-CR.gif

Maintenant, pour ceux qui doutent, la base est en ligne ici, et je mets le code pour commencer à jouer avec…

M=read.table("http://freakonometrics.blog.free.fr/public/data/PILEFACE.txt")
nl=ncol(M)
nj=nrow(M)
M1=matrix(M=="PILE",nj,nl)
MS=M1
for(l in 2:nl){MS[,l]=
(apply(M1[,1:l],1,mean)-.5)/sqrt(.25/l)}
arret=100
COULEUR= rev(rainbow(nj))
HIST=hist(MS[,arret])
plot(1:nl,MS[1,],ylim=(c(-3,3)),col="white",
xlim=c(0,250),xlab="Nombre de lancers de pièces",
ylab="Nombre moyen de PILE, centré réduit")
for(j in 1:nj){
lines(1:arret,MS[j,1:arret],col=COULEUR[j])
}
abline(v=arret)

Statistiques, rappels de probas

Avant de commencer les travaux pratiques de vendredi, quelques rappels sur les notions de convergence.

  • Convergence en loi

Soit http://freakonometrics.free.fr/blog/latex/conver15.png une suite de variables aléatoires, de telle sorte que http://freakonometrics.free.fr/blog/latex/conver62.png soit la fonction de répartition de http://freakonometrics.blog.free.fr/public/maths/conv99.png. On dira que la suite http://freakonometrics.free.fr/blog/latex/conver15.png converge en loi vers la variable aléatoire http://freakonometrics.blog.free.fr/public/maths/montone-deriv1.png de loi http://freakonometrics.blog.free.fr/public/maths/montone2.png si

http://freakonometrics.free.fr/blog/latex/conver6%c3%b4.png

en tout point de continuité de http://freakonometrics.blog.free.fr/public/maths/montone2.png. On notera

http://freakonometrics.free.fr/blog/latex/conver6.png

La convergence en loi est la forme la plus faible de convergence: en général, elle n’implique pas les autres formes de convergence que nous verrons, alors que ces autres formes de convergence impliquent la convergence en loi.

De manière équivalente, on aura convergence en loi si et seulement si, pour toute fonction bornée

http://freakonometrics.free.fr/blog/latex/conver54.png

En particulier, on aura convergence en loi si et seulement si on a convergence des fonctions génératrice en tout point.

  • Convergence en probabilité

La suite http://freakonometrics.free.fr/blog/latex/conver15.png converge en probabilité vers la variable aléatoire http://freakonometrics.blog.free.fr/public/maths/montone-deriv1.png si

http://freakonometrics.free.fr/blog/latex/conv50.png

que l’on notera

http://freakonometrics.free.fr/blog/latex/conver4.png

Parmi les résultats classiques, on a aussi que si

http://freakonometrics.free.fr/blog/latex/conver40.png

et

http://freakonometrics.free.fr/blog/latex/conver31.png

alors http://freakonometrics.free.fr/blog/latex/conver30.png. C’est ce mode de convergence qui intervient dans la loi faible des grands nombres.
si http://freakonometrics.free.fr/blog/latex/conver26.png est une suite de variables aléatoires indépendantes et identiquement distribuées d’espérance http://freakonometrics.free.fr/blog/latex/conver27.png et de variance http://freakonometrics.free.fr/blog/latex/conver28.png et que l’on considére la moyenne empirique

http://freakonometrics.free.fr/blog/latex/conver24.png

alors  http://freakonometrics.free.fr/blog/latex/conver22.png et

http://freakonometrics.free.fr/blog/latex/conver21.png

Donc http://freakonometrics.blog.free.fr/public/maths/conv98.png.

  • Convergence presque sure

La suite http://freakonometrics.free.fr/blog/latex/conver15.png converge en probabilité vers la variable aléatoire http://freakonometrics.blog.free.fr/public/maths/montone-deriv1.png si

http://freakonometrics.free.fr/blog/latex/conver13.png

ou encore

http://freakonometrics.free.fr/blog/latex/conver12.png

que l’on notera

http://freakonometrics.free.fr/blog/latex/conver2.png

La convergence presque sûre implique une convergence en probabilité.
La convergence presque sûre est utilisée dans la loi forte des grands nombres.

  • Convergences L1 et L2

La suite http://freakonometrics.free.fr/blog/latex/conver15.png converge en norme L1 vers la variable aléatoire X si

http://freakonometrics.free.fr/blog/latex/conver11.png

et on parlera de convergence en moyenne quadratique, ou L2, si

http://freakonometrics.free.fr/blog/latex/conver10.png
  • Convergences et transformations

Comme nous le (re)verrons vendredi matin, la transformation des variables aléatoires est un point fondamental en statistique. Soit http://freakonometrics.free.fr/blog/latex/conver7.png une fonction continue en tout point d’un ensemble http://freakonometrics.free.fr/blog/latex/conver9.png tel que http://freakonometrics.free.fr/blog/latex/conver8.png, alors
si http://freakonometrics.free.fr/blog/latex/conver6.png alors http://freakonometrics.free.fr/blog/latex/conver5.png
si http://freakonometrics.free.fr<br />
/blog/latex/conver4.png alors http://freakonometrics.free.fr/blog/latex/conver3.png
si http://freakonometrics.free.fr/blog/latex/conver2.png alors http://freakonometrics.free.fr/blog/latex/converg1.png
Pour aller plus loin sur les notions de convergence, je renvois au livre de Adrianus Willem van der Vaart, Asymptotic Statistics, paru chez Cambridge University Press.

  • Transformations

Sinon sur les transformations, je devrais présenter (en cours) la formule classique de transformation d’une variable à densité par une fonction bijective et dérivable, http://freakonometrics.blog.free.fr/public/perso/changevar999.png. La formule est de la forme

http://freakonometrics.blog.free.fr/public/perso/changementvar1.png

En fait, si on se souvient que http://freakonometrics.blog.free.fr/public/perso/changevar3.png, ou http://freakonometrics.blog.free.fr/public/perso/changevar4.png, on peut aussi écrire

http://freakonometrics.blog.free.fr/public/perso/changevar2.png

L’intérêt de cette relation est qu’elle peut se généraliser en dimension quelconque. Si http://freakonometrics.blog.free.fr/public/perso/changevar99.png (toujours avec une transformation bijective possédant un Jacobien http://freakonometrics.blog.free.fr/public/perso/changecar8.png strictement positif), on a

http://freakonometrics.blog.free.fr/public/perso/changevar5.png

avec http://freakonometrics.blog.free.fr/public/perso/changevar7.png. Je peux renvoyer ici pour quelques compléments (et exemples). Mais nous reviendrons plus tard sur le cas multivarié…

Propriétés des estimateurs dans une régression

Comme les rappels (je devrais plutôt dire “remise à niveau“) ont été particulièrement rapides, je vais prendre un peu de temps pour revenir sur quelques éléments du cours d’économétrie (d’autant plus que j’ai reçu quelques questions par mail).

  • sur le calcul des estimateurs

Je vais reprendre les questions que j’ai reçu, ça sera plus simple, “je comprend pas très bien quelle méthode utilise l’opération lm(Y~X…..) sur R pour trouver les coef” ou encore “je regresse un modèle simple, je sors le summary, et j’obtiens un betaREG 3 étoiles. Toujours sur R, je compare ce betaREG avec le betaCalc = (X’X)-1 * X’Y, il est différent ! Comment expliquer cela ? Est-ce possible ?“.
Commençons par la première question: la fonction lm() utilise la méthode dite des moindres carrés, ce qui revient à calculer

en adoptant une écriture matricielle, i.e. la solution est alors

Cet estimateur est celui qui minimise la somme des carrés des erreurs. Et on peut vérifier numériquement que cet estimateur est bien celui calculé par R,

>data(cars)
>X=cbind(rep(1,nrow(cars)),cars$speed); Y=cars$dist
> X[1:5,]
      [,1] [,2]
 [1,]    1    4
 [2,]    1    4
 [3,]    1    7
 [4,]    1    7
 [5,]    1    8
> solve(t(X)%*%X)
            [,1]         [,2]
[1,]  0.19310949 -0.011240876
[2,] -0.01124088  0.000729927
> t(X)%*%Y
      [,1]
[1,]  2149
[2,] 38482
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
           [,1]
[1,] -17.579095
[2,]   3.932409

ce qui correspond très précisément à ce que calcule la fonction de R

> lm(dist~speed,cars)
Coefficients:
(Intercept)        speed 
    -17.579        3.932

Et on peut aussi vérifier qu’il minimise bien la somme des carrés des erreurs,

> b0=seq(-30,10,by=3);
> b1=seq(1,7,by=.5)
> SC=matrix(NA,length(b0),length(b1))
> for(i in 1:length(b0)){
> for(j in 1:length(b1)){
 + SC[i,j]=sum((cars$dist-(b0[i]+b1[j]*cars$speed))^2)
}}
> contour(b0,b1,SC)

et je peux certifier que la valeur obtenue est bien le minimum de cette fonction (en fait c’est un résultat théorique d’optimisation).
Ca se généralise bien entendu en dimension plus grande. Alors attention, cette formule permet de calculer un estimateur du vecteur des paramètres, mais la formule est bien sûre fausse si on raisonne composante par composante.

> X=cars$speed; Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
         [,1]
[1,] 2.909132

Mais ça c’est un résultat (trivial) d’algèbre linéaire. En effet, le fait d’avoir

ne signifie pas que l’on ait ce genre de formule composante par composante, i.e.

Dit de manière un peu formelle, le produit matriciel n’est pas un produit terme à terme. La traduction économétrique de cette idée est qu’une régression multiple n’est pas une succession de régressions simples. En effet, dans l’équation précédante, le dernier terme correspond à la régression sans la constante, comme on peut le voir ci-dessous

> lm(dist~0+speed,cars)
Coefficients:
speed  
2.909

et le premier est tout simplement la moyenne empirique,

> lm(dist~1,cars)
Coefficients:
(Intercept)  
      42.98  
> mean(cars$dist)
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
      [,1]
[1,] 42.98
  • sur les propriétés des estimateurs

Une question était “doit-on vérifier les propriétés de l’estimateur (sans biais, efficace etc…) du summary de la Reg ? pour chaque régression ?
Bon, la réponse est non car c’est impossible ! Pour répéter tranquillement ce que j’avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d’hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l’estimateur par moindre carré est BLUE (best linear unbiased estimator).
Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d’être sans biais, par exemple.
On ne
peut pas
vérifier qu’un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l’absence de biais.
Mais comme on va le voir par la suite, le biais ce n’est pas forcément gênant, ce qui est important, c’est surtout la convergence. Comme le notait Clive Granger  “if you can’t even get a consistant estimator, you shouldn’t be in this business“.

  • et sur des points plus théoriques

sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un “rapport avec les plim X’E/n et plim (X’X)-1/ n“.
Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c’est à dire non corrélé avec la variable explicative, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-1.png

En effet, sinon

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-2.png

autrement dit l’estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce “problème” est d’utiliser un instrument (ou une variable instrumentale).
Mais surtout, sur cet exemple, on voit qu’on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement.
Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur).

  • un estimateur sans biais est généralement convergent (et de manière consistante)
  • un estimateur biaisé peut ne pas converger

Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l’ai expliqué auparavant, la distribution de l’estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d’aujourd’hui). La figure ci-dessous correspond au cas sans biais, et convergent

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-sans-biais-convergent_m.jpg

La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l’estimateur sera sans biais.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-convergent_m.jpg

Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l’estimateur sera toujours biaisé,
https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-non-convergent_m.jpg

Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-3.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-4.png

(c’est l’hypothèse d’homoscédasticité et d’indépendance des résidus). Bon, il faut aussi l’absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité.
On peut s’intéresser aux propriétés asymptotiques de l’estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-5.png converge en probabilité vers https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-6.png, parfois noté

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-7.png

si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-8.png

Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-9.png

si l’espérance des variables est finie, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-10.png

Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-11.png

On parlera alors d’estimateur convergent. On peut aussi montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-12.png

Mais formellement,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-13.png

si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-14.png

Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d’économétrie 1 qui commencera très bientôt… et surtout d’économétrie 2 qui abordera les aspects dynamiques.

Delta method avec R, ou transformation des paramètres

Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme

que l’on peut encore écrire

Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance.
La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses)

ce qui donne avec les notations précédentes

Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de

 

(ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt).

  • En faisant un peu de calcul formel

Pour les calculs symboliques, il y a la fonction D dans le package de base,

> f <- expression(sin(cos(x + y^2)))
> ( Dfx <- D(f, "x") )
-(cos(cos(x + y^2)) * sin(x + y^2))

ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique.

> x=c(0,1,2); y=1
> eval(D.sc)
[1] -0.7216061 -0.8316919 -0.0774320

On peut également utiliser cette fonction pour faire plein de jolis dessins,

> x <-seq(-5, 5, length=100)
> fx=x^2
> Df <-D(expression(x^2), "x")
> Dfx<-eval(Df)
> plot(x, fx, type = "l", las = 1,ylim=c(-2,2))
> for (i in 1:length(x)) {
+ abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i],
+ col = rainbow(100)[i])
+ }

Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse.

> reg1 = lm(dist~speed, data = cars)
> estmean = coef(reg1)
> estvar = summary(reg1)$cov.unscaled * summary(reg1)$sigma^2
> f = function(x,y){
+ -y/(x+y)
+ }
> X=f(estmean[1],estmean[2])
> g=expression(-y/(x+y))
> x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2])
> V=eval(c(eval(D(g, "x")),eval(D(g, "y"))))
> sqrt(t(V) %*% estvar %*% V)
          [,1]
[1,] 0.1063255

Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function.

  • En utilisant la fonction deltamethod de la librarie msm

Ici on commence comme avant, sauf que la phase de calcul est immédiat…

> reg2 <- lm(dist~speed, data = cars)
> estmean = coef(reg2)
> estvar = summary(reg2)$cov.unscaled * summary(reg2)$sigma^2
> deltamethod (~ -x2 / (x1 + x2), estmean, estvar) 
[1] 0.1063255

et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2

> deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar) 
[1] 0.1063255 6.3664368

En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique).

  • En utilisant la fonction delta.method de librarie alr3

Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante,

> reg3 = lm(dist~speed, data = cars)
> delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255

En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95%

> DM=delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255 
> DM$estimate+qt(.975,df=nrow(cars)-1)*DM$se
[1] 0.5018277
> DM$estimate+qt(.025,df=nrow(cars)-1)*DM$se
[1] 0.07448936

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].