Loi multinomiale et loi du chi-deux

La semaine passée, en cours, j’avais rappelé que quand https://latex.codecogs.com/gif.latex?\boldsymbol{N}=(N_{1},\cdots,N_{k}) décrivait le compte de https://latex.codecogs.com/gif.latex?n variable multinomiales prenant https://latex.codecogs.com/gif.latex?k modalités, la variable

https://latex.codecogs.com/gif.latex?\tilde{X}_i%20=%20\frac{N_i-np_i}{\sqrt{np_i\color{red}{(1-p_i)}}}

suit asymptotiquement une loi https://latex.codecogs.com/gif.latex?\mathcal{N}(0,1). Et plus généralement, on peut montrer que https://latex.codecogs.com/gif.latex?\tilde{\boldsymbol{X}}\sim\mathcal{N}(\boldsymbol{0},\star). Le soucis est que la matrice de variance covariance n’est pas la matrice identité. Pire que ça, elle n’est pas diagonale. Encore pire, elle n’est pas inversible.

On ne peut alors pas utiliser le joli résultat qui nous dit que si https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}), alors https://latex.codecogs.com/gif.latex?Y=(\boldsymbol{X}-\boldsymbol{\mu})^{\text{\sffamily%20T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{X}-\boldsymbol{\mu}) suit une loi du chi-deux à https://latex.codecogs.com/gif.latex?k degrés de liberté (dimension du vecteur https://latex.codecogs.com/gif.latex?\boldsymbol{X}). Car il faut que https://latex.codecogs.com/gif.latex?\boldsymbol{\Sigma} soit inversible. Or ici, je pense qu’on peut montrer que

https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=\prod_{i=1}^k%20p_i%20\cdot\left(\sum_{i=1}^k%20p_i%20-1\right)

(à un facteur https://latex.codecogs.com/gif.latex?n près). Et comme on a une loi de probabilité, le terme de droite est nul. Bref, https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=0 et son inverse n’existe alors pas. Cela dit, https://latex.codecogs.com/gif.latex?0 n’est valeur propre qu’une fois, d’où la dimension https://latex.codecogs.com/gif.latex?k-1 qui apparaît sans cesse (à condition que les https://latex.codecogs.com/gif.latex?p_i soient non-nuls, ce qui était une hypothèse dans la définition de la loi multinomiale : toutes les modalités existent vraiment).

L’idée est alors de poser

https://latex.codecogs.com/gif.latex?X_i%20=%20\frac{N_i-np_i}{\sqrt{np_i}}

suit asymptotiquement une loi https://latex.codecogs.com/gif.latex?\mathcal{N}(0,\star). Et plus généralement, là encore, https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{0},\star).

Ce qui est amusant, c’est que dans le cas binomiale (j’avais dit qu’il était intéressant de voir la loi multinomiale comme une généralisation de la loi binomiale),

https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=X_1^2+X_2^2

puisque

https://latex.codecogs.com/gif.latex?\frac{(N_1-np_1)^2}{np_1(1-p_1)}=\frac{(N_1-np_1)^2}{np_1}+\frac{(N_2-np_2)^2}{np_2}

Et en fait, un peut de calculs rapides permettent de noter que

https://latex.codecogs.com/gif.latex?(N_1-np_1)^2=(N_2-np_2)^2

de telle sorte que https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=\tilde{X}_2^2 alors que https://latex.codecogs.com/gif.latex?{X}_1^2\neq%20{X}_2^2. Bref, c’est clairement sur cette dernière quantité qu’il faut travailler !

Et là on fait des calculs.

L’astuce est d’utiliser la formule de Stirling

https://latex.codecogs.com/gif.latex?n!\sim%20n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}

alors pour https://latex.codecogs.com/gif.latex?n_1+\cdots+n_k=n

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=\frac{n!}{n_1%20!%20\cdots%20n_k%20!}p_1^{n_1}\cdots%20p_k^{n_k}soit

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim\frac{n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}}{\prod_{j=1}^k%20[n_j^{n_j+\frac{1}{2}}e^{-n_j}\sqrt{2\pi}]}p_1^{n_1}\cdots%20p_k^{n_k}

et donc

https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim%20\kappa%20\prod_{j=1}^k%20\left(\frac{np_j}{n_j}\right)^{n_j+\frac{1}{2}}

