Simuler une loi uniforme sur une surface — Les-mathematiques.net The most powerful custom community solution in the world

Simuler une loi uniforme sur une surface

ArgeuhArgeuh
dans Probabilités, théorie de la mesure
Bonjour,

Tout à l'heure on m'a posé une question et je n'ai su y répondre:

Si on se donne $f$ une fonction suffisamment réguliere, comment simuler une loi uniforme sur l'ensemble $\{(x,f(x)), x\in D\}$ où $D$ est un ouvert borné ?

J'aimerais beaucoup avoir une réponse.
Merci.
«1

Réponses

  • alekkalekk
    cela se fait souvent avec la methode du rejet:
    http://fr.wikipedia.org/wiki/Méthode_du_rejet
  • alekkalekk
    un exo que tout debutant en programmaion/probabilite devrait faire:

    un etudiant veut simuler une loi uniforme sur le disque unite. Il suit le raisonnement suivant:
    1: il simule un rayon $r$ uniformement sur $[0;1]$
    2: il simule un angle $\theta$ uniformement sur $[0; 2 \pi[$
    3: bein alors $(r \cos(\theta), \sin(\theta))$ est la simulation qu'il voulait pardi ..

    Pourquoi cet etudiant a t il completement faux ..
  • egoroffegoroff
    Le problème dans le cas de Argeuh est que son ensemble est de mesure nulle dans $\R^{n+1}$, si je comprends bien l'énoncé... Je commencerais par exprimer la loi de $X$, si $(X,f(X))$ suit la loi uniforme sur $\{(x,f(x)\}$, a priori elle sera plus concentrée là où le gradient de $f$ est grand. Ensuite on essaye de simuler $X$ suivant cette loi.
  • ArgeuhArgeuh
    Et oui, tout le problème est que cet ensemble est de mesure nulle ...
  • alekkalekk
    ah oui, effectivement. Il faudrait alors définir par ce que tu entends par "uniformément" sur cet ensemble de mesure nulle ..
  • ArgeuhArgeuh
    Puisque la surface $S$ est "un objet de dimension 2", on peut dire que $X$ est uniforme sur $S$ si et seulement si $P(X \in A) = \dfrac{|A|}{|S|}$ pour tout $A\subset S$ et où $|.|$ désigne la mesure de Lebesgue en dimension 2 ? (je sens un gros problème en écrivant ces lignes, mais disons, le uniforme est au sens intuitif du terme) (et cette question a été posée en exercice dans un TD donc ça ne doit pas être très compliqué ... enfin en théorie)
  • egoroffegoroff
    Non pas de gros problème, enfin il faut savoir de quoi on parle. Il y a une arme atomique qui permet de parler d'aire d'a peu près n'importe quelle partie de $\R^3$, c'est la mesure de Hausdorff 2-dimensionnelle. Dans le cas d'une sous-variété régulière pas la peine de sortir ce bouzin, on sait définir une mesure de surface "unitaire", qui coïncide avec la mesure de Lebesgue lorsque la variété est incluse dans un sous-espace affine. Dans ton cas, la mesure unitaire sur le graphe de $f$ donne par projection une mesure sur $D$, sauf erreur elle est de densité $\sqrt{1+\left \frac{\partial f}{\partial x} \right^2} \sqrt{1+\left \frac{\partial f}{\partial y} \right^2}$ par rapport à la mesure de Lebesgue (il manque peut-être un produit croisé), et pour obtenir une proba il suffit de diviser par la masse totale ; cela dit dans la méthode du rejet pas besoin de connaître cette masse.
  • alekkalekk
    oui, je pense que cela marche bien !

    sinon juste une question: une fonction $C^1$ admet un graphe de mesure de Hausdorff $2$ dimennsionnelle localement finie, mais pas forcement pourune fonction continue. Quand est il des fonctions derivables ?
  • egoroffegoroff
    Après réflexion je dirais que la densité est plutôt $\sqrt{1+\left( \frac{\partial f}{\partial x} \right)^2 +\left( \frac{\partial f}{\partial y} \right)^2}=\sqrt{1+||\nabla f||^2}$, ça paraît plus plausible vu l'invariance par rotation.
  • egoroffegoroff
    En effet l'autre marchait bien aussi pour les exemples simples ! du genre les "cylindres" $f(x,y)=g(x)$.

    Sinon voilà une excellente question, je vais y réfléchir :)
  • egoroffegoroff
    Bon j'ai écrit un petit programme pour tester le truc, est-ce que quelqu'un sait comment dessiner un point dans l'espace sous scilab ? En 2D je connais xrect mais ça ne marche pas avec 3 coordonnées.
  • Jean-Louis LambdaJean-Louis Lambda
    Est-ce que l'on peut construire une paramétrisation injective dérivable par morceaux de D ? Dans ce cas, on pourrait compléter le jacobien par des vecteurs orthonormaux, prendre l'inverse du déterminant, et faire une sorte de normalisation. Enfin, c'est un peu flou dans ma tête mais il doit y avoir un truc comme ça. Non, désolé, après avoir mieux lu ce que vous avez dit, je crois que ce que je viens de dire ne sert à rien.
  • Remarque0Remarque0
    egoroff : param3d1
  • egoroffegoroff
    Jean-Louis : En effet pour obtenir ma densité j'ai regardé l'image par la projection sur $D$ d'un parallèlogramme infinitésimal tangent à la surface construit sur deux vecteurs tangents.

    remarque : Merci bien je savais que je pourrais compter sur le pro du Scilab :) En attendant j'avais fait le truc avec Maple (beurk), le résultat ci-joint pour $D=[-1,1]^2$ et $f(x,y)=e^x+e^y$, 100 échantillons en utilisant la méthode du rejet pour simuler selon la densité voulue. Ca serait mieux si vous pouviez faire tourner la surface mais bon... Maintenant que j'ai la commande Scilab je vais pourvir faire des essais avec un domaine $D$ un peu plus rigolo.
    10476
    10477
  • Remarque0Remarque0
    A propos, on n'est pas obligé de se limiter à des graphes. On peut faire la même chose sur une nappe paramétrée. Ca peut être utile pour tirer uniformément sur une sphère par exemple (si j'ai bien compris le principe de la méthode de rejet). Comme je ne maîtrise pas l'aspect aléatoire de la chose, je m'abstiendrai de faire du scilab, et j'attendrai qu'egoroff nous livre le résultat de ses cogitations.
  • egoroffegoroff
    Tout à fait monseigneur ; il y a sans doute des formes fondamentales ou je ne sais quoi là-dessous, c'est à creuser, j'ai bien envie de simuler des points uniformément sur des surfaces un peu tordues genre immersion d'une bouteille de Klein ou que sais-je.

    Je vous présenterai bien les résultats de la simulation en Scilab mais mon PC est sur le point de rendre l'âme, je vais le laisser se reposer et call it a day fûr die Nacht. A demain !
  • ArgeuhArgeuh
    Waouh :)-D

    Mais, il n'y a pas de justification élémentaire ? Je suis largué là ...
  • egoroffegoroff
    Pas trop de justification élémentaire je pense, car il s'agit bien de géométrie différentielle ; mais si on en connaît un peu c'est tout simplement la formule de changement de variable avec jacobien, mais pour des "espaces courbes". M'enfin peut-être que remarque va trouver comment rendre ça plus élémentaire.
  • Remarque0Remarque0
    Ben, c'est évident, non ? C'est l'élément d'aire dans un système de coordonnées curvilignes. Dans le cas d'un graphe, ça donne la formule $\sqrt{1+\|\nabla f\|^2}$. Pour une nappe paramétrée générale $\Phi\colon\omega\to\R^3$ avec $\omega\subset\R^2$, on pose classiquement $a_\alpha=\partial_\alpha \Phi$, et la densité en question, c'est juste $\|a_1\wedge a_2\|$. C'est aussi la racine carrée du déterminant de la première forme fondamentale.
  • Remarque0Remarque0
    Et on attend toujours la bouteille de Klein remplie de petits points aléatoires en scilab et en 3d...
  • egoroffegoroff
    D'accord pour la densité pour le cas d'une nappe générale. Maintenant il ne me manque plus qu'une paramétrisation raisonnable de la bouteille de Klein, parce que pour l'instant je ne trouve que des trucs horribles qui me donnent le mal de mer :-( Je vais peut-être me rabattre sur une surface un poil plus simple...
  • egoroffegoroff
    Pour patienter en attendant la bouteille de Félix...
    10505
  • Remarque0Remarque0
    Aaaah ! Tout vient à point pour qui sait attendre ! C'est beau.
  • egoroffegoroff
    Oui désolé j'ai un peu tardé :)

    Le tore est paramétré par
    $x=(2+\cos v) \cos u$
    $y=(2+\cos v) \sin u$
    $z=\sin v$
    La densité de proba qu'on en déduit sur $[0,2\pi]^2$ est de la forme $f(u,v)=c \sqrt{4+4 \cos v + \cos^2 v}$.
  • samoksamok
    Bonsoir,

    je me pose une question, résultant du fait que le tore ou la bouteille de Klein sont des espèces de carré, dont les côtés sont identifiés entre eux d'une certaine façon.
    Donc tirer des points aléatoires de façon uniforme dans un carré me semble jouable, et si l'on souffle dans le carré (ça c'est hors de portée pour moi, informatiquement parlant), est-ce qu'au final ça va être uniformément distribué sur la surface obtenue?

    Merci de faire reculer mon ignorance infinie,
    S
  • egoroffegoroff
    Bonsoir samok,

    Qu'entends-tu par "souffler" ? :)
  • samoksamok
    Des sombres histoires de roméo-morphismes, sans doute, sourire.

    Pour le tore, je pense que tu vois bien ce que je veux dire, transformer le carré en cylindre (fini) puis incurver le tout comme il faut.

    Est-ce que cela transporte bien l'uniformité, c'est un peu ma question, j'ai envie de dire que ce serait trop simple pour être vrai.

    S
  • Remarque0Remarque0
    Non justement : on ne regarde pas une surface abstraite, mais un plongement particulier de cette surface dans $\R^3$ et on veut que la distribution apparaisse uniforme sur cette surface plongée, c'est-à-dire par rapport à la mesure superficielle, ie, un truc qui est modifié via la transformation que tu préconises. Je pense qu'egoroff pourrait, si ce n'est pas abuser, illustrer graphiquement la différence grâce à son script scilab magique.
  • egoroffegoroff
    Tout le problème est là. Lorsque je fabrique un cylindre, je peux conserver l'uniformité. Mais lorsque je veux en faire un tore, mon carré de caoutchouc n'est plus uniformément tendu. Donc des points uniformes sur un carré ne le seront plus sur le tore.

    L'idée est alors de tirer de manière inhomogène sur le carré, cette inhomogénéité étant choisie de sorte que la distribution qui en résulte sur le tore soit uniforme. Pour une paramétrisation donnée $(u,v) \mapsto (x,y,z)$, on regarde quelle est l'aire (infinitésimale) de l'image sur le tore d'un petit rectangle de côtés $(du,dv)$ ; appellons-la $A(u,v)\, du \,dv$.

    Si on tire un point du carré selon la densité de probabilité $f$, la proba de tomber dans notre carré $(du,dv)$ est $f(u,v) \, du \, dv$, on a la même proba de tomber dans son image sur le tore (si la paramétrisation est injective). Or pour avoir une loi uniforme digne de ce nom cette proba doit être proportionnelle à l'aire de ce petit bout de tore : on doit donc avoir $f(u,v) \, du \, dv = cA(u,v) \, du \, dv$.

    On en déduit que $f=cA$, avec $c=(\int A)^{-1}$ de sorte qu'on ait bien une densité. Reste à calculer $A(u,v)$ ; pour cela on fait le produit vectoriel des vecteurs tangents $t_1=\left( \frac{\partial x}{\partial u}, \frac{\partial y}{\partial u}, \frac{\partial z}{\partial u} \right)$ et $t_2=\left( \frac{\partial x}{\partial v}, \frac{\partial y}{\partial v}, \frac{\partial z}{\partial v} \right)$ et on en prend la norme.

    Reste à simuler des points $(u,v)$ de densité $f$. On utilise la méthode du rejet, qui évite d'avoir à calculer $c$. Supposons que $A$ soit bornée, disons par $M$. Il s'agit de tirer des $(u,v)$ uniformes dans le carré, et un $w$ uniforme dans $[0,M]$. On garde le point $(u,v)$ dès que $w \leq A(u,v)$, sinon on recommence. On montre que ce procédé produit en sortie des points ayant la bonne distribution : intuitivement un point a d'autant plus de chance d'être gardé que $A(u,v)$ est grand.
  • egoroffegoroff
    Aucun problème, aussitôt dit aussitôt fait. Le problème, c'est que la différence n'est pas franchement flagrante.
  • Remarque0Remarque0
    Peut-être avec un tore dont le petit rayon est proche du grand rayon ?
  • egoroffegoroff
    Oups désolé j'ai oublié d'attacher le fichier. Bon en comparant les deux on voit quand même une concentration plus forte au centre et plus faible à l'extérieur.
    10509
  • egoroffegoroff
    Alors, pour un tore $(R,r)$ la densité est proportionnelle à $A(u,v)=r (R+r \cos v)$ (merci Maple), donc la densité au centre est $\frac{R-r}{R+r}$ fois celle à l'extérieur. Effectivement avec un $r$ proche de $R$ (disons $0,9R$) on accentue la différence, mais on ne dépassera jamais 2.
  • egoroffegoroff
    Euh.. n'importe quoi ! Bon je fais tout de suite un dessin illustrant la différence pour $r=0,9R$
  • egoroffegoroff
    Bon alors voilà le dessin tant attendu, pour r=0,9R, et 60000 points (contre 2000 sur les dessins précédents). Je ne suis qu'à moitié satisfait, pour deux raisons :
    - les points semblent de répartir par "bandes" : défaut du générateur de nombres aléatoires ? de mon script ? bug d'affichage ?
    - plus grave : les points qui devraient être uniformes semblent légèrement resserres au centre malgré tout ; erreur de calcul de la densité ? effet d'optique ? je suis parano ?

    Assez parlé, voici les preuves :
    10516
  • egoroffegoroff
    En taille réelle, retourné pour ne pas exploser la largeur du forum :
    10517
  • Remarque0Remarque0
    Oui, les bandes sont bizarres... balance le script, pour voir (s'il n'est pas couvert par le secret industriel).
  • egoroffegoroff
    J'ai toujours un peu honte de montrer ce que je code à quelqu'un d'autre, parce que c'est rarement très élégant (voire un peu crade parfois) mais bon...

    Il faut changer l'extension, de .txt à .sci ; les deux dernieres fonctions tirages_uni et tirages_samok sont celles qui marchent ; la valeur de retour est toujours 0. Si tu as des améliorations à proposer tu es le bienvenu !

    Disclaimer : The Egoroff Foundation should not be held liable for any problems caused by this script.
  • Remarque0Remarque0
    Il y a définitivement des bandes, et le plus curieux c'est qu'elles semblent correspondre aux lignes de coordonnées sur le tore... j'essaie d'y voir plus clair.
  • Remarque0Remarque0
    Ah, pourrait-ce être un problème de discrétisation de la paramétrisation du tore ? (je n'ai pas encore mis le nez dans ton admirable script (c) egoroff).
  • Remarque0Remarque0
    Bof, non ça n'a pas l'air d'être ça...:-(
  • Remarque0Remarque0
    C'est quand même marrant, parce que quand on tire les points à plat, dans le plan u,v, on ne voit pas les bandes... Bon, je vais laisser refroidir un peu.
  • Remarque0Remarque0
    OK, je crois que c'est un problème graphique : le tracé du tore cache une partie des points :

    10522
  • egoroffegoroff
    ..et pour aller dans le même sens, voilà l'image sur le tore d'un maillage régulier 200x200 du carré.
    10523
    reg.jpg 194.6K
  • egoroffegoroff
    Ah ben voilà ! Je savais que tu trouverais la solution. Tu nous mets la même en version samok ?
  • Remarque0Remarque0
    Et la comparaison uniforme/samok :

    uniforme :

    10524
    10525
  • egoroffegoroff
    (tu)
  • Remarque0Remarque0
    (tu)(tu):)-D
  • egoroffegoroff
    Prochaine étape : tirer les points de sorte qu'on ait une impression d'uniformité :
    - lorsque le tore est vu du dessus, ou de côté, ou
    - lorsqu'il est vu d'un angle $\theta$ par rapport au plan horizontal (à distance infinie, ou à distance finie $D$).
    Quelqu'un est volontaire ?
  • Remarque0Remarque0
    Non. Mais n'est-ce pas une application triviale de ce qui précède en composant par une autre application ?
Connectez-vous ou Inscrivez-vous pour répondre.
Success message!