Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

Loi des grands nombres et théorème central limite en assurance

Dans le dernier numéro de la revue Risques figure un petit article sur la loi des grands nombres et le théorème central limite en assurance (ici) Le but était de remettre un peu les pendules à l’heure sur la différence entre les implications de ces deux résultats, et surtout d’amorcer une discussion sur l’importance des hypothèses (en particulier celle de risques indépendants).

Valeur de l’euro de rente

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrageévaluation du préjudice corporel, de Max le Roy. Comme l’indique le descriptif du livre, j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs, “La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

Avec la mise à jour des tables de mortalité en France, j’ai du remonter mes manches et reprogrammer tous les calculs. Sauf que cette fois, j’ai proposé un petit code R pour générer à nouveau ces tables. Le but est de calculer le prix de mille euros de rente viagère ou temporaire limitée. Une petite option permet de dire si on veut la nouvelle table, ou l’ancienne (ce qui m’a permis de vérifier que je ne disais pas des bêtises).

LF=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
header=TRUE,sep=";")
LH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
header=TRUE,sep=";") 
PRIX=function(annee=2011,age,sex="HOM",taux=0.04,
              duree,C=1000,old=FALSE){
if(old==FALSE){
file  = paste("http://freakonometrics.blog.free.fr/public/data/",
              sex,"-table-SPLx.csv",sep="")
B     = read.table(file,header=TRUE,sep=",")
an    = annee-age; if(an>2005){an=2005}
nom   = paste("X",an,sep="")
L     = B[,nom]}
if(old==TRUE){
if(sex=="HOM"){L = c(LH$Lx,rep(0,20))}
if(sex=="FEM"){L = c(LF$Lx,rep(0,20))}}
Q     = L[(age+1):length(L)]/L[(age+1)]
actualisation = (1+taux)^(0:min(duree,120-age))
prixsup = sum(Q[2:(min(duree,120-age)+1)]/
          actualisation[2:(min(duree,120-age)+1)] )
prixinf = sum(Q[1:(min(duree,120-age))]/
          actualisation[1:(min(duree,120-age))] )
return(C*c(prixsup,prixinf))}

Les tables TGH-05 et TGF-05 sont issues d’une population de rentiers et on peut consulter les tables ici (pour les 38 pages du Journal Officiel). Ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici (merci d’ailleurs à Fred qui a approuvé mon utilisation de ces tables).

Par exemple, on peut retrouver la valeur figurant dans les anciennes Annexes (ici),

> PRIX(age=50,sex="HOM",duree=150,taux=.05,old=TRUE)
[1] 13452.15 14452.15

En bouclant, on peut générer les tables données en Annexe.

AGES=c(150,65,60,55,25,20,18,16)
tableH=matrix(0,106,length(AGES))
for(i in 1:length(AGES)){
for(a in 0:min(c(AGES[i]-1,105))){
tableH[a+1,i]=PRIX(age=a,sex="FEM",duree=
        AGES[i]-a,taux=.04,old=FALSE)[1]
}}
write(tableH,"table-femme-nouvelle-04.csv",sep=";")

(on utilise la première valeur var on veut des rentes versées en fin d’année, en cas de vie). Avec un taux de 4%, on a pour les hommes la table ici, et pour les femmes la table . Si quelqu’un voit une coquille, merci de me signaler avant que le livre ne parte à l’éditeur…

Petit calcul de probabilités

En mars 1996, Parade Magazine publiait une question intéressante relative à un calcul de probabilité (la question est reproduite ici, ainsi qu’une discussion autour de l’existence de ce professeur… mais comme le note le site , il semble qu’il s’agisse d’une légende urbaine). De manière un peu surprenante, la réaction de la personne est que – si les deux élèves ont menti – la probabilité que deux personnes désigne le même pneu d’une voiture serait de 1/16. Pourtant, si on note A, B, C et D les quatre pneus de la voiture, un dénombrement rapide donne 16 situations possibles,

(A,A) (A,B) (A,C) (A,D)
(B,A) (B,B) (B,C) (B,D)
(C,A) (C,B) (C,C) (C,D)
(D,A) (D,B) (D,C) (D,D)

mais 4 paires concordantes. Soit un chance sur quatre que les étudiants indiquent le même pneu. Intuitivement, on demande à un étudiant quel pneu a crevé (jusque là pas de calculs de probabilités) http://freakonometrics.blog.free.fr/public/perso3/101_3.jpgpuis on veut voir si l’autre indique le même pneu (s’il ment, il a alors une chance sur quatre de donner le même que son voisin). Cette intuition que la probabilité serait 1/16 est alors un peu surprenante. Cela correspondrait au cas où le professeur sait quel pneu avait éclaté, mais pas les étudiants (on retrouve le mythe que l’enseignant sait toujours tout). Ou au cas où trois étudiants étaient dans la voiture et que l’on demande que les trois répondent la même chose.

Que signifie “utiliser de l’information” ?

Considérons un simple jeu de pile ou face. On suppose que l’on tombe sur pile avec une probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif (histoire de formaliser un peu). Le but du jeu est dedeviner sur quelle face la pièce va tomber: si on devine bien, on gagne, sinon on perd. Et le but de mon billet est de comprendre quelle est la stratégie qu’il convient de mettre en œuvrer si on connaît, ou si on estime http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif.

  • Connaissance parfaite des probabilités

Supposons que l’on connaisse parfaitement http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif. A priori deux stratégies sont possibles:

  • jouer au hasard, en utilisant les différentes probabilités: on prédit pile avec probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif et face avec probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-21.gif.
  • ne pas jouer au hasard, et prédire la valeur la plus probable, autrement dit si http://freakonometrics.blog.free.fr/public/perso3/pf-22.gif, on prédit pile, sinon on prédit face (à chaque lancer si on fait plusieurs lancers consécutifs).

Notons http://freakonometrics.blog.free.fr/public/perso3/pf-24.gif la variable correspondant à la valeur prise au lancer, disons 1 si on tomber sur pile et 0 sur face. Et notons http://freakonometrics.blog.free.fr/public/perso3/pf-23.gif la valeur prédite. On gagne si http://freakonometrics.blog.free.fr/public/perso3/pf-25.gif, et la probabilité de gagner est alors

http://freakonometrics.blog.free.fr/public/perso3/pf-26.gif

i.e. pour la première stratégie

http://freakonometrics.blog.free.fr/public/perso3/pf-27.gif

En revanche pour la seconde, la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-28.gif

Graphiquement, on peut représenter en ligne la prédiction que l’on fait, et en colonne la face obtenue. En bleu, on gagne, et en rouge on perd,

Supposons que http://freakonometrics.blog.free.fr/public/perso3/pf-29.gif. Considérons le polynôme http://freakonometrics.blog.free.fr/public/perso3/pf-30.gif, dont les racines sont http://freakonometrics.blog.free.fr/public/perso3/pf-31.gif et http://freakonometrics.blog.free.fr/public/perso3/pf-32.gif, i.e sur http://freakonometrics.blog.free.fr/public/perso3/pf-33.gif,
http://freakonometrics.blog.free.fr/public/perso3/pf-34.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso3/pf-35.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/pf-36.gif

(avec une inégalité stricte si http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif est strictement inférieur à 1). Bref, il est toujoursplus intéressant de parier toujours la même chose (et de ne pas essayer de jouer au plus malin).

  • Connaissance imparfaite des probabilités

Supposons maintenant que la probabilité http://freakonometrics.blog.free.fr/public/perso3/pf-20.gif soit inconnue, mais estimable. Autrement dit, on a assisté à http://freakonometrics.blog.free.fr/public/perso3/pf-03.gif lancers, on a un échantillon http://freakonometrics.blog.free.fr/public/perso3/pf-01.gif et on estime la probabilité de tomber sur pile par

http://freakonometrics.blog.free.fr/public/perso3/pf-05.gif

Comme auparavant pour la première stratégie, on utilise un générateur de nombre aléatoire http://freakonometrics.blog.free.fr/public/perso3/pf-04.gif (qui tire une loi uniforme sur le segment unité), pour la prochaine valeur, on prédit

http://freakonometrics.blog.free.fr/public/perso3/pf-06.gif

Pour la seconde stratégie

http://freakonometrics.blog.free.fr/public/perso3/pf-07.gif

Autrement dit, si on estime que la probabilité de tomber sur pile dépasse 50%, on prédit pile, et sinon on prédit face.
On gagne si http://freakonometrics.blog.free.fr/public/perso3/pf-08.gif, et la probabilité de gagner est alors

