# Extracting datasets from excel files in a zipped folder

The title of the post is a bit long, but that’s the problem I was facing this morning: importing dataset from files, online. I mean, it was not a “problem” (since I can always download, and extract manually the files), more a challenge (I should be able to do it in R, directly). The files are located on ressources-actuarielles.net, in a zip file. Those are mortality tables used in French speaking African countries, and I guess that one problem came from special characters, such as “é” or “è”… When you open the zip file, you see a folder

and in that folder, several files that I would like to import

# Log-transform kernel density estimation of income distribution

Our paper Log-transform kernel density estimationof income distribution, written with Emmanuel Flachaire is now available on http://papers.ssrn.com/id=2514882,

Standard kernel density estimation methods are very often used in practice to estimate density function. It works well in numerous cases. However, it is known not to work so well with skewed, multimodal and heavy-tailed distributions. Such features are usual with income distributions, defined over the positive support. We first show that a preliminary logarithmic transformation of the data, combined with standard kernel density estimation methods, can provide a much better fit of the overall density estimation. Then, we show that the fit of the bottom of the distribution may not be satisfactory, even if a better fit of the upper tail can be obtained in general.

# Tests basés sur la vraisemblance – score

Une autre grandeur intéressant est le score, qui est la dérivée de la vraisemblance. Intuitivement (c’est l’idée de la condition du premier ordre),  $\widehat\theta_n$ et $\theta_0$ seront proches si les dérivées en ces points sont proches. En $\widehat\theta_n$ la dérivée est nulle, donc on va se demander ici, tout simplement, si la dérivée en $\theta_0$ est proche de 0. Ou pas.

# Tests basés sur la vraisemblance – Rapport de Vraisemblance

Chose promise, chose due. J’avais dit qu’on parlerait du test de rapport de vraisemblance. L’idée – visuelle – est d’avoir une lecture dans l’autre sens : au lieu de se demander si $\widehat\theta_n$ et $\theta_0$ sont proches, on va se demander si $\log\mathcal{L}(\widehat\theta_n)$  et $\log\mathcal{L}(\theta_0)$ sont proches.

Si la fonction de vraisemblance est suffisamment régulière, on se pose la même question.

# Tests basés sur la vraisemblance – Wald

Petit rappel préliminaire. Si on dispose d’un échantillon $X_i$, i.id de loi $F_{\boldsymbol{\theta}_\star}$, où le paramètre $\boldsymbol{\theta}_\star$ est inconnu, alors on peut calculer le maximum de la vraisemblance, ce qui nous donnera un estimateur intéressant (cf. cours de stat). Illustrons avec un jeu de pile ou face. Reprenons l’échantillon de mon précédent billet,

> X=c(0, 0, 1, 1, 0, 1, 1, 1, 1, 0,
0, 0, 1, 0, 1, 0, 1, 1, 0, 1)

On peut ici tracer la fonction de vraisemblance

$\theta\mapsto \mathcal{L}(\theta,\mathcal{X})=\theta^{\sum X_i}[1-\theta]^{n-\sum X_i}$

ou (un peu plus malin) la fonction de log-vraisemblance

$\theta\mapsto \log\mathcal{L}(\theta,\mathcal{X})$

> p=seq(0,1,by=.01)
> logL=function(p)
+ {sum(log(dbinom(X,size=1,prob=p)))}
> plot(p,Vectorize(logL)(p),
+ type="l",col="red",lwd=2)
> p0=.5
> abline(v=p0,col="blue")

La valeur correspondant au trait bleu correspondant au cas que l’on va chercher à tester, de pièce non pipée, soit $\theta_0$ (de manière très générale)

# Les tests et la logique, modus tollens

Quand on apprend la logique, on apprend la notion de modus tollens, correspondant à une notion de contraposition. Si on a une proposition du genre $A \Rightarrow B$, alors la proposition contraposée est $\text{non }B \Rightarrow \text{non }A$.  Et on apprend que les deux propositions sont équivalentes (je fais de la logique classique). Par exemple si $A$ correspond à “feu” et $B$ à “fumée“, $A \Rightarrow B$ signifie que tout feu fait de la fumée. Si cette affirmation est vrai, alors il n’y a pas de fumée sans feu, et donc $\text{non }B \Rightarrow \text{non }A$, autrement dit, s’il n’y a pas de fumée, c’est qu’il n’y a pas de feu.

