International audienceEtant donné une probabilité $\mu$ sur $\R^d$ ($d$ grand), on note $X$ un vecteur aléatoire générique de loi $\mu$ et $\Phi:\R^d\rightarrow\R$ une application ``boîte noire''. Un réel $q$ étant fixé, le but est de générer un échantillon i.i.d. $(X_1,\dots,X_N)$ tel que pour tout $i$ : $X_i\sim{\cal L}(X|\Phi(X)>q)$. Lorsque $q$ est grand comparé aux valeurs typiques de la variable $\Phi(X)$, la méthode Monte-Carlo classique devient trop coûteuse. Dans ce travail nous présentons et analysons une nouvelle approche pour ce problème. Celle-ci procède en plusieurs étapes, s'inspirant de l'algorithme de Metropolis-Hastings et des méthodes dites multi-niveaux en estimation d'événements rares. Deux problèmes peuvent être traités très facilement via cette nouvelle méthode : estimation de quantiles extrêmes et estimation d'événements rares. Les idées présentées seront illustrées sur un problème de tatouage numérique

Guyader, Arnaud

Hengartner, Nicolas

Matzner-Løber, Eric

INRIA a CCSD electronic archive server

HAL Id: inria-00494788https://hal.inria.fr/inria-00494788Submitted on 24 Jun 2010HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Estimation de quantiles extrêmes et de probabilitésd’événements raresArnaud Guyader, Nicolas Hengartner, Eric Matzner-LøberTo cite this version:Arnaud Guyader, Nicolas Hengartner, Eric Matzner-Løber. Estimation de quantiles extrêmes et deprobabilités d’événements rares. 42èmes Journées de Statistique, 2010, Marseille, France, France.￿inria-00494788￿Estimation de quantiles extrêmes et deprobabilités d’événements rares1Arnaud GUYADER a , Nicolas W. HENGARTNER b et Eric MATZNER-LØBER aa Université Rennes 2 – Haute BretagneCampus VillejeanPlace du Recteur Henri Le Moal, CS 2430735043 Rennes Cedex, Francearnaud.guyader@uhb.fr,eml@uhb.frb Information Sciences Group, MS B256Los Alamos National LaboratoryNM 87545, USAnickh@lanl.govMots-clés Extrêmes, Méthodes bayésiennes.Résumé Etant donné une probabilité µ surRd (d grand), on note X un vecteur aléatoiregénérique de loi µ et Φ : Rd → R une application “bôıte noire”. Un réel q étant fixé, le butest de générer un échantillon i.i.d. (X1, . . . , XN) tel que pour tout i : Xi ∼ L(X|Φ(X) > q).Lorsque q est grand comparé aux valeurs typiques de la variable Φ(X), la méthode Monte-Carlo classique devient trop coûteuse. Dans ce travail nous présentons et analysons unenouvelle approche pour ce problème. Celle-ci procède en plusieurs étapes, s’inspirant del’algorithme de Metropolis-Hastings et des méthodes dites multi-niveaux en estimationd’événements rares. Deux problèmes peuvent être traités très facilement via cette nouvelleméthode : estimation de quantiles extrêmes et estimation d’événements rares. Les idéesprésentées seront illustrées sur un problème de tatouage numérique.Abstract Let X denote a generic random vector with probability distribution µ onRd, and let Φ be a mapping from Rd to R. That mapping can be a black box, e.g., theresult from some computer experiments for which no analytical expression is available.Our goal is to generate an i.i.d. sample (X1, . . . , XN) with Xi ∼ L(X|Φ(X) > q) togetherwith the probability P[Φ(X) > q] for any arbitrary real number q. When q lies far out inthe right-hand tail of the distribution of the random variable Φ(X), a naive Monte-Carlo1Travail effectué dans le cadre des projets Nebbiano, ANR-06-SETI-009, et LANL LDRD 20080391ER.1simulation becomes computationally intractable. In this article we present and analyse anovel simulation algorithm for this problem. It proceeds by successive elementary steps,each one being based on Metropolis-Hastings algorithm. Our technique is useful for bothestimating the probability of events and estimating extreme quantiles. We demonstratethe practical usefulness of our method by applying it to a problem in watermarking.1 Rappels sur les méthodes Monte-CarloNotons X un vecteur aléatoire générique de dimension d et de loi µ sur l’espace probabilisé(Ω,F ,P). On suppose pouvoir simuler des réalisations i.i.d. selon µ. On considère alorsune fonction « bôıte noire » Φ : Rd → R, c’est-à-dire que l’on sait seulement évaluer Φ(x)en tout point x de Rd. Le dernier paramètre est un nombre réel q. Dans ce contexte, notrebut est de simuler un échantillon i.i.d. (X1, . . . , XN) de N particules distribuées suivantla restriction de la loi de X à l’événement R = {Φ(X) > q}. La difficulté vient du faitqu’on s’intéresse typiquement au cas où la probabilité p = P(R) = P(Φ(X) > q) est trèspetite, typiquement inférieure à 10−9.Face à cette situation, la méthode Monte-Carlo classique, ou näıve, consiste à simu-ler un échantillon i.i.d. (X1, . . . , Xn) selon la loi µ et à ne garder que les N éléments,notés (X1, . . . , XN) pour simplifier, tel que ∀i = 1, . . . , N : Φ(Xi) > q. Il est clair quel’échantillon (X1, . . . , XN) a la propriété voulue, à savoir que Xi ∼ L(X|Φ(X) > q).Néanmoins, la loi forte des grands nombres implique que la proportion d’intérêt N/ntend presque sûrement vers p lorsque n tend vers l’infini. Cette méthode nécessite doncde simuler un échantillon initial de taille n ≈ N/p pour aboutir à un échantillon de tailleN . De fait, lorsque p est très petit, ceci devient numériquement déraisonnable.Supposons maintenant qu’on veuille estimer p, avec toujours p très petit. La méthodeMonte-Carlo classique consiste à simuler un échantillon i.i.d. (X1, . . . , XN) selon la loi µ,et à considérer l’estimateurp̂mc =N̂mcN=Card{i ∈ {1, . . . , N} : Φ(Xi) > q}N.Puisque N̂mc suit une loi binomiale de paramètres (N, p), il est clair que p̂mc est non biaiséet de variance relative :V (p̂mc)p2=1 − pNp≈ 1Np.Par conséquent, la taille N de l’échantillon à simuler doit être au moins de l’ordre de 1/ppour atteindre une précision acceptable, ce qui devient rédhibitoire lorsque l’événementd’intérêt est très rare.2Face à cette situation, une astuce classique consiste à effectuer un échantillonnage préférentiel(ou Importance Sampling), lequel revient à simuler l’échantillon (X1, . . . , XN) selon uneautre loi ν, pour laquelle l’événement R = {Φ(X) > q} n’est plus rare. On estime alorsla probabilité comme ci-dessus et on la corrige grâce au lien entre ν et µ. Ceci supposenéanmoins que l’on ait de l’information a priori sur le modèle afin de proposer une nou-velle loi pertinente ν d’échantillonnage. Mais puisqu’on suppose ne rien connâıtre sur Φ(qui est une bôıte noire), cette méthode est ici exclue.Supposons maintenant le problème inverse : l’estimation d’un quantile extrême. On sup-pose p très petit fixé le but est de trouver le quantile q tel que P(Φ(X) > q) = p.La méthode Monte-Carlo classique marche comme suit : simuler un échantillon i.i.d.(X1, . . . , XN) selon la loi µ, les trier de sorte que Φ(X(1)) < · · · < Φ(X(N)) et considérerl’estimateur :q̂mc = Φ(X(⌊(1−p)N⌋)),où ⌊y⌋ désigne la partie entière de y. Si la fontion de répartition F (y) = P(Φ(X) ≤ y)admet une dérivée non nulle au point q, alors on peut montrer (voir par exemple Arnoldet al. [1], p.128) que le biais est en O(1/N) et l’écart quadratique de la forme suivante :Var(q̂mc) =p(1 − p)N + 2· 1f 2(q)+ o(1/N) ≈ 1N· pf 2(q).2 L’algorithmePlutôt que la méthode Monte-Carlo classique, nous proposons l’algorithme suivant :– Simuler un échantillon i.i.d. (X1, X2, . . . , XN) selon µ et initialiserX11 = X1, . . . , X1N = XN .– Pour m = 1, 2, . . ., noterLm = min(Φ(Xm1 ), . . . , Φ(XmN )),et définirXm+1j ={Xmj si Φ(Xmj ) > LmX⋆ ∼ L(X|Φ(X) > Lm) si Φ(Xmj ) = Lm.– En notantM = max{m : Lm ≤ q},ce maximum étant par convention égal à 0 si L1 > q, nous verrons que LM ≈ q avecune très bonne précision. Ainsi l’échantillon(X1, . . . , XN) = (XM+11 , . . . , XM+1N )est i.i.d. de loi approximativement égale à L(X|Φ(X) > q).3Le point délicat de l’algorithme ci-dessus est de réussir à simuler une nouvelle particuleX⋆ distribuée selon la loi L(X|Φ(X) > Lm). L’idée est d’utiliser un noyau de transitionK qui soit µ−réversible, c’est-à-dire tel que :∀(x, x′) ∈ Rd ×Rd µ(dx)K(x, dx′) = µ(dx′)K(x′, dx).En particulier, la loi µ est alors stationnaire pour le noyau K. Partant d’un noyau instru-mental q supposé µ−irréductible, l’algorithme de Metropolis-Hastings (voir [4], [3] et [5])permet précisément d’en déduire un noyau K ayant cette propriété.A chaque étape m, notons alors Am = {x ∈ Rd : Φ(x) > Lm}, µm la restrictionnormalisée de µ à Am, c’est-à-dire µm(dx) =1µ(Am)1Am(x) µ(dx), et considérons le noyaude transition Km défini par :Km(x, dx′) = 1Acm(x) δx(dx′) + 1Am(x)(K(x, dx′)1Am(x′) + K(x, Acm) δx(dx′)).Le principe de construction de ce noyau est très simple : partant de x > Lm, une tran-sition x  x′ est proposée par le noyau K. Si Φ(x′) > Lm, la transition est acceptée,sinon elle est rejetée et x reste à la même place. Il est alors facile de voir que la loi µm eststationnaire pour le noyau de transition Km.Ainsi, à l’étape m, le but est de simuler X⋆ ∼ L(X|Φ(X) > Lm), avec X⋆ indépendant des(N − 1) particules Xmj vérifiant déjà Φ(Xmj ) > Lm. L’idée est tout simplement d’appelerV0 un vecteur aléatoire correspondant à l’une de ces particules (tirée au hasard) et deconsidérer la châıne de Markov (Vn)n≥0 de noyau de transition Km. Pour n « assez grand »,la particule X⋆ = Vn est distribuée suivant la loi L(X|Φ(X) > Lm) et est indépendantedes autres particules.3 Applications3.1 Estimation d’événement rareAvec les mêmes notations que précédemment, supposons qu’on veuille estimer la proba-bilité p = P(Φ(X) > q). Il suffit alors de considérer l’estimateur :p̂ =(1 − 1N)M.On montre que la variable aléatoire M suit une loi de Poisson, ce qui permet d’en déduirela loi de p̂ ainsi qu’une expression simple pour tous ses moments.4Proposition 3.1 La variable aléatoire M suit une loi de Poisson de paramètre −N log p.Il s’ensuit que pour tout k ∈ N :E[p̂k] = pN(1−(1−1/N)k).Ainsi p̂ est un estimateur sans biais de p, de variance :Var(p̂) = p2(p−1N − 1).Lorsque p est très petit, p̂ est donc un estimateur bien plus précis que celui obtenu parl’algorithme de Monte-carlo näıf :V (p̂)p2≈ − log pN≪ 1 − pNp=Var(p̂mc)p2.D’autre part, puisque le nombre N de particules est au moins de l’ordre de la centaine,l’approximation de la loi de Poisson P(−N log p) par la loi normale N (−N log p,−N log p)est légitime. Un intervalle de confiance de niveau 95% pour p est ainsi donné par :p̂e−2√− log p̂N ≤ p ≤ p̂e+2√− log p̂N .3.2 Estimation de quantile extrêmeRéciproquement, une probabilité p très petite étant fixée, supposons que l’on cherche àdéterminer le quantile q tel que P(Φ(X) > q) = p. Avec l’algorithme proposé, il suffit denoter :m =⌈log(p)log(1 − N−1)⌉≈ ⌈−N log p⌉,et d’estimer le quantile q par q̂ = Lm. On montre que le biais de cet estimateur est enO(1/N) et que sa variance est contrôlable.Proposition 3.2 Si la variable Φ(X) admet une densité f(q) non nulle en q, alors l’er-reur quadratique est asymptotiquement bornée comme suit :p2(− log p − 2)f(q)2≤ limN→∞N ·E [(q̂ − q)2] ≤ p2(− log p + 2)f(q)2.Ici encore, l’estimation proposée par q̂ est bien plus précise que celle offerte par l’estimateurMonte-Carlo classique q̂mc, puisque :Var(q̂) ≈ 1N· −p2 log pf(q)2≪ Var(q̂mc) ≈1N· pf 2(q).5Comme pour la comparaison entre p̂ et p̂mc, ceci montre que pour atteindre la mêmeprécision, l’estimateur Monte-Carlo classique nécessite environ (−p log p)−1 fois plus departicules que la méthode proposée dans ce papier.Notons au passage que l’écart quadratique de notre estimateur fait intervenir la valeurde la densité f au point q, laquelle est bien sûr inconnue. Ceci n’empêche néanmoinsnullement la construction d’intervalles de confiance pour q̂.Proposition 3.3 Notons :{m+ =⌈−N log p + 2√−N log p⌉m− =⌊−N log p − 2√−N log p⌋Alors l’intervalle I1−α(q) = [Lm− , Lm+ ] est un intervalle de confiance à 95% pour le quan-tile q.4 Comparaison et illustrationLors de la présentation, la méthode proposée ici sera comparée à l’état de l’art en estima-tion d’événements rares (notamment au récent travail de Cérou et al. [2]) et illustrée surune application en tatouage numérique (watermarking).Références[1] B.C. Arnold, N. Balakrishnan, and H.N. Nagaraja. A first course in order statistics.Wiley Series in Probability and Mathematical Statistics : Probability and Mathe-matical Statistics. John Wiley & Sons Inc., New York, 1992. A Wiley-IntersciencePublication.[2] F. Cérou, P. Del Moral, T. Furon, and A. Guyader. Rare event simulation for a staticdistribution. Technical Report RR-6792, Inria, January 2009. Submitted.[3] W.K. Hastings. Monte Carlo sampling methods using Markov chains and their appli-cations. Biometrika, 57(1) :97–109, April 1970.[4] N. Metropolis and S. Ulam. The Monte Carlo method. J. Amer. Statist. Assoc.,44 :335–341, 1949.[5] C.P. Robert and G. Casella. Monte Carlo statistical methods. Springer Texts inStatistics. Springer-Verlag, New York, 1999.6

