Pratique du calcul bayésien (Statistique et probabilités appliquées) [1st Edition.] 2287996664, 9782287996665 [PDF]

Zitiervorschau

Pratique du calcul bayesien

Springer Paris Berlin Heidelberg New York Hong Kong Londres Milan Tokyo

Jean-Jacques Boreux " Eric Parent Jacques Bernier

Pratique du calcul bayesien

~ Springer

Jean-Jacques Boreux

Universite de Liege (ULg) Departement des Sciences et Gestion de l' environnement 185, avenue de Longwy 6700 Arlon Belgique

Eric Parent

AgroParisTech 16, rue Claude-Bernard 75231 Paris Cedex 05

Jacques Bernier

Le Pech-de-Biaud 24250 Saint-Martial-de-Nabirat

ISBN-13 : 978-2-287-99666-5 Springer Paris Berlin Heidelberg New York © Springer-Verlag France, Paris, 2010

Imprime en France

Springer- Verlag France est membre du groupe Springer Science + Business Media Cet ouvrage est soumis au copyright. Tous droits reserves, notamment la reproduction et la representation, la traduction, la reimpression, I'expose, la reproduction des illustrations et des tableaux, la transmission par voie d'enregistrement sonore ou visuel, la reproduction par microfilm ou tout autre moyen ainsi que la conservation des banques de donnees. La loi francaise sur le copyright du 9 septembre 1965 dans la version en vigueur n' autorise une reproduction integrale ou partielle que dans certains cas, et en principe moyennant le paiement de droits. Toute representation, reproduction, contrefacon ou conservation dans une banque de donnees par quelque precede que ce soit est sanctionnee par Ia loi penale sur Ie copyright. L'utilisation dans cet ouvrage de designations, denominations commerciales, marques de fabrique, etc. meme sans specification ne signifie pas que ces termes soient libres de la legislation sur les marques de fabrique et Ia protection des marques et qu' ils puissent etre utilises par chacun. La maison d' edition decline toute responsabilite quant a I'exactitude des indications de dosage et des modes d' emploi, Dans chaque cas, il incombe a I' usager de verifier les informations donnees par comparaison ala litterature existante.

Maquette de couverture : Jean-Francois Montmarche

Collection

Statistique et probabilites appliquees dlrigee par Yadolah Dodge Professeur Honoraire Universite de Neuchatel Suisse [email protected]

Comlte editorial : Christian Genest

Stephan Morgenthaler

Departement de Mathematiques et de statistique UniversiteLaval Quebec GIK 7P4 Canada

Ecole Polytechnique Federale de Lausanne Departement des Mathematiques 1015 Lausanne Suisse

Marc Hallin

Gilbert Saporta

Universite libre de Bruxelles Campus de la Plaine CP 210 1050 Bruxelles Belgique

Conservatoire national des arts et metiers 292, rue Saint-Martin 75141 Paris Cedex 3 France

Ludovic Lebart

Telecom-Paris'Iech 46, rue Barrault 75634 Paris Cedex 13 France

Dans la meme collection : - Statistique. La theorie et ses applications

Michel Lejeune, avril 2004 - Optimisation appliquee

Yadolah Dodge, octobre 2004 - Le choix bayesien. Principes et pratique

Christian P. Robert, novembre 2005 - Maitriser l' aleatoire. Exercices resolus de probabilites et statistique

Eva Cantoni, Philippe Huber, Elvezio Ronchetti, novembre 2006

- Regression. Theorie et applications Pierre-Andre Cornillon, Eric Matzner-Lober, janvier 2007 - Le raisonnementbayesien. Modelisation et inference Eric Parent, Jacques Bernier,juillet 2007 - Premiers pas en simulation Yadolah Dodge, Giuseppe Melfi, juin 2008 - Genetique statistique Stephan Morgenthaler, juillet 2008 - Maitriser l'aleatoire. Exercices resolusde probabiliteet statistique, deuxieme edition Eva Cantoni, Philippe Huber, Elvezio Ronchetti, septembre 2009

Preface Le troisicme millenaire sera, dit-on, celui de l'information. Aussi la statistique y sera-t-elle appelee a jouer un role important et le paradigme bayesien plus que tout autre, puisqu'il offre un cadre de raisonnement bien adapte a I'integration des opinions et des faits de toutes provenances qui interviennent dans la gestion des risques et la prise de decision en contexte d'incertitude. De la collecte de donnees a la prevision, l'analyse statistique pose plusieurs defis. L'elaboration du modele rcprcsentc sans doute la phase la plus delicate de l'exercice, car elle doit repondre a un double imperatif de realisme et de parcimonie. Hormis quelques cas de figure, une demarche bayesienne n'est envisageable qu'a charge de disposer d'outils efficaces pour la quantification et la mise a jour de l'information. Jouissant d'une expertise considerable dans le dornaine, les auteurs avaient deja brosse un tableau du Traitement bauesien de l'incertitude en sciences de l'environnement dans un ouvrage paru en 2000. Six ans plus tard, Christian Robert publiait chez Springer Le choix bayesien - Principes et pratique, expose des fondements de la theorie qu'Eric Parent et Jacques Bernier completaient plaisamment en 2007 avec Le raisonnement btnjesien - Modelisatiori et inference, paru dans la meme collection. Aujourd'hui, pour notre plus grand plaisir, Jean-Jacques Boreux, Eric Parent et Jacques Bernier joignent a nouveau leurs forces pour nous instruire dans la Pratique du calcul bayesien. A l'aide d'exemples concrets, nombreux et varies, ils nous initient a la construction de modeles bayesiens et au maniement de l'imposant arsenal de calcul necessaire a leur mise en oeuvre. Au passage, ils s'efforcent aussi d'aiguiser notre esprit critique! De l'halieutique a l'hydrometeorologie, en passant par la mesure des risques d 'avalanche, de pneumoconiose ou de pollution en milieu clos, les auteurs decortiquent et analysent pour nous divers jeux de donnees issus de la pratique. Partant de series temporelles, de valeurs extremes ou d'effectifs de capturerecapture, ils nous montrent tantot comment decrire des relations entre plusieurs variables au moyen de graphes acycliques orientes, tant6t comment batir ou affiner des modeles lineaires, generalises ou hierarchiques definis par conditionnements successifs. A l'occasion, ils font aussi appel au logiciel WinBUGS pour illustrer le calcul de lois a posteriori au moyen de l'algorithme de Metropolis-Hastings ou de techniques particulaires dernier cri.

Vlll

Pratique du calcul bayesien

Dans un souci didactique evident, les auteurs ont menage une gradation dans le degre de complexite des problemes etudies, Les premiers chapitres abordent des cas relativement simples, faciles a resoudre et bien adaptes a l'apprentissage des rudiments; les enseignants s'en inspireront avec bonheur. Les applications « grandeur nature» presentees en seconde partie font quant a elles un abondant usage de structures hierarchiques, de variables latentes et autres savantes constructions; le savoir-faire statistique et le genie du calcul numerique y apparaissent ici dans toute leur splendeur. Pour reprendre l'aimable locution des auteurs, le lecteur est ainsi progressivement amene « de la plume a la souris» et il en ressort ebloui et grandi. Gageons que specialistes et utilisateurs de la statistique s'approprieront rapidement ce beau livre et qu'ils reconnaitront en lui un guide sur et accessible des principes modernes du calcul bayesian. Bonne lecture! Christian Genest, professeur Universite Laval, Quebec President sortant de la Societe statistique du Canada et de l' Association des statisticiennes et statisticiens du Quebec

Avant-propos L'anticipation est une composante essentielle des capacites d'adaptation d 'une societe et la statistique peut etre definie comme « l' art de raisonner de facon quantitative en avenir incertain ». Elle intervient dans toutes les disciplines scientifiques OU se melent savoir et donnees. Elle est done utilisee par les physiciens, les economistes, les ingenieurs, les geographes, les biologistes, les assureurs, les psychologues, les metcorologues, les gestionnaires d'entreprises, etc., bref, par tous les praticiens soucieux de batir sur des fondations solides un pont entre theorie et donnees experimentales. Comme dans toutes les disciplines scientifiques, il faut d'emblee fixer le niveau qu'on se propose d'atteindre. II nous semble que quatre niveaux suffisent a preciser les compctcnces. - Comme son qualificatif l'indique, le niveau elementaire est une prise de contact avec la discipline en question. S'agissant de la statistique, l'etudiant saisit le sens general de la modelisation probabiliste, connait les distributions de base et est autonome dans des situations simples. - Le niveau suivant vise une qualification operationnelle, Ici I'etudiant est capable de construire un modele qui repond a un questionnement. II manie les outils modernes de l'inference statistique, interprete et critique lcs resultats obtenus. - Le niveau suivant est la maitrise des concepts mathematiqucs qui justifient les procedures utilisees, A ce niveau, le statisticien fait preuve d'une tres grande creativite, comprend pourquoi une procedure faillit et sait y remedier. - Enfin, le quatrieme niveau est celui de la recherche fondamentale qui, par definition, introduit des nouvelles idees et./ou generalise des concepts existants sans avoir neccssaircmcnt de visec operationnelle au moment des travaux. Bien sur, il n'existe pas de separation nette entre ces niveaux mais, pour cet ouvrage, notre ambition est clairement une qualification operaiumnelle en statistique baueeiemic avec, peut-etre, quelques incursions au niveau maitrise. Avant de preciser cette ambition, il nous semble utile de remonter aux origines de ce livre. Construire un modele statistique paromeirique a des fins decisionncllcs, c 'est oser avoir tort en maximisant ses chances d'avoir raison! Cette repartie

x

Pratique du calcul bayesien

vint un jour a l'esprit du premier auteur face a des etudiants en sciences de l'environnement, inquiets et perplexes. II faut bien le reconnaitre, la plupart des etudiants redoutent le cours de stat, notamment a cause du langage mathematique qui Ie sous-tend. II est done tentant dadherer a tout courant de pensec qui relativise sa portee, surtout si on confond la statistique avec les statistiques qui incluent les algorithmes de calcul et les techniques d'analyses de donnees. C'est pourquoi il nous semble indispensable de bien distinguer la phase creairice de la phase calculatoire. La premiere, la modelisaiion: consiste essentiellement a imaginer un mecanisme probabiliste susceptible de produire les donnees ou observations d'interet. Bien sur, dans cette phase, le modelisateur ne s'interdit pas d'avoir aussi recours aux techniques eprouvees d'analyses exploratoires des donnees, afin de mettre en evidence rapidement les traits saillants de l'echantillon. La seconde, l'inference, a pour objet de preciser les parametres du modele probabiliste retenu en remontant des effets (les observations) vers les causes (les parametres). C'est l'inference statistique, dont la mise en oeuvre implique un savoir-faire technique, qui permet l'aide a la decision sous incertitude. En effet, des qu'on dispose de la meilleure connaissance possible des quantites incertaines - il faut pour cela mobiliser toute l'information disponible, ce qui justifie le choix bayesien adopte dans cet ouvrage - on peut donner la distribution de probabilite de toute grandeur interessante pour Ie dccideur. II faut cependant bien reconnaitre que, hormis les maitrises et doctorats en statistique, la plupart des cursus scientifiques se contentent d'une formation assez basique - le niveau elemeniairc n' est pas toujours atteint! - ce qui est paradoxal si on admet son implication dans toutes les sciences experimentales. II faut en rechercher la raison dans le passe. Naguere, le statisticien devait brider sa creativite tout en etant tres bon en math. On peut donc comprendre la reticence des etudiants et des scientifiques non matheux, forces de retenir leur imagination et de s'exprimer dans une langue qu'ils ne maitrisaicnt pas. En particulier, le paradigme bayesien, grand consommateur de calcul integral, ri'etait accessible qu'a une elite, assez peu en prise avec les problemes rencontres par les praticiens des sciences experimentales : les premiers avaient les idees, les seconds les donnees. Les PC rapides ont modifie la donne puisqu'ils ont permis I'emergence des techniques de Monte-Carlo, lesquelles, reduisant fortement les difficultes calculatoires, liberent la creativite du modelisateur. Aujourd'hui, un modele statistique parametrique bayesien est efficacement represente par un assemblage de noeuds relies par des fleches indiquant des relations de cause a effet. Les reseaux bayesiens associent la theorie des graphes, pour la complexite, a la theorie bayesienne, pour la quantification des incertitudes. Les probabilites conditionnelles sont le ciment de ces assemblages. Une fois le roseau construit, l'inference bayesienne precise la distribution de probabilite des parametres (c'est-a-dire les causes) a partir de deux sources d'information : les observations (c'est-a-dire les effets) et l'expertise. Associes aux techniques de Monte-Carlo, les reseaux bayesiens favorisent le dialogue interdisciplinaire et, par la, des modeles innovants

Avant-propos

Xl

et utiles, Notre ambition est que ce livre apporte aux etudiants et aux praticiens synthese et savoir-faire. Pour les fondements plus theoriques, nous renvoyons Ie lecteur a des ouvrages specialises, notamment celui d'Eric Parent et de Jacques Bernier, Le raisonnement bayesien - Modelisaiion et inference et a celui de Christian Robert, Le choix bayesien - Principes et pratique, tous deux publics dans cette meme collection. Pratique du calcul bayesien suit un fil conducteur qui pourrait etre resume par la locution De la plume. . . a la souris. La premiere partie, De la plume, decrit des cas reels relativement simples pour lesquels l'approche bayesienne peut etre monee a la main, sans recours a l'ordinateur. La seconde partie, a la souris, presente des modeles statistiques parametriques plus elabores, impliquant souvent des variables latentes dans une structure hierarchique, Ici, l'inference bayesienne est difficile, voire impossible, sans recours a l'ordinateur. Les reseaux bayesiens et les techniques les plus utiles de Monte-Carlo (avec dependance ou independance) font lc lien entre ces deux parties. Remerciements L'idee de cet ouvrage est nee de l'experience acquise au cours de Statistique Pratique de la collecte et du traitement de l'information environnementale : traitement bayesien de l'incertitude dispense au departement des Sciences et Gestion de 1'Environnement de l'universite de Liege, site d'Arlon (ex-FUL). Sur cette base vinrent s'appuyer les cas reels d'etudes provenant de stagiaires, dingenieurs ou de candidats au doctorat de nos institutions. Un merci tout special a Etienne Prevost (INRA) et a Etienne Rivot (Agrocampus Rennes) qui nous ont permis d'utiliser leurs donnees et travaux pour la realisation des chapitres huit et douze Iondes sur la vie des saumons. Dans Ie me me etat d'esprit, le chapitre dix doit beaucoup aux investigations de M. Philippe Girard, aujourd'hui en poste chez Nestle. Cependant, sans l'appui de nos institutions respectives, I'universite de Liege et l'Ecole nationale du genie rural des eaux et des forets (aujourd'hui AgroParisTech), nous n' aurions pu mener cette tache a bien. Nous tenons a les en remercier. Enfin, Mme Germaine Gazano no us a permis de nous isoler dans son petit paradis Corse, a l'abri des derangements de toute sorte, pour le sprint final ayant construit cet ouvrage dans sa version definitive. Mme Catherine Heyman, secretaire au departement des Sciences et Gestion de l'Environnement de l'universite de Liege, a bien voulu assumer la lourde tache de relire ce livre en no us indiquant les fautes que nous ne voyions plus. M. Jean-Yves Catheland a peint le tableau reproduit en couverture. Nous pensons que l'Art non figuratif illustre bien l'abstraction des concepts mathematiques qui, a l'image des traits et des couleurs, conduisent a une certaine comprehension du monde qui no us entoure. Que toutes ces personnes veuillent bien trouver ici un ternoignage de notre reconnaissance et de notre amitie. Arlon, septembre 2009, Jean-Jacques Boreux, Eric Parent et Jacques Bernier

Sommaire vii

Preface

ix

Avant-propos

xix

Table des illustrations

xxiii

Liste des tableaux

I

De la plume...

1

1 La Statistique : son objet, ses outils 1.1 Le travail du statisticien . 1.2 Deux eccles pour l'inference statistique . 1.2.1 L'ecole classique 1.2.2 L'ecole bayesienne . 1.3 L'analyse statistique bayesienne 1.3.1 La regle de Bayes . . . . . 1.3.2 La distribution predictive a posteriori 1.3.3 Application numerique . . 1.3.4 Retour sur Ie prior 1.4 Le choix bayesien . 1.4.1 Un precede contestable? . 1.4.2 Avantages .

3 3

5 7 9 11

12 12

15 16 16

17 18

2 Decision en avenir incertain : l'avalanche de Montroc 2.1 L'avalanche de Montroc . 2.1.1 Les faits . 2.1.2 Mise en situation . 2.1.3 Un probleme de decision. 2.1.4 Quel(s) modelels) d'echantillounage ? 2.2 Imaginer un mecanisme generateur des observations 2.2.1

Le processus de Bernoulli

2.2.2

Le processus ponctuel de Poisson

.

.

21 21 21 22 22 23 24 24

25

Pratique du calcul bayesien

XIV

2.3

Inference bayesienne . 2.3.1 Le modele beta-binomial 2.3.2 Le modele gamma-Poisson

.

27 27 30

3 Introduction a la modelisation graphique . 3.1 Introduction 3.1.1 Une courte digression . 3.2 Principe de la modelisation graphique . . 3.2.1 L'independance conditionnelle .. 3.2.2 Du reseau bayesien a la loi conjointe 3.2.3 DAG et variables latentes 3.3 Le modele de capture-recapture 3.3.1 Mise en situation . 3.3.2 La modelisation . 3.3.3 Applications . . .

33 33 34

4 Calcul des lois a posteriori 4.1 Introduction . 4.2 Quand la vraisemblance fait Ie posterior. . . . . . . 4.2.1 Approximation asymptotique de la densite

49 49 52

a posteriori . Fondements de ces approximations . . . . . Estimation asymptotique des parametres d'une population gamma . . . . . . . . . . . . . . . . . . . . . 4.2.4 Estimation asymptotique des parametres d'une regression Iineaire . . . . . . . . . . . . . ..... 4.2.5 On retiendra . . . . . .. Methodes de Monte-Carlo par chaines de Markov . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Mise en contexte . 4.3.2 Algorithme (general) de Metropolis-Hastings (MH) 4.3.3 Echantillonnage de Gibbs . . . . . . . . . . . . Methodes de Monte-Carlo. . . . . . . . . . . . . . . . . . . 4.4.1 Simulation par la methode d'acceptation-rejet . . . . 4.4.2 L'echantillonnage et le re-echantillonnage ponderes . 4.4.3 Vers les methodes particulaires . . . . . . . . . . . .

4.2.2 4.2.3

4.3

4.4

5

Le cardinal sort du rang 5 .1 Introduction........... 5.2 Modelisation hierarchique . . . . 5.2.1 Le probleme du tramway 5.2.2 Le probleme des rangs de naissance ..

36 36 38

40 41 41 41 45

53 57 59

61 65

66 66 66 69

72 73 76 81 85 85 87 87 88

Sommaire

6 Les modeles GEV et POT . 6.1 Introduction . 6.2 Le modele GEV 6.2.1 La valeur de projet . 6.2.2 Sensibilite du modele GEV aux hypotheses 6.3 Le modele POT . . . . . . . . . . . . . . . 6.3.1 La distribution de Pareto generalisee . 6.3.2 Le modele POT. . . . . . . . . . 6.4 Du modele POT au modele GEV . . . 6.5 Inference bayesienne sur les parametres d'un modele GEV . 6.5.1 La distribution conjointe a posteriori. . . . . . . . . 6.5.2 Algorithme MH sequentiel applique au modele GEV 6.6 Inference bayesienne sur les parametres d'un modele POT . . . . . . . . . . . . . . . . . . . . 6.6.1 Distribution conjointe a posteriori et inference 6.6.2 Echantillonnage de Gibbs . . . . . . . . . . . 6.7 Trois applications numeriques reelles , . . . . . . . 6.7.1 Le niveau de la mer a Port Pirie (Australie) . 6.7.2 La vitesse du vent a Tunis (Tunisie) 6.7.3 La lame d'eau a Uccle (Belgique) 7 Construire Ie prior 7.1 Introduction........ 7.1.1 Prior non informatif 7.1.2 La conjugaison . 7.1.3 L'analogie . 7.1.4 La methode par introspections successives 7.1.5 L'incertitude n'est pas l'ignorance et la subjectivite n'est pas I'absurdite . 7.2 Definition constructive d'une probabilite subjective. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Caler un prior beta sur deux quantiles elicites du parametrc d'un modele d'observable binomial . 7.3.1 L'expert donne la valeur moyenne de 1r et une incertitude sur celle-ci. . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.2 L'expert donne deux quantiles de 1r . 7.4 Caler un prior conjugue sur deux quantiles elicites des parametres d'un modele d'observable normal . 7.4.1 Dialogue avec l'expert . 7.4.2 Le parametre a elicitor est unidimensionnel 7.4.3 Le parametre a eliciter est bidimensionnel ..

xv

97

98 100 103 104 105 106 108 108 110 110 111 112 113 115 115 116 118 121

127 127 128 130 131 131 132 132 134 134 135 136 136 136 139

XVI

II

Pratique du calcul bayesien

...

a la souris

145

8 Modele de capture-recapture: application au cas des saumons147 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 8.2 Presentation du probleme . . . . . . . . . . . . . . . . . . . . . 148 8.2.1 Les trois dernieres ctapes du cycle de vie du saumon... 148 8.2.2 Variables observees . . . . . . . . . . . . . . . . . . . . . 150 8.2.3 Expertise a priori sur le comportement du saumon . . . 150 8.2.4 Les variables latentes decrivent le phenomena biologique 153 8.3 Inference bayesienne . . . . . . . . . . . . . . 155 8.3.1 Echantillonnage de Gibbs . 156 8.3.2 DAG, nceuds parents, nceuds enfants . 157 8.3.3 Actualisation bayesienne par l'echantillonnage de Gibbs 157 8.4 Resultats numeriques . . . . . . 161 161 . 8.4.1 Annee 1995 163 8.4.2 Cinq annees de donnees 164 . 8.5 Discussion 164 8.5.1 Le role du prior . 165 8.5.2 Le choix du modele. . . 165 8.5.3 Confusion des effets et importance du prior 9 Le modele lineaire generalise 9.1 Introduction . . . . . . . . . . . . . . . 9.2 Retour sur le modele Iineaire classique. . . . . . 9.3 Le modele Iineaire generalise . . . . . . . 9.3.1 Le GLM repond a ces limitations. 9.3.2 D'un point de vue pratique 9.4 La regression logistique . . . . . 9.4.1 La transformation logit . . . 9.4.2 La regression logistique .. 9.4.3 Les prothesistes dentaires seraient-ils particulierement exposes aux pneumoconioses? . . . . . . . . . . . . . . . . 9.4.4 Evaluation de l'action conjointe de deux produits . . . . 9.4.5 Regression logistique avec Ie modele de Finney (1971) .

169

10 Assemblage de modules fonctionnels normaux 10.1 Introduction . . . . . . . . . . . . . . . . . . . 10.2 Construire un modele comme on joue au Lego 10.2.1 Les moyens a mettre en ceuvre . . . . . 10.2.2 Les modeles, leur definition, leurs liens .. 10.3 Regression lineaire (M 1). . . . . . . 10.3.1 Formulation du modele M1 .. 10.3.2 Les conditionnelles completes 10.3.3 Complements sur Ie prior . . . 10.4 Un AR1 pour representor la dependance temporelle (M2)

185

169 170 172 173 175 176 176 177 178 181 182

186 188 189 189 191 192 192 193 193

Sommaire 10.4.1 Formulation du modele M2 . . . . . . 10.4.2 Les conditionnelles completes . . . . . . . . . . . . . 10.5 Modele lineaire a residus autocorreles (M3) 10.5.1 Formulation du modele M3 . . . . . . 10.5.2 Prior des parametres du modele M3 . . . . . 10.5.3 Conditionnelles completes du modele M3 10.5.4 Specification des priors du modele M3 . . . . . . . . .. 10.5.5 Applications 10.6 Modele avec erreur sur variables explicatives (M4) 10.6.1 Formulation du modele M4 . . . . . . . . . . . . . 10.6.2 Specification du parametre ¢ . . . . . . . . 10.6.3 Influence de l'erreur sur la temperature . . . 10.7 Une brique de LEGO supplementaire d'expression multinomiale 10.7.1 Formulation du modele M5 . . . . . . . . . . . . . . 10.7.2 Conditionnelles completes du modele probit (M5) . . 10.7.3 Application du modele multinomial probit (M5)

XVll

194 194 195 195 196 196 197 198 200 200 202 202 202 203 206 207

211 11 Evaluation de la pollution indoor 11.1 Introduction . . . . . . . . . . . 212 11.2 Experimentation et approche classique 212 11.2.1 Modelisation du taux d'emission . . . . . . . . . . 213 11.2.2 Modelisation du changement de masse du polluant 213 11.2.3 Breve etude critique du travail public 214 215 11.2.4 Discussion. . . . . . . . . . . . . . . . . . . . . . 216 11.3 Bruiter Ie modele deterministe . . . . . . . . . . . . . . 11.3.1 Une strategic de modelisation des incertitudes. . 216 11.3.2 Application de la regle de Bayes . . . . . . .. 217 11.3.3 Hesultats . . . . . . . . . . . . . . . .. ..... 218 12 Les avantages de la modelisation hierarchique 12.1 Donnees. . . . . . . . . . . . . . . . . . . . . . . . . . . 12.2 Modele de capture-marquage-recapture . . . . . . . . . 12.2.1 Modele Bernoulli d'alea pour la premiere phase 12.2.2 Modele Bernoulli d'alea pour la seconde phase 12.3 Modele bayesien hierarchique echangeable . . . . . . . . 12.4 Modele bayesien annuel . . . . . . . . . . . . . . . . . . 12.5 Choix des distributions a priori et analyse de sensibilite . . 12.5.1 Priors du modele avec independance annuelle 12.5.2 Priors a deux etagcs du modele hierarchique . 12.6 Resultats . . . . . . . . . . . . . . . . . . . . . . . .

221 222 222 223 224 225 228 229 229 230 231

13 Modeles de changements caches 13.1 Introduction . 13.1.1 Trois exemples hydrometeorologiques .. 13.2 La modelisation des changements .

237 238 239 240

xviii

13.3

13.4 13.5 13.6

13.7

13.8

Pratique du calcul bayesien 13.2.1 Modele M 1 : 1 seule rupture. . . . . . . 13.2.2 Modele M k : k ruptures . . . . . . . . . . . . . . 13.2.3 Modele M a (autoregressif, k ruptures) Representation des distributions a priori 13.3.1 Prior pour les dates 13.3.2 Prior pour les autres parametres Etude du modele M k Methode d'inference . . . . . . . . . Choix de k : ou selection bayesienne de modeles . . . . . . . . . . . . . . 13.6.1 Le facteur de Bayes 13.6.2 Facteur de Bayes et rapport de vraisemblance . 13.6.3 Choix de modele . . . . . . . . . . . . . . . . . 13.6.4 Note sur Ie choix de modele . . . . . . . . . . . . 13.6.5 Avantages et inconvenients des facteurs de Bayes Applications . . . . . . . . . . . . . . . . . . . . . . . . 13.7.1 Application aux modules annuels du Senegal . . . 13.7.2 Application aux apports energetiqucs annuels du SaintLaurent (1943-2000) . . . . . . . . . . . . . . . . . 13.7.3 Application du modele M a au Saint-Laurent . . . . 13.7.4 Debits maximaux annuels de la Dordogne a Cenac Discussion

14 Conclusion

240 241 243 243 244 246 247 249 250 250 250 251 251 252 253 253 254 256 258 260