http://freakonometrics.blog.free.fr/public/perso3/pf-09.gif

Dans le premier cas, la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-11.gif
dont l’espérance est http://freakonometrics.blog.free.fr/public/perso3/pf-12.gif (comme auparavant). En revanche, dans le second cas, le calcul de la probabilité de gagner est

http://freakonometrics.blog.free.fr/public/perso3/pf-13.gif

Graphiquement, on obtient, cette fois,

Si on se place au bout de 20 lancers, on obtient la courbe rouge ci-dessous pour le cas d’information imparfaite, et verte dans le cas de l’information parfaite, pour le ratio de la probabilité de gagner avec le second cas divisé par le premier cas (nous avions vu auparavant que le ratio été inférieur à 1)

à partir du code suivant

> p=seq(.5,1,by=.01)
> q=(p^2+(1-p)^2)/(p*(1-pbinom(10,20,p))+
+ (1-p)*pbinom(10,20,p))
> plot(p,q)
  • Répétons le jeu plusieurs fois et apprentissage des probabilités

Supposons que l’on répète le jeu un certain nombre de fois, avec un apprentissage au fur et à mesure que les tirages avancent, i.e. l’algorithme s’écrit

> n=10
> p=.64
> X=sample(c("a","b"),size=100,prob=c(p,1-p),replace=TRUE)
> S1=S2=rep(NA,100)
> for(s in (n+1):100){
+ P=c(mean(X[1:(s-1)]=="a"),mean(X[1:(s-1)]=="b"))
+ S1[s]=(sample(c("a","b"),size=1,prob=P)==X[s])
+ S2[s]=(names(sort(table(X[1:s])))[2]==X[s])
+ }
> sum(S1[(n+1):100])/sum(S2[(n+1):100])

Graphiquement, si on suppose que l’on apprend sur 10 tirages (et que l’on calcule le nombre moyen de victoires sur les 90 lancers suivants), on obtient un graphique assez proche de ce que l’on avait auparavant,

I.e. comme le ratio est inférieur à 1, la seconde stratégie l’emporte toujours.
Moralité, j’ai un peu l’impression qu’il est souvent inutile d’élaborer une stratégie trop complexe, car les stratégies les plus simples permettent de gagner le plus souvent…

Mal calculer les probabilités tue les probabilités

Suite à un tweet d’Olivier, j’ai découvert avec stupéfaction l’article paru dans Libé sur le risque nucléaire (ici).

On apprend que la probabilité d’accident majeur est d’environ 0.0003 0,0003 par an pour chaque réacteur. Or l’Europe possède 143 réacteurs et donc, sur 30 ans, la probabilité d’un incident majeur est

http://freakonometrics.hypotheses.org/files/2015/12/nucle-05.gif

et les auteurs de l’article concluent que “la probabilité d’occurrence d’un accident majeur sur ces parcs serait donc […] de plus de 100% pour l’Union européenne”

Etienne Ghys a déjà – brillamment – discuté ce chiffre avant hier (ici). Le soucis est que si http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif est la proba annuelle d’incident par an et par réacteur, en supposant les centrales indépendantes, la probabilité d’avoir un incident sur http://freakonometrics.hypotheses.org/files/2015/12/nucle-3.gif réacteurs sur une durée de http://freakonometrics.hypotheses.org/files/2015/12/nucle-04.gif années est

http://freakonometrics.hypotheses.org/files/2015/12/nucle-01.gif

Comme http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif est supposé petit, un développement limité donnerait

http://freakonometrics.hypotheses.org/files/2015/12/nucle-06.gif

Malheureusement, dans les développements limités de ce genre, s’il faut que http://freakonometrics.hypotheses.org/files/2015/12/nucle-02.gif soit petit, encore faut-il que la puissance ne soit pas trop grande….
En fait, si on retient que pour un réacteur nucléaire à eau pressurisée (REP) tels ceux exploités en Europe de l’Ouest, le risque de fusion du cœur est estimé à 0.00005 par centrale et par an (en reprenant les estimateurs mentionnés sur wikipedia, ici), la vraie probabilité d’avoir un incident majeur en Europe serait plutôt de l’ordre de

http://freakonometrics.hypotheses.org/files/2015/12/nuclea-07.gif

(avec une probabilité 10 fois supérieure, Etienne arrivait à 72%). C’est marrant car cette erreur commise en croyant qu’un développement limité serait légitime quand on utilise des probabilité annuelle sur des horizons plus longs était évoqué ici il y a un an jour pour jour, sur des problèmes d’actualisation….

Trop de probabilités tue les probabilités

Hier matin, au moment de quitter la maison, j’ai jeté un œil à la météo (car à Montréal le temps est très changeant, pas seulement d’un jour sur l’autre, comme je l’avis mentionné en arrivant en septembre, ici). En voyant les gros nuages gris dans le ciel, en me réveillant, je me suis demandé s’il fallait mettre un manteau et des bottes en caoutchouc aux enfants (pour aller l’école). Sur le site de météo, on a en effet la probabilité de pluie, i.e. de P.D.P.

Si on suppose que les probabilités sont indépendantes d’heure en heure, la probabilité qu’il va pleuvoir dans la matinée (dans les 5 prochaines heures pour faire simple, ou plus concret), va s’écrire comme le complémentaire de la probabilité qu’il va faire beau toute la matinée, autrement dit, avec des probabilités de pluie de l’ordre de 25%, par heure, la probabilité qu’il pleuve pendant les prochaines heures est de 75% environ

> 1-((1-.25))^5
[1] 0.7626953

Bref, on est loin de 25%. Essayons d’aller un peu plus loin pour mieux comprendre…
Notons http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie03.gif la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie04.gif. Supposons qu’il existe un processus latent sous-jacent, http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie05.gif tel que

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie06.gif

 

A une date http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie04.gif, on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie07.gifqui peut encore s’écrire

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie08.gifou, ce qui sera pleut être plus utile ensuite

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie09.gifAvec notre modèle latent, cette probabilité peut s’écrire

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie10.gifc’est à dire en notant http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie13.gif la fonction de répartition (jointe) du vecteur aléatoire http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie11.gif cette probabilité s’écrit
http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie12.gifOr Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout https://latex.codecogs.com/gif.latex?\boldsymbol{x}\in\mathbb{R}^d,
http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie14.gif

http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie15.gif
et
http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie16.gif
(en notant que ces bornes sont atteignables). Autrement dit, ici, la probabilité qu’il fasse beau 5 heures de suite est encadrée par la probabilité qu’il pleuve sur une heure… et 1.
Bref, on n’a pas grand chose à dire*.

Une possibilité est de supposer un modèle Gaussien pourhttp://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie05.gif. Supposer les lois marginales n’est en rien contraignant. En revanche, supposer que la dynamique sous-jacente est Gaussienne est une hypothèse très forte. Et pour faire simple, supposons même que la dynamique soit auto-régressive (à l’ordre 1), i.e.
http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie18.gifoù https://latex.codecogs.com/gif.latex?(\eta_t) est un bruit blanc, indépendant du processus latent. L’intérêt est que l’on a spécifié ainsi la structure de la matrice de variance covariance. En particulier, la matrice de corrélation du vecteur https://latex.codecogs.com/gif.latex?\boldsymbol{X}^\star est alors
http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie01.gif
On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de https://latex.codecogs.com/gif.latex?\rho,
http://freakonometrics.hypotheses.org/files/2015/12/animationcorrelationpluie.gif
Avec ce formalisme, on peut alors en déduire http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie12.gif car on connaît très bien la fonction de répartition d’un vecteur gaussien multivarié via

On peut alors étudier l’impact de http://freakonometrics.blog.f ree.fr/public/perso3/slatex.gif (qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation http://freakonometrics.hypotheses.org/files/2015/12/correlation-pluie21.gif (qui est associée à la dynamique du processus). Pour simplifier, si on prend http://freakonometrics.hypotheses.org/files/2015/12/slatex.gif tel que la probabilité de pluie soit de l’ordre de 25%, on a la probabilité d’avoir de la pluie suivante,

> s = .25
> r = .85
> i = 1:n
> I = matrix(rep(i,each=n),n,n)
> J = matrix(rep(i,n),n,n) 
> C = r^(abs(I-J))
> 1-pmnorm(rep(qnorm(1-s),n),mean=rep(0,n),varcov=C)[1]
[1] 0.458074

i.e. en faisant une boucle sur la valeur de la corrélation