Estimation de quantiles extrêmes et de probabilités d'événements rares

HAL-Rennes 1

Estimation de quantiles extreˆmes et de probabilite´sd’e´ve´nements raresArnaud Guyader, Nicolas Hengartner, Eric Matzner-LoberTo cite this version:Arnaud Guyader, Nicolas Hengartner, Eric Matzner-Lober. Estimation de quantiles extreˆmeset de probabilite´s d’e´ve´nements rares. 42e`mes Journe´es de Statistique, 2010, Marseille, France,France. 2010. <inria-00494788>HAL Id: inria-00494788https://hal.inria.fr/inria-00494788Submitted on 24 Jun 2010HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Estimation de quantiles extreˆmes et deprobabilite´s d’e´ve´nements rares1Arnaud GUYADER a , Nicolas W. HENGARTNER b et Eric MATZNER-LØBER aa Universite´ Rennes 2 – Haute BretagneCampus VillejeanPlace du Recteur Henri Le Moal, CS 2430735043 Rennes Cedex, Francearnaud.guyader@uhb.fr,eml@uhb.frb Information Sciences Group, MS B256Los Alamos National LaboratoryNM 87545, USAnickh@lanl.govMots-cle´s Extreˆmes, Me´thodes baye´siennes.Re´sume´ Etant donne´ une probabilite´ µ surRd (d grand), on noteX un vecteur ale´atoirege´ne´rique de loi µ et Φ : Rd → R une application “boˆıte noire”. Un re´el q e´tant fixe´, le butest de ge´ne´rer un e´chantillon i.i.d. (X1, . . . , XN) tel que pour tout i :Xi ∼ L(X|Φ(X) > q).Lorsque q est grand compare´ aux valeurs typiques de la variable Φ(X), la me´thode Monte-Carlo classique devient trop couˆteuse. Dans ce travail nous pre´sentons et analysons unenouvelle approche pour ce proble`me. Celle-ci proce`de en plusieurs e´tapes, s’inspirant del’algorithme de Metropolis-Hastings et des me´thodes dites multi-niveaux en estimationd’e´ve´nements rares. Deux proble`mes peuvent eˆtre traite´s tre`s facilement via cette nouvelleme´thode : estimation de quantiles extreˆmes et estimation d’e´ve´nements rares. Les ide´espre´sente´es seront illustre´es sur un proble`me de tatouage nume´rique.Abstract Let X denote a generic random vector with probability distribution µ onRd, and let Φ be a mapping from Rd to R. That mapping can be a black box, e.g., theresult from some computer experiments for which no analytical expression is available.Our goal is to generate an i.i.d. sample (X1, . . . , XN) with Xi ∼ L(X|Φ(X) > q) togetherwith the probability P[Φ(X) > q] for any arbitrary real number q. When q lies far out inthe right-hand tail of the distribution of the random variable Φ(X), a naive Monte-Carlo1Travail effectue´ dans le cadre des projets Nebbiano, ANR-06-SETI-009, et LANL LDRD 20080391ER.1simulation becomes computationally intractable. In this article we present and analyse anovel simulation algorithm for this problem. It proceeds by successive elementary steps,each one being based on Metropolis-Hastings algorithm. Our technique is useful for bothestimating the probability of events and estimating extreme quantiles. We demonstratethe practical usefulness of our method by applying it to a problem in watermarking.1 Rappels sur les me´thodes Monte-CarloNotons X un vecteur ale´atoire ge´ne´rique de dimension d et de loi µ sur l’espace probabilise´(Ω,F ,P). On suppose pouvoir simuler des re´alisations i.i.d. selon µ. On conside`re alorsune fonction « boˆıte noire » Φ : Rd → R, c’est-a`-dire que l’on sait seulement e´valuer Φ(x)en tout point x de Rd. Le dernier parame`tre est un nombre re´el q. Dans ce contexte, notrebut est de simuler un e´chantillon i.i.d. (X1, . . . , XN) de N particules distribue´es suivantla restriction de la loi de X a` l’e´ve´nement R = {Φ(X) > q}. La difficulte´ vient du faitqu’on s’inte´resse typiquement au cas ou` la probabilite´ p = P(R) = P(Φ(X) > q) est tre`spetite, typiquement infe´rieure a` 10−9.Face a` cette situation, la me´thode Monte-Carlo classique, ou na¨ıve, consiste a` simu-ler un e´chantillon i.i.d. (X1, . . . , Xn) selon la loi µ et a` ne garder que les N e´le´ments,note´s (X1, . . . , XN) pour simplifier, tel que ∀i = 1, . . . , N : Φ(Xi) > q. Il est clair quel’e´chantillon (X1, . . . , XN) a la proprie´te´ voulue, a` savoir que Xi ∼ L(X|Φ(X) > q).Ne´anmoins, la loi forte des grands nombres implique que la proportion d’inte´reˆt N/ntend presque suˆrement vers p lorsque n tend vers l’infini. Cette me´thode ne´cessite doncde simuler un e´chantillon initial de taille n ≈ N/p pour aboutir a` un e´chantillon de tailleN . De fait, lorsque p est tre`s petit, ceci devient nume´riquement de´raisonnable.Supposons maintenant qu’on veuille estimer p, avec toujours p tre`s petit. La me´thodeMonte-Carlo classique consiste a` simuler un e´chantillon i.i.d. (X1, . . . , XN) selon la loi µ,et a` conside´rer l’estimateurpˆmc =NˆmcN=Card{i ∈ {1, . . . , N} : Φ(Xi) > q}N.Puisque Nˆmc suit une loi binomiale de parame`tres (N, p), il est clair que pˆmc est non biaise´et de variance relative :V (pˆmc)p2=1− pNp≈ 1Np.Par conse´quent, la taille N de l’e´chantillon a` simuler doit eˆtre au moins de l’ordre de 1/ppour atteindre une pre´cision acceptable, ce qui devient re´dhibitoire lorsque l’e´ve´nementd’inte´reˆt est tre`s rare.2Face a` cette situation, une astuce classique consiste a` effectuer un e´chantillonnage pre´fe´rentiel(ou Importance Sampling), lequel revient a` simuler l’e´chantillon (X1, . . . , XN) selon uneautre loi ν, pour laquelle l’e´ve´nement R = {Φ(X) > q} n’est plus rare. On estime alorsla probabilite´ comme ci-dessus et on la corrige graˆce au lien entre ν et µ. Ceci supposene´anmoins que l’on ait de l’information a priori sur le mode`le afin de proposer une nou-velle loi pertinente ν d’e´chantillonnage. Mais puisqu’on suppose ne rien connaˆıtre sur Φ(qui est une boˆıte noire), cette me´thode est ici exclue.Supposons maintenant le proble`me inverse : l’estimation d’un quantile extreˆme. On sup-pose p tre`s petit fixe´ le but est de trouver le quantile q tel que P(Φ(X) > q) = p.La me´thode Monte-Carlo classique marche comme suit : simuler un e´chantillon i.i.d.(X1, . . . , XN) selon la loi µ, les trier de sorte que Φ(X(1)) < · · · < Φ(X(N)) et conside´rerl’estimateur :qˆmc = Φ(X(⌊(1−p)N⌋)),ou` ⌊y⌋ de´signe la partie entie`re de y. Si la fontion de re´partition F (y) = P(Φ(X) ≤ y)admet une de´rive´e non nulle au point q, alors on peut montrer (voir par exemple Arnoldet al. [1], p.128) que le biais est en O(1/N) et l’e´cart quadratique de la forme suivante :Var(qˆmc) =p(1− p)N + 2· 1f 2(q)+ o(1/N) ≈ 1N· pf 2(q).2 L’algorithmePlutoˆt que la me´thode Monte-Carlo classique, nous proposons l’algorithme suivant :– Simuler un e´chantillon i.i.d. (X1, X2, . . . , XN) selon µ et initialiserX11 = X1, . . . , X1N = XN .– Pour m = 1, 2, . . ., noterLm = min(Φ(Xm1 ), . . . ,Φ(XmN )),et de´finirXm+1j ={Xmj si Φ(Xmj ) > LmX⋆ ∼ L(X|Φ(X) > Lm) si Φ(Xmj ) = Lm.– En notantM = max{m : Lm ≤ q},ce maximum e´tant par convention e´gal a` 0 si L1 > q, nous verrons que LM ≈ q avecune tre`s bonne pre´cision. Ainsi l’e´chantillon(X1, . . . , XN) = (XM+11 , . . . , XM+1N )est i.i.d. de loi approximativement e´gale a` L(X|Φ(X) > q).3Le point de´licat de l’algorithme ci-dessus est de re´ussir a` simuler une nouvelle particuleX⋆ distribue´e selon la loi L(X|Φ(X) > Lm). L’ide´e est d’utiliser un noyau de transitionK qui soit µ−re´versible, c’est-a`-dire tel que :∀(x, x′) ∈ Rd ×Rd µ(dx)K(x, dx′) = µ(dx′)K(x′, dx).En particulier, la loi µ est alors stationnaire pour le noyau K. Partant d’un noyau instru-mental q suppose´ µ−irre´ductible, l’algorithme de Metropolis-Hastings (voir [4], [3] et [5])permet pre´cise´ment d’en de´duire un noyau K ayant cette proprie´te´.A chaque e´tape m, notons alors Am = {x ∈ Rd : Φ(x) > Lm}, µm la restrictionnormalise´e de µ a` Am, c’est-a`-dire µm(dx) =1µ(Am)1Am(x)µ(dx), et conside´rons le noyaude transition Km de´fini par :Km(x, dx′) = 1Acm(x) δx(dx′) + 1Am(x)(K(x, dx′)1Am(x′) +K(x,Acm) δx(dx′)).Le principe de construction de ce noyau est tre`s simple : partant de x > Lm, une tran-sition x  x′ est propose´e par le noyau K. Si Φ(x′) > Lm, la transition est accepte´e,sinon elle est rejete´e et x reste a` la meˆme place. Il est alors facile de voir que la loi µm eststationnaire pour le noyau de transition Km.Ainsi, a` l’e´tapem, le but est de simuler X⋆ ∼ L(X|Φ(X) > Lm), avec X⋆ inde´pendant des(N − 1) particules Xmj ve´rifiant de´ja` Φ(Xmj ) > Lm. L’ide´e est tout simplement d’appelerV0 un vecteur ale´atoire correspondant a` l’une de ces particules (tire´e au hasard) et deconside´rer la chaˆıne de Markov (Vn)n≥0 de noyau de transitionKm. Pour n « assez grand »,la particule X⋆ = Vn est distribue´e suivant la loi L(X|Φ(X) > Lm) et est inde´pendantedes autres particules.3 Applications3.1 Estimation d’e´ve´nement rareAvec les meˆmes notations que pre´ce´demment, supposons qu’on veuille estimer la proba-bilite´ p = P(Φ(X) > q). Il suffit alors de conside´rer l’estimateur :pˆ =(1− 1N)M.On montre que la variable ale´atoireM suit une loi de Poisson, ce qui permet d’en de´duirela loi de pˆ ainsi qu’une expression simple pour tous ses moments.4Proposition 3.1 La variable ale´atoire M suit une loi de Poisson de parame`tre −N log p.Il s’ensuit que pour tout k ∈ N :E[pˆk] = pN(1−(1−1/N)k).Ainsi pˆ est un estimateur sans biais de p, de variance :Var(pˆ) = p2(p−1N − 1).Lorsque p est tre`s petit, pˆ est donc un estimateur bien plus pre´cis que celui obtenu parl’algorithme de Monte-carlo na¨ıf :V (pˆ)p2≈ − log pN≪ 1− pNp=Var(pˆmc)p2.D’autre part, puisque le nombre N de particules est au moins de l’ordre de la centaine,l’approximation de la loi de Poisson P(−N log p) par la loi normaleN (−N log p,−N log p)est le´gitime. Un intervalle de confiance de niveau 95% pour p est ainsi donne´ par :pˆe−2√− log pˆN ≤ p ≤ pˆe+2√− log pˆN .3.2 Estimation de quantile extreˆmeRe´ciproquement, une probabilite´ p tre`s petite e´tant fixe´e, supposons que l’on cherche a`de´terminer le quantile q tel que P(Φ(X) > q) = p. Avec l’algorithme propose´, il suffit denoter :m =⌈log(p)log(1−N−1)⌉≈ ⌈−N log p⌉,et d’estimer le quantile q par qˆ = Lm. On montre que le biais de cet estimateur est enO(1/N) et que sa variance est controˆlable.Proposition 3.2 Si la variable Φ(X) admet une densite´ f(q) non nulle en q, alors l’er-reur quadratique est asymptotiquement borne´e comme suit :p2(− log p− 2)f(q)2≤ limN→∞N ·E [(qˆ − q)2] ≤ p2(− log p+ 2)f(q)2.Ici encore, l’estimation propose´e par qˆ est bien plus pre´cise que celle offerte par l’estimateurMonte-Carlo classique qˆmc, puisque :Var(qˆ) ≈ 1N· −p2 log pf(q)2≪ Var(qˆmc) ≈ 1N· pf 2(q).5Comme pour la comparaison entre pˆ et pˆmc, ceci montre que pour atteindre la meˆmepre´cision, l’estimateur Monte-Carlo classique ne´cessite environ (−p log p)−1 fois plus departicules que la me´thode propose´e dans ce papier.Notons au passage que l’e´cart quadratique de notre estimateur fait intervenir la valeurde la densite´ f au point q, laquelle est bien suˆr inconnue. Ceci n’empeˆche ne´anmoinsnullement la construction d’intervalles de confiance pour qˆ.Proposition 3.3 Notons :{m+ =⌈−N log p+ 2√−N log p⌉m− =⌊−N log p− 2√−N log p⌋Alors l’intervalle I1−α(q) = [Lm− , Lm+ ] est un intervalle de confiance a` 95% pour le quan-tile q.4 Comparaison et illustrationLors de la pre´sentation, la me´thode propose´e ici sera compare´e a` l’e´tat de l’art en estima-tion d’e´ve´nements rares (notamment au re´cent travail de Ce´rou et al. [2]) et illustre´e surune application en tatouage nume´rique (watermarking).Re´fe´rences[1] B.C. Arnold, N. Balakrishnan, and H.N. Nagaraja. A first course in order statistics.Wiley Series in Probability and Mathematical Statistics : Probability and Mathe-matical Statistics. John Wiley & Sons Inc., New York, 1992. A Wiley-IntersciencePublication.[2] F. Ce´rou, P. Del Moral, T. Furon, and A. Guyader. Rare event simulation for a staticdistribution. Technical Report RR-6792, Inria, January 2009. Submitted.[3] W.K. Hastings. Monte Carlo sampling methods using Markov chains and their appli-cations. Biometrika, 57(1) :97–109, April 1970.[4] N. Metropolis and S. Ulam. The Monte Carlo method. J. Amer. Statist. Assoc.,44 :335–341, 1949.[5] C.P. Robert and G. Casella. Monte Carlo statistical methods. Springer Texts inStatistics. Springer-Verlag, New York, 1999.6

https://hal.inria.fr/inria-00494788/file/p137.pdf

Estimation de quantiles extrêmes et de probabilités d'événements rares

Abstract

Similar works

Full text

Available Versions

INRIA a CCSD electronic archive server

HAL-Rennes 1