263

Annexes

265

A Annexe du chapitre 1

267

B Annexe du chapitre 2

273

C Annexe du chapitre 6

279

D Annexe du chapitre 9

287

E Annexe du chapitre 10

293

F Annexe du chapitre 11

305

G Annexe du chapitre 12

307

H Annexe du chapitre 13

313

Bibliographie

325

Index

331

Table des illustrations 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

Taille d'un gar O. Les hyperparametres m , k, a et b sont a determiner de telle facon que le prior conjoint reflete le savoir de l'expert , ici l'OMS . La figure 1.6 represente le modele normal sous la forme d'un reseau bayesien ou DA G (directed acyclic graph) . Nous aurons l'occasion de preciser cette notion dans les prochains chapitres. Les donnees apparaissent dans des carres (ou rectangles) alors que les quantites incertaines (parametres, observables) dans des cercles (ou ellipses). L'empilement de feuilles symbolise l'echantillon : une valeur par feuille. Dans un reseau bayesien, les fieches indiquent des relations causales. Le mecanisme generateur des donnees est done celui-ci : - fournir a > 0, b > 0, k E ]0, 1[ et mE lR (expertise) ; - tirer T dans Ga (a, b) puis tirer f-l dans N (m , kT) ; - pour i allant de 1 a n , tirer Zi dans N (f-l , T) .

~ k

[TI/

8

Figure 1.6 - DAG du model e normal.

D'un point de vue analytique, les calculs de la loi a posteriori sont relativement simples (voir annexe A, p. 267 et suiv.). En effet, a partir du modele (1.10), la vraisemblance d' un n-echantillon iid normal est triviale et l'app lication de la regle de Bayes conduit aux resultats suivants OU Z et s2 representent respectivement la moyenne et la variance empirique des observations. RappeIons l'essentiel des result at s donnes en ann exe A : - le posterior marginal de Test une loi gamma, de par ametre de forme a' =

1. La Statistique : son objet, ses out ils

~ (n + 2a) et de parametre d' inverse echelle b' = ~

[n s 2 + 2b +

15

n'tk

(z - m)2 ] ; le post erior marginal de J.L est une loi de Student a u' = 2a' degres de liberte, localisee sur m' = (nz + km) / (n + k) et de parametre d'echelle (J' = J b'/((n +k)a') ; - la distribution predictive a posteriori de l'observable est une loi de Student a v' degres de liberte, localisee sur m' et de para metre d'eche lle (J" (J\ /n + k + 1.

~

1.3.3

Application numerique

Revenant au probleme du radon (exemples 6 et 7), les choix m = 3.7, k = 0.5, a = 2, b = 1 donne nt un prior conjoint compatible avec les donnees de l'OMS (une concentration moyenne de l'ordre de 40 Bq/rnd , tres variable selon la nature du sol). Ce prior et toutes les donnees du tableau 1.1 conduisent aux resultats presentee a la figure 1.7 (l'axe des abscisses est en coordonnees logarithrn iques). En particulier, on voit que la probabilite de depasser la valeur guide de 400 Bq /rnd d'air est egale a 0.24. Un prior non informatif porte cette probabilite a 0.25. La difference est minime et ne modifie pas Ie risque de cancer . Cependant, si on refait le merne exercice avec seulement les deux donnees du rez-de-chaussee, on trouve respectivement 0.25 et 0.17 ! Dans ce cas, ignorer l'expertise revient a sous-estimer largement le risque.

2 -

-

Predictive : [z I y]

1\

c::::=J Pr(z > Zo I y) = 0.24

1.8

J\ J \

- - - Posterior : [Il l yJ

1.6

) l

J I

I I

1.4

}

I I f I I I I

1.2

:;

a0.8

J J J

0.6 0.4 0.2 0

.J 0

2

3

! I ! f

45

I

\ I I l I I I

1 I l

) \ \ 6

7

Log de la conc entra tion en rado n (Bq/m 3)

Figure 1.7 - Conc entration du radon en Minn esota.

8

16

Pratique du calcul bayesien

1.3.4

Retour sur Ie prior

Nous avons vu que la roue de la fortune permet d'illustrer le travail d'elicitation d'un prior. Sur de nombreux modeles et exemples, le chapitre 7 developpera les outils d'encodage du savoir de l'expert sous forme d'une distribution de probabilite. Les deux proprietes statistiques suivantes sont regulierement invoquees pour en faciliter la mise en ceuvre.

La conjugaison. L'analyste regarde la forme de la fonction de vraisemblance et choisit une famille de lois qui se « marie bien» avec elle. Par exemple, la structure de la vraisemblance d'un n echantillon iid selon une distribution exponentielle de parametre d'echelle p > 0 est en p" exp (-nyp). Le prior conjugue est une loi gamma dont la forme fonctionnelle s'ecrit pa-l exp (-bp). La distribution a posteriori de p suit immediaternent : pia, b,n, y rv gamma (a + n, b + ny). Dans le probleme du radon, nous avons utilise les proprietes de la conjugaison pour construire le prior conjoint. Application sequentielle de la regie de Bayes. Un prior peu ou non informat if sur un jeu de donnees fournit une distribution a posteriori qui peut servir de prior pour un autre jeu de donnees. Par exemple, on peut appliquer le modele developpe ci-dessus aux concentrations en radon relevces dans un comte voisin de celui de Goodhue (Minnesota, Etats-Unis). Le posterior obtenu sur ce jeu de donnees est un prior credible pour l'analyse des donnees du tableau 1.1.

Quel que soit le moyen utilise pour construire le prior, il doit etre interprete comme une succession de paris sur les valeurs du parametre, bien sur sans mobiliser les donnees impliquees dans la vraisemblance.

1.4

Le choix bayesian

La figure 1.8 synthetise Ie paradigme bayesien. Deux modeles doivent etre specifies. Le modele d'echantillonnage et le prior. C'est pourquoi les statisticiens bayesicns designent leurs modeles par des expressions du type priorvraisemblance (du moins quand le prior peut etre decrit par une distribution standard). Ainsi, on parlera des modeles beta-binomial, gamma-Poisson, normal-gamma-normal, etc. Qu'on soit classique ou bauesien, le choix du modele dechant.illonnage est decisif. II n'y a pas de recette, mais l'experience de l'analyste compte. La representation de la connaissance a priori est tout aussi delicate. L'expert passe en revue toutes les valeurs possibles du parametre et parie sur chacune d'entre elles. Ensuite, cette connaissance a priori est mise a jour par les donnees via la regle de Bayes.

1. La Statistique : son obj et , ses outils

17

Model. SIal d 'Occ urren ces

[yle] Connaissance a priori

Connaissan ce mi se it jo ur

(Exp erti •• ) Formule d. Bay••

[e]

~ L_"'-

-"-_.

e

[B] [Y IB] [BIY ] = J[B][Y IB]dB e

....

I'~I

~

A

~

Mei lleure prec ision sur les phenome nes mconnus

donnees Experime nt ales , Y~

{Y I' Y, • . . . Y k }

Figure 1.8 ~ Le paradigme bayesien : resume.

1.4.1

Un pro cede contestable?

L'analyse bayesienn e repose sur les donnees - c'e st la composant e dite objective - et sur les idees du chercheur - c'est la composante dite subject ive. Dans l' excellent livre d 'Alfred Renyi Calcul des probabilites (edition originale en langue allemande, © 1962, VEB Deutscher Verlag der Wiss enschaften, Berlin ; reimpression aut orisee de la traduction fran caise, © 1966, Dunod, Paris) on t rouve le comment aire suivant (p. 77) . [.. . ]Le theorem e de Bayes est parjaitem ent demontre, personne ne m et en dout e sa justesse; c 'est seulem ent de ses application s pratiques qu 'on dispute (sic) . [... ] Si on connait les probabilites dit es a priori , on peut appliquer le theorem e de B ayes et calculer les probabilites a posteriori . Cependant, les probabilii es a priori sont souvent inconnues et on leur att ribue qeneralem eni des valeurs arbitrai res ; c 'est ce precede qui est veritobleme nt contestable. Qu ar ante ans se sont ecoules depuis ce comment aire. Non ce n'est pas contest abl e d 'attribuer des probab ilite s a priori a des evenement s. Evaluer des chances sur la base de son exp erience est une activite int ellect uelle recurrente partagee par la majorite des et res pensan ts. Croire qu e seules les donnees garant issent l'objectivite du verdict est une erreur, car les donnees resul tent de choix, souvent impli cit es! Ainsi l'echelle europee nne de risque d 'aval an che comporte cinq indic es classes sa ns ambiguite suivant l'importan ce

18

Pratique du calcul bayesien

du risque auquel s'expose l'usager. Chaque niveau de risque est defini par une evaluation de la stabilite du manteau neigeux fondee sur une seric de criteres et des consequences a assumer en cas d'avalanches. L'expert peut attribuer a priori une probabilite Pk a l'indice k. Ce n'est pas plus arbitraire que de combiner des informations pour construire une telle echelle et la faire accepter par les pays concemcs ; ce n'est pas plus arbitraire que de selectionner quelques indicateurs parmi les dizaines qui auraient pu etre choisis. L' activite scientifique ne nie pas la subjectivite, mais elle vise son controls. Par consequent, tous les resultats generes par une demarche scientifique sont toujours conditionnels aux differents choix qui ont ete faits, qu'ils soient d'ailleurs peu ou prou justifies. Les statisticiens bayesiens se distinguent par leur volonte de les decrire clairement.

1.4.2

Avantages

Nous avons vu que le statisticien classique iuterprete son intervalle de confiance en se referant a une collection d'echant.illons qu'il aurait pu observer s'il avait reproduit son experience dans les memes conditions. Le statisticien bayesian ne rencontre pas cette difficulte. L'intervalle de confiance, qu'il appelle intervalle de credibilite pour le distinguer de son homologue classique, a une interpretation naturelle qui porte directement sur la valeur inconnue qu'il cherche a cerner. Pour un risque Q fixe, les limites de l'intervalle de credibilite sont les percentiles ()a/2 et ()1-a/2 du posterior tels que Pr (()a/2 ~ () ~ ()1-a/2) == 1- Q. Ce n'est pas seulement une question philosophique. Le concept de repetition d'cxperiences dans les memes conditions peut n'avoir aucun sens. La probabilite qu'un meteorite detruise la Terre dans les mille pro chaines annees ne peut etre fondee sur la notion de repetition. Outre les difficultes d'interpretation. le paradigme classique n'offre pas l'equivalent de la distribution predictive a posteriori. Or c'est bien le futur sachant le passe qui interesse le decideur. Par exemple, le conseil municipal de Chamonix (voir chap. 2) aurait pu se poser la question suivante : Quelle est la probabilite que le site de Montroc subisse au rnoins une avalanche dans les vingt prochaines omnees sachant qu'on y en a observe six depuis 1843? Pour repondre a cette question, il faut faire des hypotheses, postuler un modele dechantillonnage, questionner les experts, bref construire un modele statistique parametrique 7 . Cependant, le dccideur, par exemple le conseil municipal de Chamonix, n'a rien a faire des parametres du modele! Merrie s'il en ignore le nom, c'est la distribution predictive a posteriori qui I'interesse, distribution obtenue en integrant, par rapport au pararnetre, Ie produit de la probabilite de l'observable par la distribution a posteriori du parametre. Le statisticien classique ne peut pas realiser cette operation puisque, pour lui, Ie parametre ne varie pas! Enfin, les petits echantillons sont par definition peu informatifs et le theoreme central limite ne tient plus! Le statisticien classique est peu arme pour traiter ces cas difficiles. Le statisticien bayesien, lui, peut palier un manque de 7

la statistique bayesienne non parametrique n'est pas l'objet de ce livre.

1. La Statistique : son objet, ses outils

19

donnees en introduisant de l'expertise dans le modele. Ces situations ne sont pas rares en sciences et particulierement en sciences de l'environnement. Ainsi, on a vu que pour mesurer la concentration en radon dans une maison, il faut laisser le detccteur sur place (par exemple dans la piece la plus frequentee) durant une periode de 2 a 3 mois. Trois mois, une mesure ! Pourquoi se priver d'une seconde source d'information qu'est l'avis de l'expert (p. ex. l'OMS) ? La controverse philosophique ecole classique versus ecole bayesienne est finalement peu interessante. II faut faire un choix et le notre est clair: c'est le paradigme bayesien, La suite de ce livre est une collection de modeles utiles, car de portee assez generale. Chacun d'entre eux constitue un chapitre. Les difficultes calculatoires sont mises en evidence et une solution est proposee. II est possible que d'autres solutions, plus elegantes, existent. Tous ces modeles sont illustres avec des exemples concrets (donnees reelles). Nous postulons que le lecteur a une culture generale en mathematique du niveau du baccalaureat es sciences. La connaissance des distributions standards est indispensable. Elles sont reprises dans l'appendice A de (Gelman et al., 2004).

Epilogue Ce premier chapitre a introduit l'idee que construire un modele statistique parametrique revient a imaginer un mecanisme probabiliste susceptible de reproduire les observations. L'observable est une variable aleatoire pour laquelle on postule une distribution de probabilite souvent nommee modele de connaissance. II s'agit en fait d'une famille de lois de probabilite indexee par un parametre inconnu de dimension finie. On le notera souvent B. Le choix d'un modele de connaissanee est done une affaire dexperience matinee d'imagination et d'audace. Sous le paradigme bayesien, () est incertain mais prend ses valeurs dans un espace de dimension finie, 8, appele ensemble des ciats de la nature. Avant de disposer de l'echantillon de donnees, un specialiste du probleme etudie pourra souvent dire quelque chose sur (). II pariera plus volontiers sur telle plage de valeurs que sur telle autre. Ainsi, l'incertitude sur () peut etre decrite par une distribution de probobilite a priori ou prior. La regle de Bayes reactualise cette expertise en multipliant le prior par la vraisemblance de l'echantillon. Apres normalisation, le resultat obtenu est la distribution a posteriori de B (ou posterior). Toute utilisation ulterieure, notamment l' aide a la decision, sera fondee sur la distribution a posteriori de B. Le fil conducteur de cet ouvrage pourrait se resumer par l'aphorisme de la plume ala souris. En effet, naguere le modelisateur ne disposait que de ses idees, d'un porte-plume et d'une feuille blanche. Aujourd'hui, l'ordinateur personnel a remplace le porte-plume et demultiplie les capacites de traitement. Toutefois qu'on ne s'y trompe pas! L'imagination et la creativite constituent toujours les pierres angulaires du raisonnement conditionnel bayesien. Sans modele, le stockage des donnees dans un ordinateur, meme performant, est improductif!

20

Pratique du calcul bayesien

A contrario, l'art de la construction de modeles probabilistes ressemble a l'apprentissage de la musique : il faut commencer par le solfege. La premiere partie de cet ouvrage propose l'etude des gammes, la seconde nous entraine vers des partitions plus evoluees. Le chapitre 2 presente un probleme decisionnel complet, fonde sur un fait reel et tragique : I'avalanche de Montroc. Les modeles sous-jacents - Ie modele beta-binomial et le modele gamma-Poisson - sont calculables « a la plume »,

Chapitre 2

Decision en avenir incertain I'avalanche de Montroc

• •

Prologue Quand on s'interesse a une experience aleatoire dichotomique, l'hypothese que les observations successives constituent un processus de Bernoulli peut etre justifiee, tantot par la nature de l'experience aleatoire (p. ex. jeu de pile ou face), tantot constituer une hypothese de pure commodite «pour voir». Dans tous les cas, elle conduit au modele beta-binomial et, quand l'evenement d'interet est rare, au modele gamma-Poisson. Ces modeles, tres simples, nous permettent de construire un probleme fictif d'aide a la decision, fonde sur un drame reel ayant fait la une des journaux : A urions-nous pu eoiier La catastrophe de Montroc?

2.1 2.1.1

L'avalanche de Montroc Les faits

Le 9 fevrier 1999, une avalanche meurtriere (douze deces) a detruit une partie du hameau de Montroc pres de Chamonix. Cette coulee de neige a englouti vingt-trois chalets dans une zone declaree « constructible», car consideree comme hors d'atteinte d'apres la cartographie des risques etablie en 1992. En fait, avant la date fatidique, la derniere avalanche sur ce site avait ete observee en 1945. Toutefois, selon Le Dauphine libere, cinq avalanches survenues entre 1843 et 1945 n'auraient pas ete prises en compte", 1

Le maire de Chamonix a ete condamne

a 3 mois

de prison avec sursis le 17 juillet 2003.

22

Pratique du calcul bayesien

2.1.2

Mise en situation

Nous sommes en 1992 et le conseil municipal de Chamonix attend votre etude pour prendre sa decision, c'est-a-dire accepter ou refuser de declarer la zone constructible. Vous savez que la derniere avalanche a ete observee en 1945 et cinq autres avalanches ont affecte Ie site entre 1843 et 1945 (vous ignorez les annees}, II est clair que vous ne pouvez pas dire et faire n'importe quoi. Votre horoscope ou une analyse statistique naive ne constituent pas une methode devaluation des risques conforme aux regles de la demarche scientifique. En particulier, les chances qu'a un evenement de se produire demain ne sont pas insignifiantes simplement parce qu'il ne s'est pas produit depuis quarante-sept ans! L'avenir est incertain et le calcul des probabilites est l'instrument de mesure de toute incertitude. Avant de passer a l'action, il peut etre utile de se rememorer quelques pensees. - «La verite, ce n'est pas le certain et l'incertain, ce n'est pas l'ignorance », (Ilya Prigogine (1917-2003), prix Nobel de chimie (1977))

- «II est bon de suivre sa pente, pourvu que ce soit en montant », (Andre Gide (1869-1951), prix Nobel de Litteraturc (1947), Les faux-monnayeurs (1925))

- «All models are false, some are useful». (Bernardo & Smith, (Bernardo et Smith, 1994))

2.1.3

Un probleme de decision

En 1992, ce qui interesse le decideur, ici Ie conseil municipal de Chamonix, c'est le risque associe aux deux decisions qu'il peut prendre: d1

== declarer la zone constructible et perdre C 1 M

d2

== refuser lc projet et perdre C2 M

EUROS si le site subit au moins une avalanche grave dans les h pro chaines annees (indemnisation des victimes) ; EUROS si le site ne subit aucune avalanche dans les h prochaines annees (les non- recettes).

Le decideur doit donc fixer un horizon de prevision h et votre travail est d'evalucr la probabilite p (h) d'observer au moins une avalanche destructrice sur cette periode, Le tableau 2.1 resume ce probleme de decision en termes de perte associec a chaque decision selon que I'evenement redoute se realise dans les h annees, avec la probabilite p (h), ou ne se realise pas avec la probabilite complementaire. Selon la theorie de la decision, voir p. ex. (Bernier et al., 2000), une reqle de decision coherente consiste a opter pour la decision qui minimise la valeur attendue de la perte totale Ct E (Ctld 1 , h) E (C t ld2 , h)

p(h) X C 1 (1 - p (h))

(2.1) X

C2

(2.2)

2. Decision en avenir incertain : l'avalanche de Montroc Couts

d1 d2

23

Etat de la nature ()

p (h) 01 0

1-p(h) 0

C1

Tableau 2.1 - Montroc : pertes associees aux decisions. Par consequent le rapport (2.3) fournit une regle de decision rationnelle (2.4)

Remarque 2.1 II n'est pas necessaire d'estimer ces couts avec une grande precision. D'une part, le bon sens permet de soutenir que la destruction d'un site habite coute plus cher que les non-recettes : C1 > C 2 . D'autre part, il est recornmande de batir divers scenarii C 1/C2 et de considerer divers horizons de prevision

h.

Pour chaque couple

(h, g~)

Ie calcul du rapport r indique la

decision qui est rationnelle (voir fig. 2.3). Imaginons que le decideur fixe h == 30 ans et estime que 0 1/02

~

10. Si,

a l'issue d'un raisonnement coherent, vous trouviez p (30) ~ 0.08 alors vous devriez recommander la decision d 1 , car r ~ 0.87 < 1 (eq. 2.4). Et si l'an-

nee suivante une coulee de neige rasait le site, auriez-vous pour autant mal travaille ? La reponse est categorique : non, car la probabilite est un concept previsionnel, ante evenement. Si l'evenement rare se realise, vous n'avez tout simplement pas eu de chance et il faut l'accepter. De telles situations se presentent dans la vie de tous les jours. Par exemple, la perte des quatre moteurs d'un avion est un evcnement qui a une probabilite tres faible, mais cet evenement s'est produit et des gens sont morts. Bien entendu, le taux d'echec a l'issue de demarches folkloriques est incomparablement plus eleve.

Remarque 2.2 Bien sur, il est possible de discuter la valeur du rapport C 1/C2 , car r augmente avec lui. Ainsi, sous les memes hypotheses, des que le rapport des couts vaut 12 il faut recommander d 2 . On peut d'ailleurs faire une analyse de sensibilite sur ce rapport.

2.1.4

Quel(s) modelers] d'echantillonnage?

Convenons qu'une annec quelconque est « noire» (code 1) si on y observe au mains une avalanche importante sur le site d'interet. Elle est « blanche»

24

Pratique du calcul bayesien

(code 0) dans le cas contraire. A Montroc, on a releve six annees « noires» sur la periods 1843-1992. Le choix d'un modele d'echam.illonnagc (on dit aussi modele de population) fait partie des hypotheses de la modelisation, Entrent dans les raisons de ce choix des considerations de cornmodite mathematique, de realisme et de parcimonie des parametres. Tous les resultats obtenus sont necessairement conditionnels a l'adoption de ce modele. Ce chapitre se limite aux modeles de connaissance suivants : Ie modele binomial et le modele de Poisson.

Remarque 2.3 D'autres modeles d'echantillonnage sont possibles. En effet, il ne faut pas confondre l'absence d'information avec l'absence reelle d'avalanche, car on peut tres bien imaginer que des coulees de neige n'aient pas ete enregistrees. La modelisation de ce modele de donnees manquantes sort du cadre de cet expose.

2.2 2.2.1

Imaginer un mecanisme generateur des observations Le processus de Bernoulli

A chaque annee t, on associe une variable aleatoire de Bernoulli, disons Yt, qui prend la valeur 1 avec la probabilite 7rt si le site de Montroc subit au moins une avalanche grave et la valeur 0 avec la probabilite complementaire dans Ie cas contraire. Si on postule que ces variables aleatoires sont ituiependanies et identiquement disiribuees (Vt : 7rt == 7r), la suite {Yt} constitue un processus de Bernoulli:

Remarque 2.4 Conceder que le modele d'echantillonnage est un processus de Bernoulli est d'abord un choix de commodite, En effet, on sait que le climat a change depuis Ie milieu du XIX e sieclc (c'est-a-dire la composition de l'urne a change) et il est meme possible qu'il y ait de la memoire dans Ie systeme. Cependant, en acceptant ce modele, au moins pour un temps, on va pouvoir quantifier Ie risque associe a chaque decision. Ensuite, il faudra discuter les resultats a I' aulne des hypotheses qui y ont conduit. Le modele binomial Puisque chaque annee est representee par une variable aleatoirc de Bernoulli, leur somme n

x ==

LYt t=l

2. Decision en avenir incertain : l'avalanche de Montroc

25

est une variable aleatoire binomiale, de parametres n, 1T, dont la densite s'ecrit : (2.5)

ou

n! (~)==-(n - x)!x!

2.2.2

Le processus ponctuel de Poisson

Le processus ponctuel de Poisson (voir annexe C) est un modele, un processus sans memoire, qui interdit les simultaneites et qui considere que les occurrences apparaissant dans des intervalles de temps disjoints sont independantes.

Comrnencons par preciser ce qu'est un evenemcnt ponctuel Sur la periode d'interet de longueur finie l, on divise l'axe du temps en n periodes elementaires de duree constante ~l : l == n.Sl, Des lors, n ----t 00 comme ~l ----t O. Mais l fini et ~l ----t 0 signifient que l'evenement d'interet est un point sur l'axe du temps, c'est-a-dire un evenement ponctuel.

Exemple 2.1 On observe un carrefour pendant 5 ans. Un jour est «rouge» si on y constate au moins un accident avec lesions corporelles. ~l == 1/365 ~ n ~ 1.8 x 103 jours. • La distribution de Poisson est un cas limite de la distribution binomiale Numerotons les periodcs elementaires dans l'ordre de succession depuis 1 periode elementaire, on peut associer une variable aleatoire de Bernoulli qui prend la valeur 1 avec la probabilite invariante 1T si l'evenement d'interet se realise. Si x periodes elementaires parmi les n ont vu l'evenement d'interet se realiser, on a un processus de Bernoulli:

a l. A chaque

[xln,1T]

,

