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).
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 là. Si quelqu’un voit une coquille, merci de me signaler avant que le livre ne parte à l’éditeur…
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 là, 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,
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) 
puis 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.
Considérons un simple jeu de pile ou face. On suppose que l’on tombe sur pile avec une probabilité 
(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 .
Supposons que l’on connaisse parfaitement . A priori deux stratégies sont possibles:
et face avec probabilité 
.
, on prédit pile, sinon on prédit face (à chaque lancer si on fait plusieurs lancers consécutifs).Notons 
la variable correspondant à la valeur prise au lancer, disons 1 si on tomber sur pile et 0 sur face. Et notons 
la valeur prédite. On gagne si 
, et la probabilité de gagner est alors
i.e. pour la première stratégie
En revanche pour la seconde, la probabilité de gagner est
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,
 |
 |
. Considérons le polynôme 
, dont les racines sont 
et 
, i.e sur 
,
ou encore
i.e.
(avec une inégalité stricte si 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).
Supposons maintenant que la probabilité soit inconnue, mais estimable. Autrement dit, on a assisté à 
lancers, on a un échantillon 
et on estime la probabilité de tomber sur pile par
(qui tire une loi uniforme sur le segment unité), pour la prochaine valeur, on prédit
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 
, et la probabilité de gagner est alors
Dans le premier cas, la probabilité de gagner est
(comme auparavant). En revanche, dans le second cas, le calcul de la probabilité de gagner est
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
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…
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
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 
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 
réacteurs sur une durée de 
années est
Comme est supposé petit, un développement limité donnerait
Malheureusement, dans les développements limités de ce genre, s’il faut que 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
(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….
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 
la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant 
. Supposons qu’il existe un processus latent sous-jacent, 
tel que
A une date , on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut
qui peut encore s’écrire
ou, ce qui sera pleut être plus utile ensuite
Avec notre modèle latent, cette probabilité peut s’écrire
c’est à dire en notant 
la fonction de répartition (jointe) du vecteur aléatoire 
cette probabilité s’écrit
Or Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout ,
où
et
(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 pour. 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.
où 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
est alors
On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de ,
Avec ce formalisme, on peut alors en déduire 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 
(qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation 
(qui est associée à la dynamique du processus). Pour simplifier, si on prend 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 à
qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,
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
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ù 
1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si -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.
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 
de services possibles. A son premier service, il choisit le service 
. La première difficulté est que le service soit correct (ou passe le filet pour faire simple), et on note 
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 
. Bien que ne connaissant pas vraiment le tennis, il n’est pas absurde de supposer qu’il y ait une forme d’arbitrage,
petit), mais s’il passe le filet le serveur sera à son avantage, et la probabilité de gagner le point sera élevée ( grand) grand), mais la balle sera alors facile pour l’adversaire, et la probabilité de gagner le point sera faible ( petit)A un service 
donné, la probabilité de gagner le point (sur ce premier service) est alors 
.
Supposons que soit “ordonné” au sens où 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
soit à valeur dans 
, avec une fonction croissante pour 
, puis décroissante pour 
.
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 
correspondant à son premier et second service, respectivement. La probabilité de gagner le point, à stratégie 
donnée, est alors
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, à 
donné
la condition du premier ordre est
i.e
Autrement dit, le second service doit être celui qui maximise sa probabilité de gagner le point, i.e. 
.
Mais si on sait qu’il est optimal de jouer , dans ce cas, pour le premier service, le joueur cherche à résoudre
dont la condition du premier ordre est
i.e.
Le point s_1 qui vérifie cette condition peut se visualiser sur la figure ci-dessous, car 
s’interprète comme 
, i.e. la tangente de la fonction rouge, et on veut que cette tangente soit de pente 
, 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…
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)
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é
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é
et encore une fois, on a un impact presque nul du changement d’heure,
 |
 |
J’avais commencé plusieurs analyse (ici ou là) 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…