Et plus généralement si on croise le seuil (i.e. la probabilité horaire de pleuvoir) et la probabilité, on a l’abaque suivante

On pourrait toutefois objecter que l’important n’est pas forcément qu’il ne pleuve pas de la matinée, mais qu’il ne pleuve pas quand ils sont dehors, i.e. (pour simplifier), à 8 heures, midi et 16 heures (c’est à dire toutes les 4 heures). On s’intéresse alors à
http://freakonometrics.hypotheses.org/files/2015/12/pluiecorrel-31.gif
qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,
http://freakonometrics.hypotheses.org/files/2015/12/ppluie-coreel.gif
Dans ce cas, la probabilité de ne pas avoir de plus à ces trois moments de la journée est alors de la forme suivante

La corrélation joue alors un rôle relativement faible car il intervient à la puissance 4 dans la matrice d’auto-corrélation.
Dans ce cas, on notera que l’hypothèse d’indépendance n’est peut être pas absurde. Et avec une P.D.P. de 30%, sur toute la journée, on peut probablement se dire que la probabilité d’avoir de la pluie soit le matin, soit le midi, soit le soir sera de l’ordre de

> 1-(1-.3)^3
[1] 0.657

A tester !
* En fait si, on pourrait dire beaucoup de choses (et écrire une dizaine de billets) en particulier sur la borne supérieure, ici 1, dès lors que la probabilité horaire d’avoir de la pluie dépasse 20% – car la borne supérieure est

http://freakonometrics.hypotheses.org/files/2015/12/prooooobacommentaire.gif

En fait, la borne inférieur n’est pas une fonction de répartition de manière générale. Toutefois, la borne est atteignable pour une certaine structure de dépendance. Ce qui explique que par la suite, avec un modèle Gaussien, on sera relativement loin de cette borne supérieure. La borne inférieure (i.e. la borne supérieure sur les fonction de répartition) est elle atteinte dès lors que les risques sont variables sont comonotones, c’est à dire dans le cas Gaussien si la corrélation entre deux instants (et entre tous les instants) vaut 1. La petite animation sur la matrice de corrélation permet de mieux comprendre la difficulté de comprendre le cas où la corrélation rho est négative: dans le cas où http://freakonometrics.hypotheses.org/files/2015/12/rhofleche.gif1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si http://freakonometrics.hypotheses.org/files/2015/12/rhofleche.gif -1, on a des corrélations proches de -1 dans la moitié (environ) des cellules de la matrice, et +1 dans l’autre moitié. La borne supérieure sur la probabilité savoir de la pluie dans matinée correspond donc à une structure plus complexe que la borne inférieure, car on ne peut pas avoir une corrélation négative forte avec tous les lags.

Pourquoi frapper aussi fort la balle au service ?

Au tennis, il n’est pas rare de voir les joueurs servir très fort sur leur première balle, puis de se calmer sur la seconde. Est-ce pour se reposer un peu après le premier effort ou bien y-a-t-il une réelle logique derrière ?