_no), ,7fx (1 _ 7f )n-x (n x .x. n(n-1)···(n-x+1) X( )n-x - - - - - - - - - - 1 T 1 - 1T x!

(1Tn)X x!

(1 _ ~) ... (1 _ ~) (1 n

n

1T)n

(1 - 7f)X

L'evenement d'interct est un evenement rare si x « n, c'est-a-dire En posant A == 7fn > 0 OU n ----t 00 et 1T ----t 0

on obtient la distribution de probabilite de Poisson

1T ----t

O.

26

Pratique du calcul bayesien

AX

[XIA] == ,exp (-A)

(2.6)

x.

En effet : lim

n-+oo

(1 - ~) ... (1 - ~) == 1 n

n

lim (1 - 7f) X == 1

7r-+O

lim

n-+oo

(1- ~)n n

=exp(-A)

Remarque 2.5 La distribution de probabilite de Poisson (eq. 2.6) est definie sur l'ensemble des entiers naturels N : 1 == exp (- A)

AX

AX

L , {:} x=o L ,x. x=o x. 00

00

== exp (A)

Dans Ie chapitre 5, nous utiliserons une variable aleatoire de Poisson prenant ses valeurs dans No == N\ {O}. On dit que la distribution de Poisson est tronquee sur No. Des lors 00 AX == exp (A) - 1

L ,x.

x=l

et la distribution de probabilite de Poisson de parametre A > 0 tronquee sur

No s'ecrit :

[xIA] ==

AX

1

-;y -ex-p-(A-)---1

Le processus ponctuel de Poisson Soit X j le temps qui separe deux occurrences successives de l'evenement d'intcret. Si on postule que les durees X j sont iid selon une loi exponentielle de parametre d'echelle 1/A, alors la distribution du nombre d'occurrences, disons Y, sur une periode de l unites est donnee par la loi de Poisson de parametre Al. La reciproque est vraie : si la distribution du nombre d'occurrences sur une periode de longueur lest donnee par la loi de Poisson de parametre Al, alors les durees sont iid selon une loi exponentielle de parametre d'echelle 1/ A.

x, IA ~d dgamma (xiI, A-I) {:} YIA, l

r-;»

dpois (yIAl)

(2.7)

~~

Le parametre A est la cadence des occurrences, c'est-a-dire leur nombre sur la periode de reference.

2. Decision en avenir incertain : l'avalanche de Montroc

2.3

27

Inference bayesienne

Ayant imagine un processus susceptible de generer les observations, il faut maintenant estimer son parametre caracteristique - qui peut avoir plusieurs composantes auquel cas c'est un vecteur - et quantifier l'incertitude afferente, Comme on l'a vu au chapitre 1, le paradigme bayesian offre un cadre de raisonnement particulierement fiable et fecond. La vraisemblance est conditionnelle au parametre et la distribution a priori du parametre, ou prior, decrit l'incertitude de l'expert sur celui-ci. La regle de Bayes dit comment reactualiser cette expertise disposant des donnees: il suffit de multiplier la vraisemblance par le prior. La distribution a posteriori du parametre, ou posterior, implique la normalisation de ce produit. Cette operation peut se reveler compliquee, voire impossible, sans le concours de methodes speciales. Ce ne sera pas le cas ci-dessous. A condition de connaitrc les fonctions eulerietuies de premiere et seconde espece, respectivemment la fonction gamma (symbole f) et la fonction beta (symbole B), tous les calculs peuvent etre faits a la plume. Ces fonctions ne doivent pas etre confondues avec les distributions de probabilite gamma et beta auxquelles elles ont d'ailleurs donne leur nom (voir annexes B et B).

2.3.1

Le modele beta-binomial

La vraisemblance Rappelons que la vraisemblance mesure les chances d'observer I'echantillon conditionnellement au parametrc. Pour l'avalanche de Montroc, le modele d'observation est la loi binomiale (eq, 2.5). La vraisemblance est donc immediate

(2.8) Choix du prior et application de la regie de Bayes Quand on regarde la vraisemblance (eq. 2.8), on reconnait immediatement la signature fonctionncllc/ d'une densite de probabilite beta. On dit qu'un prior beta est conjugue a une vraisemblance binomiale. La conjugaison a deja ete abordee au chapitre 1 (p. 16) et sera davantage explicitee au chapitre 7. Bien sur, il faut preciser les parametres du prior beta, disons a > 0 et b > 0 :

[Bla, b]

ex ga-l

(1 _

(})b-l

(2.9)

Les parametres des lois a priori decoulent de l'expertise reconnue et sont ponctuels, c'est-a-dire sans incertitude. La litterature scientifique les designe souvent sous Ie nom d 'tujperparametres. 2 L'expression «signature fonctionnelle » traduit I'idee que la relation rnathematique en main constitue la partie essentielle d'une densite de probabilite. II ne reste plus qu'a la normaliser.

28

Pratique du calcul bayesien

Comment determiner les hyperpararnetres d 'un prior beta? RappeIons que le parametre 8 (eq. 2.8) represente la probabilite qu'une annee calendaire, choisie au hasard, voit au moins une avalanche debouler sur le site de Montroc . Ces annees « noires» sont plut6t rares, sinon Ie probleme de decision n'aurait aucun sens. Pour l'exemple , imaginons qu'un specialiste des avalanches accorde une chance sur dix a 8 de depasser la valeur 0.05 et cinq chances sur cent, d'etre inferieure a la valeur 0.01. Ces paris lui sont propres et temoignent de son savoir. Pour l'analyste, l'expert a fourni les quantiles 8go ~ 0.05 et 85 ~ 0.01. A partir de ceux-ci, une methode numerique lui permet de determiner les hyperparametres a et b: a ~ 3.82 et b ~ 124.1 (voir chapitre 7, p. 135). Ces valeurs ne varient que si l'expert change d'avis, ce qui est son droit. Tant qu'il ne Ie fait pas, elles sont connues sans incertitude. La regle de Bayes reactualise cette expertise en tenant compte des donnees : x = 6 pour n = 150. La distribution a posteriori de 8 est encore une densite beta (interet de la conjugaison) , dont les parametres integrent l'expertise et les observations, c'est-a-dire toute l'information disponible : 81k, n, a, b rv dbeta (81x + a, n - x

+ b)

(2.10)

La figure 2.1 montre le prior et le posterior ainsi obtenus.

40,---

-,----

---,---

----,-

0.01

0.02

0.Q3

----,-

-

,------

-,----

---,---

----,-

0.06

0.07

0.08

----,-

---.,

35

30

~

:B

25

~

~ 20 2

.~

c

~ 15

10

0.04

0.05

e

0.09

0.1

Figure 2.1 - Avalanche de Montroc : Ie mod ele beta-binomial.

2. Decision en avenir incertain : l'avalanche de Montroc

29

La distribution predictive a posteriori Le probleme decisionnel, autoriser ou refuser de declarer la zone constructible, doit etre pose dans une perspective predictive. Cette decision, dont les consequences concernent les h annees futures, est fondee sur les informations du passe, ce qui justifie les calculs presentes ci-apres. On pourrait alors s'interroger sur l'apparente contradiction entre l'hypothese iid et la pretention de prevoir l'avenir en se servant du passe. En fait, le lien entre le futur et le passe s'appuie sur la connaissance du parametre (). L'idee est la suivante. On s'interroge sur la probabilite d'observer Y annees « noires» dans les h prochaines annees sachant que, dans le passe, on en a reellement observe x en n annees. Y est une variable aleatoire discrete prenant ses valeurs dans l'ensemble n (Y) == {O, 1, ... ,h}. La distribution predictive a posteriori donne les chances de chacune des occurrences yEn (Y) en impliquant les acquis (a,b, n, x) et l'horizon de prevision envisage (h). On l'obtient en integrant la distribution jointe de Y et () sur toutes les valeurs possibles de () (voir chapitre 1, eq, 1.9) :

[Y = ylh, a, b, n, x] = ([yIO, h] [Olx, ti, a, b] dO

Je

(2.11)

- [yl(), h] est la probabilite de y donnee par la loi binomiale (eq, 2.5), - [()Ix, n] est lc posterior beta obtenu ci-dessus (eq, 2.10) L'integration (eq, 2.11) ne pose aucun probleme, La distribution predictive a posteriori du modele beta-binomial est la distribution de Polya :

[Y == y Ih ,a, b,n, x ]==(h)B(y+x+a,h-y+n-x+b) B (x + a,n - x + b) y

(2.12)

Remarque 2.6 II est important de noter que la calcul de la distribution predictive a posteriori est realise en integrant un produit de distributions de probabilite, En d'autres mots, il faut tenir compte des constantes de normalisation.

La figure 2.2 montre la distribution de Polya pour quatre horizons de prevision. La probabilite d'observer au moins une annee « noire» a l'horizon h est le complement de n'en observer aucune

p (h)

[Y 2 1Ih,a,b,n,x] 1 - [Y == Olh, a, b, n, x] 1- B(x+a,h+n-x+b) B(x+a,n-x+b)

(2.13)

Maintenant nous sommes en mesure d'appliquer notre regle de decision (eqs, 2.2 et 2.4) avec differents scenarii 0 1/02 (fig. 2.3). On voit que refuser de rendre le site de Montroc constructible (decision d2 ) est une decision rationnelle des que l'on envisage un horizon de prevision compatible avec un projet de lotissement (p. ex. 20 ans ou plus).

30

Pratique du calcul bayesien

h = 5 ans

08

0.6

~ 0.6 :c m

~

0.4

0.4 0.2

0.2

4 1 2 3 Nombre d'arneesnoires :y

0

h = 20 ans

0.5 0.4

9 10

h =30ans

0.3

~

'" z 0.3

~

m

02

02

o,

0.1

0.1 0

1 2 3 4 5 6 7 8 Nombre d'arneesnoires : y

0.4

.QJ

~

h = 10 ans

0.8

0

5 10 15 Nombre d'aonees noires : y

20

5 10 15 20 25 Nombre d'annees noires: y

30

Figur e 2.2 - Avalanche de Montr oc : dist ribut ion de Polya pour quatre horizons de prevision.

D iscussion

L'hypoth ese ii d, fond atrice du raisonnement , est discutable. Comment en effet soute nir que le processus est sans memoire et que e est invariant sur la periode 1843 --> 1992 + h ? Le modele beta-binomial est done critiquable. Cependant, da ns l'etat act uel des connaissances, l'hypoth ese ii d n'est pas plus discutable que son contraire et c'est peut-etre la seule qui soit compatible avec la « pauvrete » de l'information conte nue dans l'enonce du problerne. Dans Ie cas OU cet te hypothese serait rejetee, il faudr ait alors developp er un modele beaucoup plus sophist ique. Neanmoins, nous sommes persuades que si le conseil municipal de Chamonix avait pu beneficier de l'inform ation generee par ce modele, il aurait refuser de prend re le risque d'un drame humain .

2.3. 2

Le modele gamma-P oisson

La distribution a posteriori

Si on considere que la periode 1843-1992 est la periode unite (150 ans), alors l = 1 (eq. 2.7) et la vraisemblance s'ecrit (eq. 2.6) :

[X = XIA] ex AXexp (- A)

(2.14)

2. Decision en avenir incertain : l'avalanche de Mont roc

20.----

-

- .-

-

-

-

,...---

-

-

- .-

-

-

-

,...---

-

-

31

-,

18 16

14

10

15

20

25

Figure 2.3 - Avalanche de Montroc : regie de decision. On voit que la forme fonctionnelle d'une distri bution de P oisson est la merne que celle d'une densite gamma. Ceci suggere de decrire l'in certit ude a priori sur le par ametre de Poisson, ici A, a l'aide d'une densite de probabilite gamma dont il fau t fixer Ie par metre de forme, a > 0, et Ie par am etre d 'echelle, b » 0 :

[Ala,b] ex Aa - 1 exp (- bA)

(2.15)

Remarque 2.7 Dans cette formulat ion de la distribution gamma, E (Al a, b) = al b et V (Ala, b) = alb 2 . La regie de Bayes fournit la distribution a posteriori du par ametre de Poisson : (2.16) AIH ", dgamma (Alx + a, 1 + b) ou la let tre H repr esent e to utes les hypo th eses, notamment les hyperparametres a et b, et les donn ees, ici le nombre d'ann ees « noires » , x .

La distribution predictive a posteriori Soit Y, la variable aleatoire « nombre d'annees noires a Mont roc » dans les h prochaines annees. On sait qu'on en a observe x = 6 sur une periode I de 150 ans. Dans l'ann exe B, on montre que la distribu tion predict ive a posteriori est

32

Pratique du calcul bayesien

une loi binomiale negative dont la distribution de probabilite s'ecrit

(2.17)

ou 7r

== h~1~b'

r == x

+a

(2.18)

Dans le cas de Montroc, meme avec des priors non informatifs, les modeles beta-binomial et gamma-Poisson produisent une aide a la decision vraiment similaire a celIe montree a la figure 2.3.

Epilogue Dans un contexte decisionnel, lorsque les enjeux sont importants, la quantification du risque attachee a chacune des decisions en competition est une etape obligatoire. Dans cette perspective, la modelisation statistique bayesienne mobilise les donnees disponibles et l'expertise reconnue pour fournir une information utile au decideur. La credibilite du paradigme bayesien reside dans sa transparence et dans la rigueur de la demarche. Les hypotheses sont sur la table et la regle de Bayes assure la coherence du raisonnement. La puissance de cette approche est renforcee par la distribution predictive a posteriori qui n'a pas d'equivalent classique. Ainsi, la distribution de Polya est la distribution predictive a posteriori du modele beta-binomial. Ce dernier permet de traiter des problemes OU l'observable est une variable aleatoire dichotomique : I'evenement d'interet se realise ou ne se realise pas. Ce modele est approprie quand la succession des observations constitue un processus de Bernoulli, le nombre d'essais etant fixe. La loi de Poisson etant un cas limite de la loi binomiale, le modele gamma-Poisson s'applique quand l'evenement dichotomique d'interet est rare. Sa distribution predictive a posteriori est la loi binomiale negative. Meme si l'hypothese « processus de Bernoulli» n'est pas toujours facile a justifier, ces modeles simples (mais pas simplistes!) sont utiles. Ainsi, la tragedie de Montroc nous a permis de batir un contexte decisionnel, certes fictif, mais riche d'enseignements. Bien que critiquable, la quantification du risque realisee ci-dessus a du sens. En tout cas, elle aurait pu alimenter les debats et influer sur la decision finale. Une decision est rationnelle s'il est clairement etabli qu'elle participe a la satisfaction de l'objectif declare en respectant un certain nombre de principes juges essentiels. Ainsi, la clarte du dialogue entre l'analyste et le decideur ; la pertinence des informations et le respect du cahier des charges sont des exigences qui nous semblent incontournables (Bernier et al., 2000). L'acceptation de la methode par toutes les parties n'est pas la moindre des difficultes, Elle repose en partie sur la comprehension qu'elles en ont et la representation graphique du modele va dans ce sens, C'est ainsi que Ie chapitre 3 precise la notion de reseau bayesian et introduit les variables latentes et la modelisation hierarchique.

Chapitre 3

Introduction a la modelisation graphique Ie modele de capture-recapture

• •

Prologue Les modeles graphiques associent la theorie des graphes, qui modelise des reseaux, a la theorie des probabilites, qui quantifie l'incertitude. L'idee fondamentale est la modularite : un modele complexe est construit en combinant des modeles simples. Les modeles graphiques eclairent parfaitement la notion dindependance conditionnelle. Le modele dit de capture-marquage-recapture constitue un exemple pedagogique d'autant plus intercssant qu'il recoit de nombreuses applications pratiques dans les sciences naturelles et humaines.

3.1

Introduction

Sans en formaliser la presentation, nous avons deja montre des modeles graphiques dans les deux premiers chapitres. Ainsi, a la page 11 de la section 1.3, nous avons rcprcscnte de trois facons differentes la probabilite conjointe d'une observable Y et d'un parametre () prenant ses valeurs dans l'ensemble des etats de la nature de dimension finie, 8. Nous avons retenu qu'un reseau bayesien, ou DAG, represente un modele statistique parametrique a l'aide de nceuds relies par des flechcs indiquant les liens de dependance entre des quantites incertaines. Le DAG lc plus simple relie un parametre a une observable (fig. 3.1). Le parametrc, (), et l'observable, Y, sont des nceuds stochastiques representee par

34

Pratique du calcul bayesien

des cercles ou des ellipses. La fleche indique une relation de filiation . En vocabul air e graphique, 8 est le nceud parent et Y est le nceud enf ant. Des que l'on fixe 8 on peut generer des valeurs y de l'observable Y . C'est en ce sens que 8 joue le role d'une « cause » et que l'observation y joue Ie role d 'un « effet », Apprendre quelque chose sur 8 revient a cherch er la loi conditionnelle de 8 sachant l'observation y en mobilisant eventuellernent une certaine expertise sur 8. Cela revient a invers er Ie sens de la fleche pui squ 'on remonte de l'effet vers la cause (fig. 3.1).

Parame trc

Modele

Observ able

Inference

Observ ation

Figure 3.1 - Le DAG Ie plus simple.

Exemple 3 .1 Si Y est le t emps qui s'ecoule entre deux manifestations d'un evenement dommageable, on peut le modeliser a l'aide d'une distribution exponentielle de parametre 8 :

[y I8] = 8exp( -8y)

'* E(Y) =

1/8

(3.1)

Si on ne disposait qu e d 'une seule observation (a deconseiller) , ce mod ele serait represents par Ie DAG de la figure 3.1. • La figure 3.2 montre un DAG plus sophistique. Les fleches doubles indiquent des operations logiques. Par exemple, p~ = P2 - m 2, r/ = T"J - ml . Les valeurs fixees son t representees par des carte s ou des rect angles. Ainsi , P2 est une const ante.

3.1.1

Vne courte digression

Dans la section 1.4 nous avons justifie Ie choix bayesi en , mais nous ne nous sommes pas encore vraiment interreges sur la pertinence de l' approche probabilist e. En fait , si la st at ist ique permet d 'interpret er un phenomene naturel, elle ne l'explique pas (Robert , 2006)! L'exemple suivant va nous permettre d'illustrer le propos.

3. Introduction

a la modelisation graphique

35

Figure 3.2 - Un reseau bayesien plus soph istique. Exemple 3.2 La troisieme loi de Kepler (1571-1630) , decouverte en 1618, repose sur l'analyse des donnees de Tycho Brahe (1546-1601) : quelle que soit la planete, le carre de sa periode de revolution, T, divisee par le cube de son demi grand axe, a, est un e constante. A l'epoque, Kepler disposait des donnees pour six plan etes (fig. 3.3) . L'alignement (en coordonnees logarithmiques) est remarqu able . •

100

satcme

+

Jupiter .

10

ars +

Venus . Mercu re

10

+

Demi grand axe (Terre = 1)

Figure 3.3 - La troisierne loi de Kepler. Imaginon s qu 'un st atisticien d'aujourd'hui ignorant tout de l'astronomie remonte le temps. Il propose aux conte mporains de Kepler une modelisation probabiliste du phenomena observe (fig. 3.3). Personn e ne sait que la course des planetes aut our du Soleil est det errninee par la loi de la gravitat ion universelle de Newton. Alors Ie statisticien propose un modele d'echantillonnage sense pouvoir reproduire les observation s. Il opt e pour le modele norm al deja rencontre au chapit re 1 (section 1.1) et davantage explicite au chapitre 9 (section 9.2).

36

Pratique du calcul bayesien

ln c, == In17+jJlnTi +ci,

e,

rv

iid

dnorm(O,a)

Ce modele probabiliste et le traitement statistique bayesien qui en decoule (Ie prior est non informatif) conduisent aux resultats suivants (tableau 3.1).

ic« jJ

17 a

2.5 0.650 0.976 0.014

50 0.665 1.003 0.025

97.5 0.681 1.032 0.065

Tableau 3.1 - La troisiemc loi de Kepler.

Pour aboutir a la troisicme loi de Kepler, le statisticien devrait « oublier » les incertitudes et decreter que 17 == 1 et jJ == 2/3. Le statisticien ne se le permettra pas. Ainsi, comme le fait remarquer (Robert, 2006), apposer un modele probabiliste sur un phenomene inexplique peut paraitre tres reducteur. II est vrai que quand on connait la mecanique newtonienne, notre modele de regression semble d'autant plus demuni qu'illui est impossible de reconstruire la loi de la gravitation universelle (meme pas la troisieme loi de Kepler) a partir des observations. En d'autres termes, un modele probabiliste n'explique jamais le phenomena reel d'interet ! II se contente d'en fournir une representation a des fins operationnelles, Uranus u'etait pas connue a l'epoque de Kepler, mais, en mode predictif, notre modele calcule sa distance au Soleil a partir de sa periode de revolution. L'erreur relative mediane est inferieure a 2 %! Ca ne vaut pas la loi de Kepler, mais ce n'est pas si mal si on n'en dispose pas. Bien entendu, l'approche probabiliste exige que le modele d'echantillonnage choisi « convienne » au probleme etudie. Ce choix est capital! Le metier et le bon sens sont ici des atouts precieux.

3.2

Principe de la modelisation graphique

Le lecteur interesse trouvera dans (Cowell, 1998) une excellente introduction a la modelisation graphique, notamment l'exemple 3.3 dont sont issues les figures 3.7 a 3.9.

3.2.1

L'independance conditionnelle

La figure 3.4 illustre la notion tres importante d'independance conditionnelle. Pour apprendre quelque chose sur Z, il n'est pas necessaire de considerer Y si on dispose de X. On notera

(ZIX == x) .L (YIX == x) ou, plus simplement

(3.2)

3. Introduction it la modelisation graphique

(Z 1. Y)

IX

37

(3.3)

qui se lit : «Z est conditionnellement independant de Y relativement it l'information X == x », En d'autres mots, disposant de l'information X == x, un apport d'information sur Y, soit Y == y, ne modifie pas l'incertitude sur Z

[ZIY == y, X == x] == [ZIX == x]

(3.4)

Figure 3.4 - V n heritage : Ie nceud Zest conditionnellement independant du nceud Y sachant le nceud X.

A contrario, la figure 3.5 montre que pour apprendre quelque chose sur Z il faut considerer les noeuds X et Y.

Figure 3.5 - Une naissance : le nceud Z depend des nceuds X et Y.

La figure 3.6 montre que Zest independant de X sachant Y.

Pratique du calcul bayesien

38

Figure 3.6 - Une chaine : le nreud Zest independant du nceud X conditionnellement au nceud Y.

3.2.2

Du reseau bayesien

a la

loi conjointe

Un reseau bayesien (DAG) a une structure definie comme suit: a chaque nceud X est associee une distribution de probabilite conditionnelle dont Ie conditionnement porte uniquement sur les parents du nceud Pr (Xlpa (X))

(3.5)

La distribution jointe d'un ensemble de nceuds, disons U, est le produit de

toutes ses distributions conditionnelles (fig. 3.7) :

Pr(U)

==

IIpr(Xlpa(X)) x

Figure 3.7 - Distribution jointe d'un reseau bayesian.

(3.6)

3. Introduction

a la

modelisation graphique

39

Exemple 3.3 Soit U l'ensemble (conjonction) des nceuds : U == {A,B,C,D,E,F,G,H,I}. Pr (U) == Pr (A) Pr (B) Pr (C) x Pr (DIA) Pr (EIA, B) Pr (FIB, G) x Pr (GIA, D, E) Pr (HIB, E, F) Pr (JIG, F)

(3.7)

II est interessant de noter que la marginalisation sur un nceud sans descendant revient a enlever ce nceud du reseau ainsi que tous les liens y aboutissant. Par exemple en marginalisant sur le nceud H (fig. 3.8) : Pr (A, B, G, D, E, F, G, I) ==

:L Pr (U)

(3.8)

H

Figure 3.8 - Marginalisation sur un nceud.

On peut toujours ecrire un roseau bayesicn en placant les lettres de telle sorte que les parents d'un nceud le precedent dans la liste. Un tel arrangement est une typologie. Pour un DAG donne, il y a de multiples typologies. Ainsi, par rapport a la figure 3.7, (A, B, C, D, E, F, G, I), (B, A, C, F, E, D, I, G) et (C, A, B, E, D, F, G, I) conviennent. •

Deux proprietes markoviennes 1. Independonce conditionnelle. Un nceud est conditionnellement independant de ses non-descendants etant donne ses parents :

(E -.L nd (E)) Ipa (E)

(3.9)

En d'autres mots, disposant de l'information pa (E), un apport d'information sur nd (E) ne modifie pas l'incertitude sur E.

40

P ratique du calcul bayesien 2. Modula rit e (fig. 3.9). La loi d 'u n nceud sa chant le reste du resea u ne depend qu e de ses par ents, de ses enfants et des copar ents de ses enfant s. PI' (EIA, B , C, D , F, G, 1)

= PI' (EIA, B , D , G)

(3.10)

00 0

r!J t5'(j K ~

\8

Figure 3.9 - Modularite : loi d'un nceud sachant Ie reste du reseau. Dan s le chapit re 8, nous appliquons ces prop rietes - en det aillant les operat ions - a I'exemple des sa umons (voir p. 155 et suivantes) .

3 .2.3

DAG et variables latentes

Un modele st atisti que bayesien est utilement represent e par un DAG . Les qu ant ites incert ain es constituent des noeuds stochas tiques. Les parametres du mo dele sont des nceuds san s par ent et les observables sont des nceuds sa ns enfant . Tout nceud stochastique qui n 'est ni un par am etre ni une obser vable est une vari able lat ente. Ainsi, dan s la sect ion suivante, nous verr ons qu e le cardinal! d 'un ensemble qu 'on ne peu t recenser apparait comme une varia ble latent e dan s le mod ele dit de copture-m arquaqe-recopture. Le plu s souvent, l'int roduct ion de ce ty pe de variable dans le modele est justifiee pa r le souci de prendr e en compte des influences cachees qui affectent l'observable (voir chap. 8). Quan d on le peut (c'est une question de clarte du DAG ), les paramet res du mod ele forment la couche super ieure du DAG et les observa bles, la couche inferieure, Les variables latentes constit uent une couche interrnedi aire, pri se en sa ndwich ent re les par am etres et les observables, qui confere au modele une st ructure hierarchique. 1

Le ca rd ina l d ' un ensemble fin i E des ign e Ie nombre d 'element s de E .

3. Introduction a la modelisation graphique

3.3

41

Le modele de capture-recapture

Dans son application la plus courante, il s'agit d'estimer la taille d'une population statistique hors recensement.

Remarque 3.1 Bien que nous soyons encore formellement dans la premiere partie de cet ouvrage, nous devrons utiliser l'ordinateur pour resoudre le modele de capture-marquage-recapture. Que le lecteur veuille bien ne pas trouver la une incoherence de notre part. La locution « de la plume a la souris» doit etre comprise comme un cheminement et non comme une separation nette. Le recours a l'ordinateur est done preponderant dans la seconde partie de cet ouvrage sans etre completement exclu de la premiere.

3.3.1

Mise en situation

Le recensement est une operation statistique de denombrement d'une population generalement realise a des fins decisionnelles. Les premiers recensements connus ont eu lieu des l'Antiquite, notamment a Rome, dans le but de connaitrc la richesse du pays, afin de repartir l'impot. Mais une telle operation exige du temps et consomme des moyens importants quand elle n'est pas tout simplement impraticable. Uno alternative au recensement consiste a estimer la taille de la population d'interet a partir d'un double echantillonnage. On preleve au hasard, c'est-a-dire on peche'', un certain nombre d'individus que l'on remet dans leur milieu apres les avoir marques d'une manierc quelconque. Apres brassage, un second echantillonnage fournit un lot d'individus dont certains sont marques - ils sont recaptures - d'ou la denomination du modele. Sous certaines conditions, les effectifs des deux peches et les recaptures suffisent pour obtenir la distribution a posteriori de la taille de la population rl'interet. Ce

modele trouve de nombreuses applications pratiques dans les sciences naturelles et humaines.

3.3.2

La modelisation

Soit a estimer la taille, 1], d'une population donnee. II peut s'agir du nombre de poissons dans un lineaire de riviere, du nombre de sans-abri dans une ville, du nombre de chenes dans une foret, du nombre de declarations suspectes dans le ressort d'un percepteur, etc. En d'autres mots, 1] est le cardinal inconnu d'un ensemble bien defini qu'on veut inferer. Une premiere «peche» fournit un certain nombre d'individus que l'on marque d'une maniere quelconque avant de les « relacher » dans leur milieu. Soit ml ce nombre. Lors d'une seconde « peche » on prend P2 individus dont m2 sont marques, c'est-a-dire recaptures. 2 Ce modele est tres utilise en pisciculture, notamment pour contr6ler des peuplements ou s'assurer de I'efficacite des mesures de repeuplement.

42

Pratique du calcul bayesien

Hypotheses

A chaque individu du milieu (indice i), on associe une variable aleatoire de Bernoulli, disons Yik, qui prend la valeur 1 avec la probabilite 7Tik, s'il est capture a la k-ieme peche (k == 1,2), et la valeur 0 avec la probabilite complementaire s'il ne l'est pas. 1. Ces variables indicatrices sont independantes et identiquement distribuees :

Vi,Vj # i,Vk: Pr(Yik == IIYjk) == Pr(Yik == 1) ==

7T

(3.11)

2. Les deux peches sont independantes : (3.12) 3. II n'y a ni source, ni puits, ni emigration, ni immigration, c'est-a-dire TJ est invariant, au moins pendant la duree des operations. Un modele probabiliste de connaissance pour Y est fonde sur une loi de Poisson de paramctre A > 0, lui-meme tire dans une loi gamma dhyperparametres a > 0 et (3 > O. Cela revient a dire que la distribution de TJ est une binomiale negative. Un prior non informatif est obtenu en posant a == (3 ---+ 0 =} [A] ex A-1. La loi binomiale negative impropre de ce prior non informatif a pour

esperance 1 et a une variance infinie. Un premier modele dechantdllonnage mime la collecte des donnees Sous ces hypotheses : - les ml individus marques lors de la premiere peche sont les succes obtenus a l'issue d'une sequence de TJ cprcuvcs de Bernoulli a TJ fixe; la distribution de m, est binomiale, de probabilite 7T et d'ordre TJ (3.13) - les m2 individus recaptures lors de la seconde peche sont les succes obtenus dans une sequence de m.; cpreuvcs de Bernoulli a ml fixe (3.14) - les individus non marques et captures lors de la seconde peche, soit P~ == P2 - m2, sont les succes obtenus lors d'une sequence de TJ' epreuves de Bernoulli OU TJ' == TJ - ml est fixe (3.15)

a la

3. Introduction

modelisation graphique

43

Le DAG montre a la figure 3.2 (p. 35) representait, sans le dire, ce modele dans lequel TJ I == TJ - m1 et P2' == P2 - m2· Puisque la distribution conjointe d'un reseau bayesien est egale au produit des distributions de chaque nceud stochastique sachant ses noeuds parents, on a:

La distribution a posteriori des parametres suit (regle de Bayes) :

En posant" (3.18)

la vraisemblance s'ecrit : \] TJ! S ( [S, C I1f, A, TJ ex: (TJ _ c)! 1f 1 -

1f

)2TJ -