de telle sorte que

https://latex.codecogs.com/gif.latex?\log%20\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=k+\sum_{j=1}^k%20\left(n_j+\frac{1}{2}\right)%20\log%20%20\left(\frac{np_j}{n_j}\right)On peut alors faire des développements limités, en faisant intervenir

https://latex.codecogs.com/gif.latex?x_i=\frac{n_i-np_i}{\sqrt{np_i}}

et non pas la version avec https://latex.codecogs.com/gif.latex?\sqrt{np_i(1-p_i)} au dénominateur. On peut alors réécrire la somme en faisant intervenir https://latex.codecogs.com/gif.latex?k-1 termes (et non pas https://latex.codecogs.com/gif.latex?k), exactement comme dans les calculs faits dans le premier cours, avec les calculs de variances. Bref, un peu de calcul supplémentaire permet de montrer que

https://latex.codecogs.com/gif.latex?\sum_{i=1}^k%20X_i^2=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}\overset{\mathcal{L}}{\rightarrow}%20\chi^2(k-1)

C’est donc

https://latex.codecogs.com/gif.latex?Q=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}

et non pas

https://latex.codecogs.com/gif.latex?\tilde{Q}=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i\color{red}{(1-p_i)}}

qui sera utilisée comme statistique de test par la suite.

Pour utiliser le test du chi-deux comme test d’ajustement de loi, je pourrais reprendre l’exemple des buts lors des coupes du monde de Football

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
soccer1998=read.table("http://freakonometrics.free.fr/soccer1998")
S98=(soccer1998$V1+soccer1998$V2)
soccer2010=read.table("http://freakonometrics.free.fr/soccer2010")
S10=(soccer2010$V1+soccer2010$V3)

La Science Economique, et son maudit Nobel

Je voulais revenir rapidement sur le prix Nobel décerné il y a moins de 48 heures à Jean Tirole, suite à plusieurs billets que j’ai pu lire, ici et là sur internet. Il y a quelques années, dans un chapitre de le modèle et le récit (publié en 2001), Pierre André Chiaporri racontait la petite anecdote suivante

quand on lui demandait pourquoi il faisait de l’économie “mathématique” plutôt que “littéraire”, Paul Samuelson répondait qu’il aurait souhaité pouvoir faire de l’économie littéraire – mais que c’était trop difficile pour lui. Ce sentiment est, je crois, partagé par beaucoup d’économistes de nos jours : de nombreux problèmes en économie théorique sont trop complexes pour être étudiés sans l’aide d’un minimum de formalisation.

(je recommande d’ailleurs la lecture de l’intégralité de l’article, qui discute l’importance de la modélisation en science économique, avec en toile de fond la théorie du consommateur, reprenant ainsi des idées que l’on trouvait déjà dans la théorie du consommateur est-elle réfutable ? ). En lisant ce paragraphe, je n’ai pu m’empêcher de penser au petit livre Economic Fables (publié par openbookpublishers). J’avais déjà évoqué ce livre dans un précédant billet. Pour ceux qui ne connaisse pas son auteur, Ariel Rubinstein (אריאל רובינשטיין) est un économiste, professeur à Tel Aviv et à NYU. Il a publié des livres de référence – dont une très bonne introduction à la théorie des jeux . Ce qui est amusant, ce que ce livre est un livre de fables, ce qui peut surprendre pour quelqu’un qui m’avait habitué à

Continue reading La Science Economique, et son maudit Nobel

Somewhere else, part 173

Some posts and articles worth reading

Continue reading Somewhere else, part 173

Somewhere else, part 172

Some interesting posts and articles, here and there