Modélisons un peu tout ça (en reprenant une modélisation présentée par David Gale en 1971). Un joueur a un ensemble http://freakonometrics.blog.free.fr/public/perso3/tennisserv01.gif de services possibles. A son premier service, il choisit le service http://freakonometrics.blog.free.fr/public/perso3/tennisserv02.gif. La première difficulté est que le service soit correct (ou passe le filet pour faire simple), et on note http://freakonometrics.blog.free.fr/public/perso3/tennisserv03.gif la probabilité qu’un tel service passe le filet. Si le service est passé, il a une certaine probabilité de gagner le point, que l’on notera http://freakonometrics.blog.free.fr/public/perso3/tennisserv05.gif. Bien que ne connaissant pas vraiment le tennis, il n’est pas absurde de supposer qu’il y ait une forme d’arbitrage,

  • un service puissant a peu de chance de passer (http://freakonometrics.blog.free.fr/public/perso3/tennisserv06.gif petit), mais s’il passe le filet le serveur sera à son avantage, et la probabilité de gagner le point sera élevée (http://freakonometrics.blog.free.fr/public/perso3/tennisserv05.gif grand)
  • un service faible a de grandes chances de passer (http://freakonometrics.blog.free.fr/public/perso3/tennisserv06.gif grand), mais la balle sera alors facile pour l’adversaire, et la probabilité de gagner le point sera faible (http://freakonometrics.blog.free.fr/public/perso3/tennisserv05.gif petit)

A un service http://freakonometrics.blog.free.fr/public/perso3/tennisserv07.gif donné, la probabilité de gagner le point (sur ce premier service) est alors http://freakonometrics.blog.free.fr/public/perso3/tennisserv08.gif.

Supposons que http://freakonometrics.blog.free.fr/public/perso3/tennisserv01.gif soit “ordonné” au sens où http://freakonometrics.blog.free.fr/public/perso3/tennisserv06.gif peut être supposée strictement croissante.
Compte tenu de la remarque que nous avions faite initialement, il n’est pas absurde de supposer que la fonction

http://freakonometrics.blog.free.fr/public/perso3/tennisserv09.gif

soit à valeur dans http://freakonometrics.blog.free.fr/public/perso3/tennisserv10.gif, avec une fonction croissante pour http://freakonometrics.blog.free.fr/public/perso3/tennisserv.gif, puis décroissante pour http://freakonometrics.blog.free.fr/public/perso3/tennisserv12.gif.

Le joueur cherche à maximiser sa probabilité de gagner le point.

Mais au tennis, si le premier service ne passe pas, le serveur a le droit a une deuxième chance. Le serveur choisi alors en réalité un couple http://freakonometrics.blog.free.fr/public/perso3/tennisserv15.gif correspondant à son premier et second service, respectivement. La probabilité de gagner le point, à stratégie http://freakonometrics.blog.free.fr/public/perso3/tennisserv16.gif donnée, est alors

http://freakonometrics.blog.free.fr/public/perso3/tennisserv17.gif

On peut noter que la stratégie optimale pour la seconde balle est indépendante de ce qui a été fait pour la première. En effet, si on cherche à résoudre, à http://freakonometrics.blog.free.fr/public/perso3/tennisserv18.gif donné

http://freakonometrics.blog.free.fr/public/perso3/tennisserv20.gif

la condition du premier ordre est

http://freakonometrics.blog.free.fr/public/perso3/tennisser21.gif

i.e

http://freakonometrics.blog.free.fr/public/perso3/tennisserv22.gif

Autrement dit, le second service doit être celui qui maximise sa probabilité de gagner le point, i.e. http://freakonometrics.blog.free.fr/public/perso3/tennisserv23.gif.
Mais si on sait qu’il est optimal de jouer http://freakonometrics.blog.free.fr/public/perso3/tennisserv23.gif, dans ce cas, pour le premier service, le joueur cherche à résoudre

http://freakonometrics.blog.free.fr/public/perso3/tennissserv24.gif

dont la condition du premier ordre est

http://freakonometrics.blog.free.fr/public/perso3/tennisserv25.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/tennisserv26.gif

Le point s_1 qui vérifie cette condition peut se visualiser sur la figure ci-dessous, car http://freakonometrics.blog.free.fr/public/perso3/tennisserv27.gif s’interprète comme http://freakonometrics.blog.free.fr/public/perso3/tennisserv28.gif, i.e. la tangente de la fonction rouge, et on veut que cette tangente soit de pente http://freakonometrics.blog.free.fr/public/perso3/tennisserv29.gif, qui est la pente de la courbe bleue (comme l’avait noté David Gale),

Autrement dit, il est optimal de servir plus fort la première fois que la seconde… Comme quoi finalement tout peut se démontrer avec un joli modèle…

Minuit, l’heure du crime ? heure d’hiver ou d’été ?

Suite à mon précédant billet où j’essayais d’étudier la dynamique des appels aux 911 à Montréal, j’avais noté qu’un seul type d’appels avait lieu aux alentours de minuit, à savoir les trouble à la paix. Mais minuit est une notion quelque peu relative. En tous les cas, c’est la réflexion que je me faisais alors que les jours commencent à rallonger… car l’hiver à Montréal, minuit ça peut presque 8 heures après le coucher du soleil, et l’été à peine 4 heures…C’est aussi ce que me demandait @simon dans un commentaire.

En fait, de manière assez étrange, les cycles sont relativement stables en fonction des saisons (avec dans le cercle au centre l’ensoleillement de la journée – avec en bleu la nuit)

http://freakonometrics.blog.free.fr/public/perso3/911-mois.gif

Je mettrais peut être cela sur le compte que Montréal est relativement proche de l’heure solaire, au sens où midi est assez souvent le milieu de la journée et que minuit est le milieu de la nuit.

Dans le même état d’esprit, on peut se demander s’il y a un impact du changement d’heure. Par exemple si on considère le changement d’heure qui survient au printemps. Le graphique de gauche représente les appels le week-end qui précédait le changement d’heure (vendredi, samedi et dimanche qui précédaient) et à droite le week-end qui a suivi (même trois jours),

Il y a peut être un léger désaxage vers la droite, mais à peine…

Et si on regarde d’autres types d’appels, comme les coups de feu, on observe le début de la période d’appel (un peu plus intense) suit le commencement de la nuit: s’ils finissent vers 4 heures du matin, les appels arrivent plus tôt l’hiver que l’été

http://freakonometrics.blog.free.fr/public/perso3/911-mois-409.gif

en revanche, on a un impact presque nul du changement d’heure,

On peut aussi regarder les conflits, qui présentaient une dynamique relativement différente, on observe ici une relative stabilité

http://freakonometrics.blog.free.fr/public/perso3/911-mois-610.gif

et encore une fois, on a un impact presque nul du changement d’heure,

ma que calor !

J’avais commencé plusieurs analyse (ici ou ) sur les appels passé au 911, à la police, à Montréal, en fonction de l’heure et des jours de la semaine afin de détecter d’éventuels cycles. Mais ces derniers jours, maintenant qu’il commence à faire beau et chaud, je me suis demandé s’il n’était pas possible de croiser avec la température… Car j’aurais l’idée naïve qu’on ne va pas commettre un cambriolage quand il fait -25°C, ou que les conflits sont en hausse s’il fait +35°C…

Si on regarde assez rapidement, en comptant le nombre d’appels en fonction de la température, et en faisant une régression splines, on obtient

autrement dit, entre -20°C et +30°C, il y a environ 50% d’appels en plus. Ce qui est loin d’être négligeable ! Mais on peut aussi croiser la température avec l’heure d’appel, avec l’heure en abscisse et la température en ordonnée

Aussi, la région ou le nombre d’appels est le plus élevé est indiquée ci-dessous: on note que l’on a plus d’appels quand il faut chaud, et plutôt en fin de journée,

Pour les conflits, par exemple, on a davantage de conflits en fin d’après midi comme on l’avait relevé, mais surtout, ces derniers sont très liés à la température: il y a d’autant plus d’appels qu’il faut chaud.

Les appels sont d’autant plus importants dans la région bleue ci-dessous,

Si l’on se concentre uniquement sur la température et les conflits, on a l’évolution suivante

après presque deux fois plus d’appels s’il fait +30°C par rapport à -20°C. De même pour les agressions sexuelles, qui augmentent très fortement avec la température (environ +25%),

En revanche les holds ups sont indépendants de la température,

En revanche les alarmes cambrioleurs semblent générer plus d’appels s’il fait chaud, mais aussi s’il fait froid…

Oui, je suis fâché avec les chiffres

On a tous des secrets plus ou moins difficiles à révéler… Enfin, quand je dis secret c’est un bien grand mot (c’est un peu au sens que pourrait utiliser ma fille quand elle me dit “papa j’ai un secret: ce midi, j’ai mangé deux fois du dessert“). Par exemple, quand un collègue me dit “tu pourrais m’envoyer les transparents de ta présentation“, et que je dois lui répondre “tu sais, ils sont sur mon blog“: je dévoile un faux secret… Oui, tous les profs qui tiennent un blog en parlent un jour ou l’autre: avouer à des collègues qu’on tient un blog, c’est un secret que l’on révèle encore péniblement… Péniblement n’est peut être pas le mot car si une fois sur deux cela provoque une espèce de réaction condescendante, une fois sur deux cela permet au contraire d’entamer des discussions intéressantes car les interlocuteurs peuvent être curieux…

Non, en l’occurrence, je pensais plutôt aux conversations que j’ai au moins une fois par semaine, quand il faut récupérer un enfant chez un(e) ami(e) à lui (ou à elle). On discute alors avec les parents, et invariablement arrive la question “et vous faites quoi dans la vie ?”. J’aimerais bien pouvoir dire “je suis astronaute” ou “je suis constructeur de robots” ou “je suis dompteur de lions”, enfin un truc qui fait rêver. Mais non, je dois me contenter d’un “je suis prof de maths“. “Pardon?”. “Je suis prof de maths“. “Ahhhhhhh“. Oui, il y a toujours une espèce de déception que je peux lire dans le regard qui commence soudain à se dérober, comme gêné (genre il aurait espéré “je suis chirurgien cardiaque, et j’ai encore sauvé deux enfants de la mort aujourd’hui“, qui aurait donné lieu à un autre type de “ahhhhhhh”). Et là, arrive inéluctablement une phrase qui me place dans l’embarrât: au choix “moi j’ai jamais rien compris aux chiffres” ou “ben vous devez aimer les chiffres alors“…

Et là j’avoue que je ne sais jamais quoi dire… Je pourrais le prendre à la rigolade genre “oui, mais pas tous. Genre 37 je ne l’aime pas du tout“, mais le fait est que je n’aime pas le chiffres ! C’est d’ailleurs troublant que mathématique soit associé dans l’imaginaire populaire à calculs ou arithmétique. Au XVIIIème siècle, les gens qui faisaient des mathématiques étaient des géomètres. Et, en fait je préférerais qu’on me dise “ben vous devez aimer les dessins alors”… Ou à la rigueur des lettres, car si on regarde deux minutes le blog, mes papiers ou mes notes de cours, on voit des nombres d’accidents qui s’appellent http://freakonometrics.blog.free.fr/public/perso3/mm02.gif, des coûts qui s’appellent http://freakonometrics.blog.free.fr/public/perso3/mm01.gif, des tailles des gens qui s’appellent http://freakonometrics.blog.free.fr/public/perso3/mm04.gif, des choses inconnues qui s’appellent http://freakonometrics.blog.free.fr/public/perso3/mm03.gif ou http://freakonometrics.blog.free.fr/public/perso3/mm05.gif suivant le problème… bref, j’aime les lettres mais surtout pas les chiffres (sauf peut être quelques uns, je ne tiens pas à me brouiller avec tout le monde) !

Je me faisais cette réflexion ce matin, en lisant l’histoire d’une autre monsieur qui était accusé d’être fâché avec les chiffres, ici

Car il était annoncé que «  les deux tiers des échecs scolaires, c’est l’échec d’enfants d’immigrés ». Loin de moi l’idée de commenter l’utilisation de ces chiffres, surtout que je pense que ce monsieur et moi n’avons en commun que le fait d’être tous deux fâchés avec les chiffres. A la lecture de l’article, nous explique ensuite que  “la France compte 13% d’élèves de familles immigrées, soit 87% d’élèves de familles autochtones” selon des sources OCDE. On nous montre alors un graphique expliquant que “18% (des 87%) des enfants d’autochtones et 38% des enfants d’immigrés sont sous le niveau 2 qui correspond, selon le ministre, à une situation d’échec scolaire.”

Vient alors une explication “la France compte donc 15,7% (18% de 87%) d’enfants de 15 ans issus de familles autochtones sous le niveau 2 et 4,9% d’enfants de 15 ans issus de familles immigrés (38% de 13%) sous le niveau 2. Donc sur les 20,6% de « sous-niveau 2 » nous n’avons pas deux tiers mais plutôt un quart d’enfants d’immigrés.

Je sais pas vous, mais c’était au réveil, j’avais encore mon café dans la main, et je relisais cette phrase sans rien comprendre… et j’avais beau relire, impossible de comprendre qui disait n’importe quoi… Car personnellement, tant qu’on n’utilise pas un peu de formalisme, je suis perdu…
On peut par exemple faire un dessin… j’ai toute ma population schématisée par le grand rectangle. Dedans j’ai mes “enfant d’immigré” en rouge et “d’autochtones” en vert, et les enfants en échec en bleu, et ceux qui ne le sont pas en violet. Et la question c’est: au sein de la population bleue, quelle est la proportion qui sontbleu et rouge (qui est en jaune en bas à droite) ?

Si on note http://freakonometrics.blog.free.fr/public/perso3/mm06.gif le fait d’être en échec scolaire (population bleue) et http://freakonometrics.blog.free.fr/public/perso3/mm07.gif le fait d’être “enfant d’immigré” (population rouge, je ne discuterais pas la définition, mais les calculs), et que l’on suppose que http://freakonometrics.blog.free.fr/public/perso3/mm08.gif ce sont les enfants “d’autochtones” (population verte), alors on cherche http://freakonometrics.blog.free.fr/public/perso3/mm09.gif, et le monsieur prétend que ça doit être 2/3.
Bon, on sait que http://freakonometrics.blog.free.fr/public/perso3/mm10.gif et http://freakonometrics.blog.free.fr/public/perso3/mm11.gif, et que http://freakonometrics.blog.free.fr/public/perso3/mm12.gifet http://freakonometrics.blog.free.fr/public/perso3/mm13.gif.
Lançons nous dans les calculs: d’après la formule de Bayes,

http://freakonometrics.blog.free.fr/public/perso3/mm14.gif

L’astuce est de noter que

http://freakonometrics.blog.free.fr/public/perso3/mm15.gif

d’où la formule

http://freakonometrics.blog.free.fr/public/perso3/mm16.gif

et ici tout est connu… on a alors

http://freakonometrics.blog.free.fr/public/perso3/mm17.gif

Aussi on retrouve les 24% évoqués dans l’article…
Ce qui est bien quand on formalise, c’est que l’on est certain des chiffres qu’on avance… mais aussi, cela permettrait d’aller plus loin (conditionnellement aux chiffres à notre disposition…). Par exemple, on peut faire un tableau de contingence. On connaît les probabilités marginales

http://freakonometrics.blog.free.fr/public/perso3/mm18.gif
http://freakonometrics.blog.free.fr/public/perso3/mm19.gif

(on a fait le calcul au dessus pour écrire le dénominateur)
contenant les probabilités jointes,

http://freakonometrics.blog.free.fr/public/perso3/mm20.gif
http://freakonometrics.blog.free.fr/public/perso3/mm21.gif
http://freakonometrics.blog.free.fr/public/perso3/mm22.gif
http://freakonometrics.blog.free.fr/public/perso3/mm23.gif

Bref, on pourrait alors faire des tests du chi-deux si on savait sur combien d’élèves les tests ont été fait pour avoir de tels chiffres…

> M=matrix(c(.38*.13,.18*.87,
+ (1-.38)*.13,(1-.18)*.87),2,2)
> n=1000
> f=function(n){chisq.test(n*M)$p.value}
> F=Vectorize(f)
> P=seq(10,1000,by=10)
> plot(P,F(P),type="l",ylim=c(0,.4))
> abline(h=.05,lty=2)

i.e. on a graphiquement en fonction de la taille de l’échantillon (en abscisse) la p-value (en ordonnée) du test d’indépendance du chi-deux.

Autrement dit, sauf si ces statistiques sont tirés d’un échantillon de 200 élèves, on rejette l’hypothèse d’indépendance entre les résultats scolaires et l’origine des élèves. Mais j’avais dit que je ne commenterais pas…

Bref, mon point est qu’il me semble qu’en formalisant un peu, au lieu de s’embrouiller dans un discours rempli de chiffres, on fait moins d’erreurs.

Il y a des jours avec, et des jours sans

Je ne sais pas si ça vous est déjà arrivé, mais quand on fait des crêpes, la probabilité d’en réussir une est proportionnelle à la probabilité que l’on vient de réussir les dernières. Autrement dit si j’ai réussi les premières, je suis dans un bon jour, et je réussi davantage celles qui suivent. Alors que si je commence par rater les premières, je suis dans un mauvais jour, et je les rate toutes (pour les sportifs, il est aussi possible de noter que, par exemple, au tennis, la probabilité de réussir son service est d’autant plus forte que j’ai réussi les derniers: c’est pareil).

  • Une première formalisation du problème

Supposons que sur les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais (qui servent d’étalonnage), http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif fois j’ai réussi ma crêpe, et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif fois, je l’ai ratée (trop de trous, trop irrégulière, trop cuite, que sais-je encore). Mon a priori sur la probabilité de ne pas avoir de chance (et de rater ma crêpe) est http://freakonometrics.blog.free.fr/public/perso3/crepes-04.gif. Supposons que la probabilité de réussir à l’instant http://freakonometrics.blog.free.fr/public/perso3/crepes-06.gif (on commence à compter après les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais) soit la fréquence empirique observée sur le passé. Autrement dit

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est le nombre de succès. En particulier http://freakonometrics.blog.free.fr/public/perso3/crepes-09.gif. On retrouve ici le fait que plus j’en ai raté, plus j’ai de chances de rater les prochaines… c’est la mauvaise main en quelque sorte.
Afin de mieux comprendre, supposons http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif. Autrement dit, sur mes deux premiers essais (qui ne servent qu’à étalonner), j’ai raté et réussi une fois. Juste avant l’instant 1, j’ai deux possibilité: réussir ou rater, avec la probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-11.gif. I.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-12.gif

Si on se place une crêpe plus loin, et que l’on regarde les probabilités conditionnelles

http://freakonometrics.blog.free.fr/public/perso3/crepes-13b.gif

avec quelque-chose de similaire si on conditionne par l’autre événement,

http://freakonometrics.blog.free.fr/public/perso3/crepes-14.gif

Et donc la probabilité non conditionnelle s’écrit

http://freakonometrics.blog.free.fr/public/perso3/crepes-15.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-16.gif
http://freakonometrics.blog.free.fr/public/perso3/crepes-18.gif

et

http://freakonometrics.blog.free.fr/public/perso3/crepes-19.gif

Moralité

http://freakonometrics.blog.free.fr/public/perso3/crepes-20.gif

Autrement dit, on poursuivant la récurrence, on a que http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est uniformément distribué sur http://freakonometrics.blog.free.fr/public/perso3/crepes-21.gif
Étonnant, non ? Autrement dit, si j’ai commencé par un échec et une réussite, la probabilité d’avoir 50 réussites sur 100 essais sera très exactement http://freakonometrics.blog.free.fr/public/perso3/crepes-23.gif. Qui est la même probabilité que d’avoir 1 réussite, ou 99 !
On a aussi que

http://freakonometrics.blog.free.fr/public/perso3/crepes-22.gif

Mais que se passe-t-il si http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif ne sont plus égaux à 1 ?

  • Lien avec l’approche bayésienne et la théorie de la crédibilité

Ce qui est amusant, c’est que la probabilité

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

est bien connue par les personnes qui connaissent l’assurance et la théorie de la crédibilité. Supposons que je note http://freakonometrics.blog.free.fr/public/perso3/crepes-24.gif la variable indicatrice correspondant à une réussite ou un succès au http://freakonometrics.blog.free.fr/public/perso3/crepes-26.gifème essai (avec probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-25.gif). D’un point de vu bayésien, la loi conjuguée de la loi binomiale (on va compter les réussites sur n essais) est une loi Beta. Rappelons que la loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif a pour densité

http://freakonometrics.blog.free.fr/public/perso3/crepes-28.gif

sur http://freakonometrics.blog.free.fr/public/perso3/crepes-29.gif, et que l’espérance de cette variable est http://freakonometrics.blog.free.fr/public/perso3/crepes-30.gif, i.e. http://freakonometrics.blog.free.fr/public/perso3/crepes-31.gif. Donc cette modélisation devrait parfaitement convenir à notre problème. Ce qui est intéressant, c’est que comme on travaille avec la loi conjuguée, la loi ex post est encore une loi Beta. Autrement dit, ici, si je suppose que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit a priori une loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif, et que la loi de http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif sachant http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit une loi binomiale http://freakonometrics.blog.free.fr/public/perso3/crepes-32.gif, alors la loi a posteriori de http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif sachant http://freakonometrics.blog.free.fr/public/perso3/crepes-33.gif est une loi http://freakonometrics.blog.free.fr/public/perso3/crepes-34.gif (c’est ce que l’on avait vu ici). Et en particulier

http://freakonometrics.blog.free.fr/public/perso3/crepes-35.gif

Autrement dit on obtient une combinaison convexe entre la moyenne de la loi a priori (à gauche) et la fréquence empirique (à droite).
On retrouve ici ce que nous avions noté auparavant, en notant que http://freakonometrics.blog.free.fr/public/perso3/crepes-36.gif, et http://freakonometrics.blog.free.fr/public/perso3/crepes-37.gif.
On peut toutefois se demander si le résultat que l’on avait obtenu sur la loi de S_n reste valide ici…
Notons que

http://freakonometrics.blog.free.fr/public/perso3/crepes-38.gif

i.e. (moyennant quelques calculs)

http://freakonometrics.blog.free.fr/public/perso3/crepes-40.gif

Dans le cas où http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif, en réécrivant aussi le premier terme, on a

http://freakonometrics.blog.free.fr/public/perso3/crepes-41.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-42.gif

On retrouve bien la formule obtenue auparavant… Ce qui est intéressant aussi, c’est que si on note que

http://freakonometrics.blog.free.fr/public/perso3/crepes-43.gif

or

http://freakonometrics.blog.free.fr/public/perso3/crepes-44.gif

et donc

http://freakonometrics.blog.free.fr/public/perso3/crepes-45.gif

et par récurrence, on obtient que

http://freakonometrics.blog.free.fr/public/perso3/crepes-46.gif

et que de plus

http://freakonometrics.blog.free.fr/public/perso3/crepes-47.gif

Autrement dit, effectivement, si j’en ai raté beaucoup au début, ma probabilité d’en rater beaucoup est grande (et elle est d’ailleurs constante). Certains parleraient de no-free lunch.

Découpage équitable de la pizza avec des olives

mais non papa, tu comprends jamais rien ! ok, c’est bien si on a autant de pizza tous les deux… mais on veut aussi avoir autant d’olives tous les deux ! et autant de morceaux d’ananas !”. Ouais, bon, on n’est pas couché là !

On va quand même essayer de la couper cette pizza (avant qu’elle ne refroidisse) !

  • Partage de la pizza et des olives

Commençons par supposer que seules les olives soient importantes. La première idée que j’aurais serait de faire une régression médiane, qui va découper le nuage de points (les olives) en deux…

Sous R (oui, j’ai installé R sur la table de la cuisine), c’est assez facile à mettre en œuvre… Pour des raisons de symétrie, je vais faire les régressions de l’abscisse sur l’ordonnée, et réciproquement…

> library(quantreg)
> lq=rq(y ~ x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 22
> sum(y0<=predict(lq))
[1] 22
>
> lq=rq(x ~ y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients[2],
+ a=-lq$coefficients[1]/lq$coefficients[2],lwd=5,col="red")
> sum(x0>predict(lq))
[1] 22
> sum(x0<=predict(lq))
[1] 22

Le soucis est que si la régression n’est pas contrainte, il y a une répartition 50-50 des olives, mais les parts ne sont pas équitables,

> library(quantreg)
> lq=rq(y ~ 0+x, tau = 0.5, data=D)
> abline(lq,lwd=5,col="red")
> sum(y0>predict(lq))
[1] 20
> sum(y0<=predict(lq))
[1] 24
>
> lq=rq(x ~ 0+y, tau = 0.5, data=D)
> abline(b=1/lq$coefficients,a=0,lwd=5,col="red")
> sum(x0>predict(lq))
[1] 25
> sum(x0<=predict(lq))
[1] 19

et si je supprime la constante les parts sont équitables (la régression passe par le centre la pizza), mais le partage des olives n’est pas équilibré…

On peut aussi calculer pour tous les angles possibles ce que donnerait des partages de la pizza (qui passent par le centre). C’est ce que représente la courbe rouge ci-dessous: l’extérieur du disque le ratio des olives dépasse 1, à l’intérieur il est inférieur à 1. Les découpages optimaux correspondent au cas où le ratio vaut précisément 1… ce qui est possible si on coupe comme indiqué ci-dessous, à droite,

Je pense qu’on doit pouvoir montrer qu’il est toujours possible de couper la pizza en deux, de manière équitable, avec autant d’olives de part et d’autre (s’il y a un nombre pair d’olives).
oui mais papa… tu nous avais dit que tu savais pas couper en passant par le centre…”. Ah oui, c’est vrai ! Damned…. bon, utilisons la méthode décrite dans le billet précédant (ici), a priori ça doit être possible de trouver un sommet pour lequel le partage des olives est équitable. En fait, on pourrait aussi tirer au hasard des points dans la pizza, et compter les olives en faisant la répartition de manière alternée,

> sommet=matrix(NA,500000,3)
> for(s in 1:nrow(sommet)){
+ a0=runif(10)*2-1; b0=runif(10)*2-1
+ if(a0^2+b0^2)<=1){
+ slope1=(y0-b0)/(x0-a0)
+ I1=((slope1>0)&(slope1<1))|((slope1<(-1)))
+ sommet[s,]=c(sum(I1),a0,b0)
}}

On a alors en rouge l’ensemble des points à partir desquels on peut faire des découpages en 8, avec des pointes à 45° (en commençant par couper horizontalement)

Et ça semble marcher !

oui mais papa… tu as oublié les ananas…”.

  • Partage de la pizza et des olives et des ananas

En utilisant la méthode précédente, rien de plus simple ! En plus de compter la répartition des olives, on va compter la répartition des ananas…! Sur le dessin à gauche on a les olives et les ananas. A droite, les points rouges correspondent à des répartitions équitables des olives, et les bleus des répartitions équitables des ananas… Reste à trouver des points à la fois rouges et bleus,

ce sont finalement les points violets ci-dessous ! On a alors une méthode pour couper la pizza en deux parts de même taille, avec autant d’olives et autant d’ananas de part et d’autres !

Découpage équitable de la pizza

Bon, maintenant que l’année scolaire est passée, je vais arrêter un peu les billets sérieux pour revenir pour mes mésaventures à la maison. Par exemple, hier soir, c’était pizza. Et mes enfants aiment partager à deux une pizza (sur laquelle on a rajouté plein de choses qu’ils aiment, genre des olives, ou des morceaux d’ananas). Mais il faut un découpage équitable. Les enfants ne veulent pas un partage juste, mais un partage moitié-moitié.

Classiquement, je coupe la pizza en deux (cf figure ci-contre, avec un côté en bleu pour mon fils, et un en rouge pour ma fille, histoire de simplifier un peu). Mais comme me l’a fait noté mon fils, si mon couteau ne passe pas par le centre de la pizza (on va dire le centre du cercle parce qu’il va falloir faire des maths… pour les recettes de cuisine, il faut changer de site, genre ici).
Donc pour faire simple, la proba d’avoir deux parts rigoureusement égale est très très faible… car il faut couper suivant une droite qui passe exactement au centre (sinon, il y a quelqu’un qui est défavorisé, alors on tire au sort… a prioric’est juste juste, mais ex-post, il y en a toujours un ou une qui trouve le système injuste et qui veut qu’on retire au hasard…)
Fort heureusement, il existe une technique simple, et rigoureusement juste, permettant d’avoir des parts équitables…
L’idée est la suivante: on commence par couper horizontalement, n’importe où. Ensuite, on coupe verticalement (à angle droit, en prenant son équerre) de manière à récupérer 4 parts. Au point d’intersection,  on coupe à 45° direction nord-est, et 45° direction nord-ouest…. Une fois le découpage effectué, on réparti les parts, à tour de rôle, les parts bleues pour mon fils, et celles en rouges pour ma fille,

Car le truc génial est que le découpage – en plus d’être fort élégant – est parfaitement équitable ! Comme les enfants sont toujours sceptiques, il faut bien sûr leur prouver que oui, ils en auront exactement pareil. Pour ça on reprend le couteau…

Si on recoupe les parts suivant les traits-ci-dessus, chaque part bleue pourra être superposée à une part rouge (la preuve mathématique repose sur des simple notions de symétrie…). Génial non ? On a même le théorème du mardi soir, à savoir que “si on coupe une pizza en 8 parts, avec des pointes à 45° partir d’un point choisi arbitrairement dans la pizza, alors la somme alternée des parts sera égale” (cf ici). En 1967, Upton a posé ce problème,

Larry Carter et Stan Wagon affirmaient que ce résultat resterait vrai si on coupait en 12 ou 16 parts (avec des angles égaux pour les pointes), mais pas pour 6, 10, 14 ou 18 par exemple. Le blog du chou romanesco (ici et ) prétend que c’est vrai pour les découpages en un nombres pairs. En fait, Rick Mabry et Paul Deiermann ont un théorème assez général ( pour la version complète) et qui généralise aussi en dimension 3 (pour une calzone en quelque sorte). Mais je sens que ça n’intéresse plus trop les enfants, qui ont déjà fini la pizza….

Les plus curieux peuvent aussi consulter le court article ci-dessous,

Comment parle un président français le 31 décembre ?

Il y a pas mal d’agitation présidentielle sur les blogs (français) cette semaine, que ce soit au sujet des 30 ans de l’élection de François Mitterrand, ou l’annonce de la date de la prochaine élection présidentielle (et probablement d’autres raisons). Un peu naïvement je me suis demandé si les présidents étaient vraiment si différents que ça ? Ou plus précisément, est-ce qu’une fois investis par la fonction, tous ces hommes ne se comporteraient-ils pas – finalement – tous pareils ?

Pour cela je me suis simplement demandé si on pouvait distinguer les discours des différents Présidents. J’ai donc pris quelques discours entendus les 31 décembre, lors de vœux retransmis à la télévision (en ligne ici, en excluant les discours aux forces armées, qui me semblaient biaiser un peu mon échantillon). On commence par récupérer la liste de tous les mots prononcés (les fichiers csv sont en ligne ici)

setwd("/home/arthur/Documents/r-data/voeux/")
FILES=list.files()
library(tau)
BASE=scan(FILES[1],what="character")
WORD=fixEncoding(BASE,latin1=TRUE)
 for(i in 2:length(FILES)){
BASE=scan(FILES[i],what="character")
 WORD=c(WORD,fixEncoding(BASE,latin1=TRUE))
 }

Ensuite on regarde le top 500 (j’exclue les mots communs genre le ou la, mais je garde les pronoms, comme mon ou notre),

> WORD=toupper(WORD)
> LW=rev(sort(table(WORD)))
> LISTE=names(LW[1:500])
> LISTE=LISTE[-c(1,2,3,4,5,6,7,8,9,10,11,14,15,18,19,20,21,22,23,25,
               27,43,40,36,41,60)]
> LISTE[1:25]
 [1] "NOUS"   "JE"     "VOUS"   "NOTRE"  "PLUS"   "FRANCE" "PAS"    "NOS"   
 [9] "EST"    "AVEC"   "SE"     "MES"    "CEUX"   "MAIS"   "SONT"   "OU"    
[17] "CETTE"  "AUSSI"  "TOUS"   "SUR"    "SI"     "PAYS"   "LEUR"   "FAIRE" 
[25] "VOTRE"

et créer une base de décompte, par discours,

ANNEE=as.numeric(substr(FILES,6,7))
AN=1900+ANNEE*(ANNEE>15)+(100+ANNEE)*(ANNEE<15)
PRES=c("CHIRAC","CHIRAC","CHIRAC","SARKOZY","SARKOZY",
"SARKOZY","VGE","VGE","VGE","MITTERRAND","MITTERRAND",
"MITTERRAND","MITTERRAND","MITTERRAND")
NOM=paste(PRES,AN,sep="-")
COMPTE=matrix(NA,length(FILES),length(LISTE))
for(i in 1:length(FILES)){
BASE=scan(FILES[i],what="character")
WORD=toupper(fixEncoding(BASE,latin1=TRUE))
for(j in 1:length(LISTE)){
COMPTE[i,j]=sum(WORD==LISTE[j])
}}
colnames(COMPTE)=LISTE
rownames(COMPTE)=NOM

On peut alors faire tourner des procédures d’analyse de données, comme une ACP

library("FactoMineR")
res.pca = PCA(COMPTE, scale.unit=TRUE, ncp=5, graph=FALSE)
plot.PCA(res.pca, axes=c(1, 2), choix="ind")

Les projections des individus (les discours) sont les suivantes

alors que celles pour les variables (les mots) sont, toujours sur les deux premiers axes factoriels,

On peut jouer des heures à se demander qui dit quoi, mais ce qui m’intéressait était de voir si les présidents étaient proches, en terme de discours. On peut alors faire des classifications hiérarchiques,

DISTANCE = dist(COMPTE)
h = hclust(DISTANCE) 
plot(as.dendrogram(h), horiz = TRUE)

ou ici avec un autre algorithme de classification

On notera les regroupements dépendent peu des distances et des méthodes de constitution de classes… et que les regroupements se font par président: tous les discours de François Mitterrand semblent proches, et tous ceux de Jacques Chirac aussi (en tous les cas ceux que j’ai téléchargé, mais l’analyse n’est pas exhaustive). On notera aussi la proximité entre Valéry Giscard d’Estaing l’année de son élection, et Nicolas Sarkozy l’année de son élection….

C’est finalement un peu rassurant de voir que les Présidents ne sont pas tous les mêmes, et qu’on peut les reconnaitre en fonction de leur discours… c’est tout du moins la prochaine étape, ou on va voir si on peut affecter des nouveaux discours (que je n’ai pas utilisés pour l’analyse) dans des classes, et repérer qui l’a écrit…

Multivariate probit regression using (direct) maximum likelihood estimators

Consider a random pair http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-01.gif of binary responses, i.e. http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-02.gif with http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-03.gif taking values 1 or 2. Assume that probability http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-04.gif can be function of some covariates.

  • The Gaussian vector latent structure

A standard model is based a latent Gaussian structure, i.e. there exists some random vector http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-06.gif such that http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-07.gif if http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-08.gif is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that

http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-09.gif

where we can assume that http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-10.gif is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-11.gif.

> logV=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2)
+     a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
+ for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
+     a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))}
+     a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     1*log(a01) +
+     2*log(a10) +
+     3*log(a00) )
+ }
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par