(3.19)

S

Un prior non informatif pour nest uniforme sur [0,1] et un prior non informatif pour A est proportionnel a A-1. Enfin, TJ est tire dans une loi de Poisson de parametre A. En substituant dans 3.17, il vient :

[TJ, n , AIs, ] c

ex:

A17-1exp(-A) S(

(TJ _ c)!

1f

1-

()

)217-s 1f

1'12C

TJ

(3.20)

Une double integration par rapport aux parametres n et A fournit la marginale a posteriori de TJ : 1

f(TJ)

[TJls, c] = K x (TJ _ c)! B (s + 1, 2TJ - s + 1) 1'12 C (TJ)

(3.21 )

ou la constante de normalisation, K, peut etre definie sur une grille de valeurs de TJ. La marginale a posteriori de nest facilement obtenue via l'algorithme suivant:

1. fixer N et compteur < -1 2. tant que compteur est inferieur

a N,

repeter :

- tirer une valeur de TJ dans [TJls, c] ; - tirer une valeur de nlTJ dans dbeta( nls - compteur < -compteur + 1. 3

A l'issue des operations, on sait qu'il

+ 1, 2TJ - s + 1) ;

y a au mains c individus dans le milieu.

44

Pratique du calcul bayesien

Remarque 3.2 La fonction factorielle rend l'infini pour les grands entiers naturels (sur mon ordinateur, x! E N {:} x ~ 170). Un changement d'unite est possible (par exemple, travailler en dizaine d'individus) a condition de remplacer la fonction factorielle par l'integrale d'Euler. - La constante de normalisation s'ecrit :

~

K = B (8 + 1, 2c _ 8 + 1) +

LJ

1]=c+l

B (8 + 1, 21] - 8 + 1) (17 - c) B (c, 17 - c)

(3.22)

- La marginale a posteriori de 17 suit :

[1]18, c] =

{

-k B (s + 1, 2c - s + 1) {:} 17 == K(~-c) B (8 + 1, 21] -

C

8 + 1) / B (c, 1] - c) {:} 1]

>c

(3.23)

Capture et recapture par un echantillonnage multinomial

A l'issue des deux peches, conditionnellement a 17, un individu quelconque est necessairement dans un des quatre etats possibles : capture-capture (cc), capture-manque (cm), manque-capture (mc), manque-manque (mm). Les effectifs de ces quatre etats sont donnes dans le tableau 3.2 OU le nombre d'individus jamais captures est inconnu. 1\11

c

m

c

m2

ml-m2

ill

P2 -m2

Total

P2

Total

17 - ml - P2 + m2 17 - P2

ml

17 - m.; 17

Tableau 3.2 - Une truite est capturee (c) ou manquee (m).

Le tableau 3.3 donne les probabilites associees l\ll c m

c 1r'2

(1 - 1r) 1r

a chacun de ces quatre etats,

m 1r (1 - 1r) (1-1r)'2

Tableau 3.3 - Probabilites des etats,

Des lors, les effectifs du tableau 3.2 sont vus comme le resultat de 17 tirages independants dans une loi multinomiale de parametre

La figure 3.10 montre Ie DAG dans lequel y represente Ie vecteur des effectifs. Clairement, 1r .L A. Le prior de 1r est une distribution beta de parametres a et b et celui de A est une distribution gamma de parametrc p et q. Ils sont non informatifs en posant p == q == 0 et a == b == 1.

3. Introduction

a la modelisation graphique

45

Figure 3.10 - Le modele de capture-recapture: echantillonnage multinomial.

La vraisemblance s'ecrit :

[yl1T,7]J

ex:

7]! 1T s (1_1T)2 rJ(1] - s + m2)!

S

(3.24)

et ... c'est le meme modele que ci-dessus (eq, 3.19).

3.3.3

Applications

Estimation de l'incidence de la tuberculose pediatr'ique en BasseNormandie Les objectifs de cette etude (Brouard et al., 1995) etaient de verifier la pratique de la declaration obligatoire (DO) et le respect des critercs de declaration, d'estimer l'incidence de la tuberculose pediatrique en Basse-Normandie par la methode de capture-recapture et ainsi I'exhaustivite de la DO. Deux sources de donnees ont ete explorccs :

1. les DO enregistrees dans les directions departementales des affaires sanitaires et sociales (DDASS) des trois departements de la Basse-Normandie; 2. les enregistrements d'isolement de Mycobacterium tuberculosis (MT) sur l'ensemble des laboratoires d'analyse medicale (LAM) de Basse-Normandie. Cette enquete, du type retrospectif, est fondee sur les cas pediatriques identifies entre le 1er janvier 1992 et le 30 juin 1993. Les resultats sont les suivants (tableau 3.4) OU : - R == 6 est le nombre de cas diagnostiques selon la source 1 (DO); - S == 8 est le nombre de cas diagnostiques selon la source 2 (LAM) ; - C == 4 est Ie nombre de doublons. Pour estimer le nombre de cas, N, les auteurs utilisent les formules elaborees par Chapman et Seber en 1949 :

46

Pratique du calcul bayesien

LAM LAM Total

DO

DO

Total

C

N2

8

N1

R

N

Tableau 3.4 - Tuberculose pediatrique en Basse-Normandie.

ic...;

N ± ZI-a/2VVar (N)

N

(8+1)(R+1) -1 C+ 1 (8 + 1) (R + 1) N 1N2 (C+1)2(C+2)

Var (N)

(3.25) (3.26) (3.27)

Les resultats sont les suivants: N == 11.6, Var (N) == 3.36 et les bornes d'un intervalle de confiance a 90 % sont respectivement 8.6 et 14.6 cas. Notons que les auteurs les presentent de facon un peu plus optimiste puisqu'ils concluent : « par la methode capture-recapture, le chiffre des tuberculoses pediatriques est cstime a 11, le calcul de la variance donne un ecart de ce chiffre de plus ou moins 3 (11 ± 3) » . Avec prior non informatif, le modele de capture-recapture developpe donne les resultats suivants : a == 0.1 N 1T

Pa/2 10 0.30

Pso

13 0.55

Pl-a/2 20 0.75

Tableau 3.5 - Estimation bayesienne de N.

Force est de constater que l'estimateur classique sous-estime le nombre moyen de cas et sa dispersion!

Evaluation de l'incidence du paludisme dans les arrnees francaises en

1994 Cette etude (Deparis et al., 1997) est fondee sur deux systemes reglementaires de surveillance epidemiologique : - le recueil et l'exploitation des donnees epidemiologique« des arrnees (REDEA); - la surveillance epidemiologique specifique du pal udisme (SESP). Ces deux sources de donnees sont supposees independantes (tableau 3.6). Pour estimer Ie nombre de cas, X, les auteurs utilisent les formules elaborees par Chapman et Seber en 1949. Soit a le nombre de cas declares dans les

3. Introduction

REDEA REDEA Total

a la

modelisation graphique

SESP

SESP

Total

238 186 424

242

480

47

X

Tableau 3.6 - Paludisme dans les armees francaises (1994).

deux systemes (a == 238) ; b, le nombre de cas declares uniquement a la SESP (b == 186); c, le nombre de cas declares uniquement au REDEA (c == 242).

(a+b+1)(a+c+1) -1 a+1 (a + b + +1) (a + c + 1) bc (a+1)2(a+2) Un intervalle de confiance mateur X) :

a 95 % suit

(3.28) (3.29)

(hypothese de normalite sur l'esti-

le95 == x ± 1.96sx Sur cette base, l'incidence annuelle du paludisme s'eleve intervalle de confiance a 95% egal a [803, 905].

X 7r

q5

Q50

Q95

750 0.42

860 0.52

1030 0.63

(3.30)

a 853 cas

avec un

Tableau 3.7 - Paludisme : estimation bayesienne du nombre de cas X.

Ici aussi, I'estimateur classique sous-estime le nombre moyen de cas et sa dispersion!

Epilogue Un modele statistique bayesien mime la nature en ce sens qu'il vise a genercr des donnees similaires aux observations reelles. Ce faisant, il permet d'interpret.er le phenomene d'interet, souvent dans une perspective decisionnelle. Un roseau bayesien ou DAG est une representation graphique astucieuse du modele. D'une part, il aide a sa conception: d'autre part, il favorise la multidisciplinarite, car le dessin est un langage accessible a tous. Dans le DAG, un parametre est un nceud stochastique sans parent et une observable, un nceud stochastique sans enfant. Les variables latentes sont des quantites incertaines, qui ne sont ni l'un ni l'autre. Elles constituent une couche interrnediaire, prise en sandwich, entre les parametres et les observables. Une Heche indique un lien causal entre

48

Pratique du calcul bayesien

deux noeuds stochastiques : l'etat du nceud recepteur est conditionnel a celui du nceud emetteur, L'inference bayesienne consiste a inverser le sens des fleches, c'est-a-dire a remonter vers les parametres (causes) en partant des observations (effets), en tenant bien sur compte de l'expertise (priors). Le DAG met bien en evidence les notions d'independancc conditionnelle et de modularite, La distribution conjointe de tous les nceuds stochastiques est simplement egale au produit de chaque nceud connaissant ses nceuds parents. Elle s'exprime donc en termes de distributions conditionnelles et marginales. Pour apprendre quelque chose sur un nceud stochastique, il suffit de connaitre ses parents, ses enfants et les coparents de ses enfants. C'est modularite est mise a profit dans les logiciels comme WinBUGS. Nous avons illustre les avantages du reseau bayesian en nous appuyant sur le modele de capture-marquage-recapture. C'est un modele tres utile pour les sciences naturelles et humaines puisqu'il permet d'inferer la taille d'une population statistique inaccessible par recensement. S'il est conceptuellement facile a comprendre, il necessite deja un recours a l'ordinateur, car une solution analytique complete implique d'integrer la relation (3.20) par rapport a 1r et 7], ce que personne ne sait faire. D'une maniere tres generale, les modeles realistes n'ont pas de solution analytique. Les reseaux bayesiens sont done indissociables des methodes modernes de calcul sur ordinateur. Dans le chapitre 4, nous presentons quelques methodes de reference du calcul numerique stochastique. Pour cela, nous nous appuierons sur le modele luieoire et quelques modeles lineaires generalises.

Chapitre 4

Pratique du calcul des lois a posteriori Prologue Qui dit modeles rcalistes, dit aussi difficultes calculatoires. Le but de ce quatrieme chapitre est de donner un apercu des principales familles de methodes d'approximation des distributions a posteriori. Dans lc cas tres particulier OU le prior est non informatif et que la taille de l'echantillon est grande, la densite a posteriori peut etre approchee par une loi normale multidimensionnelle. Cette approximation asymptotique repose sur les proprietes des estimateurs du maximum de vraisemblance (section 4.2). A l'erc des ordinateurs personnels puissants, cette approximation - fondee sur des hypotheses assez restrictives est avantageusement abandonnee au profit des methodes numeriques stochastiques. Ce sont d'abord les methodes de Monte-Carlo par chaines de Markov (MCMC). Ces techniques de simulation avec dependance sont presentees en section 4.3, notamment l'algorithme general de Metropolis-Hastings et Vechantillonnage de Gibbs. Ces deux algorithmes sont d'ailleurs implantes dans Ie logiciel WinBUGS. Les techniques classiques de simulation avec independance ou methodes de Monte-Carlo (MC), issues de l' echantillonnage potidere, avec ou sans re-echantillonnage, ont eu plus rccemment des developpements importants sous le nom generique de methodes des particules (section 4.4).

4.1

Introduction

Un modele de connaissance - on dit aussi modele d'echomiillonnaqe - est une famille de lois de probabilite parametree par () E 8 OU 8, souvent appcle ensemble des eiais de la nature, est de dimension finie : dim e == d E No. Par consequent, le parametre () est tantot un scalaire (d == 1), tantot un vecteur (d > 1). La notation est la meme, et c'est le contexte qui fait la difference.

50

Pratique du calcul bayesien

Disposant d'un modele de connaissance et d'une loi a priori pour (), la reactualisation du savoir sur () associee a une information y est donnee par la regle de Bayes (chap. 1, p. 12). Cette distribution a posteriori est Ie socle sur lequel repose l'aide a la decision en avenir incertain. Ainsi, la distribution predictive a posteriori (chap. 1, p. 12) quantifie les chances d'observer une future valeur fj quand on dispose de l'information y :

Wly]

=

l WIB]

[Bly] dB

(4.1)

D'une maniere plus generale, le statisticien bayesien est amenc des integralcs de la forme E (h (B) Iy) =

r

Je h (B) [Bly] dB =

Ie h (()) [yl()] [()]d() Ie [yIB] [B]dB

a calculer (4.2)

ou h (()) est une fonction reelle, Pratiquement, des solutions analytiques n'existent que pour des modeles particuliers, les structures non hierarchiques de la famille exponentielle (Parent et Bernier, 2007). Certes, certains logiciels offrent des algorithmes dintegration numeriqucs de bonne qualite, Mais, on l'a vu, les modeles bayesiens realistes prcsentent souvent une structure hierarchique impliquant des variables latentes. Du point de vue calculatoire, les variables latentes peuvent etre considerees comme des parametres supplementaires (::::} dim 8 » 1). Or l'imprecision des methodes dintegration numerique croft dramatiquement avec la dimension de 8 (Robert, 2006). Des lors que Ie nombre des quantites incertaines (parametres + variables latentes) excede quelques unites, les methodes dintegration numeriques sont supplantees par les methodes numeriques stochastiques d'approximation, MCMC et particules. Le principe de base des methodes numeriques stochastiques d'approximation est simple. Soit a rcsoudre l'integrale (4.2). Si on considere une suite de variables aleatoires independantes (()1, ... .U": ... ) et distribuecs selon la loi a posteriori de (), on obtient un echantillon de nombres reels en tirant au hasard une valeur dans chacune d'entre elles. La moyenne arithmetique de leur image par la fonction h, soit G

~ Lh (B i ) i=l

1

converge (presque surement ) vers la cible quand G nombres), ce qui justifie l'approximation : G

~ ~ h (B i ) ~ E (h (B) Iy) =

l

---+ 00

(loi des grands

h (B) [Bly] de

(4.3)

1 La convergence presque sure est analogue a la convergence simple de l'analyse mathematique, sauf en quelques points. Elle entraine la convergence en loi.

4. Calcul des lois a posteriori

51

De plus si la variance a posteriori de h (()) est finie, disons a 2 > 0, le theoreme central limite nous dit que cette moyenne arithmetique est distribuee selon une loi normale, de variance a 2 /G et que l'ordre de grandeur de l'erreur relative est 1/V'G, ce qui permet de calculer des intervalles de confiance sur I'integrale. Le principe de base des methodes de Monte-Carlo par chaines de Markov est analogue si ce n'est que, cette fois, la chaine (()l, ... .B": ... ) est generee par un noyau de transition [OJ I()j-l] dont on considere la moyenne (4.4) Cette moyenne converge vers la cible quand G ---+ 00 pour autant que Go soit assez grand et que la chaine de Markov possede la propriete dite d'ergodiciu: (Robert et Casella, 1999), une propriete generiquement verifies sous des conditions peu strictes pour les chaines de Markov homogenes. L'echantillon prealable de i == 1 a i == Go, laisse de cote, s' appelle echantillon de chauffe.

Remarque 4.1 Les processus stochastiques sont des modeles permettant d'etudier les phenomenes aleatoires evoluant au cours du temps. Parmi ceux-ci, les chaines de Markov sont les modeles (a temps discret) les plus simples, lorsqu'on abandonne l'hypothese dindependance. Pour plus de details, on consultera avec profit (Foata et Fuchs, 1998). Remarque 4.2 II importe de remarquer que la formule 4.3 s'applique tout aussi bien au calcul d'une probabilite P(A) == Prob(() E A) par Ie biais d'une c . . di . h(O) == I A (0) == 01 si . tonction In icatrice si e eE t/:. A A ' puisque : Prob(B E A)

=

L

IA(B) [Bly] dB

(4.5)

Certes, le praticien des sciences experimentales est souvent plus interesse par les sorties des modeles que par les mathematiques qui les soutiennent et c'est certainement une des raisons du succes planetaire du logiciel WinBUGS (Spiegelhalter et al., 2003). Ce logiciel, gracieusement mis a la disposition de la communaute scientifique, permet un apprentissage rapide du raisonnement conditionnel bayesien, II distingue clairement la partie creative, c'est-a-dire l'elaboration du DAG, element de l'interface graphique Doodle, de la partie calculatoire. Pour son utilisateur, l'estimation des quantites incertaines est transparente. II lui suffit de savoir qu'une marche aleatoire dans l'ensemble des etats de la nature, 8, genere une chaine de Markov ()l, ... ,()j, . . . ,()N et que, hormis une periode dite « de chauffe », a ecarter puisque cette chaine part d'un point arbitraire, la repartition des () ainsi generes converge en distribution

52

Pratique du calcul bayesien

vers sa cible. Par exemple, l'histogramme marginal normalise de chaque composante de () approche d'aussi pres que l'on veut (en augmentant Ie nombre de simulations N) la loi marginale a posteriori de cette composante. WinBUGS est donc un excellent outil pedagogique qui peut certainement resoudre pas mal de vrais problemes. Mais, comme ses concepteurs, nous defendons l'idee que l'emploi intelligent de la souris demande un minimum de comprehension des methodes sous-jacentes. Elles sont indispensables a l'etudiantjchercheur qui souhaite ecrire ses propres codes, par exemple en R (R Development Core Team, 2009). Le lecteur interesse par lcs fondements theoriques et les subtilites des methodes de calcul bayesien consultera avec profit les ouvrages specialises, notamment : (Tanner, 1996), (Robert et Casella, 1999), (Chen et al., 2000), (Gelman et al., 2004), (Robert, 2006), (Parent et Bernier, 2007).

4.2

Quand Ia vraisernblance fait Ie posterior

L'inference bayesienne mobilise deux sources d'information : d'une part les donnees, via la vraisemblance, et d'autre part le savoir de l'expert via la distribution a priori sur les parametres et via les hypotheses structurelles sur lesquelles repose le modele utilise. Un exemple tres simple va nous montrer que quand la taille de I'echantillon est grande (n - t (X)) ou quand Ie prior est tres vague, il y a un lien lineaire approche entre la log-densite a posteriori et la logvraisemblance. Regardons les poids respectifs du prior et de la vraisemblance sous deux configurations. 1. Quand la taille de l'echantillon est grande (n ---t (X)), l'influence du prior s'estompe et c'est la vraisemblance qui fait le posterior. Exemple 4.1 Soit y un rz-echantillon iid issu d'une distribution exponentielle parametree par () :

[yl()]

== ()n exp ( -ny())

Le prior conjugue est une distribution gamma

La reglc de Bayes donne la distribution a posteriori de () :

[()Iy, a, b] ex

()n+a-l exp (-

(ny + b)B)

On reconnait la forme analytique d'une nouvelle distribution gamma. Si la taille de l'echantillon est telle que n » a et n » b (cette condition est verifiee lorsque n ---t (X)), alors n+a-l ~ n et ny+b ~ ny. Dans ce cas, le posterior et la vraisemblance ont la meme forme analytique. Comme

4. Calcul des lois a posteriori

53

la vraisemblance est dcfinie a un facteur de proportionnalite pres, on voit apparaitre un lien lineaire approche entre la log-densite a posteriori et la log-vraisemblance

[OIY, a, b] ex [yIO]

=?

In [OIY, a, b] == In [yIO]

+ cte

• 2. La vraisemblance fait aussi le posterior quand l'etat de connaissance sur la problematique en main n'autorise qu'un prior tres vague, c'est-a-dire quand le prior est peu informatif.

Exemple 4.2 Dans l'exemple 4.1, un prior vague est obtenu en faisant tendre les parametres a et b vers o. La forme analytique de la densite a posteriori devient [Oly] ex on-l exp (-nyO) Quand nest assez grand, alors n et n - 1 possedent le meme ordre de grandeur : on retrouve alors Ie lien Iineaire approximatif entre la log• densite a posteriori et la log-vraisemblance.

4.2.1

Approximation asymptotique de la densite

a posteriori Dans un probleme realiste, la quantite de donnees n'est jamais infinie. Dire que la taille n de I'echantillon est grande est une assertion reposant sur les proprietes asymptotiques du modele en main. Ces proprietes ont ete utilisees tres tot en Statistique. Le premier utilisateur en a ete Laplace qui, independamrnent de son auteur historique, a retrouve la formule de Bayes et utilise les principes dinference bayesiens (Sivia, 1996). L'interet de ces proprietes asymptotiques est d'etablir une relation lineaire approchee entre la Iog-densitc a posteriori et la log-vraisemblance. Cette relation Iineairc est fondee sur des hypotheses generales concernant cette vraisemblance lorsque n est grand. Nous en faisons une presentation heuristique avant de preciser les conditions theoriques qui en assurent la validite pratique. Pour simplifier, nous ne traiterons ici que le cas d'un echantillon d'observations suppose iid. Nous le noterons y == (Yl,··· ,Yn).

Le parametre du modele d'echantillonnage est unidimensionnel Soit un modele statistique bayesien caracterise par un parametre 0 unidimensionnel. La densite a posteriori de 0 est donnee par la regle de Bayes :

- [yIO] est la vraisemblance de I'echantillon y ; - [0] est la densite a priori du parametre 0;

54

Pratique du calcul bayesien

- [y] est la constante de normalisation :

[y] =

L

[yle] [e] de

En prenant Ie logarithme des deux membres : In [ely] == In [yle]

+ In [e] + cte

(4.6)

Supposons que la densite a posteriori de eest unimodale. Un developpement de Taylor au voisinage de son mode, disons donne/ :

e;,

e;

Pour Ie - I petit, on peut negliger les termes d'ordre superieur et en tenant compte de la regle de Bayes (eq. 4.6) :

In[ely]

~

In [yle~] + (e -

e~) 8lna~le] + a~e[e] lo=oz 2

~ (e _ e*) 2 8 In [yle]

+2

-l-In

8e

y

2

2

8 In [e] I + 8e 2 ()=::()Z

[e;] + cte

"-v--' cte

Si la densite a priori de e est « plate» au voisinage du mode premiere est nulle en ce point (les derivees superieures aussi) :

e;, sa derivee

Dans le cadre asymptotique (n ---t 00), l'influence du prior sur le densite a posteriori est tres faible et le mode de la vraisemblance, disons By, se confond avec le mode du posterior,

e; :

2

1 ( e- e Ay ) 2 8 In[ely] ~ In [yleAy ] + 2 ae 2 In[yle] Io=iJ + cte y

Le terme impliquant la derivee premiere a disparu, car By est le mode de la

vraisemblance

(to In[yIell o=fjy =

0) .

En posant

2

On suppose que toutes les derivees existent.

4. Calcul des lois a posteriori

55

on obtient

(4.7) Remarque 4.3 Le terme P script en le notant

(By)

ne depend que des donnees. On simplifie le

Py .

Pour revenir ala densite a posteriori de (), il suffit de prendre l'exponentielle des deux membres :

[ely]

~ cte x exp ( - ; (e _By) 2)

(4.8)

On reconnait la signature fonctionnelle d'une densite normale, localisee sur

By et de precision P (By) :

II est important de se rappeler que cette approximation n'est valable que sous les hypotheses: n grand, un seul mode et prior plat au voisinage de celui-ci.

Exemple 4.3 On verifie sans peine que l'estimateur du maximum de vraisemblance d'un echantillonnage exponentiel iid est {) ~ l/Y. Un developpement de Taylor jusqu'a l'ordre 2 de la log-vraisemblance autour de {) s'ecrit

In [yle] =In OU

0

[YIB] - 2~2

(e-ef +o(h)

(h) rcprcsentc un infiniment petit par rapport a h == () - () : lim o(h) ~ 0

h---+O

h

(4.9)

Puisque Ie premier terme du developpement est constant, on a done obtenu le resultat approche suivant

[ely] c::: exp ( - 2~2 (e - Bf)

(4.10)

On reconnait le terme caracteristique d'une loi normale unidimensionnelle, de moyenne {) ~ 1/Y et de precision T ~ n / {)2. A partir de donnees simulees par n tirages aleatoires indcpendants dans une loi exponentielle de moyenne J-L ~ 0.2, la figure 4.1 montre l'influence de la taille de l'echantillon sur la qualite de l'approximation asymptotique du « vrai » posterior gamma par une loi normale.

•

56

P ratique du calcul bayesian

15 r---,-------.------.---,--,r== -

~

/

10

II

c::

\

= = =======n

-

1-

-

Posterior gamma

I

- Approximation normale

\~

5 b /

0.05 3 ,----,-

0.1 ------.-

2 (r)

II

c::

/

0 ---

o

/

/

/

/

0.15 -----.--

0.2

0.25

0.3

0.35

0.4

0.45

- ,-

- ,-

-,--

--,--

---,---

---,------,

/ - --- ./

e

0.5

.........

---"-=~-

- - - Posterior gamma - - Approximation normale

/

0.05

01

0.15

0.2

0.25

e

0.3

035

0.4

0.45

0.5

Figure 4.1 - Approximation asymptotique du posterior d'un echantillonnage exponentiel (n = 30 et n = 3.)

Le parametre du modele d'echantillonnage est multidimensionnel Soit () = (()1, ' " , ()d)T E 8 Ie par am etre d 'un mod ele statist ique bayesien (dim 8 = d) et soit y un n-echa nt illon ii d. Si Ie pr ior , [()], est non inform atif, to ute l'inform ation disponible pour quant ifier l'incer ti tude sur () est , encore ici, portee par l'echantillon et la regle de Bayes s'ecrit : rely] ex [yl()] =:;. In [()Iy] = In [yl()]

+ cte

L'id ee est encore d 'approcher la distributio n a posterio ri de () par un e loi normale, cette fois multivariee, localisee sur Ie mode de la vr aisembl ance, ()y , et de mat rice de precision P . Un developpement de Taylor de la log-vraisembl ance autour de son mod e ()y j usqu'a l'ordre 2 donne :

(4.11) Dans cet te expression, P est une matrice sym etrique definie positive, dite matrice de precision . Son terme genera l s'ecrit : P' 'J· · -- _ aaoIn[y ao,OI I 2

i

j

o=O y

, Z. ,J. = 1, 2, " ' , d

(4.12)

En revenant au post erior , on a I'approximation

(4.13)

4. Calcul des lois a posteriori

57

On reconnait la forme structurelle canonique d'une loi normale multivariee, localisec sur le mode Oy et de matrice de precision P == ~-l ou ~ est la matrice de variance-covariance de cette distribution multinormale, soit exactement (4.14)

au Ip I est

le determinant de Ia matrice

p.

Exemple 4.4 Soit un n-echantillon iid selon une loi normale de parametrc

o== (J-L, T). En mobilisant la moyenne et la variance empiriques, la log-vraisemblance s'ecrit

n In T In [yIO] == 2"

nr 2

-

2+ (fj - J-L) 2)

( Sy

Le calcul de toutes les derivecs qui nous interessent donne al~:le] = nr (y _ p,) ; 81n[yI8] 8T

Le mode

By est

==

.!!:- _

2T

!! 2

(s2 + (-Y _ J-L )2)., Y

solution du systeme

La matrice de precision

Py

suit

L'approximation asymptotique de la densite a posteriori de 0 est donc une loi normale bidimensionnelle, localisee sur By et de precision Py . •

4.2.2

Fondements de ces approximations

Preambule - La recherche des extrema est fondamentale en statistique.

Definition 4.1 Soit f une application de classe C2 definie sur JRd a valeurs dans]R et soit x == (Xl, ... ,Xd) E ]Rd. On appelle matrice hessienne de f en x la matrice des derivees partielles secondes, c 'est-a-dire la matrice H (x) de terme general : (H (x )) ij --