Why is so much writing so bad? Why is it so hard to understand a government form, or an academic article or the instructions for setting up a wireless home network? The most popular explanation is that opaque prose is a deliberate choice. Bureaucrats insist on gibberish to cover their anatomy. Plaid-clad tech writers get their revenge on the jocks who kicked sand in their faces and the girls who turned them down for dates. Pseudo-intellectuals spout obscure verbiage to hide the fact that they have nothing to say, hoping to bamboozle their audiences with highfalutin gobbledygook.  But the bamboozlement theory makes it too easy to demonize other people while letting ourselves off the hook. In explaining any human shortcoming, the first tool I reach for is Hanlon’s Razor: Never attribute to malice that which is adequately explained by stupidity. The kind of stupidity I have in mind has nothing to do with ignorance or low IQ; in fact, it’s often the brightest and best informed who suffer the most from it. I once attended a lecture on biology addressed to a large general audience at a conference on technology, entertainment and design. The lecture was also being filmed for distribution over the Internet to millions of other laypeople. The speaker was an eminent biologist who had been invited to explain his recent breakthrough in the structure of DNA. He launched into a jargon-packed technical presentation that was geared to his fellow molecular biologists, and it was immediately apparent to everyone in the room that none of them understood a word and he was wasting their time. Apparent to everyone, that is, except the eminent biologist. When the host interrupted and asked him to explain the work more clearly, he seemed genuinely surprised and not a little annoyed. This is the kind of stupidity I am talking about. []

Continue reading Somewhere else, part 172

Introduction aux tests nonparamétriques

Aujourd’hui, on va commencer la partie du cours sur les tests, avec une introduction par les tests non-paramétriques. On parlera (dans la continuité de ce qu’on a fait jusqu’à présent) de tests d’ajustement de lois, mais aussi de tests plus divers, comme des tests d’indépendance (entre observations successives).

Le test le plus classique est le test des signes. L’idée est de regarder puis de compter le nombre de fois où https://latex.codecogs.com/gif.latex?X_{i-1}\leq%20X_i (différence positive) et le nombre de fois où https://latex.codecogs.com/gif.latex?X_{i-1}\leq%20X_i (différence négative).

> (X=rnorm(20))
 [1]  0.81352495  0.78076709 -0.58654359  0.52805014 -0.31101282 -0.65022034 -0.93619729 -0.16483729 -1.97396753 -0.22506985  1.14959851 -0.38895448 -0.09298421
> (S=c("N","P")[1+1*(diff(X)>0)])
 [1] "N" "N" "P" "N" "N" "N" "P" "N" "P" "P" "N" "P" "P" "N" "P" "N" "P" "P" "N"

On verra ce test en cours. On peut par exemple compter le nombre de fois où la différence était positive,

> sim=function(a=0){
+ n=100
+ epsilon=X=rnorm(n)
+ for(t in 2:n) X[t]=a*X[t-1]+epsilon[t]
+ S=c("N","P")[1+1*(diff(X)>0)]
+ return(sum(S=="P"))}

(je fais ici une fonction un peu plus compliqué pour voir ce qui se passe si les données ne sont pas indépendantes, justement). Si on génère plein d’échantillons, on peut voir que le nombre de variations positives suit une loi normale, si les observations sont effectivement indépendantes,

> library(ggplot2)
> set.seed(1)
> vectS=unlist(lapply(rep(0,1000),sim))
> ggplot(data=data.frame(vectS = vectS), 
+ aes(x = vectS)) +
+ geom_histogram(aes(y = ..density..),
+ binwidth = 1, fill = "light blue", 
+ col = "black") +
+ stat_function(fun=function(x) dnorm(x,
+ (n-1)/2+.5,sqrt((n+1)/12)),
+ col = "red", size = 1) + ylab("Densite") +
+ ggtitle("Histogramme")

Continue reading Introduction aux tests nonparamétriques

Somewhere else, part 171

Some posts are articles worth reading

Continue reading Somewhere else, part 171

Exposés de Méthodologie Statistique (M2)

Tout au long de l’année, je vais animer un cours d’exposés de méthodologie statistique, pour les étudiants de dernière année de Master. L’objectif de ce travail est d’apprendre à présenter des outils ou techniques statistiques avancées. Je reviendrais jeudi sur l’idée du cours, et sur ce qui est demandé (cf doc).

Le travail sera réalisé par groupe de 2 étudiants. Des créneaux seront réservés sur l’emploi du temps pour le suivi des projets. Des soutenances seront organisées  auxquelles tous les étudiants devront assister.