(the code is a bit long since I had trouble working properly with matrices – or more precisely to vectorize my functions – so I used loops… I am sure it is possible to write a better code).
It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators,

> library(mnormt)
> set.seed(1)
> n=1000
> r=0.5
> X1=runif(n)
> X2=rnorm(n)
> Y1S=1+5*X1
> Y2S=8-5*X1
> RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2))
> YS=cbind(Y1S,Y2S)+RES
> Y1=(YS[,1]>quantile(YS[,1],.5))*1
> Y2=(YS[,2]>quantile(YS[,2],.5))*1
> base=data.frame(i,Y1,Y2,X1,X2,YS)
> head(base)
  i Y1 Y2        X1          X2      Y1S      Y2S
1 1  0  0 0.2655087  0.07730312 3.177587 5.533884
2 2  0  0 0.3721239 -0.29686864 1.935307 5.089524
3 3  1  0 0.5728534 -1.18324224 4.757848 5.172584
4 4  1  0 0.9082078  0.01129269 4.600029 3.878225
5 5  0  1 0.2016819  0.99160104 2.547362 6.743714
6 6  1  0 0.8983897  1.59396745 5.309974 4.421523

(the two columns on the right are latent observations, that cannot be used since theoretically they are unobservable). Note that it is a simple regression, one of the component is here only to bring some noise. First of all, let us look at marginal probit regression