2

8 f(x)

8Xi8xj'

.. -

'l,

J -

1 2 d , ,"',

C'est la matrice d'une forme quadratique sumetrique.

58

Pratique du calcul bayesien - Soit une experience alcatoire qui fournit l'observation y. Fisher (1925) propose de mesurer l'information apportee par cette observation au parametre e du modele statistique choisi pour mimer cette experience, [yle] . Definition 4.2 Si le domaine de l'observable Y ne depend pas du parala quatuiie d'information de Fisher apportee par l'information metre Y == y sur e est une matrice sumetrique definie positive, dite matrice d'information de Fisher, reliec a la matrice hessienne comme suit:

e,

I (e) == -E {H (e)} oii l' esperance est prise par rapport

a l' observable Y.

Fondements Les formules heuristiques precedentes sont fondees sur la theorie asymptotique du maximum de vraisemblance (Berger, 1985).

Theoreme 4.1 Sous un certain nombre d'lurpotheses generales concernant le modele de vraisemblance [yIO], si un cchaniillon iid, soit y, est tire du modele particulier [yleo], alors il existe une solution By de l'equation

8 In [yle] 8e == 0 qui converge en probabilit« vers 00 quand n ---+ 00. De plus, le vecteur By-Oo converge en loi vers une distribution normale multivariee, localisee sur 0 et de matrice de precision eqale a la matrice d'information de Fisher I (eo) .

Complements - En suivant (Berger, 1985), la matrice de precision P (eq, 4.12) est appelee matrice d'information empirique. Elle peut etre utilisee comme approximation asymptotique de la matrice de Fisher dans une expression approchce de la distribution a posteriori des parametres. - L'information de Fisher d'un n-echantillon iid est simplement (4.15)

Cette matrice intervient dans la construction de priors non informatifs. C'est ainsi que (Jeffreys, 1939) a propose la construction d'un prior vague a partir du determinant de 1(0) :

[0] == Jdet I (e) == II (0)1

1 2 /

(4.16)

Ce prior impropre a comme propriete de fournir une inference insensible reparametrisation du modele de vraisemblance.

a une

4. Calcul des lois a posteriori

59

- Parmi toutes les conditions theoriques fondant ces proprietes (Cramer, 1946), la plupart sont des hypotheses de regularite mathematique des convergences assez generales. II en est toutefois une qui est tres imp ortante sur le plan pratique: le domaine des observables ne doit pas dependre du parametre (). - Maintenant, que veut dire pratiquement l'expression n grand? II n'y a pas de regie, seulement des cas d'cspcce, II faut un peu d'experience, C'est pourquoi sur le plan bayesian il peut etrc prudent de conforter ces calculs approches par les resultats d'algorithmes de Monte-Carlo qui ne sont pas tributaires de I'hypothese n grand.

Exemple 4.5 (Exemple 4.4 continue). Calculons la matrice d'information de Fisher et Ie prior vague de Jeffrey. La matrice hessienne Hey est HCIL,r) =

-n

T

[

fL - Y

j.L-Y] 2~2

Puisque E (y) == u, l'information de Fisher et le prior de Jeffreys sont immediats I(tL,7)

== n

[~

1], 27 2

[j.L, T]

ex: Jdet I(tL,7) ex: 1/

VT

On remarquera que ce n'est pas Ie prior obtenu quand on fait tendre les parametres d'une loi gamma vers zero. En fait, les composantes du vecteur () apparaissant dans la matrice de Fisher sont inconnues. Aussi, on peut remplacer () par son mode O.

e= (fJ, l/s~)

Ie

=

n [1/;~

SD2]

La matrice de Fisher I et la matrice d'information empirique P sont parfois identiques.

•

4.2.3

Estimation asymptotique des parametres d'une population gamma

Exemple 4.6 Operation Sources", Soit un n-echantillon iid issu d'une population gamma, de parametre de forme a > 0 et de parametre d'echelle inverse (3 > o. II s'agit ici de la concentration en nitrates relevee dans n == 94 points d'eau repartis sur le territoire belgc en mars 1994. Les moyennes geornetrique et arithmetique des observations etaient respectivement 9 == 25.4 mg/I, fj == 39.3 mg/I. On demande d'estimer lc posterior par une loi normale bivariee ainsi que la probabilite 1r qu'une nouvelle observation depassc la norme Yo == 50 mg/I (fig. 4.2). L'operation « Sources », initiee par le Pr Louis De Backer (DeL), consiste a evaluer la qualite des eaux souterraines belges par un test colorimetrique realise par les enfants des eccles primaires. La fiabilite du test est controlee par un titrage des nitrates au laboratoire. 3

60

P ratique du calcul bayesien

r>.

~

Figur e 4.2 - Op eration Sources : un mod ele gamma pour la concent ration en nitrates.

•

Soit y , un n-echa nt illon iid issu d 'une population ga mma de par am etre de form e 0: > 0 et de param etre d 'echelle fJ > o. Posons e = (o:, fJ). La vraisemblance s'ecrit :

[yle]

gr n

fJD:

(0:) yf -l exp (- fJYi)

au y et g represent ant respectivement la moyenne arithmetique et geornetriqu e des observations. La log-vraisemb lan ce suit : L

(e ) == In [yle] = no: In fJ + n (0: - 1) ln g - n y fJ - n In r (0:)

Le calcul des derivees donne

a;~Ii) = n ln fJ + n ln g - mp (0:); = n o:fJ- l - ny ;

a;~ )

4. Calcul des lois a posteriori

61

ou rljJ (a) et rljJ' (a) sont respectivement les fonctions digamma et trigamma (disponibles dans R). Rappelons qu'elles sont definies comme :

(4.17) La matrice hessienne s' ecrit :

