[SON] Bienvenue dans cette séance où nous allons introduire la méthode du rejet. Dans une séance précédente, nous avons vu comment simuler une v.a. réelle en inversant sa fonction de répartition, et malheureusement, il y a de nombreuses situations où l'on ne peut pas inverser la fonction de répartition, et en fait, même la répartition de répartition en elle-même est trop compliquée pour être utilisée. La méthode du rejet précisément, est très importante parce qu'elle évite cet inconvénient, elle évite d'avoir à utiliser la fonction de répartition de la loi qui nous intéresse. Je vous présenterai cette méthode en deux temps. Dans un premier temps on verra cette méthode dans un contexte simple, et ensuite on verra la version générale de la méthode du rejet. Quelle est la question que l'on se pose à la base? On suppose qu'on a une v.a. X, qui a une densité, F(x), on suppose qu'elle est connue. Je vous rappelle qu'une densité, c'est une fonction qui prend des valeurs positives ou nulles, et dont l'intégrale vaut 1. Autrement dit, l'aire sous la courbe de f vaut 1. Donc la version la plus simple de cette méthode, qui aide à comprendre comment elle fonctionne, c'est la suivante : on va faire deux hypothèses. La première c'est que f est à support compact, cela signifie que la fonction f est nulle en dehors d'un intervalle [a, b] et on va supposer qu'il existe un majorant, on va l'appeler M de f(x). Autrement dit, si l'on prend n'importe quel point dans l'intervalle [a, b], f(x) va être plus petit ou égal à M. On peut faire un petit dessin pour résumer cette situation. Donc j'ai pris un exemple comme ça, arbitraire, où l'on a une fonction f, on a son graphe, on a a et b, donc en-dehors, cette fonction est nulle, et M est là, cela correspond à l'ordonnée du sommet. Le maximum du graphe. On a la surface du rectangle a b c d qui est M fois (b- a) la surface sous la courbe, comme on l'a dit c'est une densité, c'est la surface qui est sous cette courbe, là, ça vaut 1, et la zone hachurée là qui est en grisé, c'est exactement F(x0), où F est la fonction de répartition. Donc je vous remets juste en-dessous. Si j'intègre de- l'infini à x0 f l'intégrale commence en fait à a, puisque f est nul en desors de l'intervalle [a, b] et ça donne F(x0) moins F(a) et F(a) est nulle. Donc on trouve bien que la surface ici c'est F(x0). Quelle est la méthode, donc, dans ce cadre particulier? On tire un point que je vais appeler A de coordonnées (x, y) et ça va être un point uniformément distribué sur le rectangle a b c d. Donc qu'est-ce qu'on fait, plus précisément? On tire une v.a. uniforme U1 dans l'intervalle [0, 1], et X on le construit, il faut juste décaller et multiplier U1 par ce qu'il faut, donc c'est (b- a) U1 et on ajoute a. Et on tire indépendamment de U1, une autre v.a. uniforme sur [0, 1] et y ça va être M fois U2. Ensuite, on teste, on vérifie que Y est plus petit ou égal à f(X), et si c'est le cas on garde X, on l'accepte. Sinon, on le rejette et on tire un autre point. Et biensûr, si l'on n'a pas réussi, à cette étape, à trouver, on recommence jusqu'à la première fois que l'on trouve que cette condition est satisfaite. Ce que l'on peut démontrer, que l'on va faire, qui est assez simple, c'est que la variable X qu'on a ainsi obtenue, donc la première fois qu'on a satisfait à la condition en opérant cet algorithme, cette v.a. a bien pour fonction de répartition F. Alors on va le démontrer. Je vous refais le petit dessin ici, c'est pratique. Donc on a la courbe qui est comme ça, on va de a à b, M est ici. Nous avons notre rectangle a b c d, et la zone hachurée que l'on avait avant, elle est ici. Ici c'est le point x0. Donc ça c'est le graphe de F. Ce qui nous intéresse donc, c'est calculer la probabilité que X soit plus petit ou égal à x0 et on va voir que c'est bien F(x0), que c'est bien la fonction de répartition qui nous intéresse. Qui est donc bien que l'aire ici hachurée. Alors, la première remarque, c'est que par définition du point A qu'on tire dans ce rectangle, la probabilité que X soit plus petit que x0, sachant que l'on garde A, donc le point de coordonnées (X, Y) ces deux probabilités sont égales, puisque le fait que l'on garde le point A on a tiré indépendamment X et Y, donc le fait que Y soit en dessus ou en-dessous, n'influe pas sur le fait que X soit plus petit que x0. Donc on a bien égalité de cette probabilité, et de cette probabilité conditionnelle. L'idée maintenant, j'écris juste, l'algorithme vous dit que si on le garde c'est que Y est plus petit que f(X) et maintenant j'écris juste que c'est la probabilité conditionnelle c'est la probabilité que X soit plus petit que x0 et que Y soit plus petit que f(X), divisé par la probabilité que Y soit plus petit que f(X). Là j'ai juste utilisé la définition de la probabilité conditionnelle. Et je vais calculer séparément le numérateur et le dénominateur. Donc le numérateur, [AUDIO_VIDE] n'oubliez pas qu'on tire des variables uniformes, donc on a tout simplement que c'est la surface hachurée, qui est exactement sur le dessin, divisée par la surface du rectangle a b c d. [AUDIO_VIDE] Et donc on trouve la surface de la zone hachurée, c'est F(x0) et on divise par la surface du rectangle, qui est M (b- a). Calculons le dénominateur. Cette fois,-ci on doit diviser la surface sous la courbe, par la surface du rectangle. [AUDIO_VIDE] Et on trouve donc, la surface sous la courbe n'oublions pas que c'est une densité, donc elle vaut 1, et la surface du rectangle vaut M (b- a). Le résultat c'est que, P (X inf à x0) c'est égal au produit donc, de F(x0) / M(b- a) et donc de l'inverse de cette quantité, c'est-à-dire M(b- a) / 1 et on trouve bien F(x0). On a donc montré que, en suivant cet algorithme, donc on tire un point au hasard uniformément dans le rectangle, et si on tombe sous la courbe, on a bien une v.a. X, qui a pour fonction de répartition F. Donc le théorème qui est ici est bien établi. Je fais une petite remarque, sur cet algorithme, donc, un peu quantitative, la première chose qu'il faut remarquer, c'est que cet algorithme est aléatoire, puisque le temps de premier succès, qui représente le fait que la condition soit satisfaite, c'est-à-dire que l'on tombe bien sous la courbe, c'est une v.a. qu'on peut appeler N, et en fait si l'on pousse un peu plus le travail, on pourrait montrer que N est une loi géométrique, de paramètre 1/M. Nous allons maintenant généraliser la méthode que je viens de vous expliquer, il y a plusieurs raisons à cela, la première étant que si l'on utilise une variable uniforme, on va donc tirer dans un rectangle, et il y a des chances qu'on tombe très souvent entre le graphe du haut du rectangle et le graphe de la densité qui nous intéresse, donc il y a beaucoup de points qui vont être rejetés, ça c'est la première chose ; la deuxième chose c'est que, si la densité a un support qui n'est pas compacte, il est clair qu'on ne peut pas utiliser un tirage uniforme dans un rectangle, il faudrait un rectangle de longueur infinie, ce qui n'a absolument pas de sens, en fait, on va utiliser une idée assez jolie, qui est d'utiliser une v.a. auxiliaire, qu'on sait, elle, simuler. Alors plus précisément on va supposer qu'il existe une densité g et une constante a positive, ça on va se les donner nous-mêmes, tel que pour tout x, f(x) est plus petit que ag(x). Et g, qu'est-ce que ça va être? Cela va être la densité d'une v.a. que l'on sait cibler. Et typiquement, on utilise une v.a. qu'on sait simuler en inversant sa fonction de répartition. La première petite remarque que je voudrais faire c'est que, si vous intégrez cette inégalité, vous voyez que a doit être être plus grand ou égal à 1, puisque f et g sont des densités et leur intégrale doit valoir 1. Ce sont deux fonctions positives d'intégrale 1. Et en fait, vous voyez même que a > 1. En effet, si a = 1, cela veut dire que les deux densités coïncident. Ce qui fait que ce n'est pas du tout ce qui nous intéresse. On a envie justement d'utiliser cette fonction auxiliaire, cette densité auxiliaire, alors qu'on ne sait pas trop quoi faire de la densité f. Alors, voici l'algorithme du rejet, que je vais appeler algorithme du rejet généralisé. On va tirer un point aléatoire de la manière suivante. A qui qui aura pour coordonnées (X, Y). X va suivre une loi de répartition G, et que je l'ai dit, typiquement, va inverser la fonction de répartition en utilisant une variable aléatoire uniforme sur l'intervalle [0, 1]. Quant à Y, cela va être a * g (X) * U, où U c'est une variable aléatoire uniforme sur l'intervalle [0, 1], indépendante de U tilde. Maintenant, on teste si f (X) est plus grand ou égal à Y. On garde X si c'est le cas, sinon on tire un autre point, de la même manière. On va le faire bien sûr, à chaque fois, indépendamment de ce qu'on a fait précédemment. L'idée, c'est qu'on réitère jusqu'à ce que cette condition, f (X) plus grand ou égal que Y, soit remplie. C'est vraiment la généralisation de ce qu'on a fait précédemment, en utilisant cette densité auxiliaire. Tout de suite, ce qu'on peut remarquer, c'est que le temps qu'il faut pour satisfaire cette condition, donc le nombre de fois qu'on va faire la boucle de cet algorithme, c'est lui-même une variable aléatoire. Et on verra plus tard ce qu'on peut en dire. Je vais vous illustrer l'algorithme et la méthode avec la loi normale. La loi normale, elle a une fonction de répartition qui est compliquée et qui n'est pas connue explicitement. C'est une formule qui n'est pas pratique du tout. Et, encore moins son inverse, sa fonction réciproque. Alors, je vous rappelle, si X suit, on va prendre une loi normale, centrée et réduite, cela veut dire qu'on a une densité, qui est donnée par f (x) = exp (- x 2 / 2) / racine de (2 * pi). Ce que je dis, c'est qu'en fait on va pouvoir majorer cette densité par exp ( 1 / 2- valeur absolue de x ) / racine de (2 * pi). Cela, cela va marcher pour tout x dans R. Alors pourquoi? Le fait qu'on ait cette inégalité, c'est équivalent au fait qu'on peut simplifier les racines de (2 * pi). On peut mettre tout dans une exponentielle d'un côté. Et cela revient à dire qu'une exponentielle, exponentielle de cette quantité est plus grande ou égale à 1. Autrement dit, que ce qu'il y a dans l'exponentielle est plus grand ou égal à 0. Et là, on reconnaît le développement d'un carré, qui est manifestement positif ou nul. Autrement dit, on a bien cette inégalité. Alors, on prend ce terme, et on voit qu'on peut le mettre sous la forme racine de ( 2 * e / pi ) * 1 / 2 * e (- valeur absolue de x). On peut encore l'arranger, en remarquant que cette quantité n'est rien d'autre que la densité de la loi de Laplace standard, qu'on a vu dans une séance précédente. Autrement dit, on a bien la densité de la loi normale centrée réduite, qui est dominée, qui est majorée par a * g (x), où g (x), c'est la densité d'une loi de Laplace standard, et a vaut racine de ( 2 * e / pi ), ce qui est à peu près 1,32. Et dans une séance précédente, on a vu que si on prenait une variable uniforme dans l'intervalle [-0,5, 0,5], et qu'on calculait moins son signe fois le logarithme de (1- 2 * valeur absolue de U), on avait comme cela, une variable aléatoire qui suit une loi de Laplace standard. Voilà comment on peut par exemple appliquer la méthode du rejet, au cas de la loi normale. Donc, la constante a, on a envie de la prendre la plus petite possible. L'idée étant que plus elle est petite, moins on va tirer de points entre le graphe de la densité g et celui de la densité f. Donc, là j'ai mis le graphe de Laplace multiplié par à peu près 1,8. Dessous, il y a le graphe de la densité gaussienne. Et, si on abaisse la courbe rouge, qui revient à diminuer a, la valeur qu'on a trouvé tout à l'heure, en fait elle est optimale. On peut le montrer, dans le sens où si on descend plus bas, on se retrouve à voir la courbe rouge qui passe par endroits sous la courbe verte, sous la courbe de la loi de la densité gaussienne. Donc, on verra dans une séance ultérieure comment on peut manipuler cela avec une expérience numérique interactive. Pour ceux parmi vous qui sont les plus courageux, on peut se demander si l'algorithme qu'on a introduit précédemment est bien justifié mathématiquement, si on peut vraiment montrer qu'il a un sens. Et pour ce faire, on va formaliser un peu plus ce que j'ai dit précédemment. On va supposer qu'on a une suite de variables aléatoires X n. Chaque X n va être l'inverse par la fonction de répartition G d'une variable aléatoire uniforme U tilde n. Y n va être de la forme a * g (X n) * U n. n va représenter l'itération, l'itérée n de notre boucle dans notre algorithme. Que sont ces variables U tilde n et U n? D'une part, cela c'est une suite de variables aléatoires indépendantes identiquement distribuées, de loi uniforme sur [0, 1]. Pareil pour U n, et en plus, on va supposer que ces deux suites de variables uniformes sont indépendantes entre elles. Maintenant, il faut introduire cette variable aléatoire grand N, qui va être le premier instant où f (X n) va être plus grand ou égal à Y n, ce qui est la condition qu'on se donnait dans l'algorithme de rejet. Et que va être la variable aléatoire acceptée, celle qui va effectivement nous donner la bonne densité, d'après l'algorithme, si on le croit? On va l'appeler X chapeau. C'est égal à X n, donc c'est égal à, si grand N = petit n. Donc, on regarde la première fois que la condition est satisfaite, en on obtient pour ce n là, un X n qui va donner la variable aléatoire qui a la bonne densité. Cela, c'est ce que prétend l'algorithme. En fait, on a le théorème suivant. La première chose c'est que la variable aléatoire N suit une loi géométrique de paramètre 1 / a, en particulier donc l'espérance vaut a. Et la variable aléatoire X chapeau suit la, effectivement, la loi de densité f. Et en plus, on a que ces deux variables aléatoires sont indépendantes. Nous allons maintenant démontrer ce théorème. On va commencer par introduire des événements E n, que je vais définir de la façon suivante. Ce sont les événements qui sont la condition que f (X n) est plus grand ou égal à Y n. Et n va valoir 1, 2, etc. Alors, ces événements, par construction même, ils sont indépendants. et ils ont la même probabilité. Et quelle est-elle? Je vais l'appeler alpha. Pour l'instant, on ne connaît pas ce nombre, et on va le déterminer. A cause du fait qu'on a que des variables uniquement distribuées, alpha, ce nombre donne la probabilité de tous les ensembles, les événements E n. alpha, on va le déterminer tout à l'heure. Et en particulier, on constate que la probabilité que grand N, le premier temps auquel cette condition est satisfaite, cette variable grand N > n, c'est (1- alpha) puissance n. Donc, on reconnaît là une loi géométrique. Cela, c'est la première étape. La deuxième étape, c'est de considérer une fonction continue bornée quelconque. Et on va utiliser le théorème de, qui permet de calculer, de caractériser une loi. C'est le théorème qui se trouve dans la séance 6 du cours 3. On va exprimer l'espérance de h (X chapeau), de telle manière qu'on puisse interpréter sa loi. Donc, là on est est en train de calculer cette espérance. L'idée, c'est de la décomposer sur toutes les valeurs possibles qui sont données par les événements E n. On peut écrire cette espérance comme la somme de n = 1 à l'infini, de l'espérance de h (X n) * l'indicatrice de E n * l'indicatrice de (N > n- 1). Là, on a simplement décomposé ce qui se passe pour X chapeau. Étant donné que a priori, grand N peut prendre n'importe quelle valeur, et que la première fois qu'il atteint une valeur n, cela correspond à un événement E n, qui est le fait que la condition f (X n) plus grand que Y n est satisfaite. L'observation clé, c'est que cette variable aléatoire est indépendante de cette variable aléatoire là. Le fait que la condition soit satisfaite à un n donné, qui est cette inégalité f (X n) plus grand que Y n, ne dépend pas du moment auquel cela a lieu. C'est dû aux hypothèses d'indépendance qu'on a faites. Donc, on peut factoriser l'espérance. On se retrouve avec le produit de l'espérance, d'espérance de X n * indicatrice de E n * l'espérance de l'indicatrice que N > (n- 1). Et ici, vu ce qu'on a fait au-dessus, cela ce n'est rien d'autre que (1- alpha) puissance (n- 1), puisque l'espérance de cette indicatrice, c'est la probabilité que N > n- 1, et c'est donc (1- alpha) puissance (n- 1). Maintenant, calculons l'espérance de h (X n) * l'espérance de l'indicatrice de E n. Donc, cela s'écrit comme l'intégrale de h (z) * g (z), g (z), donc c'est la densité auxiliaire, fois l'intégrale de 0 à 1, indicatrice que f (z) est plus grand ou égal à a * g (z) * x. On intègre par rapport à x, et on intègre. Donc ici, on intègre par rapport à x qui suit une loi uniforme sur 0, 1. On intègre l'indicatrice de cet événement, qui, puisque tous les événements E n, leurs probabilités, leurs espérances ne dépendent pas de n, à cause du fait qu'on utilise des variables identiquement distribuées. Et ensuite, ce résultat dépend lui-même encore de z, et il faut encore intégrer contre la densité g. Donc, maintenant, cela se calcule explicitement, cela. Ici, vous voyez qu'on a intégré x de 0 jusqu'à f (z) / g (z). Ce qui donne exactement, f (z) / (a * g (z)) d z. Il y a donc une simplification. Et on obtient 1 / a, qui est ici, * intégrale de h (z) * f (z) d z. On remarque en particulier que si h = 1, parce que cela c'est vrai pour tout fonction continue ou bornée, on peut prendre la fonction qui vaut constamment 1, alors alpha c'est donc l'espérance de l'indicatrice de E n. C'est donc la probabilité de E n. C'est 1 / a, d'après la formule qu'on vient d'avoir, * f (z) d z, et comme f (z) d z, l'intégrale cela vaut 1, on en déduit que alpha = 1 / a. On a donc identifié alpha, le paramètre de la loi géométrique pour n. En particulier, notez que 0 < alpha < 1. Il nous reste la dernière étape. On reprend l'identité précédente : espérance de h (X chapeau). On a vu que cela s'écrivait somme de n = 1 jusqu'à l'infini, espérance de h (X n) * indicatrice de E n * (1- alpha) puissance (n- 1), c'est ce qu'on a obtenu il y a un petit instant. On vient de voir que cette partie est égale à alpha * l'intégrale de h (z) * f (z) d z, qu'on peut d'ailleurs sortir de la série, parce qu'elle ne dépend pas de n, cette expression. [AUDIO_VIDE] Et là, on voit qu'on obtient h (z) * f (z) d z et alpha / (1- (1- alpha)). Tout cela donne exactement 1. Donc, en vertu du théorème d'identification des lois d'une variable aléatoire, qui est le théorème qui était dans la séance 6 du cours 3, X chapeau a pour densité f. Je vous laisse vérifier pourquoi X chapeau et n sont bien des variables aléatoires indépendantes, et nous avons terminé donc la démonstration du théorème. Et donc, cela termine cette introduction à la méthode du rejet.