>  reg1=glm(Y1~X1+X2,data=base,family=binomial)
>  reg2=glm(Y2~X1+X2,data=base,family=binomial)
> summary(reg1)
 
Call:
glm(formula = Y1 ~ X1 + X2, family = binomial, data = base)
 
Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.90570  -0.50126  -0.00266   0.49162   2.78256
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.291725   0.267149 -16.065   <2e-16 
X1           8.656836   0.510153  16.969   <2e-16 ***
X2           0.007375   0.090530   0.081    0.935
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386.29  on 999  degrees of freedom
Residual deviance:  726.48  on 997  degrees of freedom
AIC: 732.48

Number of Fisher Scoring iterations: 5
> summary(reg2)
Call:
glm(formula = Y2 ~ X1 + X2, family = binomial, data = base)
Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.74682  -0.51814  -0.00001   0.57969   2.58565
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.91709    0.24399  16.054   <2e-16 ***
X1          -7.89703    0.46277 -17.065   <2e-16 ***
X2           0.18360    0.08758   2.096    0.036 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1386.29  on 999  degrees of freedom
Residual deviance:  777.61  on 997  degrees of freedom
AIC: 783.61
Number of Fisher Scoring iterations: 5

Here, the optimization yields,

> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
> OPT[1]
[1] 0.5261382
> matrix(OPT[2:7],2,3)
          [,1]      [,2]       [,3]