He == -n ( rljJ' (a) _(3-1 Le mode {} =

(&, S) est solution du systeme rljJ ( a) - In (3 { a(3-1 == fj

== In g

En substituant la seconde equation dans la premiere on obtient une equation en a qu'on peut facilement resoudre numeriquement dans le logiciel R : fj

rljJ (a) -Ina + In g

~

== 0 =? & ~ 1.2867 =? (3 ~ 0.0327

La matrice d'information empirique (4.12) arrive en substituant ces valeurs dans la matrice hessienne : P

3 e ~ 10

(

0,11 -2.87

-2.87) 112.83

La figure 4.3 illustre les resultats. A gauche, on a le mode de

0=

(a, i3) et

quelques isodensites (vue en plan). A droite, on montre la densite de probabilite de la probabilite 7f qu'une source non encore observee ait une concentration en nitrates depassant la norme Yo : 7f

4.2.4

== Pr (Y > yoln, fj, g) et Pr (0.22 S

7f

S 0.36) == 0.95

Estimation asymptotique des parametres d'une regression lineaire

Ce modele archiconnu ne pose aucun probleme de calcul. Nous l'avons d'ailleurs mentionnc des le chapitre 1 et renvoye le lecteur au chapitre 9, page 170, pour plus de details. Cependant, le traiter par voie asymptotique est un bon exercice de maniement. Dans sa version la plus simple, ce modele postule que la reponse reellc, Y, a un stimulus reel, x, a une distribution normale, localisee sur a + (3x et de precision T (Fig. 4.4).

62

Pratique du calcu l bayesien

0.045

14 12

0.04 10 0.D35

8

c::l.

6

0.03

4 0025 2 002 0.8

1.2

1.4

0 01

1.6

02

a.

0.3 0.4 Pr(y > Yo)

0.5

Figure 4.3 - Operation Sources : estimation asymptotique.

Un prior non informatif courant est [e] = [a , ,8, 7 ] ex 7 - 1 . Cependant, si on travaille avec le parametre In 7 defini sur JR., [e] ex cte, c'est-a-dire que In [ely, x] = In [yle, x] + cte. La log-vraisemblance s'ecrit : L (e) == In [yle, x] =

"2n In 7

n

"27 ,,", ~ (Yi

-

2

- a - ,8Xi )

+ cte

i= l

Les definitions suivantes apparaitront dans les developpe ments" .

~~X2

x2

n

t

~ (Yi - a - ,8xi)2 1 n (Xi - x)2

sce( a , ,8)

:; : L

s;

i=l

1

-n L (Xi - x) (Yi - Y)

S xy

n

i= l

Le calcul des derivees premi eres 0~(1I)

= 7~ (Yi - a - ,8xi )

o~~) = 7~ (Yi

0£(1I) oinT 4

=

- a - ,8Xi) Xi 2 sce

!l _ I.

2

see signifi e somme des carres des ecaris.

4. Calcul des lois a post eriori

63

Figure 4.4 - Le modele lineaire simple. ent raine le mode de la vraisemblance

o, =

( " 0:,

(3, T' ) =

( _

' _ Sx y Y - (3x ,- 2 ' Sx

n "

)

see(o:,(3)

La matrice hessienne -nTX -nTx 2

TE (Yi -

(3Xi ) Xi

0: -

ent raine Ia quantite d'information de Fisher

1(0) = - E (H o) = tir

(X = :2 o o

~)

1/2T

car E (see) = EE (Yi - 0: - (3Xi )2 = Evar (Yi) = n l r , La mat rice d'information empirique P suit :

P_

n2 - see(& , S)

(1

X

0

xx2 0

'Clog(u) Counter=Counter+1; 1ambda (i) =cand; else lambda(i)=lambda(i-1); end· end. ' Rat~=counter/4000

LAMBDA=lambda(100l:4000); % Predicti ve PRED=~oissrnd(LAMBDA)+l;

IC90= [prcti le(LAMBDA, 5), prctile(LAMBDA, 50),prctile(LAMBDA,95)]

Figure 5.2 - Algorithme MH pour les rangs de naissance.

Resultats Apres 4000 iterations dont 1000 pour la periode de chauffe, avec un taux d'acceptation de 0.45, on obtient les resultats suivants (tableau 5.3 , fig. 5.3). Inference bayesienne sous WinBUGS La relation (5.8) donne la contribution du rang de naissance rj a la vraisemblance. Si n; etudiants declarent le rang de naissance r, sous l'hypothese

92

Pratique du calcul bayesien

1090 A Taille

50 2.57 3

5 2.47 1

95 2.68 6

Tableau 5.3 - Rang de naissance : IC90.

dindependance, la contribution du rang r (qui arrive n r fois) s'ecrit :

a la vraisemblance

(5.12) Cette distribution n'est pas disponible dans WinBUGS mais l'astuce suivante permet de s'en sortir.

Le zero trick

Soit [yIO] la contribution de l'observation y a la vraisemblance pour un modele d'observable parametre par O. On sait que si une variable aleatoire x est distribuee selon une loi de Poisson de paramctre a > 0, la probabilite qu'elle prenne la valeur zero est exp (-a). Maintenant, si on identifie a a l'oppose du logarithme de la vraisemblance, on a :

a == -In [yIB] > a ::::} [x == ala] == exp (-a)

== exp (- (-In [yIB])) == [yIO]

(5.13)

Ainsi la contribution de l'observation y a la vraisemblance d'un echantillon issu d'un modele d'observable parametre par B, est identique a la contribution d'un zero a la vraisemblance d'un echantillon de zeros, issu d'une loi de Poisson parametree par a == -In [yIO]. Mais attention! Rien ne garantit que -In [yIO] > o. Aussi, la vraisemblance etant definie a une constante pres, on doit ajouter une constante 0 »> 1 a la log-vraisemblance de telle sorte que l'on soit certain que -In [yIO] + 0 > o.

Application du zero trick au probleme des rangs de naissance A partir de la relation (5.12), on a 00

In [riA, nr]

AZ

= -nrA + n; In ~ r (z + 2) 1{1,... ,z+l} (r)

II suffit donc de tirer des zeros dans une loi de Poisson de parametre TJr == -In [riA, n r ] + 0

5. Le cardinal sort du ra ng

93

3.5 K

3 2.5

2 L..-.............~~.L...-~~~'-'--~~~.LI...----'- -'--'-~u.J 10°

10

1

2

3

10

10

4

10

March e aleatoire

0.4 Q)

0.3

.~

tl

"D

-

0.2

'~

0...

0.1 0

I

-

2

I---

I--

3

4

5

Il---t 6

7

8

9

10

Nombre d'enfants

Figure 5.3 - Le problems des rangs de naissan ce. Profi l d 'une mar che aleatoire et distribution predictive a posteriori de la t aille de la fratrie type.

Remarque 5.1 Sous WinBUGS, la fonction step permet de coder facilement la condition 1{1 ,... ,z + l } (r)

1{1 ,... ,z + l } (r)

1¢?r:::; z + 1

= step (z + 1 - r) = { 0 ¢? r > z + 1

(5.14)

La fonctio n loggam calcule In r (z + 2). Par consequent, exp (loggam (z + 2))

= r (z + 2) = (z + 1) z!

(5.15)

Les figures (5.4) et(5 .5) mont rent respectivement le DAG et Ie code WinBUGS. Apres 4000 iterat ions dont 1000 pour la periode de chauffe, le tableau (5.4) donne un intervalle de credibilite a 90 % pour >. et T . On retrouve (evidemment) les memes resultat s que ceux obt enus sous R.

94

Pratique du calcul bayesien

A Ta ille

5 2.47 1

50 2.57

3

95 2.69 7

Tableau 5.4 - Rang de naissance : IC90 (WinB UGS) .

eta [r]

Figure 5.4 - Le problerne des rangs de naissance. Representat ion du modele hierarchique par un DAG sous WinBUGS.

Epilogue Inferer le cardinal type d'une collect ion d'ensembles ordonnes a partir de la seule connaissan ce du rang d'un de leurs elements est un probleme generique qui a des applicati ons pra tiques. Le modele du tramway voit le nombre de tramways circulant en ville comme un par ametre et c'est pourquoi on peut postul er un prior , en l'occurrence un prior non informatif de la forme [N] ex. N - 1 , car Nest vu comme un parame tre rl'echelle. Ce modele simple n'introduit pas de variable latente. Il exploit e direct ement tou te l'information disponible : conditionnellement a N , le ra ng du tramway observe est vu comme un tirage aleat oire dans une loi discrete uniforrne prenant ses valeurs dans n = {I , 2, ' " ,N}. On pourrait etre tente de s'en servir pour le probl eme des ra ngs de naissance : une fratri e = une ville et le rang de naissance de l'etudiant = Ie numero du tramway. Mais la generalisat ion de ce modele, pris tel quel, a plus d'une fratri e n'est pas simple. Le modele des rangs de naissance impliqu e des variables

5. Le cardinal sort du rang

95

latentes dans une structure hierarchique. La variable latente, Zj, represente Ie nombre de freres et soeurs de l'etudiant j. On a postule pour ces variables une distribution de Poisson de parametre A. Le rang de naissance est alors distribue uniformement sur l'ensemble f2j == {l,··· ,Zj + l}. L'inference sur A est realisee via un algorithme MH programme dans R ou dans WinBUGS via le zero trick. Dans les deux cas, chaque valeur de la chaine AIOOI, ... ,Aj,··· A4000 gencre un nombre de freres et soeurs, c'est-a-dire une valeur Zj, via un tirage aleatoire dans une loi de Poisson de parametre Aj. La predictive a posteriori est obtenue en ajoutant 1 a Zj puisque l'etudiant appartient a sa fratrie. A ce stade, nous avons decouvert et manipule des outils puissants pour resoudre des problemes de plus en plus interessante et utiles. Le chapitre suivant introduit la modelisation des evenements extremes via les modeles GEV et POT. Le defi est reel car, par definition, ces evenements sont rares alors que les enjeux sont importants. II y a donc peu de donnees et l'expertise est reduite. Pourtant, les fondements de ces modeles sont solides et leur utilisation rationnelle permet de mettre en place des protections qui fonctionnent. '# Ou rangde Ilaissance d"un .tudlant ala taille de sa fratrie (Sam,n Data Sets 119t

Utilisation du "zerotrick" Lavariable latente. 2, representele nombre defreres et sceurs (hors I'etudiant) ~ Elleasttiree dans une loidePois$onparametree parlambda :> 0 L'observable asile rang denaissance, r, deretudiant ~ II esttiredans une loidiscrete unifarme. denoie sur1,2, z+1 Enpredictif, N:: zet lataille delafratne, T,estdone egale aN+1 Notana que t ::: Z + 1

modet ( lambda .... dgamma(a,b) for(rin 1: m ) ( zero(r] 0, telles que le maximum normalise, Zn, converge en loi vers une v. a. r. Z lorsque n tend vers l'infini

La theorie donne une reponsc affirmative a cette question et precise la distribution de Z. Le comportement asymptotique de la loi du maximum M n depend de la fonction de repartition initiale F. (Fisher et Tippett, 1928) ont etabli qu'il n'y a que trois types de lois limites possibles: Frechet, Weibull 4 et Gumbel. La majorite des lois de probabilite usuelles appartiennent a l'un des trois domaines dattraction''. Par exemple, les distributions gamma et log-normale appartiennent au domaine d'attraction de Gumbel regroupant la majorite des distributions a queue fine; les distributions de Pareto, log-gamma et de Student appartiennent au domaine d'attraction de Frechet regroupant la majorite des distributions a queue lour de ; la distribution uniforme appartient au domaine d'attraction de Weibull regroupant la majorite des distributions sans queue. On suppose que la limite existe. pas confondre avec la loi de Weibull utilisee dans Ie domaine de la fiabilite. 5 On appelle domaine d'attraction d'une loi H l'ensemble des lois F pour lesquelles Ie maximum d'un echantillon, M n , converge en loi vers la loi des extremes du type H. 3

4

A ne

100

Pratique du calcul bayesien

En fait, on peut caracteriser ces trois types de distribution par une distribution unique, la loi qeneralieee des valeurs extremes ou modele GEV (generalized extremes values) (Gnedenko, 1943), (Jenkinson, 1955). Le modele GEV est coherent avec lc modele POT (peak over threshold) qui voit les valeurs extremes d'une observable comme les depassements d'un seuil fixe assez haute Ces depassements constituent un processus de Poisson marque, les excedents etant distribues selon une loi de Pareto qeneralisec qui n'est rien d'autre que l'oppose du logarithme du modele GEV. Ainsi, les modeles GEV et POT sont en quelque sorte les deux faces d'une meme medaille. Ils sont d'application dans les situations OU il est raisonnable de postuler que les evenements extremes sont independants. Dans le cas contraire, des modeles plus sophistiques existent (Drees, 2008). Les modeles GEV et POT sont caracterises par un parametre tridimenConduire une inference baycsienne sur implique de recourir aux sionnel methodes speciales du chapitre 4. Pour le modele GEV, aucune des trois conditionnelles completes n'est standard, mais un algorithme de Metropolis-Hastings sequentiel est relativement facile a regler. Pour le modele POT, deux des trois conditionnelles completes sont standards et l'utilisation d'une grille pour la troisieme permet de programmer facilement un echantillonnage de Gibbs. Le lecteur interesse trouvera dans (Coles, 2001) un excellent ouvrage d'introduction a la modelisation statistique des valeurs extremes traitee essentiellement sous le paradigme classique (Coles donne un exemple d'inference bayesienne dans la premiere section de son dernier chapitre).

e.

e

Note 6.1 Le statisticien bayesien raisonne toujours conditionnellement aux parametres. Cependant, pour allegcr les ecritures, il arrivera que le conditionnement soit implicite, notamment dans les developpernents.

6.2

Le modele GEV

Soit {X t } un processus stochastique a temps discret". Soit Xl, ... ,Xn une serie de n v. a. r. iid de fonction de repartition F. On peut ordonner cet echantillon par ordre croissant: X(l) < X(2) < ... < X(n). Intuitivement, on comprend que le maximum Mn == X(n) est une valeur extreme si nest assez grand. La probabilite que ce maximum soit inferieur a une valeur z don nee est triviale Pr (Mn < z) == (F (z))n Lorsque n tend vers l'infini, cette distribution est nulle en tout point z < z., ou z., est la plus petite valeur de la v. a. r. M n pour laquelle F == 1. On dit d'une telle distribution qu'elle est degeneree. L'idee est d'appliquer une 6 Sous Ie nom de processus stochastique it temps discret, on entend un modele permettant de decrire un phenornene aleatoire evoluant au cours du temps, OU les observations sont realisees en des instants t ETC Z.

6. Les modeles GEV et POT

101

transformation Iineaire au maximum M n telle que, lorsque n tend vers l'infini, la distribution limite, G, soit non degeneree, Les deux theoremes suivants fondent la theorie des valeurs extremes.

Theorems 6.1 (Fisher et Tippett, 1928). S'il existe des suites normalisantes {an} et {b n > O} telles que Pr (Zn = M nb- an n


O a

(6.2)

avec

(6.3) Remarque 6.1 La difficulte posee par la determination des coefficients an et bn > 0 n'est qu'apparente car

entraine

Pr (M n ~ an + bnz) ~ G (an

+ bnz) == G* (z)

ou G* est un autre membre de la famille GEV. Par consequent - comme on doit conduire une inference bayesienne sur les parametres pour identifier le membre de la famille GEV en adequation avec les donnees et l'expertise disponibles en pratique on ne se preoccupe pas de ces coefficients et il est licite d'ecrire, a n grand fixe :

102

Pratique du calcul bayesien

La loi generalisee des valeurs extremes postule que le maximum normalise Zn converge en loi vers la v. a. r. Z de fonction de repartition G (eq, 6.1) lorsque n tend vers l'infini. La v. a. r. Zest donc bien une valeur extreme. Le signe du parametre de forme ~ (prononcer xi) est capital. - Si ~ > 0, la loi de la valeur extreme Z est un membre de la famine des lois de Frechei (lois a queue lourde). - Si ~ < 0, la loi de la valeur extreme Z est un membre de la famine des lois de Weibull (lois bornees superieurement, donc sans queue). - Le cas ~ == 0 doit etre interprets comme la limite du modele (eq. 6.1) lorsque ~ ---+ 0, ce qui conduit a la famine des lois de Gumbel. Proposition 6.1 La limite du modele (eq. 6.1) lorsque ~ ---+ 0 conduit a la famille des lois de Gumbel definies sur ffi. par la fonction de repartition suivante G (ZIJLl a) = exp (- exp ( _ z:

A partir de [eq.

(6.5)

6.1), en raison de la condition (eq. 6.2) on a

M)-l/e

Z

(

JL) )

1+~-a-

(1 (

Z

M))

=exp -~ln l+~-a-

Le passage it la limite conduit it une indetermination (0/0) levee en appliquant la reqle de l'Hospital

. 1 (

M) = = -M -

Z lim - In 1 + ~-a

e~o~

Z -

a

Par consequent lim

e~o

M) -lie

1 + ~-a Z -

(

== exp

(z - M) --a

ce qui enirainc le resuliai (eq. 6.5). En pratique, le statisticien bayesien pose le modele (6.1) et c'est la distribution a posteriori de parametre ~ qui lui revele le domaine d'attraction de l'observable. La fonction de densite de probabilite du modele GEV s'obtient en differenciant (6.1) par rapport a Z :

(6.6)

6. Les modeles GEV et POT

6.2.1

103

La valeur de projet

La modelisation des valeurs extremes est du plus grand interet pour les sciences appliquees, notamment pour dimensionner les ouvrages de protection (digues, reseaux devacuat.ion des eaux de ruissellement, barrieres antiavalanche, etc.). En general, les dommages seront une fonction croissante de la difference positive entre I'intensite de I'evenement redoute et le niveau de protection. On appelle valeur de projet la valeur zp qui ala probabilite p d'etrc depassee

p == Pr (Z

> zplB)

(6.7)

La quantite T == »:' definit la periode de retour de l'evenement Z > zp. Elle est nommee ainsi car elle represente l'intervalle de temps moyen, par exemple en annees calendaires, separant deux occurrences successives de cet evenement. Ainsi, un evenement de periode de retour de T annees a la probabilite p == T- 1 de survenir chaque annee, Posons Xp

= -In (1 -

p)

~~

si p est petit

(6.8)

En general, la probabilite pest fixce par le decideur qui veut, par exemple, se proteger contre une crue qui revient tous les 100 ans, c'est-a-dirc qui a la probabilite p == 0.01 de se produire chaque annee, On deduit la valeur de projet zp associee a p en distinguant le cas OU ~ i=- 0 du cas OU ~ == o. Apres quelques manipulations elernentaires, on trouve :

a ==> zp = fJ, ~ ~ (1 - X;~)

o =?

zp

== J-L -

a

In x p

(6.9)

(6.10)

Dans un repere cartesien, les couples (zp, -In x p) dessinent une droite si ~ < 0 (Weibull) ou concave si ~ > 0 (Frechet) 7 . On peut en effet montrer que Ie ratio

~

== 0 (Gumbel), une courbe convexe si ZlO-3 -

zlO-2

zlO-2 -

ZlO-l

c'est-a-dire, Ie rapport de l'accroissement des quantiles du centenal au millenal sur l'accroissement des quantiles du decennal au centenal, est plus grand que 1 (comportement explosif si ~ > 0) tandis que les accroissements relatifs entre chaque ordre de grandeur de la periode de retour decroissent (atteinte d'une borne sup si ~ < 0). 7 Si l'axe des abscisses est en coordonnee logarithmique, on arrive aux memes conclusions avec les couples (zp, xp).

104

Pratique du calcu l bayesien

C>

'..L (1

+ ~z ~ c) -liE]

(6.21)

Du modele POT au modele GEV

La ressemblance des modeles (eq. 6.1) et (6.21) est frappante. En fait, la fonction de repartition de l'intensite maximale sur une fenetre unitaire (L == 1)

6. Les modeles GEV et POT

109

est eiroitement rcliee a la fonction de repartition du maximum des valeurs elementaires sur cette meme periodc. En posant L == 1, la relation (eq. 6.21) s'ecrit :

Pr (Z

< z) =

~ ~ c) -1/ E]

exp [_ A ( 1 + z

(6.22)

L'experience montre qu'un reparametrage des deux modeles facilite les demonstrations et l'ecriture des programmes informatiques. Plus important encore, un tel reparametrage permet de simplifier l'echantillonnage de Gibbs dans le cas du modele POT (Parent et Bernier, 2003). Pour bien distinguer les developpements, nous affecterons les parametres du modele GEV de l'indice o.

Po == a-I> 0,

f30 == -Po~o,

P ==

1]-1

> 0, f3 ==

-P~

(6.23)

Le seuil c etant convenablement fixe, les modeles POT (eq. 6.22) et GEV (eq. 6.1) deviennent respectivement

< zl,8, A, p) =

POT:

Pr (Z

GEV:

Pr (Z ::::; zl,8o, u, Po)

exp

[-A (1 - ,8 (z - c))P/,6]

= exp

[- (1 - ,80 (z - J-t))p0/,6o]

(6.24) (6.25)

La similitude des deux modeles est evidente, Avec ce reparametrage, les familles des lois de Frechet, Gumbel et Weibull correspondent respectivement a (3 < 0, (3 == 0 et (3 > O. Remarque 6.4 La loi du maximum selon le modele POT differe de la loi du maximum selon le modele GEV. - La variable GEV a une limite inferieure necessaire pour que Pr (Z ~ z) soit definie quand (3 < 0 (Frechet) : Z

> Zmin == J-l + /30- 1

Cette limite inferieure tend vers -00 dans le cas Gumbel. Au-dela de cette limite technique, la v. a. r. Z peut prendre n'importe quelle valeur superieurc. - Dans le cas du modele POT, la loi du maximum est une distribution censuree dans le sens OU elle depend d'un seuil c. Au-dela de ce seuil, les observations sont marquees (depassetnents), en deca de ce seuil, les observations n'interviennent que par le processus de Poisson Pr(X ~ ciA, L == 1) == exp (-A)

110

Pratique du calcul bayesien

Cependant pour les grandes valeurs de Z, au-dela de seuils c realistes, les deux modelcs devront donner des calculs de Pr(Z > c) tres voisins pour autant que les observations, differcntes dans chaque cas, et la validite des hypotheses le permettent. II ne faut pas oublier qu'en fonction de l'information disponible les estimations des parametres f3 et P de la distribution de Pareto generalisee peuvent varier selon le seuil. Cependant, leur homogeneite theorique est essentielle comme on I' a vu dans la discussion sur le choix du seuil. Par consequent, on peut aussi obtenir la valeur de projet a partir d'un modele POT. A partir du modele POT (eq. 6.24) et de la definition d'un quantile d'ordre 1 - p on obtient successivement :

1 - p = Pr (Z

< zp1(3, .A, p) =

exp [-

>. (1 - (3 (zp - c))P/ 13]

1( (1

zp=c+j3 1- ->:In(l-p)

)(3/P)

(6.26)

Avec le parametrage initial, compte tenu de la relation (eq. 6.8), on a aussi

(6.27) On comparera ce resultat avec la relation (eq. 6.9) rappelee ei-dessous

6.5

Inference bayesienne sur les parametres d 'un modele G EV

Le processus stochastique a temps discret est divise en k blocs generant une scrie de maxima Zl,'" ,Zk. Puisque les populations sous-jacentes sont independantes (hypothese iid), ces maxima le sont aussi et, pourvu que la taille des blocs soit assez grande, on peut considerer qu'ils sont identiquement distribues selon le modele GEV.

6.5.1

La distribution conjointe a posteriori

Le modele GEV (eq, 6.25) est done caracterise par Ie parametre () == (!3o, f-L, Po) et Ia densite de probabilite eorrespondante s'ecrit

[zIB] == Po (1 - f30 (z - f-L) ou Po E

lRt, !3o

E

v:

lR o, f-L E lR et f30 (z - f-L) < 1.

(30-

1

G (zIB)

(6.28)

6. Les modeles GEV et POT

111

L'hypothese iid entraine la vraisemblance d'un k echantillon de maxima:

[Zl,· .. ,zkIB]

=

P~

IT {[1 - f30 (Zi - JL)]pol/10-1 G (ziI B) } k

i=l

Pour le prior, on postulera l'independance des composantes du vecteur 0

et un prior non informatif simple a la forme suivante : 1 [fL] [Po] ex Po Pour le construire, nous avons pris 130 et u uniformes sur un domaine assez grand. Pour le parametrc d'echelle Po > 0, le prior habituel est une distribution gamma dont les parametres tendent vers zero

[0]

== [130]

[pola, b] ex

pg- 1 exp (-bpo)

---t

1 Po

-

a,b---+O

L'application de la regle de Bayes nous donne le posterior non normalise

[Blz1, ... ,Zk] ex p~-l

IT {[1 - f30 (Zi - JL )]pol k

/10-

1

G (Zi IB) }

(6.29)

i=l

La normalisation par calcul integral n'est pas possible et aucune conditionnelle n'est standard. L'inference peut se faire via un algorithme de MetropolisHastings.

6.5.2

Algorithme MH sequentiel applique au modele GEV

II sera commode de poser

f (f3o, JL, Po)

=

p~-l

II {[1 - f30 (Zi k

JL)]p01/10-1 G (ziIB)}

i=l

Puisque Po > 0, l'algorithme MH est plus facile changement de parametre suivant

a mettre en oeuvre avec le

¢ == In Po {:} Po == e ¢ La transformation logarithmique donne

In f (f3o, JL, Po)

= (k - 1) ¢ + (;: - 1)

t

In [1 - f30 (Zi - JL)]

k

- L [1- f30 (Zi - JLW'"I/1o ou

i=l

130 (Zi - fL)

< 1;

i == 1, . . . ,k

112

Pratique du calcul bayesien

L'algorithme Soit une marche aleatoire realisee dans JR3

(138,JLo,¢0) .

a partir d'un point

initial

()o ==

Pour loi instrumentale, nous avons choisi le produit de trois densites normales unidimensionnelles independantes :

130

r-v

dnorm (13~-l,vf3o);

JL*

r-;»

dnorm (JLi-1,vM)

;

¢* r-v dnorm (¢i-l,V 0 permet de creer un lien entre l'incertitude de l'expert et l'ecart-type de la loi beta.

7. Construire Ie prior Si on voit que E (IT) (1 - E (IT)) ==

rs

(r + s)

135

2

la seconde equation devient kc 2

==

m(l-m) r

+s+1

{=}

r

+ s ==

m(l-m) kc 2

- 1>0

(7.5)

On a donc obtenu la somme r + s a partir de la connaissance de m et de c (expertise) et en fixant une valeur k (on commence par exemple avec k == 1). Bien entendu, cette somme doit etre strictement positive, c'est-a-dire que les valeurs m, c et k doivent respecter l'inegalite 7.5. Sous cette condition, les hyperparametres recherches sont

r==m(r+s),

s==(l-m)(r+s)

(7.6)

Pour terminer, on presente la densite beta ainsi obtenue a l'expert en lui demandant si elle reflete bien son savoir. S'il n'est pas tres satisfait, il faut d'abord jouer avec k et, en cas dechecs repetes, recommencer les loteries.

7.3.2

L'expert donne deux quantiles de

Par construction, on a obtenu les quantiles IT q et IT p

7f :

Avec les notations de R, pbeta est la fonction de repartition de la distribution beta. Les hyperparametres recherches sont les solutions du systeme suivant

p - pbeta (ITp , r, s) == 0 { q - pbeta (ITp , r, s) == 0

(7.7)

II faut disposer d'un solveur numerique,

Exemple 7.3 Une machine de production est en cours de reglage. Le parametre IT est la probabilite qu'une piece choisie au hasard soit conforme au cahier des charges. Selon l'operateur, il y a 95 chances sur 100 que IT excede 0.5 et 10 chances sur 100 qu'il excede 0.9. • La resolution numerique du systeme 7.7 avec p == 0.9, q == 0.05, ITp == 0.9 et ~ 7.51 et s ~ 2.62. Encore une fois, on presente ce resultat a l'expert et on remet l'ouvrage sur le metier si necessaire,

IT q == 0.5 donne r

136

Pratique du calcul bayesien

7.4

Caler un prior conjugue sur deux quantiles elicites des parametres d'un modele d'observable normal

Note 7.1 Un prior depend d'hyperparametres. Par exemple, le prior pour le parametre de localisation fL d'un modele d'observable normal, X N (fL, T), peut dependre de la precision T et s'ecrire fLIT N (m, kT). En toute rigueur, on devrait mettre les hyperparametres dans le conditionnement et done ecrirc fLIT, m, k N (m, kT). Pour allegcr les ecritures on ne le fera pas (on ne conditionne pas sur les hyperparametres). t"'V

t"'V

t"'V

7.4.1

Dialogue avec l'expert

Nous considererons ici le cas d 'un parametre () reel, interpretable dans le cadre d'un modele donne. II peut s'agir, par exemple, du parametre d'un modele d'observable exponentiel (dim () == 1) ou d'un modele d'observable normal (dim () == 2). II est evident que la difficulte augmente avec la dimension de (). En effet, prenons l'exemple d'un parametrc tridimensionnel : () == (()l, ()2, ()3). Certes on peut toujours decomposer la distribution conjointe comme suit :

(7.8) On comprend que pour l'expert, il soit plus aise de donner un avis sur un quantile marginal, par exemple la medians de (J2, que sur un quantile conditionnel, par exemple la mediane de ()2 quand il dispose de l'information ()3' Nous reviendrons bientot sur cette difficulte en illustrant la procedure avec le modele d'observable normal. En general, l'expert n'est pas convie a proceder a une introspection detaillec pour elicitor toutes les caracteristiques d'un prior. Comme indique ci-dessous, il est en effet beaucoup plus courant de lc limiter a fournir quelques valeurs typiques : mediane (J50, quartile (()75 ou (J25), decile (B go ou (JIO), etc. Ces caracteristiques peuvent suffire a caler des distributions de probabilite de forme analytique connue a un nombre de parametres indetermines pres si ce nombre est egal au nombre de caracteristiques elicitees, C'est la methode dite des quan-

tiles. Remarque 7.3 Les parametres du prior sont souvent appeles hyperparametres pour les distinguer des parametres du modele d'observable.

7.4.2

Le parametre

a eliciter est

unidimensionnel

C'est, par exemple, le parametre de localisation du modele normal de variance unitaire : E (Y) == B. Les premieres questions a poser a l'expert doivent concerner le support de B, c'est-a-dire l'etendue de l'intervalle [()min, ()sup]. Bien souvent l'expert sera

7. Construire le prior

137

dans l'incapacite d'evaluer precisement ces limites, auquel cas il est preferable d'utiliser des distributions a priori dont les bornes sont mathematiquement infinies et de lui soumettre la tache d'eliciter des quantiles, grandeurs statistiques plus aisernent interpretables. C'est le cas notamment de la mediane (}50 de (). Si meme l'intervalle [(}min, (}sup] est indeterrnine, l'expert peut etre capable de repondre a la question suivante : Quelle est pour vous la valeur M telle que Pr(() < M) == Pr((} 2:: M) ? La valeur M qu'il donne est la mediane (}50. Ensuite, on peut lui poser la question suivante : Quelle est maintenant, selon uous, la valeur Q de () telle que Pr(M ~ () ~ Q) == Pr((} 2:: Q) ? Puisque M est la mediane, Q est necessairement le troisicme quartile, c'est-a-dire (}75 == Q. En poursuivant ces questions sur des segmentations d'intervalles en probabilites egales on peut atteindre toute proximite d'un quantile (}p quelconque. Certaines de ces questions peuvent etre un controle de coherence. Ainsi apres une premiere elicitation du troisieme quartile (}75, l'expert peut etre amene a repondre a : Quelle est la valeur Q telle que Pr( Q ~ () ~ (}75) == Pr(() ~ Q) ? Si Q est differente de la mediane M trouvee precedemment, alors l'expert doit etre confronte avec cette incoherence et doit la resoudre. Si la notion de quantile devient plus precise dans l'esprit de l'expert, on peut lui demander de repondrc a des questions plus elaborees concernant des fonctions simples, comme des ecarts ou des rapports de quantiles : - Quelle est la valeur la plus probable de X90 - X50 d'une grandeur oleatoire X eiudiee ? - Quelle est la valeur la plus probable de X90 / X50 ?

!

Remarque 7.4 Si l'expert est plus familier des distributions de probabilites et de leurs caracteristiques statistiques, on peut lui demander d'exercer son introspection directe pour eliciter des statistiques plus synthetiques comme des esperanees mathematiques ou des variances de grandeurs d'interet. Nous connaissons des specialistcs de la geophysique, tres exerces en analyse des donnees de leur domaine de recherche, pour qui des statistiques comme un coefficient de variation (ecart-type exprime en unite de moyenne) ont des significations physiques parfaitement quantifiables a priori. Bref, il est possible d'obtenir de l'expert quelques valeurs typiques de (), souvent la mediane M et un quantile (}p qu'il juge avoir une probabilite de depassement egale a 1 - p. Le travail de l'expert s'arrete la et le statisticien peut alors lui soumettre un modele de distribution a deux parametres (loi normale ou loi gamma par exemple) qui pourra etre cale sur les deux informations fournies.

Calage d'un prior normal indexelicitationsuelicitation !d'un prior beta@d'un prior beta L'expert ayant fourni la mediane (}50 et un percentile (}p, le calage d'une loi normale ncccssite de determiner sa moyenne, Me, et son ecart-type, O"e > o. Cette tache est particulierement facile car

138

Pratique du calcul bayesien

!-Le

(7.9)

== B50

B

B

-

50 ae == -p - -

(7.10)

zp

OU zp designs Ie p-ieme quantile de la loi normale standard (p. ex. p == 90 zp ~ 1.28).

=}

Calage d'un prior gamma

indexelicitationsuelicitation !d'un prior gamma Note 7.2 Nous avons pris l'habitude d'ecrire la densite de probabilite gamma comme suit

[Ola, (3]

=

~:) 0

0

-

1

exp (-(30)

Ainsi, le parametre d'echelle, {3 > 0, s'exprime dans les unites inverses de la variable aleatoire B. C'est pourquoi, nous l'appelons souvent paramctre d'echelle inverse.

L'expert ayant fournit deux quantiles, Bp et Bq , le calage d'un prior gamma necessite de determiner deux parametres, a > 0 et {3 > o. Cette operation implique de connaitre le theoreme suivant.

¢/{3

Theoreme 7.1 Si ¢ rv dgamma (a, 1) alors B {3 > 0 est le parametre d' eclielle inverse.

rv

dgamma (a, {3) OU

Corollaire 7.1 Le quantile de B correspondant au quantile de ¢ est Bp == {3-1¢p OU ¢p == qgamma (p, a, 1)2. Il s 'en suit que le rapport des quantiles Bp / Bm est uulepetuiami de {3.

Le rapport Bp/B q etant connu, l'equation en a qgamma(p, a, 1) _ Bp qgamma(q, a, 1) Bq

== 0

(7.11)

est resolue graphiquement ou par un solveur numerique (on remarque que a ne depend pas des unites). Sachant la solution &, le parametre {3 suit par ~ = qgamma(p,&,

Bp

2

1) = qgamma(q,&, 1)

La fonction qgamma est disponible dans R.

Bq

(7.12)

7. Const ruire Ie prior

139

Remarque 7.5 Dans Ie cas particulier ou l'expert a donne Ie mode de 0, soit et un quantile d'ordre p, soit Op , on doit resoudre Ie systeme suivant :

Om ,

0m -{ (3 =

a -I j3

qgam~: (p,a , l)

On commence par resoudre l'equation en 0: Op

x (0: - 1) -

Om X

et pn termine par

qgamma (p , 0:, 1)

= 0 ~ 0:

0: - 1 (3 = A

Om

Exemple 7.4 Pour l'expert , la duree de vie mediane d'un compose electro• nique vaut 15 unites de temps et Ie nonanti erne percentile en vaut 25.

A partir des relation s 7.9 et 7.10, un prior normal sera localise sur /-lo = 15 avec un ecart-ty pe (TO ~ 7.80. Apart ir des relations 7.11 et 7.12 et du graphique (fig. 7.1) un prior gamma aura les param etr es suivants : 0: ~ 5.55 et (3 ~ 0.35. (0) ~ 6.8. L'esperance et l'ecart-type de 0 suivent : E (0) ~ 15.9 et

/v

0.7,-------,-

----,-

----,--

----,--

---,--

-,--

-,--

--.--

--.-----,

06 0.5

0.4

~

0.3

: : 0.2

0.1

-0. 1

Fig ure 7.1 - Det ermination graphique du par am etre de form e d 'un prior gamma .

7.4.3

Le parametre

a eliciter est bidimensionnel

indexelicitationsjelicitation !d' un prior beta@d'un prior beta

140

Pratique du calcul bayesien

Soit l'observable Y supposec distribuee selon une loi normale N(jJ;, T). On l'a deja dit, le prior conjoint peut toujours s'ecrire comme le produit d'une distribution conditionnelle par une distribution marginale : (7.13) L'elicitation d'un quantile d'une distribution conditionnelle comme [jJ;IT] est beaucoup moins aisee que l'elicitation d'un quantile d'une distribution marginale comme [jJ;]. Notons que l'expert peut n'avoir aucune raison de lier jJ; a T soit parce qu'il sait que ces deux parametres sont independants (jJ; 1- T), soit parce que son savoir est tellement reduit qu'il ne saurait defendre un lien et donc, par defaut, il postule leur independance : (7.14) - Le calage d'un modele gamma a partir de la mediane de T et d'un quantile d'ordre p signifiant pour l'expert (par exemple Q7,0.90) se fait selon la methode decrite ci-dessus. - Pour u, l'elicitation de la mediane M~ et d'un quantile Q~,p se fait sans reference a T. Ensuite, le calage d'une loi standard depend du lien entre jJ; et T. Si l'hypothese dindependance est retenue, le calage d'une loi normale sur jJ; est chose aisee. Dans le cas contraire, il s' agit de caler une loi de Student sur u selon une procedure un peu plus subtile.

Cas oft les deux parametres sont independants Exemple 7.5 En septembre 2000, une equipe de l'Institut national de la recherche agronomique (INRA) mesura la taille des juveniles des saumons sauvages sur le Scorff, une riviere de Bretagne. Des peches electriques ont permis de prelever des echantillons le long de la riviere en 38 sites regulierernent espaces, Une question concerne le differentiel de croissance des juveniles localises en 3 sites sur 38 (en amont, en aval et a proximite) d'une pisciculture industrielle. Les effluents de la pisciculture ont-ils une influence sur la taille des saumons sauvages? Si oui, de quel signe? • Pour le parametre u, lc chercheur a l'INRA s'est appuye sur les observations des 35 sites restants. Sur ces 35 jeux de donnees, il a calcule 35 moyennes empiriques. A la vue de leur histogramme, il a propose un prior normal, centre sur m == 100 cm avec un facteur d'echelle de 8 == 10 cm. Quant a la precision, les statistiques de dispersion empiriques ont conduit a une loi gamma de parametres a == 3.4 et b == 250 => E (T) ~ 10- 2 , V (T) ~ 5.44 X 10- 5 . Le prior conjoint suit :

1 (1

[jJ;, T] == - - exp ~8

--(jJ; - m) 2 28

2) x -bT a

r(a)

a-I

exp (-bT)

(7.15)

7. Construire Ie prior

141

Imaginons que l'expert ait donne Qp"O.90 == 115 cm en lieu et place de s. Les proprietes de la loi normale permettent imrnediatement de trouver la valeur s correspondante (eq. 7.10). On trouve s ~ 11.7 cm. L'hypothese d'independance entre J-L et T, effectivement commode pour l'elicitation, implique que la distribution conjointe a priori n'est pas un conjugue naturel du modele normal.

Cas oil les deux parametres sont dependants On se place dans la situation OU J-L est lie a T (modele 7.13). Le prior (conjoint) conjugue du modele normal est explicite en fonction d'hyperparametres a, b, m, k, c'est-a-dire (en utilisant les notations de R) :

T

rv

dgamma(a, b)

(7.16)

J-LIT

rv

dnorm(m, kT)

(7.17)

Pour le parametre T, il suffit de repeter la procedure suivie pour le calage d'une distribution gamma (fig. 7.1). Pour le parametre J-L, on l'a deja dit, un expert prefere parier sur des quantiles signifiants pour lui, et ceux associes a une distribution conditionnelle ne lui disent en general pas grand-chose.

Theoreme 7.2 A partir des relations 7.13, 7.16 et 7.17, la distribution marginale de J-L est une loi de Student, a v == 2a deqres de liberte, localisec sur m et de pararnetre d'echelle O"p, == Jb/ak. Corollaire 7.2 La variable oleaioire t == (J-L - m) /0" est distribuee selon une loi de Student standard,

a v == 2a

deqres de liberu:

Compte tenu de la symetrie de la loi de Student, on obtient Ie systeme suivant :

Mp, == m QI",P = m +

(7.18)

V(-;;:kTx .tmv(p, 2a)

(7.19)

OU tinv (p, 2a) donne le quantile d'ordre p d'une loi de Student standard a 2a degres de liberte. En resolvant ce systeme par rapport a m et k, on trouve

k ==

~ (tinv(p, 2a)) 2 a

Qp"p -

Mp,

(7.20)

Conditionnellement a la connaissance de a et b, il suffit donc que l'expert donne m.; et un quantile Qp"p pour calculer k. Or, dans l'exemple des saumons

142

Pratique du calcu l bayesien

sauvages du Scarff, l'expert a donne m Jl = 100 em, s Jl = 10 em et Ie calage d'un modele gamm a sur la median e et un quantile de T a donn e a = 3.4 et b = 250. Si la distribution mar ginale de fl etait normal e, on sait qu 'un ecart-ty pe au-d ela de la moyenne correspond pr atiqu ement au qua tr e-vingtqua tri eme percentile : m + S ~ QJl ,O.84 ' Mais la distribution marginale de fl est de Student , distribution qui est plus etalee que la distribution normal e. En consequence, m + S corres pond a un quanti le moindre (0.50 < p < 0.84). On est en train de caler une distribu tion sur les connaissances de I'expert et on peut decider que m + S corres pond au troi sieme quarti le de la loi de Student : QJl ,O.75 = m + s. Des lars 250 ( tinV(0.75,6.8)) 2 x ~0 .37 k= 3.4 10 La figure 7.2 montre Ie prior du chercheur de I'INRA en tro is dimensions.

8

o

6

Figure 7.2 - Represent ation du pr ior de !'expert en 3D.

Epilogue Nous avons pose et surtout rapp ele un certain nombre de prin cipes et de precaut ions a prendr e pour conduire Ie necessaire dialogue expert-statisticien dans cette tac he commune d 'elicitation. II s'agit d'obtenir de l'expert des evaluat ions quantitatives permet t ant de parier sur les valeurs possibles des inconnues, les par ametr es du modele.

7. Construire Ie prior

143

Nous avons souligne l'importance des modeles d'elicitation et des priors modelises - sans nier qu'ils peuvent etre choisis aussi pour des raisons de commodite mathematique - en considerant directement le modele probabiliste des observables. Les modeles classiques ont ete inventories en les accompagnant de techniques d'elicitation adequates comme les methodes des quantiles. Cependant, cet inventaire ne clot pas la liste des methodes disponibles. Nous avons deja parle de l'utilite a cet egard des modeles hierarchiques que nous verrons plus en details en progressant dans la lecture de la seconde partie de cet ouvrage. Les modeles hierarchiques permettent l'introduction rationnelle d'informations complementaires objectives pour quantifier les hyperparametres des priors de premier niveau qui apparaissent alors comme des parametres de niveau superieur. On pourrait dire que cette hierarchisation repousse le probleme d'elicitation des priors ace niveau superieur. Neanmoins, la sensibilite des resultats de l'analyse statistique finale aux incertitudes des priors diminue alors considerablement. Pour conclure, la panoplie des methodes pratiques delicitation est maintenant assez large pour permettre une application complete de toute la chaine des raisonnements bayesiens et en garantir l'efficacite. Avec le prochain chapitre, nous entrons formellement dans la seconde partie de cet ouvrage en traitant un probleme reel d'halieutique.

Deuxieme partie

... a la souris

Chapitre 8

Modele de capture-recapture par assemblage de modules fonctionnels binomiaux : application au cas des saumons Prologue La modelisation statistique bayesienne revient a imaginer un modele probabiliste, susceptible de reproduire les observations (chap. 1), souvent pour fournir une aide a la decision en avenir incertain (chap. 2). Ce modele est avantageusement represente par un DAG (chap. 3). D'un point de vue operationnel, il faut eliciter le prior (chap. 7) et inferer ses parametres par application de la regle de Bayes (chap. 1). La determination des distributions a posteriori peut ncccssiter un recours aux methodes de Monte-Carlo (chap. 4). C'est notamment le cas pour les modeles realistes, lesquels impliquent souvent des variables latentes (chap. 3). Leur DAG montre une structure hierarchique et modulaire. Cette modularite confere une grande souplesse au modele comme le montre ce chapitre dedie a l'evaluation des stocks de saumons.

8.1

Introduction

Le modele d'evaluation des stocks de saumons presente ici ne repose que sur des equations de bilans, des tirages binomiaux et des priors sous forme de lois beta. Le DAG est une representation conceptuelle des differents evenernents

148

Pratique du calcul bayesien

qui peuvent se produire dans une population de saumons qui remontent la rimere Scorff, utilisee comme cas ri'etudc (Parent et Prevost, 2003). Ces donnees reelles proviennent d'un projet commun entre l'Institut de recherche agronomique (INRA), le Conseil superieur de la peche, et la Federation de peche et de protection des ccosystemes aquatiques du Morbihan". Les scientifiques et les gestionnaires de la rivierc ont besoin non seulement de l'estimation de la taille de la population de saumons (la valeur la plus probable), mais aussi de l'estimation de l'incertitude la concernant (Clobert et Pradel, 1993). Trois types de quantites incertaines apparaitront dans le DAG : les observables (notees Yindice) , les variables latentes (notees Xindice) - ou variables phenomenologiques auxiliaires non observees - et les parametres (dcsignes par des lettres grecques). Les lois a priori seront construites a dire d'expert (chap. 7). Enfin, nous realiserons l'inference par echantillonnage de Gibbs (chap. 4).

Remarque 8.1 Dans un souci pcdagogique, nous faisons ici une exception a notre parti pris de depart et nous utiliserons done une lettre latine majuscule, par exemple Y, pour designer une observable ou une variable latente, et la minuscule correspondante, soit y, pour designer une valeur particuliere. On ne peut evidemment pas respecter une telle convention pour les parametres (c'est ce qui ad'ailleurs justifie notre convention initiale).

8.2 8.2.1

Presentation du probleme Les trois dernieres etapes du cycle de vie du saumon : remonter la r ivlerc, echappcr aux pecheurs a la ligne et survivre jusqu'a la saison du frai

Les saumons atlantiques (SaZmo saZar), qui reviennent adultes dans les rivieres de Bretagne (France), sont repartis en deux categories : le saumon de printemps qui a passe deux annees en mer (exceptionnellement trois) et les castillons qui reviennent dans leur riviere natale l'annee qui suit leur migration vers la mer. Les castillons constituent l'essentiel des adultes (r-v 90 %) qui reviennent dans la riviere, principalement de la fin du printemps a la premiere moitie de l'ete, Sur la riviere Scorff, un dispositif experimental de controle des migrations a ete installe, et les adultes de retour sont denombres par la technique du marquage-recapture. Le rnarquage est opere dans un dispositif de piegeage situe a l'embouchure de la riviere. L'efficacite du piege varie selon le debit de la riviere, L'etude de cas presentee ici ne traite que du retour des castillons. La figure 8.1 decrit le sort d'un saumon rentrant dans sa riviere d'origine apres 1 La collecte des informations sur le terrain a ete effectuee par les techniciens de la station experimentale du Moulin des Princes, Nicolas Jeannot et Francois Burban, aides de Jean-Yves Moelo.

8. Modele de capture-recapture: application au cas des saumons

149

son voyage dans l'Atlantique. Trois evenernents principaux peuvent arriver au candid at repro ducteur. 1.

A l'entree dans le Scorff, le saumon peut etre capture, marque et relache, C'est la premiere etape de la procedure d'estimation du stock.

2. Ensuite, une certaine quantite d 'individus - marques ou non - sera prelevee par les pecheurs a la ligne. La loi francaise exige que la prise de saumon soit officiellement declares, mais cette obligation legale n'est pas toujours respectee . Une et ude locale supp lementaire permet de completer ces renseignements. Ces deux sources permettent d'obtenir une premiere evaluation du nombre de saumons « reellement captures», et un certain nombre de saumons preleves est apporte aux techniciens de l'INRA pour identification du marquage. 3. Enfin, le poisson qui a echappe a la peche a la ligne devra survivre jusqu'a la saison de reproduction. Pendant le frai hivernal, les chercheurs se rendent sur les sites de reproduction et completent les et udes statistiques par une phase de recapture.

Environnement naturel

Non Marque recapture Y6 Marque Recapture

Y,

Pieges et Marques

Y1 x",f

Xur

Libres

Marque Vu pour sur

Marque et peche

X

..••..•..•.••..... •.••.•

(Pecheurs

COIlUO/e des Non Marque

Non Mar que Vu pour sur ~qll es

Non marque et pechex"c

Declare r

Iarque Declare

Y

3

Y~

Y4

')

;:::::.:..~ ::::::

Figure 8.1 - Le destin d 'un saumon qui revient remonter Ie Scorff.

.

150

Pratique du calcul bayesien

8.2.2

Variables observees

Les donnees du tableau 8.1 concernent six variables (en colonnes) suivies pendant six annees consccutives. Les donnees de la premiere annee (1994) sont exclues de l'etude, car elles sont significativement differentes des autres. La procedure u'etait pas completement rodee et l'efficacite du piege et la recapture au moment du frai ont ete moins bonnes. Les variables observees portent les informations suivantes : - Y1 : nombre d'individus captures, marques et relaches : - Y2 : nombre de poissons marques, peches a la ligne et rapportes par les pecheurs pour la detection du marquage; - Y 3 : idem pour les poissons non marques; - Y4 : total des poissons provenant d'observations sur les sites de peche (Y4 > Y2 + Y3 ) ; - Ys : nombre de poissons marques, recaptures pendant ou apres Ie frai ; - Y 6 : idem pour les poissons non marques.

Annee 1994 1995 1996 1997 1998

Y1 156 500 502 320 442

1999

167

Y2 3 39 25 17 50 16

Y3 14 10 8 7 5

Y4 42 75 87 33 66

4

24

Ys 4 31 45 19 56 16

Y6 14 28 14 9 13 11

Tableau 8.1 - Donnees du Scorff.

Expertise a priori sur Ie comportement du saumon

8.2.3

Des parametres techniques, inconnus, mais supposes stationnaires Des caracteristiques stochastiques regissent le comportement individuel d'un saumon. Ces quantites sont censees rester identiques d'un poisson a l'autre. Les sept parametres techniques suivants, inconnus et incertains, sont conceptuellement essentiels pour les biologistes : 1. s.: nombre de castillons qui remontent le courant; 2. () : probabilite qu'un castillon soit capture et marque, au passage du piege ; 3.

probabilite qu'un saumon non peche survive jusqu'a la periode de reproduction; Q

:

4. (3: probabilite qu'un castillon soit preleve par les pecheurs ;

5.

T : probabilite qu'un saumon peche soit enregistre comme «prise certaine » ;

8. Modele de capture-recapture: application au cas des saumons

151

6. 6 : probabilite qu'un saumon peche et enregistre soit declare et Ie marquage verifie par les techniciens; 7.

probabilite qu'un castillon soit recapture apres la periode de reproduction.

1r :

Encoder l'expertise a priori Pour le cas du Scorff, la connaissance a priori (resumee par H dans le raisonnement conditionnel) peut etre synthetisee comme suit. 1. Etant donne la taille de la riviere, les donnees anterieures sur la production juvenile dans la riviere (Bagliniere et Champigneulle, 1986) et le nombre de survivants apres le sejour en mer (Potter et Crozier, 2000), les experts sont prets a parier a 9 contre 1 que le nombre de saumons rentrant dans le Scorff, n, se situe dans l'intervalle [100,3000] avec une valeur hautement probable autour de 700 individus. 2. On ne connait guere la probabilite de capture, (), au piege place pres de l'embouchure du Scorff mais on peut imaginer une repartition symetrique avec 0.5 comme moyenne et seulement 10 % de chances pour cette probabilite d'etre infcrieure a 0.1 ou superieure a 0.9. 3. La premiere estimation du taux de survie des saumons dans la riviere, a, est superieure a 0.9. Les experts sont pratiquement surs (avec une probabilite a priori de 0.9) que ex est superieur a 0.75. 4. Le taux d'exploitation de la peche a la ligne, {3, est sans doute situe autour de 0.1- 0.3. II semble peu credible (moins de 10 % de chance) que f3 depasse 0.7. 5. La probabilite, T, qu'un saumon attrape soit reconnu par les controles locaux comme prise certaine est superieure a 0.9 et il semble hautement improbable (5 %) qu'elle soit inferieure a 0.5. 6. On sait peu de choses sur la probabilite, 6, qu'un saumon reconnu soit prcsente au controle du marquage. Une repartition symetrique avec 0.5 comme moyenne et seulement 10 % de chances d'etre inferieure a 0.1 ou superieure a 0.9 traduirait cette meconnaissance (prior plutot vague). 7. En considerant le nombre de sites etudies et les efforts de survie durant la recapture, la probabilite de recapture, 1r, est tres vraisemblablement inferieure a 0.25, peu probablement comprise entre 0.25 et 0.5 et il est presque impossible qu'elle soit superieure a 0.5.

Remarque 8.2 Dans ce qui precede, tres vraisemblablement signifie qu'il y a neuf chances contre une, presque impossible represente moins de une chance sur cent. La probabilite restante (environ 9 %) quantifie le qualificatif improbable.

152

Pratique du calcul bayesien

Construction des lois a priori

a dire d'expert

La figure 8.2 montre une loi de probabilite discrete de forme acceptable pour representer l'expertise H sur le parametre «. Cette distribution a ete obtenue par une discretisation d'une distribution gamma avec un parametre de forme egale a 2.4 et un parametre d'echelle egale a 5002 (voir chap. 7, p. 138). Cette distribution est tronquee a l'intervalle [0,4000] en raison des ressources limitees de la rivicre. Tronquer au-dela de 4000 permet aussi un calcul d'integration plus commode, mais une analyse de sensibilite montre que c'est largement justifie. Cette distribution presente un mode aux environs de f\; ~ 700 et met 90 % de la masse de probabilite dans l'intervalle [100, 3000]. (8.1)

Les six autres parametres (), Q, {3, 7, b, 7r sont des probabilites, Leur prior est donc avantageusement represente par une distribution beta sur l'intervalle reel [0.1] (voir annexe B). Pour chacun d'entre eux, il faut donc fixer deux coefficients, an et bu, de sorte que cette distribution reflete bien l'expertise. La figure 8.3 et le tableau 8.2 montrent les resultats de l'elicitation de la loi de probabilite beta(aH,b H ) pour traduire l'expertise a propos des differents parametres techniques. Cette elicitation a ete conduite a partir de techniques presentees au chapitre 7 a partir des equations 7.3 a 7.5. Comme la connaissance a priori de chaque parametrc est etablie independamment, le prior conjoint est le produit de tous les priors.

b

Interpretation Efficacite du picge Taux de survie Taux de capture Suivi sur site Suivi techniciens

7r

Taux de recapture

() Q

(3 7

Expertise H ()0.05 == 0.1; ()0.95 == 0.9 M a ~ 0.95; QO.l == 0.75 M(3 ~ 0.2; (30.9 == 0.7 M T ~ 0.9; 70.05 ~ 0.5 b O.0 5 == 0.1; bO.95 == 0.9 M 1r ~ 0.2; 7r0.99 == 1/2 [0.25 < 7r < 0.5] == 0.09 7r0.9 == 1/4

aH 1.53 10 1.3 5.5 1.53

bH 1.53 1.5 2.2 1.5 1.53

1.6

11

Tableau 8.2 - L'expertise a priori H est encodee via des distributions beta.

Remarque 8.3 L'expertise sur 1f implique quatre conditions. La determination des parametres de la loi beta implique de resoudre un probleme d'optimisation sous contraintes. 2

E (~) == 2.4 x 500.

8. Modele de capture-recapture: application au cas des saumons

X

10

153

-4

7 ['--'-'-----,---,-----,----,--

-

-r--

-

,.---

---,-

-

--,

Mode=700

6

Intervalle de credibilre a 90 %

500

1000

1500

2000

2500

3000

3500 4000 Tattle du stock '0(

Figure 8.2 - Loi a priori pour la taille du stock , parametre K .

8.2.4

Les variables latentes decrivent le phenomene biologique

Les paramet res inconnus et les variables observees ne sont pas suffisants pour decrire les peregrinations d'un saumon. Des variables non observees, mais ayant une signification physiqu e, sont alors introduites. Elles sont utiles pour aider a comprendre les etapes int ermediaires de la modelisation condit ionnelle. Evid emment , le modele doit et re complete ment defini ce qui exige que les distribution s conditionnelles des variables latentes sacha nt les par ametres et les observabl es doivent et re precisees. Les vari ables lat entes suivant es presentent un interet particuli er pour la modelisation : - X u u == saumons non captures, par consequent non marques (indice uu pour unmarked, uncaptured) ; - X m c == individus marques peches a la ligne ; - X u c == individus non marques peches a la ligne (unmarked, captured) ; - X m j == individus marqu es rest es libres pendant la period e de peche (marked, free) ; - X uj == individus non marques testes libres pendant la period e de peche ; - X m r == individus marques enregistres comme reellement at tra pes (marked, registered) ;

154

Pratique du calcul bayesien

1.5

6

2

5

1.5

4

3 2

0.5 0

0

0.5

0

0

8

0.5

0

\, 0

0.5

~

Ct.

1.5

4

1( \

3

1 2

0

0

0.5 T

0.5

2

0

0

0

0.5

8

\ 0

0.5 11

Figure 8.3 - Loi a priori pour les parametres descriptifs de comportement.

- X u r == individus non marques enregistres comme reellement attrapes ; - X m s == casiillons marques survivants jusqu'au frai ; - X u s == casiillons non marques survivants jusqu'au frai.

Certaines combinaisons de variables latentes sont importantes pour etablir les comptes-rendus des scientifiques. A titre d'exemple, scientifiques et responsab les de la peche aimeraient connaitre le champ des valeurs credibles pour X m c + X u c , nombre total de saumons attrapes par les pecheurs a la ligne. D'autre part X m s + X u s , qui represente « T'echappement », apparait comme une valeur de pour connaitre la perennite de l'et at du stock.

Le model e st a t ist iq ue sous la forme d'u n graphe a cy cli que oriente Les equations completes du modele comprennent des equations deterrninistes de bilan et des equations stochastiques de comportement binomial. Elles s'ecrivent comme suit (notation R) :

8. Modele de capture-recapture: application au cas des saumons

Y1

~ dbinom(~,

155

0)

X uu == ~ - Y1 X mc ~ dbinom(Y1 , (3), X uc ~ dbinom(Xuu, (3) Xmj == Y1 - X mc, X uj == X uu - X uc Y4 ~ dbinom(Xuc + X mc, T), X mr ~ dbinom(Xmc, T) Y4 == X ur + X mr Y2 ~ dbinom(Xmr, 6), Y3

~

(8.2)

dbinom(Xur, 6)

X ms ~ dbinom(Xmj, a), X us ~ dbinom(Xuj, a) Y5 ~ dbinom(Xms, 1r), Y6 ~ dbinom(Xus, 1r) La figure 8.4 rcprescnte toutes les quantites par des noeuds (soit stochastiques, soit deterrninistes) sur un graphe oriente d'influence, OU les fleches penetrent dans un nccud depuis les variables qui exercent une influence directe sur celui-ci. La figure 8.5 donne le graphe acyclique oriente qui correspond au graphe d'influence de la figure 8.4 en effectuant l'elimination des noeuds deterministes : seules sont conservecs les quantites aleatoires sur lesquelles portera l'inference bayesienne, Pour la commodite du dessin, on a associe les variables (X mc, X uc) en un meme noeud. Le graphe de la figure 8.5 represente le raisonnement conditionnel sur lequel le modele est fonde : les fleches du raisonnement conditionnel descendent des parametres conceptuels jusqu'aux quantites observees, Pour realiser le fonctionnement de l'inferencc bayesienne, on peut imaginer les hyperparametres des lois a priori comme des nceuds parents pour les parametres du modele. Le modele interannuel est un empilement des modeles annuels. II s'appuie sur la tres forte hypothese de stabilite des parametres, permettant une coherence interannuelle, done un transfert d'information d'annee en annee. Partageant les valeurs communes de (), Q, {3, T, 6 et tt ; Ie DAG d'un tel modele s'obtient en « empilant » une repetition de structures annuelles identiques au DAG de la figure 8.5. Cette hypothese est particulierement discutable en ce qui concerne l'efficacite de la capture () et la probabilite de recapture 1r qui peuvent certainement varier d'une annee sur l'autre en fonction du debit de la rivierc et des conditions hydrometeorologiques,

8.3

Inference bayesienne

Toutes les etapes decrites ci-apres ne sont que des applications des principes et methodes apprises dans les chapitres precedents. Cependant, il nous semble utile de les appliquer a partir du DAG et explicitant les proprietes dindependance conditionnelle et de modularite, La densite conjointe a priori des parametres s'ecrit [~, (),

Q,

{3, T, 6,1rIH]

156

Pratique du calcul bayesien

K

Effectifd l 'enlTee thJ. sco rff

®

ProbabilitedesruviE

Probabilitede recapture

D

Variable latente

o

Inlenntidiaim dlitenRiniste

Variableob,ervie Paramem.de comportement

Figure 8.4 - La vie d 'un sau mon apres sa remont ee dan s le Scarff sous la forme d'un di a gram me d ' influence.

au la

lettre H rapp elle que l'on conditi onne sur un savoir initial et des hypotheses de const ruction. On le sait , l'inferen ce bayesienn e consist e a met tre cette loi a jour en impliquant les observations disponibles :

[1",0, a , (3, T , 15, 1r1H, y] ex [yll" , 0, a, (3, T , 15, 1r]

X

[1" ,0, a , (3, T , 15, 1r1H]

ou On s'en dou t e, la septuple int egration n'est pas possible

8.3.1

«a la

plume ».

L'echant.illonnage de Gibbs divise le probleme en plusieurs sous-problemes simples

Soit un poin t initial, arbit ra irement choisi dans l'espace des par ametres. En t irant to ur a tour dans chac une des sept cond itionnelles complete s, et en repet an t ce cycle un grand nombre de fois, on peut obtenir un echa nti llon de 7-uplet s provenant de la loi a posteriori conjointe des par ametres.

8. Modele de capture-recapture : applicat ion au cas des saumons

157

ProbabiIitededirlaration

D

Yariable /mente

o

Yariableob.en'Iie Parametres de comportement

Fi gure 8.5 - La vie d 'un sa umon a pres sa rernont ee dans Ie Sca rff sous la form e d 'un DAG.

8.3.2

Dans Ie DAG, la conditionnelle complete d'un noeud impliquent seulement ses nceuds parents, ses nceuds enfants et les nceuds coparents de ses enfants

Le tableau 8.3 donne pour chaque variable d'interet de l'inference bayesienne (c'est-a-dire par ametre ou variable latente stochastique), l'ensemble des variables condit ionna ntes associees. Ce tableau se const ruit a partir de la figure 8.5 ou chaque noeud a ete relie a ses nceuds parent s, fils ou coparents de ses enfants. Dans la sect ion suivante, nous verrons que certains nceuds ont une loi conditionnelle dont la st ruct ure est connue (par conjugaison); en revanche, pour d'autres nceuds, la forme de leur conditionnelle complete ne sera pas dans la bibliotheque des dist ribut ions de probabilite standa rds et il faudra l'expliciter.

8.3.3

Actualisation bayesienne des elements d'un DAG par I'echant.illonnage de Gibbs

On remarquera que seuls les nceuds stochastiques (a l'exception des observables qui sont des nceuds terminaux) peuvent et re mis a jour par le t heoreme

158

Pratique du calcul bayesien

Nceud a mettre a jour ()

(3 T

. et 80 (fig. 11.4) . Finalement, la distribution a posteriori du profil du taux d'emission du tapis au cours du temps peut etre resumee par les courbes des quantiles 5 %, 50 % et 95 % calculees pour chaque valeur du temps (fig. 11.5). La ligne continue montre la mediane et les lignes pointillees representent l'intervalle de credibilite a 90 % tandis que la ligne en gras montre le profil t emporel obtenu avec les estimations ponctuelles des auteurs (Hayter et Dowling, 1993).

18

1.7

1.8

1.5

.0 " 13

1.1

O'~

0.01

o_~

o.m

Q~ x

O~

0.00

0.07

QOO

Figure 11.4 - Correlation interpararnetres objectifs a posteriori.

C ommentaire On l'a dit, l'analyse bayesienne produit beaucoup plus de resu ltats que les methodes classiques, surtout en ce qui concerne la quantification des incertitudes. Or celles-ci doivent etre considerees lors de la prise de decision. Par exemple, on rappelle que le demi -temps de vie d 'un materiau emetteur est Ie temps necessaire pour qu e son activite diminue de moitie

8(t) = 80 ex p (- >.t ) } 8(t) =0.580

T _ln2 :::}

-

>.

Imaginons que le legisla teur fixe une norme a quarante-huit heures, c'esta-dire qu'il veut que le demi-temps de vie soit inferie ur a ce delai. Dans ce cas , si on se contente des estimations ponctuelles (11.2) , Ie tapis a perdu la moitie de son activite emet t rice apres 29 heures et cette norme est respectee. La pris e en compte des incertitudes aboutit a la conclusion inverse . Evidemment, comme toujours, ces resultats dependent de toutes les hypotheses sur lesquelles l'analyst e s'est appuye.

220

P ra tique du calcul bayesien

- '" 0"

1.2

'"

-H&D

o. 02

12

15 18 Temps (jour)

21

24

27

30

Figure 11.5 - Profil te mporel des taux d'emission et inter valle de credibilite

a 95 %.

Epilogue Ce cha pit re illustre les apports de I'an alyse bayesienne pour l'etude des emissions de formaldehydes d'un echant illon de tapis. Le profil du taux d'ernission au cours du temps du materiau et udie est I'obj ectif de I'experience, mais il n'est pas dir ect ement observable. On utilise un instrument approprie : un modele reduit de chambre aeree conte nant l'echantillon polluant . Les donn ees apparaissent comme des series discretes d'observations appariees, repr esent ant les niveaux de concent ration de polluant dan s la cha mbre , au cours du temps. Celles-ci sont utilisees pour modeliser Ie profil du niveau de concentration de polluan t dans la chambre au cours du temps qui est intrinsequement non lineair e. On peut ensuite I'u tiliser a son tour, pour est imer Ie profil du taux d'emission au cours du te mps de I'echantillon et udie. Sous Ie paradigme bayesien, un modele statistique simple nous a permis de quant ifier les incertitudes at tachees a une estimation pon ctuelle des parametres du modele. En utilisant un prior joint non informat if, nous avons utilise les techniqu es de Monte-Carlo par chaine de Markov pour calculer la dist ribu tion a post eriori mar ginale de chaque par ametre objectif. Prend re en compte les incertitudes permet des recomma nda tions operationnelles de prudence : par exemple, au vu des donn ees experiment ales, il est fort plausible que Ie profil du taux d 'emission au cours du temps ne soit pas nul passe 10 jours, mais on peut parier avec confiance qu 'il Ie sera au-d ela de 20 jours. De te ls resultats sont essent iels pou r la prise de decisions dans Ie domaine des normes de securite en sante publique.

Chapitre 12

Les avantages de la

modelisation hierarchique : application a la capture-marquage-recapture des saumons Prologue Voici un modele bayesien hieturchique (MBH) pour l'analyse des donnees de capture-marquage-recapture de saumons. Ce chapitre se presente comme une suite au chapitre 8 et s'appuie sur l'etude (Rivot et Prevost, 2002). Chaque annee i, ces deux chercheurs de l'INRA de Rennes veulent estimer le nombre inconnu Vi de saumons qui remontent la riviere Oir pour frayer ainsi que la probabilite de capture ()i du piege utilise pour effectuer ces mesures. Ils disposent d'une seric d'observations allant de 1984 a 2000 collectees sur le terrain par les techniciens de la station experimentale du Moulin des Princes, Nicolas Jeannot et Francois Burban, aides de Jean-Yves Moelo. Pour analyser de telles donnees, on peut vouloir, en premier lieu, faire I'hypothese d'indcpcndance complete entre les annees, c'est-a-dire imaginer que les resultats des experiences de capture-marquage-recapture d'une annee ne nous amenent aucune information quant aux resultats possibles des autres annees. A l'oppose, on peut etre tente d'ignorer la variabilite entre chaque annee en regroupant en un memo echantillon les donnees de toutes les annees comme si elles provenaient du meme modele d'observation. Le modele hierarchique realise un compromis astucieux entre ces points de vue extremes. II suppose que les annees ne sont ni completement identiques ni completcment independantes et considere que les ()i et les

222

Pratique du calcul bayesian

Vi sont issus d'une memc distribution de probabilite dont les parametres sont inconnus. Lorsqu'il y a peu de donnees, un modele qui suppose l'independance entre les annees menera a des inferences a posteriori pauvres. En effet, pour ces annees avec un faible effectif mesure, les donnees apportent peu d'information, ce qui produit des distributions a posteriori imprecises et difficilement exploitables. La superiorite du modele hierarchique vient de ce qu'il organise le transfert d'information entre les annees puisque ce sont des unites statistiques qui partagent une caracteristique commune. II pallie egalement un autre inconvenient de I'independance interannuelle qui conduit a des resultats beaucoup plus sensibles au choix des distributions a priori (( Gazey et Staley, 1986), (Chao, 1989)) que lorsqu'on impose une structure hierarchique.

12.1

Donnees

Les series de donnees, relativement longues mais peu abondantes (petite taille de I'echantillon), sont assez frequentes quand on veut estimer par des techniques de capture-marquage-recapture la taille d'une population sauvage durant plusieurs annees. Par exemple, sur la rivicre Oir, en Bretagne, les agents de l'INRA ont collecte des donnees sur les saumons adultes qui reviennent frayer, pour chaque annee i de 1984 a 2000. Les donnees du tableau 12.1 se presentent sous la forme suivante : Ci represente le nombre de saumons captures au piege a l'embouchure de la rivierc (station du Cerisel). Un nombre Xi de poissons captures ne sont pas relaches en amont, soit qu'ils meurent en cours de manipulation, soit qu'ils soient gardes pour des experiences ou pour la production d'oeufs. On appelle m, ~ c, - Xi le nombre de poissons marques et relaches. Ces poissons relaches dans la riviere sont marques individuellement avant de poursuivre leur remontee pour frayer. L'echantillonnage de recapture est rassemble pendant et apres le moment du frai. Appelons r, le total de tous ces poissons recaptures ou observes : parmi ceux-ci, on retrouve Yi poissons deja marques.

12.2

Modele de capture-rnarquage-recapture

Les inconnues du probleme sont evidemment le nombre de saumons (Vi) qui remontent la rivierc Oir pour frayer et la probabilite de capture ((OJ) du piege utilise l'annee i pour effectuer ces mesures, comme le schematise la figure 12.1. Sachant la valeur de ces parametres inconnus (Vi et OJ), la vraisemblance donnera la loi des variables aleatoires (C i, Xi, u; u; Yi). Dans la suite du chapitre, on utilisera le terme data pour designer l'ensemble de ces observations des donnees (Ci, Xi, tiu, r., Yi). Note 12.1 Encore une fois, pour des raisons pedagogiques, on distingue la variable aleatoire X de sa realisation x. Les lettres latines sont reservees aux

12. Les avantages de la modelisation hierarchique Annec 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000

c,

Xi

mi

r,

Yi

167 264 130 16 226 235 15 44 31 100 32 109 70 56 34 154 53

13 48 37 4 43 36 8 0 11 19 14 7 15 22 4 6 0

154 216 93 12 183 199 7 44 20 81 18 102 55 34 30 148 53

22 25 9 24 12 56 17 24 9 7 5 46 82 15 36 35 37

12 21 5 2 12 56 2 23 4 4 1 39 25 12 6 23 4

223

Tableau 12.1 - Donnees de capture-marquage-recapture pour les saumons au cours de leur remontee migratoire dans la rivicre Oir pour frayer. grandeurs aleatoires observables et les parametres inconnus sont dcsignes par des minuscules grecques. Comme tout modele, ce modele necessite egalement quelques hypotheses simplificatrices. - H1 : Tous les saumons de l'ensemble Vi sont supposes pouvoir etre individuellement et egalement pris dans le piege, avec la me me probabilite OJ. OJ est constante pendant la saison de migration. - H2 : La taille de la population de saumons en amont du piege reste invariable au cours de l'experience. Aucun saumon ne peut redescendre le courant apres avoir franchi le piege, II n'y a ni mortalite par suite du marquage, ni mortalite naturelle entre le moment du marquage et celui de la recapture. - H3 : II n'y a pas de pertes specifiques dues au marquage une fois le poisson relache en amont du piege, - H4 : La probabilite de recapture est la me me pour tous les poissons marques ou non.

12.2.1

Modele Bernoulli d'alea pour la premiere phase

Sous I'hypothese d'egalite des chances d'etre capture HI, on peut considerer la migration des saumons Vi comme des experiences de Bernoulli independantes

224

Pratique du calcul bayesien

t

Amont •

Recaptures y = nbde marques sur r captu res

Remise a l ' eau

m =c-x nb de m arques

Poissons evitant le pieg e

~_""'\ - . Piege de

remontee c = nb de captur es

~ ~~ Aval

Stock entrant

'Y =

p op. de saum ons adultes

Figure 12.1 - Schema du piege de remontee des saumons adultes.

avec une probabilite de succes of. En consequence , C i est le resultat d 'une exp erience binomiale avec Vi repetitions OU chaque saumon a une probabilite of d'etre pris :

[C· = ~

12.2.2

c·lv t

t ,

01.1]

= CC i (Ol) Ci (1 _ Ol)!li- ci Vi z Z

(12.1)

Modele Bernoulli d'alea pour la seconde phase

La loi de la variable Xi (perte par manipulation, proportion gardee pour la reproduction) n'a pas d 'importance en soi, en tout cas a l'egard du probleme de representation qui nous interesse, et on travaillera dans la suite de ce chapitre conditionnellement ala connaissance de Xi = X i ( X i connu). Les hypotheses H2H3 permettent de considerer que la difference Vi - Xi est le nombr e de saumons dans les frayeres au-dela du piege et qu 'il se trouve m, saumons marques parmi eux. Dans cet te seconde phase du precede, on cherche la loi de Y; sachant

M, =

ttu ,

En premiere approche, l'echantillonnage de recapture fonctionn e comme si l'echantillonnage de recapture r i et ait tire au hasard dans la population totale Vi-Xi . Plus exactement , l'echantillon de recapture conti ent Y i poissons marques provenant des m; marques durant la premiere phase et Z; = ri - Yi non marques parmi un nombre to tal Vi - X i - m; de poissons non marques en seconde phase. L'experience de recapture peut se concevoir comme le resultat r est le pararnetre de recapture de deux experi ences binomiales OU le parametre O

12. Les avantages de la modelisation hierarchique

a I'annee i

225

:

(12.2)

12.3

Modele bayesian hierarchique echangeable

Le modele bayesien hierarchique, note ci-apres MBH, impose une structure hierarchique sur l'efficacite de la capture et egalernent sur la taille de la population. D'une annee (ii) a l'autre (i 2 ) , les nombres de saumons ViI ou Vi2 qui remontent la riviere Oir pour frayer ne sont pas les memes, mais ils possedent quelque chose en commun : il s'agit de la meme population ecologique implantee sur la rivierc Oir, et les quantites ViI et Vi2 se ressembleront sans, bien sur, etre identiques. De la meme maniere, comme c'est toujours la meme procedure de piegeage que l'on met en ceuvre pour effectuer le marquage, la probabilite de capture (OTI ) du piege utilise l'annee i 1 partage quelques caracteristiques avec 0T2 , la probabilite de capture a I'annce i 12 . Elles ne sont pas egales, car il reste une certaine variabilite Iiee aux facteurs non maitrisables qui influent sur les deux experiences (debits de la riviere, conditions meteorologiques, etc.). II en va de meme pour la probabilite de recapture 02 de la seconde phase. Dans la suite, on appelle Oi == (OI, 0;) le parametre vectoriel des capturabilites des deux phases de chaque annee i. Le modele hierarchique rassemble toutes les annees a travers un niveau de repartition des parametres inconnus qui explique les similitudes entre l'experience de capture-marquage-recapture et la dependance parmi les tailles de population. Fondamentalement, l'hypothese hicrarchique se traduit par un niveau conditionnel supplementaire des distributions de probabilite (Gelman et al., 1995b) decrivant chaque experience annuelle. Considerons d'abord la mise en place d'une structure hierarchique pour les efficacites de capture et de recapture. Le MBR suppose que les Oi sont issus d'une meme distribution de population [Oi I,], conditionnelle a un vecteur , dhyperparametres inconnus. On attribuera une distribution a priori [,] aces hyperparametres. En effet, les efficacites annuelles, Oi, resultent dexperiences analogues utilisant le meme equipement et le meme protocole experimental. La structure hierarchique de la distribution de probabilite marque la dependance entre les Oi en exprimant les similitudes et l'heterogeneite des Oi. Des variations entre les Oi peuvent etre dues a des changements imprevisibles de l'environnement (niveau de la riviere, temperature) ou du comportement du poisson au cours des annees. Le MBR considere , comme la quantite inconnue d'une distribution unifiant toutes les annees, Ces hyperparametres reglent en particulier la variance et la moyenne des Oi : une variance nulle signifie qu'une meme valeur du parametrc Oi doit etre adoptee pour toutes les annees, tandis qu'une variance infinie pour la distribution des ()i signifie que chaque annee est independante. Entre ces deux extremes, la mise a jour par inference bayesienne de

226

Pratique du calcul bayesien

Priors sur Jet r

hyperparametres

parametres

annee i

Mode e d 'aleasn turels donnees

Figure 12.2 - Le modele hierarchique introduit un niveau de coherence interannuel par l'intermediaire des hyperpararnctres (r, 8)

la distribution des hyperpararnetres realise le transfert de I'information d'une annee sur ses voisines. Unc structure hierarchique est egalement possible pour assurer une certaine coherence interannuelle des tailles Vi de la population des saumons de l' Oir. Cette population est observee pendant plusieurs annees (correspondant aux indices i == 1, ... , I ). Le developpement (progression ou reduction) de la taille de la population depend des memes processus ecologiques quelle que soit I'annee. Les estimations derivees des annecs {1, ... ,i - 1, i + 1, ... ,I} apportent aussi des informations par rapport a la taille de la population d'une annec donnee i. On introduit aussi un niveau hierarchique pour la repartition des Vi via une loi [Vi 16] avec des hyperparametres inconnus 6 et leur propre distribution a priori

[6].

Le MBH decrit a la figure 12.2 traite de facon conjointe les series de chacune des annees i == 1, ... ,I. Les grandeurs (Vi,Oi) ont un statut mixte. Elles dependent du vecteur des hyperparametres ¢ == (r, 6) et sont des variables aleatoires non observables qui conditionnent les observables (Ci , Ii, R i ) elles recouvrent alors de ce fait un statut de parametres inconnus. La distribution a priori conjointe n(v,O,¢) repose sur deux hypotheses: prcmiercmcnt, I'independance entre (0, r) et (v,6) et deuxiemement I'cchangeabilitc de Oiet de Vi (Gelman et al., 1995b). La loi a priori sur tous les parametres s'ecrit finalement:

12. Les avantages de la modelisation hierarchique

227

La distribution a posteriori conjointe [v, Oldata] est obtenue par la combinaison de la distribution a priori jointe [v, 0, ¢] et de l'expression de vraisemblance [datalv, 0, ¢] puis par elimination sur les hyperparametres :

[v, Bldata]

0:

J

[v, B, 4>] [dataIN, B, 4>] d4>

(12.5)

La distribution a priori de (0, r) marque la dependance interannuelle des Oi. L'echangeabilite est un concept plus general que l'independance statistique. II est fonde sur I'hypothese qu'avant de voir les resultats de l'experimentation capture-marquage-recapture (en l'absence de donnees), il n'y a pas d'argument pour distinguer a priori Oi. En termes mathematiques, l'echangeabilite signifie que la distribution jointe des Oi ne change pas quand on permute les indices i. L'ordre dans lequel les donnees ont ete rassemblees n'a pas d'importance. Comme le suggere (Gelman et al., 1995b), la distribution cchangeablc la plus appropriee pour (0, r) considere chaque 0i comme un echantillon independant de la distribution conditionnelle de la taille de la population, parametres par r,7f(Oilr). Nous faisons la meme hypothese rl'echangeabilite pour les Vi. L' hypothese d'echangeabilite combinee avec l'independance entre (O,r) et (v,6"), conduit a la distribution a priori jointe (12.4). Le terme de vraisemblance [datalv, 0, ¢] est le produit des fonctions de vraisemblance annuelles [datailvi,Oi, ¢] note L, dans ce qui suit. L, est issu du modele stochastique qui sert de base au processus d'echantillonnage des experiences de capture-marquage-recapture. La vraisemblance ne depend pas du vecteur ¢ des hyperparametres d'ou la simplification I

[datalv, 0, ¢] == [datalv, 0]

==

I

II [datailvi,Oi] II t; ==

i=l

i=l

Les distributions dechantillonnage utilisees impliquees dans les L, sont des produits de formes binomiales (eq. 12.1) et (eq. 12.2) (Gazey et Staley, 1986). La distribution a posteriori complete conjointe de [v, 0, ¢] s'ecrit, a une constante de normalisation pres, grace a la formule de Bayes : (12.6) Pour obtenir la distribution a posteriori des quantites intercssantcs (eq. 12.5), il faut done integrer la distribution a posteriori complete conjointe, selon chaque composante du vecteur ¢ des hyperparamctrcs (considere ici comme un parametre de nuisance que l'on eliminc par integration)

228

Pratique du calcul bayesien

[v, eldata]

ex

~ (Dr lei I')'] [Vi 18] Li) [')'] [8] d')'d8

(12.7)

"

Le MBH permet un transfert d'informations entre unites statistiques echangeables dont Ie mecanisme probabiliste est decrit par les equations donnees dans l'annexe G.

12.4

Modele bayesian annuel

Le modele bayesien annuel (MBA) suppose l'independance interannuelle des efficacites des deux phases et du nombre d'adultes remontant la riviere Vi. Ce modele peut etre considere comme un cas particulier limite du modele hierarchique echangeable plus general. II suffit d'imaginer que la variabilite de la dispersion interannuelle est tellement grande que les ()i sont tires independamment les uns des autres depuis l'urne hierarchiquc virtuelle qui rassemble tous les ()i. hyperparam etres

o·

1

rt rs u

P i rs sur v

parametres

v.1

e

annee i

Mode e d 'aleas naturels donnees

Figure 12.3 - Le modele annuel avec independance considere chaque annee isolement par I'interrnediaire des hyperpararnetres ("'(1, .."'(1,61, ..61).

Le MBA schematisc a la figure 12.3 suppose donc a priori l'independance complete entre les ()i d'une part et les Vi d'autre part. Sous cette hypothese les donnees de capture-marquage-recapture de l'annee i ne sont utiles que pour estimer les ()i et Vi correspondant a la merne annce. Le MBA a de plus besoin

12. Les avantages de la modelisation hierarchique

229

de specifier une distribution a priori a chaque (Ji et Vi et - au prix d'un abus de langage qui facilitera la comparaison avec le MBH - on appelera ces hyperparametres avec les memes lettres (

I+1)/2

dz (A.24)

Cette dernicre integrale peut etre resolue numeriquement.

Remarque A.4 1. Le prior est non informatif si k, a, b -+ o. Dans ce cas, [J-L,7] ex 7- 1/2. - Le posterior marginal de 7 est une loi gamma, de parametre de forme n/2 et de parametre d'inverse echelle ns 2 /2. - Le posterior marginal de J-L est une loi de Student a n degres de liberte, localisee sur z et de parametre d'echelle s/ yTi. - La distribution predictive a posteriori de l'observable est une loi de Student a n dogres de liberte, localisee sur z et de parametre d'echelle sJ(n + 1) In. 2. Pour un prior informatif, le prior marginal de J-L est une loi de Student a v == 2a degres de liberte, localisee sur m et de parametre cl'echelle a == Jb/ (ka). On remarquera la similitude des expressions avec le posterior marginal de J-L.

Chapitre B

Annexe du chapitre 2 : les modeles discrets de base La lecture de cette partie plus technique est indispensable. II vous est fortement conseille de refaire les calculs au moins une fois.

Note B.1 Le sigle v. a. r. signifie variable aleatoire reelle. Le sigle pdf signifie fonction de densite de probabilite (probability density function). Par abus de langage, on peut l'utiliser pour decrire la distribution de probabilite d'une variable aleatoire discrete (v. a. d.).

Le processus de Bernoulli 1. Imaginons qu'on dispose d'une serie d'urnes remplies avec un tres grand

nombre de boules identiques sauf leur couleur qui est blanche (code 0) ou noire (code 1). On attribue un numero a chaque urne et la proportion de boules noires dans l'urne k est notee 1rk. En general, \:Ik, \:Ij i= k : 1rj i= 1rk, car chaque urne a une composition qui lui est propre. On extrait une boule de chaque urne. Les tirages sont indepcndants mais pas identiquement distribues :

2. Maintenant, imaginons une seule urne dans laquelle on realise des tirages avec remise mais sans la melanger. Les boules tirees puis remises ant donc plus de chances d'etre reprises. Les tirages sont dependants mais identiquement distribues, car la composition de l'urne ne change pas d'un tirage a l'autre :

274

Pratique du calcul bayesien

3. Ensuite, on considere une seule urne contenant un nombre (pas trop grand) de boules blanches et noires en proportion inconnue dans laquelle on effectue des tirages avec remise en y ajoutant chaque fois m boules de la meme couleur (tirages de Polya). Ici, les tirages sont dependants et la composition de l'urne change a chaque tirage (7i"t+l -1= 7i"t) :

4. Enfin, on considere une seule urne contenant des boules blanches et noires en proportion inconnue et on effectue des tirages avec remise en prenant bien soin de la melanger avant chaque nouveau tirage. II est clair que nous sommes dans le cas OU les tirages sont uulepeiulants et identiquement distrioue» (hypothese iid) :

(B.1) Cette derniere procedure dechantillonnage est connue sous Ie nom de processus de Bernoulli.

L'hypot.hese iid L'hypothese iid est tres importante en statistique. D'une maniere generale, supposer l'echantillon iid revient a admettre que les donnees seront tirees independamment les unes des autres dans la meme loi de probabilite, Cette hypothese est done toujours eonditionnelle au modele d'echantillonnagc ehoisi, lequel est caracterise par un parametrc () (notation generique) de dimension finie (p. ex. dim () == 2 pour une loi normale).

La distribution gamma La pdf d'une variable aleatoirc X definie sur l'intervalle reel [0,oc] est une loi gamma de parametrc de forme a > 0 et de parametre d'echelle b > 0 si et seulement si :

[xla, b] = r

(~) b x a

a

-

1

exp ( -~)

(B.2)

Le reel T (a) est defini par I'integralc d' Euler suivante :

1 00

a> 0: r(a) =

u a - 1 exp(-u) du

(B.3)

L'integrale d' Euler (eq. B.3) est dite fonction eulerienne de premiere espece.

Exercice B.I Montrez que

E (X) == ab,Var (X) == ab2

(B.4)

Annexe B

275

Notons qu'il est courant de definir la fonction de densite de probabilite gamma en utilisant un parametre d'echcllc inverse (c == lib> 0) :

(B.5) L'integrale d' Euler n'est rien d'autre qu'une generalisation de la fonction factorielle :

(B.6)

n!==f(n+l) Quel que soit le reel positif, a on a : I' (a + I) == af (a)

(B.7)

La distribution beta La densite de la distribution de probabilite d'une variable aleatoire X definie sur l'intervalle reel [0,I] suit une loi Beta de parametres r > 0 et s > 0 si et seulement si I r-1 ( )8-1 (B.8) [x I r, ] S = B (r, s) x 1- x ou Ie reel B (r, s) est defini par l'integrale d'Euler suivante dite fonction eulerienne de seconde espece : r, S

> 0 : B (r, s)

=

II

ur -

1

(1 -

ur-

1

du

(B.9)

II existe un lien entre les fonctions euleriennes gamma et beta :

B( r,s )

== f(r)f(s) f(r+s)

(B.IO)

Cette identite sera tres souvent utilisee,

Esperance d'une v. a. r. X distribuee selon une loi beta sur

[0,1]

Par definition, la valeur attendue ou esperance mathematique de X est

E (X) =

II

x [xlr, s] dx

Ce calcul est trivial

E(X)

=

B(r+1,s) B(r,s)

=

_r_ r+s

(B.II)

276

Pratique du calcul bayesien

Variance d'une v. a. r. X distribuee selon une Ioi beta sur [0,1] Par definition, la variance de X est

Var (X)

:=

E (X 2) - [E (X)]2

Par consequent

La variance suit

Var(X)

= B(r+2,s) B(r,s)

_

(_r_)2 r+s

Or

B(r+2,s) B(r,s)

----:=

f(r+2)f(s) f(r+s) x--f(r+2+s) f(r)f(s) (r+l)rf(r) f(r+s) --------- x --(r+s+l)(r+s)f(r+s) I'{r) r+l r ---x--

r+s+l

r+s

Finalement

Var(X) _ _ r_ ( r+ 1 _ _ r_) r+s r+s+l r+s := _r_ ((r+l)(r+S)-r(r+s+l)) r+s r+s+l

= r:s

C+:+l) rs

(B.12)

(r+s)(r+s+l)

Mode d'une v. a. r. definie sur [0,1] et distribuee selon une loi beta II suffit d'annuler la derivee premiere de la densite :

!(x)=x r - 1 ( 1 - x r - 1

*

df dx

=o¢}X

=

s+r#2

S

r-l

+r -

2

=.M

(B.13)

Annexe B

277

La loi de Poisson comme limite de la loi binomiale On part de la loi binomiale :

Le nombre de combinaisons que l'on peut faire en prenant n objets par paquet de x peut encore s'ecrire

n., xl (n - x)!

=

X!!!

nx

z. )

x-1 (

1 - ;,

Par consequent

Pr (X = xlO, n) =

(nf)) x

ot!!

(1 _ f)) n

--;y- (1-

x -1 (

i )

1-;,

Faisons tendre le nombre d'essais n vers l'infini, la probabilite de succes f) vers zero et leur produit vers une limite finie A E On a:

IRt.

lim TI~-=-1 (1 1,-0

n-+oo

lim (1 - f)) -

0-+0

lim

n-+oo

(1 -

x

i) == 1 n

== 1

~)n n

== exp (-A)

On obtient la loi de Poisson qui est une loi d'evenements rares : Pr (X

AX == XIA) == -, exp (-A) x.

(B.14)

La distribution binorniale negative C'est la loi du nombre d'echecs y avant d'obtenir Ie r-ieme succes (r ~ 1). Le nombre d'epreuves z avant d'obtenir le r-ieme succes decoule de la loi binomiale, car en z - 1 epreuves on a r - 1 succes :

Le nombre d'echecs est y == z - r. Par consequent

Tenant compte de

x! = f (x + 1) = z.F (x),

B (a, b) = r~(2~~~)

278

Pratique du calcul bayesien

on a: 1

( y+r- l ) ==

B(y,r)

r-l

x~ y

Finalement, la distribution de probabilite binomiale negative s'ecrit : (B.I5)

La predictive a posteriori d'un modele gamma-Poisson On sait que sur une periode de longueur l, le posterior s'ecrit

[Alx] ex Ax +a - 1 exp (-A (l + b- 1 ) )

(B.I6)

ou a et b sont respectivement le parametre de forme et le parametre d'echelle du prior gamma. QueUe est la probabilite d'observer y evenements sur la future periode h sachant que, dans le passe, on en avait observe x sur la periode l ?

[ylh, x, I] = =

1

00

[yl)" h] [),Ix, I] d)'

y (l

r

b-1)x+a

~ + y! r (x + a) Jo hY (l+b-1)x+a

11

00

),y+x+a-l

exp (-), (h + I + b- 1 ) ) d)'

f(y+x+a) (h+l+b-1)y+x+a

r(x+a)

f(y+x+a) hY (l+b-1)x+a r (y) r (x + a) (h + l + b-1)y+x+a

Y

1

B (y, x 1

(l+b-1)x+a

hY

+ a) Y

1 (

= B(y,x+a)y

(h + I + b-1 )y+x+a I + b:'

h+l+b- 1

) x+a (

h

h+l+b- 1

) y

(B.17)

Posons

(B.I8)

r==x+a l + b- 1

7r

==

h + I + b-

1

{::}

1-

7r

h

== - - - -1 h+l

+ b-

(B.I9)

II vient

[ylh, x, I]

= B (

1

) ~1fr (1 -

y,r y

1f)Y

(B.20)

En comparant cette derniere avec la distribution B.I5, on voit que x joue le role de r en l'etendant aux reels positifs.

+a

Chapitre C

Annexe du chapitre 6 : le modele des fuites et Ie modele GEV sous WinBUGS Du processus ponctuel de Poisson au modele des fuites On s'interesse a un evenement ponctuel (p. ex. un point sur un axe ou un pixel sur une surface) marque par une certaine intensite, Sur une fenetre donnee (p. ex. une periode de temps fixee, un troncon de longueur fixee, une surface d'aire fixee) , le nombre d'occurrences est note N d'intensite Z == (Zl,'" ,ZN). Quand on sait que N == n et que Z == z, on dit que l'information est complete. Cependant, il existe des situations OU on ne dispose que du cumul des intensites, Dans un tel cas, Nest une variable latente et on dit que l'information est incomplete. Le modeles des fuites est du a (Morlat, 1968) pour representer les fuites sur les conduites de gaz. II ignorait leur nombre et done leur intensite respective, mais il connaissait le cumul des pertes par la difference entre les debits d'entree et de sortie du troncon d'interet.

Le processus ponctuel de Poisson Pour en simplifier l'expose, nous nous refererons a des tops arrivant au hasard sur l'axe du temps. A chaque date t, on peut associer une variable de Bernoulli qui prend la valeur 1 avec la probabilite 7f si un top est observe a cette date. Un processus ponctuel est une suite de variables de Bernoulli indeperulomtes et identiquement distribuees (hypothese iid).

280

Pratique du calcul bayesien

Le processus ponctuel de Poisson est un modele statistique parametrique fonde sur trois hypotheses. HI. Le processus est sans memoire, c'est-a-dire que la probabilite d'observer 1 evenement sur une periode de longueur h suffisamment petite est proportionnelle a h :

[N == llA, h] == Ah + 0 (h) -

0

(h) represente un infiniment petit par rapport a h himo(h) ==0 h---+O

h

- A est l'intensite du processus, supposee invariante dans Ie temps. H2. II n'y a pas de simultaneite, c'est-a-dire que la probabilite d'observer plus d'un evenement sur une periode de longueur h est negligeable si h est petite:

[N ~ 21h] == 0 (h) H3. Les evenements qui se produisent sur des periodes disjointes, soit hI et h2 , sont independants

Sur cette base on montre que le nombre de tops, N, sur une periode unite, c'est-a-dire une fenetre dont la longueur est egale a 1 unite de temps (p. ex. le mois) , est distribue selon une loi de Poisson de parametre A == E (N), d'ou Ie nom du processus :

An

[nIA] == exp (-A) , n.

(C.l)

Le processus ponctuel de Poisson marque On ajoute une quatrieme hypothese au processus ponctuel de Poisson. H4. Les intensites Z, des occurrences sont independantes de leur nombre N, independantes entre elles et identiquement distribuecs selon une loi exponentielle de parametre p telle que E(Zilp) == lip:

Vi Vi, vi Vi Conditionnellement

=1=

z, -l N i: z, -l z, Zilp

r-;»

(C.2)

dgamma (1, p)

a n et a p, Ie cumul H ==

n

E Z; est

i=l

une loi gamma

Hlp, n

0:.

dgamma (n,

p)

distribue selon

Annexe C

281

Soit une periode unite (indice t) sur laquelle on a observe nt tops dont le cumul des intensites est ht . La vraisemblance de cet echantillon est triviale :

(C.3) Si ce processus est stationnaire sur T periodes independantes de meme longueur (L == 1), la vraisemblance de I'echantillon d == {(nt, ht ) : t == 1, 2, ... T} est simplement

(C.4)

ou Sn ==

T

L:: tu,

t=1

Sh ==

T

L:: ht

t=1

Le prior le plus simple postule l'independance de A et de p avec A dgamma (a>.., b>..) et p rv dgamma (ap , bp ) :

Alsn, a>.., b>.. plsn, Sh, a p , bp

dgamma (sn + a>.., T + b>..) dgamma (sn

+ ap , Sh + bp )

(C.5)

Ainsi, a posteriori, la valeur attendue du nombre de tops reste independante de leur intensite : A ..1 p.

Le modele de depassement Soit une observable Y qui evolue dans le temps. Un top arrive quand cette observable depassc un certain seuil u fixe. Le nombre de tops sur une periode unite, par exemple l' annee, et leur intensite respective (au-dessus du seuil) constitue un processus de Poissonmarque. Les marques au-dessus du seuil sont supposees iid selon une certaine loi. Le modele POT du chapitre 6 postule que si le seuil est assez haut, cette loi est la distribution de Pareto generalisee, Dans le cas d 'une loi exponentielle, on retrouve les resultats indiques ci-dessus.

Le modeles des fuites Reprenons le processus de Poisson marque. En cas d'information imparfaite, la seule observable est le cumul des intensites Hs, c'est-a-dire que lc nombre de tops N, est une variable latente intervenant dans la loi du cumul Hi. En posant h == (hI,··· ,h T , ) et N == (N I , · · · ,NT), la vraisemblance completee s'ecrit :

Avec les priors utilises ci-dessus, les conditionnelles completes a posteriori sont respectivement -Xlh,N,p

dgamma (SN + a>.., T + b>..)

plh,N,A

dgamma (SN + ap , Sh + bp )

(C.6)

282

Pratique du calcul bayesien

On remarquera la similitude des relations C.5 et C.6. Bien sur, dans la seconde, la somme des variables latentes SN == N 1 + ... + NT est inconnue. II faut donc ajouter un module pour realiser l'inference via un echantillonnage de Gibbs. Pour t == 1, ... ,T, la conditionnelle complete de la variable latente N, n'est pas standard mais peut etre definie sur une grille (voir chap. 4). En posant N_ t == N\ {Nt} on a :

Un algorithme de Gibbs (voir chap. 4) est facile simuler Alh, N, p

rv

dgamma

(SN

simuler plh, N, A rv dgamma (SN

a programmer:

+ a>.., T + b>..) + a p , Sh + bp )

Pour t == 1,2,··· ,T : simuler [Ntlp, A, N_ t , h] On peut voir le modele des fuites comme le modele des depassernents avec un seuil nul alors que l'information est imparfaite.

Les valeurs extremes sous WinBUGS Le modele GEV Soit N blocs de n observations (n assez grand). Sur chacun d'entre eux, on s'interesse au maximum Zk pouvant prendre la valeur Zk (k == 1, ... ,N). La contribution de l'observation Zk a la vraisemblance est donnee par la densite :

Ce n'est pas une densite standard de WinBUGS, mais on s'en sort en utilisant l' astuce « zeros trick », La densite de Poisson s'ecrit

[yIA]

== exp

AY

(-A) ,

y.

Si on ne tire que des zeros, la contribution d'une observation blance est simplement [OIA] == exp (-A)

a la vraisem-

Ainsi WinBUGS considcre un ensemble de donnees constitue de N zeros tires dans une loi de Poisson de parametre Ak == -In [zkIB] + C.

Remarque e.l La constante C assure Ak > 0 et ne pose aucun probleme puisque la vraisemblance est definie a une constante pres.

Annexe C

283

Choix du prior En general, le savoir a priori est tres reduit et il n'y a aucune raison de lier Ie parametre d'echelle, p, au parametre de forme, {3. En revanche, le parametre de localisation, /1, est lie au parametre de forme a cause de la condition

Des lors, le prior conjoint peut s' ecrire comme suit :

[0] == [{3, JL, p] == [JLI{3] [{3] [p] Le choix classique (et judicieux) pour un parametre d'echelle est une loi gamma et une loi normale pour un parametre de localisation :

p~

dgar,nr,na(f,e)

{3 ~ dnormim, t) Pour p, on obtient un prior non informatif avec f == e ---t o. Pour {3, on l'obtient avec m == 0 et une precision t ---t 00 (WinBUGS prend par defaut

f == e == 10- 3 , t == 10- 6 ) .

Remarque C.2 Par experience, on sait que 1{31 est inferieur a quelques unites. Par consequent, t == 10- 3 est suffisant (question de vitesse de convergence). Remarque C.3 On peut aussi poser ¢ == -In pet prendre ¢ ~ dnorr,n(O, 10- 6 ) (ce sera notre choix). Pour [/1113], un prior non informatif est une loi uniforme sur un intervalle reel dependant de 13: [JL I13] ~ dunif(r,s). En posant a == min {Zk} - 13- 1 et b == max { Zk} - {3-1, Ic respect de la condition entraine - 13 > 0 => JL > b (Weibull) - 13 < 0 => JL < a (Frechet) Dans WinBUGS, on tire u ~ dunif( -00, a) et v ~ dunif(b,oo) et on construit le prior sur JL comme suit: JLI{3, a, b f -

A l'issue

U

* step (-{3) + v * step (13)

de l'inference, on revient au parametrage initial du modele GEV

avec

(J

~

== exp (1)) == - 13 exp (¢)

284

Pratique du calcul bayesien

DAG associe au modele GEV Le DAG simplifie ci-dessous (fig. C.l) montre le modele GEV via le « zeros trick» de WinBUGS OU le parametre de Poisson s'ecrit :

Figure C.l - DAG du modele GEV sous WinBUGS.

On rappelle que le niveau de retour zp associe ala periode de retour T == lip est donne par la relation :

z == JL p

e#o

~ (1 - x-e) = ~

p

/3#0

J-l + ~

(3

(1 -

x/3 exP(¢)) p

OU

x p=:-ln(l-p) Le code WinBUGS est le suivant. Pour le niveau de la mer a Port Pirie, les resultats du tableau C.1 sont obtenus apres 40000 iterations dont 20000 pour la pcriode de chauffe. Deux chaines sont lancees pour controler la convergence. Ainsi le maximum annuel du niveau de la mer a Port Pirie converge en loi vers la loi des extremes de Weibull (~ < 0). Chaque annee de la periode 1923-1987, il y a une chance sur cent (p == 0.01) d'observer une hauteur d'eau superieure a 4.80 m (avec un risque d'erreur fixe a 5 %).

Annexe C

285

model; { a < -zmin-1/beta b