D’un point de vue plus pratique, il s’agit de préparer un exposé de 45 minutes présentant une méthode statistique choisie dans la liste ci-dessous. L’exposé s’adressera à un public de scientifiques et/ou d’économistes ayant déjà un bagage mathématique et statistique et devra leur permettre de comprendre la méthode, son utilisation et son application. La présentation de composera

  • d’une introduction permettent notamment de se faire une idée des applications possibles de la méthode présentée ; on situera aussi la méthode par rapport à d’autres méthodes existantes.
  • d’une partie théorique.
  • d’un exemple illustratif (les étudiants choisiront eux-mêmes le problème, avec des données associées et l’exemple sera traité avec un logiciel statistique). L’application doit bien montrer en quoi les résultats obtenus apportent des informations importantes et utilisables dans le cadre d’une étude.

De plus, il sera important de dominer la portée interprétative des résultats statistiques en relation avec la réalité étudiée et en rapport avec les raisons pour lesquelles cette réalité est étudiée comme on le fait dans un contexte professionnel.

  1. Régressions Ridge et Lasso
  2. Régression PLS
  3. Arbres de décision et méthode CART
  4. Bagging et Boosting
  5. Réseaux de neurones
  6. Perceptron et SVM
  7. Processus AR à seuil
  8. Processus AR fractionnaire
  9. Modèles de chaînes de Markov cachées
  10. Quantiles et expectiles
  11. Algorithme EM
  12. Estimation de mélanges de lois
  13. Théorie des Valeurs Extrêmes
  14. Interpolation spatiale et Krigeage
  15. ICA, Analyse en Composantes Indépendantes
  16. Fonction de profondeur et quantiles multivariés

What happens if we forget a trivial assumption ?

Last week, @dmonniaux published an interesting post entitled l’erreur n’a rien d’original  on  his blog. He was asking the following question : let https://latex.codecogs.com/gif.latex?a, https://latex.codecogs.com/gif.latex?b and https://latex.codecogs.com/gif.latex?c denote three real-valued coefficients, under which assumption on those three coefficients does https://latex.codecogs.com/gif.latex?ax^2+bx+c has a real-valued root ?

Everyone aswered https://latex.codecogs.com/gif.latex?b^2-4ac\geq%200, but no one mentioned that it is necessary to have a proper quadratic equation, first. For instance, if both https://latex.codecogs.com/gif.latex?a and https://latex.codecogs.com/gif.latex?b are null, there are no roots.

Continue reading What happens if we forget a trivial assumption ?

Somewhere else, part 170

Some articles and posts worth reading

Continue reading Somewhere else, part 170

La magie des pourcentages

Ce matin, je me suis fait interpeller sur twitter par @PierreAgeron qui m’interrogeait sur la Figure suivante

Il me demandait comment les hausses observées en 2014  par trimestre (+0.0%, +0.0%, +0.1%, et +0.1%) pouvait donner +0.4% “globalement” pour 2014.

Le principe avec les pourcentages, c’est qu’on parle de variation, par rapport à une référence, et si on change la référence, alors forcément, on perd toute cohérence. On mélange des choux et des carottes. Cela dit, ce qui me gênait le plus, ce n’était pas 2014, mais2012. En effet, on a une hausse de +0.2% suivi d’une baisse de -0.2%, et une hausse de +0.3% suivi d’une baisse de -0.3%. Dit comme ça, on devrait s’attendre à une baisse globale, d’après une égalité que l’on connaît tous,

https://latex.codecogs.com/gif.latex?(1+x)(1-x)=1-x^2

Une hausse de https://latex.codecogs.com/gif.latex?x% suivie d’une baisse de https://latex.codecogs.com/gif.latex?x%, c’est une baisse de https://latex.codecogs.com/gif.latex?x^2%. Et pourtant, en 2012, on nous dit qu’on a une hausse.

Heureusement, Damien Babet, a.k.a. @DaBab, m’a proposé un joli exemple numérique que je vais reprendre ici . Considérons les huit valeurs suivantes pour 2011 et 2012,

> v=c(1011,1010,1012,1014,1016,1014,1017,1014)
> names(v)=c(paste("2011Q",1:4,sep=""),
+ paste("2012Q",1:4,sep=""))
> barplot(v,ylim=c(1000,1020))

En fait, cet exemple correspond effectivement au cas décrit initialement, puisqu’en 2012, on a une hausse puis une baisse du même taux,

> d=diff(v)/v[1:7]
> barplot(dv)