[1,] -2.451721  4.908633 0.01600769
[2,]  2.241962 -4.539946 0.10614807

Note that the coefficients we have obtained are almost identical to the ones obtained with R standard procedure,

>  library(Zelig)
>  REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
>  summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
    rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
                 Min        1Q     Median      3Q     Max
probit(mu1) -10.5442 -0.377243  0.0041803 0.36709 8.60398
probit(mu2)  -7.8547 -0.376888  0.0083715 0.42923 5.88264
rhobit(rho) -13.8322 -0.091502 -0.0080544 0.37218 0.85101
 
Coefficients:
                  Value Std. Error   t value
(Intercept):1 -2.451699   0.135369 -18.11116
(Intercept):2  2.241964   0.125072  17.92536
(Intercept):3  1.169461   0.189771   6.16249
X1:1           4.908617   0.252683  19.42602
X1:2          -4.539951   0.233632 -19.43203
X2:1           0.015992   0.050443   0.31703
X2:2           0.106154   0.049092   2.16235
 
Number of linear predictors:  3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
&n
bsp;
Dispersion Parameter for binom2.rho family:   1
 
Residual Deviance: 1460.355 on 2993 degrees of freedom
 
Log-likelihood: -730.1774 on 2993 degrees of freedom
 
Number of Iterations: 3