# 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

# Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction 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).

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

# Loi multinomiale et loi du chi-deux

La semaine passée, en cours, j’avais rappelé que quand $\boldsymbol{N}=(N_{1},\cdots,N_{k})$ décrivait le compte de $n$ variable multinomiales prenant $k$ modalités, la variable

$\tilde{X}_i = \frac{N_i-np_i}{\sqrt{np_i\color{red}{(1-p_i)}}}$

suit asymptotiquement une loi $\mathcal{N}(0,1)$. Et plus généralement, on peut montrer que $\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 $\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, alors $Y=(\boldsymbol{X}-\boldsymbol{\mu})^{\text{\sffamily T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{X}-\boldsymbol{\mu})$ suit une loi du chi-deux à $k$ degrés de liberté (dimension du vecteur $\boldsymbol{X}$). Car il faut que $\boldsymbol{\Sigma}$ soit inversible. Or ici, je pense qu’on peut montrer que

$\det(\boldsymbol{\Sigma})=\prod_{i=1}^k p_i \cdot\left(\sum_{i=1}^k p_i -1\right)$

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

$X_i = \frac{N_i-np_i}{\sqrt{np_i}}$

suit asymptotiquement une loi $\mathcal{N}(0,\star)$. Et plus généralement, là encore, $\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),

$\tilde{X}_1^2=X_1^2+X_2^2$

puisque

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

$(N_1-np_1)^2=(N_2-np_2)^2$

de telle sorte que $\tilde{X}_1^2=\tilde{X}_2^2$ alors que ${X}_1^2\neq {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

$n!\sim n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}$

alors pour $n_1+\cdots+n_k=n$

$\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=\frac{n!}{n_1 ! \cdots n_k !}p_1^{n_1}\cdots p_k^{n_k}$soit

$\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim\frac{n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}}{\prod_{j=1}^k [n_j^{n_j+\frac{1}{2}}e^{-n_j}\sqrt{2\pi}]}p_1^{n_1}\cdots p_k^{n_k}$

et donc

$\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim \kappa \prod_{j=1}^k \left(\frac{np_j}{n_j}\right)^{n_j+\frac{1}{2}}$

de telle sorte que

$\log \mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=k+\sum_{j=1}^k \left(n_j+\frac{1}{2}\right) \log \left(\frac{np_j}{n_j}\right)$On peut alors faire des développements limités, en faisant intervenir

$x_i=\frac{n_i-np_i}{\sqrt{np_i}}$

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

$\sum_{i=1}^k X_i^2=\sum_{i=1}^k \frac{(N_i-np_i)^2}{np_i}\overset{\mathcal{L}}{\rightarrow} \chi^2(k-1)$

C’est donc

$Q=\sum_{i=1}^k \frac{(N_i-np_i)^2}{np_i}$

et non pas

$\tilde{Q}=\sum_{i=1}^k \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)
S98=(soccer1998$V1+soccer1998$V2)
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é à

# The book is out

That’s official ! Not only on amazon.com or crc.com… I did actually hold a copy of the book in my hands. I am loving it ! (if I forget about small typos we’ve seen over the past weeks…)

And the package is online on http://cas.uqam.ca/. An additional webpage will be uploaded, soon, for typos and additional material !

# 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ù $X_{i-1}\leq X_i$ (différence positive) et le nombre de fois où $X_{i-1}\geq X_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")

# 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 $a$, $b$ and $c$ denote three real-valued coefficients, under which assumption on those three coefficients does $ax^2+bx+c$ has a real-valued root ?

Everyone aswered $b^2-4ac\geq 0$, but no one mentioned that it is necessary to have a proper quadratic equation, first. For instance, if both $a$ and $b$ are null, there are no roots.