On est donc bien dans la situation du graphique introductif. Et maintenant, quand on agrège nos données trimestrielle, on a bien une hausse

> v2=c(sum(v[1:4]),sum(v[4+1:4]))
> names(v2)=c(2011,2012)
> barplot(v2,ylim=c(4000,4080))
> (v2[2]/v2[1]-1)*100
     2012 
0.3459353

(si vous n’êtes pas convaincu, mettez une grosse grève qui aurait paralysé la France pendant deux trimestre en 2011, avec 0 en Q1 et Q2…). Étonnant, non ? J’avais prévenu ! Les pourcentages, c’est de la magie !

How to import some parts of a large database

In the introduction of Computational Actuarial Science with R, there was a short paragraph on how could we import only some parts of a large database, by selecting specific variables. The trick was to use the following

> read.table.select.columns=function(datatablename,
I,sep=";"){
+ datanc=read.table(datatablename,header=TRUE,
sep=sep,skip=0,nrows=1)
+ mycols=rep("NULL",ncol(datanc))
+ names(mycols)=names(datanc)
+ mycols[I]=NA
+ datat=read.table(datatablename,header=TRUE,
sep=sep,colClasses=mycols)
+ return(datat)}

For instance, if we use the same dataset as in the introduction, we can import only two variables of interest,

> loc="http://myweb.fsu.edu/jelsner/extspace/extremedatasince1899.csv"
> dt1=read.table.select.columns(loc,c("Region",
"Wmax"),sep=",")
> head(dt1,10)
    Region      Wmax
1    Basin 105.56342
2    Basin  40.00000
3    Basin  35.41822
4    Basin  51.06743
5  Florida  87.34328
6    Basin  96.64138
7     Gulf  35.41822
8       US  35.41822
9       US  87.34328
10      US 106.35318
> dim(dt1)
[1] 2100    2

Continue reading How to import some parts of a large database

Cross Validation for Kernel Density Estimation

In a post publihed in July, I mentioned the so called the Goldilocks principle, in the context of kermel density estimation, and bandwidth selection. The bandwith should not be too small (the variance would be too large) and it should not be too large (the bias would be too large). Another standard method to select the bandwith, as mentioned this afternoon in class is the cross-validation technique (described in Chiu (1991)). Here, we would like to minimize

https://latex.codecogs.com/gif.latex?\mathbb{E}\left[\int%20[\widehat{f}_h(x)-f(x)]^2dx\right]

The integral can be writen

https://latex.codecogs.com/gif.latex?\int%20\widehat{f}_h(x)^2dx-2\int%20\widehat{f}_h(x)f(x)dx+\int%20f(x)^2dx

Since the third component is constant, we have to minimize the expected value of the sum of the first two.

The idea is to approximate it as

https://latex.codecogs.com/gif.latex?J(h)=\int%20\widehat{f}_h(x)^2dx-\frac{2}{n}\sum_{i=1}^n%20\widehat{f}_{(-i)}(X_i)

which can easily be computed. Consider here some sample, with 50 observations, from a Gaussian distribution,

> set.seed(1)
> X=rnorm(50)

From Silverman’s rule of thumb (which should be appropriate here since the sample has a Gaussian sample) the optimal bandwidth is

> 1.06*sd(X)*length(X)^(-1/5)
[1] 0.4030127

Using the cross-validation technique mentioned above, compute

> J=function(h){
+ fhat=Vectorize(function(x) density(X,from=x,to=x,n=1,bw=h)$y)
+ fhati=Vectorize(function(i) density(X[-i],from=X[i],to=X[i],n=1,bw=h)$y)
+ F=fhati(1:length(X))
+ return(integrate(function(x) fhat(x)^2,-Inf,Inf)$value-2*mean(F))
+ }
> vx=seq(.1,1,by=.01)
> vy=Vectorize(J)(vx)
> df=data.frame(vx,vy)
> library(ggplot2)
> qplot(vx,vy,geom="line",data=df)

The function has the following shape

and the optimal value is

> optimize(J,interval=c(.1,1))
$minimum
[1] 0.4687553

$objective
[1] -0.3355477

Note that, indeed, it is close to Siverman’s optimal bandwidth.

Somewhere else, part 169

Some articles and posts worth reading

Continue reading Somewhere else, part 169