eyes_hair

Test du chi-deux et indépendance

Considérons le tableau de contingence suivant

> N=margin.table(HairEyeColor, c(1,2))
> N
       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16

avec ici les comptages, que l’on peut aussi traduire sous la forme de probabilités jointes (empiriques)

> N/n
       Eye
Hair    Brown  Blue Hazel Green
  Black 0.115 0.034 0.025 0.008
  Brown 0.201 0.142 0.091 0.049
  Red   0.044 0.029 0.024 0.024
  Blond 0.012 0.159 0.017 0.027

 

Pour faire un test d’indépendance entre la couleur des yeux et la couleurs, rappelons que sous l’hypothèse d’indépendance

doit être égale à

et ce,  pour tout couple .  Ici, sous hypothèse d’indépendance, on devrait avoir les probabilités jointes suivantes

> n=sum(N)
> NHair=apply(N,1,sum)
> NEye =apply(N,2,sum)
> Nind=NHair%*%t(NEye)/n
> Nind/n
           Brown       Blue      Hazel      Green
Black 0.06779584 0.06625502 0.02865915 0.01972243
Brown 0.17953342 0.17545311 0.07589367 0.05222790
Red   0.04456949 0.04355654 0.01884074 0.01296567
Blond 0.07972288 0.07791100 0.03370104 0.02319211

ce qui devrait se traduire par les comptages suivants

> Nind
          Brown      Blue    Hazel     Green
Black  40.13514  39.22297 16.96622 11.675676
Brown 106.28378 103.86824 44.92905 30.918919
Red    26.38514  25.78547 11.15372  7.675676
Blond  47.19595  46.12331 19.95101 13.729730

La distance que l’on va utiliser entre les deux tableaux de contingence est la distance du chi-deux,

soit ici

> X=(N-Nind)^2/Nind
> Qobs=sum(X)
> Qobs
[1] 138.2898

Notons qu’il est toujours informatif de regarder les termes de la somme, qui sont les contributions au chi-deux

> X
       Eye
Hair    Brown  Blue Hazel Green
  Black 19.34  9.42  0.22  3.81
  Brown  1.52  3.80  1.83  0.11
  Red    0.00  2.99  0.72  5.21
  Blond 34.23 49.69  4.96  0.37

Rappelons que sous l’hypothèse d’indépendance,

On peut visualiser la densité de cette loi

> library(ggplot2)
> x=seq(0,140)
> ggplot(data=data.frame(x=x), aes(x = x))+
+ stat_function(fun=Vectorize(function(x)
+ dchisq(x,df=(ncol(N)-1)*(nrow(N)-1))),
+ geom="line", col = "red", size = 1) + 
+ ylab("Fonction de densité")

On voit que la valeur observée n’est pas vraisemblable pour une telle loi. La p-value est proche de 0 ici,

> 1-pchisq(Qobs,df=(ncol(N)-1)*(nrow(N)-1))

C’est en fait ce que l’on retrouve en faisant le test

> chisq.test(N)

	Pearson's Chi-squared test

data:  N
X-squared = 138.2898, df = 9, p-value < 2.2e-16

Notons que cette fonction permet d’obtenir le tableau de contingence sous hypothèse d’indépendance

> chisq.test(N)$expected
       Eye
Hair        Brown      Blue    Hazel     Green
  Black  40.13514  39.22297 16.96622 11.675676
  Brown 106.28378 103.86824 44.92905 30.918919
  Red    26.38514  25.78547 11.15372  7.675676
  Blond  47.19595  46.12331 19.95101 13.729730
Selection_719

Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction with Ewen Gallic is now online, on hal.archives-ouvertes.fr/.

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and difficult to implement. We provide a simple technique — based of properties of Gaussian kernels — to efficiently compute weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. We illustrate the use of that technique to visualize hot spots of car accidents and campsite locations, as well as location of bike thefts.

There are new applications, and new graphs, too

Most of the codes can be found on https://github.com/ripleyCorr/Kernel_density_ripley (as well as datasets).

Somewhere else, part 175

(Halloween Tree, by Glen Brogan)

Some posts and articles worth reading, here and there

A bank’s earnings are a quantum event; they are entirely probabilistic, and the answer you get depends on who’s doing the observing. You make some guesses with some degree of statistical likelihood, and then you apply one of a half-dozen accounting regimes to the guesses, and you get a number, and then you’re like, ooh, look at this number, it’s so numeric.

Continue reading

Selection_701

Removing Uncited References in a Tex File (with R)

Last week, with @3wen, we were working a the revised version of our work on smoothing densities of spatial processes (with edge correction). Usually, once you have revised the paper, some references were added, others were droped. But you need to spend some time, to check that all references are actually mentioned in the paper. For instance, consider the following compiled tex file :

Continue reading

Somewhere else, part 174

Some posts and articles worth reading

Continue reading

Loi multinomiale et loi du chi-deux

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

http://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 http://latex.codecogs.com/gif.latex?\mathcal{N}(0,1). Et plus généralement, on peut montrer que http://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 http://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}), alors http://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 à http://latex.codecogs.com/gif.latex?k degrés de liberté (dimension du vecteur http://latex.codecogs.com/gif.latex?\boldsymbol{X}). Car il faut que http://latex.codecogs.com/gif.latex?\boldsymbol{\Sigma} soit inversible. Or ici, je pense qu’on peut montrer que

http://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 http://latex.codecogs.com/gif.latex?n près). Et comme on a une loi de probabilité, le terme de droite est nul. Bref, http://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=0 et son inverse n’existe alors pas. Cela dit, http://latex.codecogs.com/gif.latex?0 n’est valeur propre qu’une fois, d’où la dimension http://latex.codecogs.com/gif.latex?k-1 qui apparaît sans cesse (à condition que les http://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

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

suit asymptotiquement une loi http://latex.codecogs.com/gif.latex?\mathcal{N}(0,\star). Et plus généralement, là encore, http://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),

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

puisque

http://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

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

de telle sorte que http://latex.codecogs.com/gif.latex?\tilde{X}_1^2=\tilde{X}_2^2 alors que http://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

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

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

http://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

http://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

http://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

http://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

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

et non pas la version avec http://latex.codecogs.com/gif.latex?\sqrt{np_i(1-p_i)} au dénominateur. On peut alors réécrire la somme en faisant intervenir http://latex.codecogs.com/gif.latex?k-1 termes (et non pas http://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

http://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

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

et non pas

http://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

Somewhere else, part 173

Some posts and articles worth reading

Continue reading

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

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ù http://latex.codecogs.com/gif.latex?X_{i-1}\leq%20X_i (différence positive) et le nombre de fois où http://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

Somewhere else, part 171

Some posts are articles worth reading

Continue reading

Selection_595

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 http://latex.codecogs.com/gif.latex?a, http://latex.codecogs.com/gif.latex?b and http://latex.codecogs.com/gif.latex?c denote three real-valued coefficients, under which assumption on those three coefficients does http://latex.codecogs.com/gif.latex?ax^2+bx+c has a real-valued root ?

Everyone aswered http://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 http://latex.codecogs.com/gif.latex?a and http://latex.codecogs.com/gif.latex?b are null, there are no roots.

Continue reading

Somewhere else, part 170

Some articles and posts worth reading

Continue reading

An Open Lab-Notebook Experiment


Research blogs