> matrix(coefficients(REG)[c(1:2,4:7)],2,3)
          [,1]      [,2]       [,3]
[1,] -2.451699  4.908617 0.01599183
[2,]  2.241964 -4.539951 0.10615443

The correlation here is also the same

> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5260951

That procedure works well an can be extended to ordinal responses (not only binary ones, or to three dimensional problems,

logV=function(beta){
BETA=matrix(beta[4:(3+ncol(Y)*ncol(X))],ncol(Y),ncol(X))
z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]),X%*%(BETA[3,]))
r12=beta[1]
r23=beta[2]
r31=beta[3]
s1=s2=s3=1
sigma=matrix(c(s1^2,r12*s1*s2,r31*s1*s3,
               r12*s1*s2,s2^2,r23*s2*s3,
               r31*s1*s3,r23*s2*s3,s3^2),3,3)
sigma1=matrix(c(s2^2,r23*s2*s3,
                r23*s2*s3,s3^2),2,2)
sigma2=matrix(c(s1^2,r31*s1*s3,
                r31*s1*s3,s3^2),2,2)
sigma3=matrix(c(s1^2,r12*s1*s2,
                r12*s1*s2,s2^2),2,2)
    a111=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a111=c(a111,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a011=pmnorm(z[1,2:3],varcov=sigma1)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a011=c(a011,pmnorm(z[i,2:3],varcov=sigma1)-pmnorm(z[i,],varcov=sigma))}
    a101=pmnorm(z[1,c(1,3)],varcov=sigma2)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a101=c(a101,pmnorm(z[i,c(1,3)],varcov=sigma2)-pmnorm(z[i,],varcov=sigma))}
    a110=pmnorm(z[1,1:2],varcov=sigma3)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a110=c(a110,pmnorm(z[i,1:2],varcov=sigma3)-pmnorm(z[i,],varcov=sigma))}
    a100=pnorm(z[1,1],sd=s1)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a100=c(a100,pnorm(z[i,1],sd=s1)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a010=pnorm(z[1,2],sd=s2)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(2,3)],varcov=sigma1)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a010=c(a010,pnorm(z[i,2],sd=s2)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(2,3)],varcov=sigma1)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a001=pnorm(z[1,3],sd=s3)-pmnorm(z[1,c(2,3)],varcov=sigma1)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a001=c(a001,pnorm(z[i,3],sd=s3)-pmnorm(z[i,c(2,3)],varcov=sigma1)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
    a000=1-a111-a011-a101-a110-a001-a010-a100
 
a111[a111<=0]=1e-50
a110[a110<=0]=1e-50
a101[a101<=0]=1e-50
a011[a011<=0]=1e-50
a100[a100<=0]=1e-50
a010[a010<=0]=1e-50
a001[a001<=0]=1e-50
a000[a000<=0]=1e-50
 
-sum(((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==0))*log(a111) +
    4*log(a011) +
    5*log(a101) +
    6*log(a110) +
    7*log(a001) +
    8*log(a010) +
    9*log(a100) +
    10*log(a000) )
}

A strong assumption in that bivariate model is that residuals have a Gaussian structure. It is possible to change that assumption

  • marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
  • in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel’s copula (also called the bivariate logistic copula), or Clayton’s

Here, the following code can be used to extend the model to non Gaussian structures,

> F=function(x,r){pmnorm(x,rep(0,length(x)),
+                 varcov=matrix(c(1,r,r,1),2,2))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
> 
> logVgen=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+     a11=F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))}
+     a10=Fx(z[1,1])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))}
+     a01=Fy(z[1,2])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     11*log(a01) +
+     12*log(a10) +
+     13*log(a00) )
+ }
>
> beta0=c(0,1,1,1,1,1,1)
> (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par)
[1]  0.52613820 -2.45172059  2.24196154  4.90863292 -4.53994592  0.01600769
[7]  0.10614807
There were 23 warnings (use warnings() to see them)

E.g.

> library(copula)
> F=function(x,r){pcopula(pnorm(x),
               claytonCopula(2, r))}
> Fx=function(x1){F(c(x1,1e40),0)
}
> Fy=function(x2){Fx(x2)}
  • An application to school tests

Consider the following dataset,

hsb2=read.table("http://freakonometrics.free.fr/hsb2.csv",
        header=TRUE, sep=",")
math_male=hsb2$math[female==0]
write_male=hsb2$write[female==0]
math_female=hsb2$math[female==1]
write_female=hsb2$write[female==1]
plot(math_female, write_female, type="p",
     pch=19,col="red",xlab="maths",ylab="writing",cex=.8)
points(math_male, write_male, cex=1.2, col="blue")

with here maths versus writing, with girls in red and boys in blue, where variables here are

  female :
    0: male
    1: female
  race :
    1: hispanic
    2: asian
    3: african-amer
    4: white
  ses :
    1: low
    2: middle
    3: high
  schtyp : type of school
    1: public
    2: private
  prog : type of program
    1: general
    2: academic
    3: vocation
  read : reading score
  write : writing score
  math : math score
  science : science score
  socst : social studies score

We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ?

Here the code is simply

> W=hsb2$write>=50
> M=hsb2$math>=50
> base=data.frame(Y1=W,Y2=M,
+             X1=hsb2$female,X2=hsb2$read)
>
> library(Zelig)
> REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
> summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
    rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
                Min        1Q  Median      3Q    Max
probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592
probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072
rhobit(rho) -4.9821  0.010481 0.13500 0.40776 2.9171
 
Coefficients:
                  Value Std. Error  t value
(Intercept):1 -5.484711   0.787101 -6.96825
(Intercept):2 -4.061384   0.633781 -6.40818
(Intercept):3  1.332187   0.322175  4.13497
X1:1           1.125924   0.233550  4.82092
X1:2           0.167258   0.202498  0.82598
X2:1           0.103997   0.014662  7.09286
X2:2           0.082739   0.012026  6.88017
 
Number of linear predictors:  3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
 
Dispersion Parameter for binom2.rho family:   1
 
Residual Deviance: 364.51 on 593 degrees of freedom
 
Log-likelihood: -182.255 on 593 degrees of freedom
 
Number of Iterations: 3
> (exp(summary(REG)@coef3[3])-1)/(exp(
summary(REG)@coef3[3])+1)
[1] 0.5824045

with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here

> beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1),
+      summary(REG)@coef3[c(1:2,4:7),1])
> beta0
              (Intercept):1 (Intercept):2          X1:1          X1:2
0.58240446   -5.48471133   -4.06138412    1.12592427    0.16725842
X2:1          X2:2
0.10399668    0.08273879
> (OPT=optim(fn=logV,par=beta0,method="BFGS")$par)
(Intercept):1 (Intercept):2          X1:1          X1:2
0.5824045    -5.4847113    -4.0613841     1.1259243     0.1672584
X2:1          X2:2
0.1039967     0.0827388

i.e. we obtain (almost) exactly the same estimators. But here I have used as starting values for the optimization procedure the estimators given by R. If we change them, hopefully we have a robust maximum likelihood estimator,

> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")$par)
              (Intercept):1 (Intercept):2          X1:1          X1:2
   0.58233360   -5.49428984   -4.06839571    1.12696594    0.16760347
         X2:1          X2:2
   0.10417767    0.08287409
There were 12 warnings (use warnings() to see them)

So once again, it is possible to optimize numerically a likelihood function, and it works.

  1. Y[,1]==0)&(Y[,2]==1 []
  2. Y[,1]==1)&(Y[,2]==0 []
  3. Y[,1]==0)&(Y[,2]==0 []
  4. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 []
  5. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 []
  6. Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 []
  7. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 []
  8. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 []
  9. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 []
  10. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 []
  11. Y[,1]==0)&(Y[,2]==1 []
  12. Y[,1]==1)&(Y[,2]==0 []
  13. Y[,1]==0)&(Y[,2]==0 []

"sendo l'intento mio scrivere cosa utile a chi la intende…"