Loi uniforme dans un tore

Saturne · July 2018

Bonjour,

Un tore est obtenu comme solide de révolution d'un cercle autour d'un axe.
Je pioche un angle au hasard, uniformément. J’obtiens le cercle correspondant. Je pioche un point dans ce cercle uniformément. J'espérais ainsi obtenir une loi uniforme dans le tore, mais ce n'est pas le cas. Auriez-vous une idée du pourquoi ?

ev · July 2018

Est-ce que ça aide ?

e.v.

Saturne · July 2018

Hi ev,

Non. Dans un tore, pas sur la surface.

aléa · July 2018

Bonsoir,

Questions:

1) révolution d'un cercle ou d'un disque ?
2) Qu'est-ce qui t'amène à penser que ta simulation ne donne pas ce que tu veux ?

mojojojo · July 2018

Pourtant c'est le même type de phénomène que pour la surface d'un tore.

En fait c'est parce qu'il faut compenser par rapport à la distance à l'axe de rotation. Par exemple imagine qu'au lieu d'avoir un tore plein tu en ais 2 concentriques de même section. Avec ta méthode tu choisis un angle uniformément, ça te donnes deux disques de même aire puis tu choisis un point uniformément dans l'union de ces deux disques. De cette façon on a autant de chance de tomber dans un tore ou l'autre, pourtant le tore extérieur a un plus grand volume et le tore intérieur un plus petit volume.

Bref en coordonnées cylindriques l'élément infinitésimal de volume est $r\mathrm dr \mathrm d z \mathrm d \theta$ et donc tu dois mettre une densité de probabilité sur ton disque en $C/r$ où $r$ est la distance à l'axe de révolution et $C$ une constante pour bien arriver à une mesure de probabilité uniforme.

Edit : joli lien EV, plein de jolis dessins.

Saturne · July 2018

Je me réveille et je suis dans le coltard.

@Alea:

1) Un disque
2) J'ai évalué une intégrale avec ces simulations et je ne trouve pas le bon résultat.

@majojojo
Faut que je sorte du coltard, attends.

aléa · July 2018

J'approuve tout à fait ce que dit mojojojo.

Une autre manière de le dire, c'est que si un point suit la loi uniforme sur un ensemble de révolution, la loi de l'angle est bien la loi uniforme sur $[0,2\pi]$, mais la loi de la position dans la coupe n'est pas uniforme.
Pourquoi ?
Si on considère un sous-ensemble $K$ de l'objet de révolution $D$, la probabilité que le point uniforme sur l'ensemble de révolution ait une position de coupe dans $K$ n'est pas $\lambda^2(K)/\lambda^2(D)$, mais
$\lambda^3(r(K))/\lambda^3(r(D))$, où $r(A)$ est l'objet de révolution induit par $A$.

C'est là qu'intervient le théorème de Guldin : la mesure du volume engendré par la révolution d’un élément de surface plane autour d’un axe situé dans son plan et ne le coupant pas est égale au produit de l’aire de la surface par la longueur de la circonférence décrite par son centre de gravité.

Cette longueur de circonférence est précisément le $2\pi r$ évoqué par mojojojo.

Saturne · July 2018

Désolé mais j'ai pris des médocs pour dormir hier et je suis encore dans les vappes. Je regarderai plus tard. Mais en un mot c'est facile ou pas ? Merci pour votre aide en tout cas.
C'est bizarre, j'ai un autre algo qui marche par acceptation/rejet, et parfois il donne le même résultat que l'intégrale correspondante, parfois non.

aléa · July 2018

Ce n'est pas très facile.

En un mot tu ne dois pas simuler la loi uniforme sur le disque, mais la loi de densité proportionnelle à
$x1_D$, où $x$ est la distance du point $(x,y)$ à l'axe de révolution.
Mettons que $D$ est le disque de centre $(0,h)$ et de rayon $r<h$.

On peut calculer sans trop de difficulté la fonction de répartition de $X$. Reste à voir si ça s'inverse bien pour éviter les méthodes de rejet. Je regarde.

aléa · July 2018

Ah, la fonction de répartition est moche. Il va falloir faire du rejet.

Saturne · July 2018

Hmm ok, merci. Je vais m'en tenir à la méthode de rejet alors (qui n'est pas de moi). Ceci dit j'aimerais comprendre pourquoi elle ne marche pas à tous les coups...
C'est dingue comme une forme élémentaire comme le tore pose des difficultés.

aléa · July 2018

Simulation de la loi uniforme sur le tore obtenu par rotation du disque d'équation
$$\{z=0;(x-h)^2+y^2\le r^2\}$$ autour de l'axe des $y$.

On commence par prendre $(X,Y,T)$ indépendants suivant respectivement les lois uniforme sur $[h-r,h+r]$, $[-r,r]$, $[0,h+r]$.

Si $(X-h)^2+Y^2>r^2$ ou $T>X$, on rejette.

Sinon, on renvoie
$M_{\theta}\begin{pmatrix}X\\Y\\0\end{pmatrix}$, où $M_{\theta}=\begin{pmatrix}\cos \theta & 0&-\sin\theta\\ 0 &1&0\\ \sin\theta&0&\cos\theta\end{pmatrix}$ et où $\theta$ est choisi (indépendamment du reste) suivant la loi uniforme sur $[0,2\pi]$.

Saturne · July 2018

Merci je vais regarder ça et comparer à l'implémentation fournie par le package R alphahull3d (qui ne donne pas toujours le même résultat que l'intégration numérique).

En fait j'ai fait un package R qui implémente des algorithmes pour simuler la loi uniforme sur/dans diverses formes géométriques : sphère, triangle sphérique, "sphere patch", cube, ellipsoïde, triangle, polygone, $\ell^p$-boule, anneau, simplexe. J'aimerais bien y ajouter le tore.

Saturne · July 2018

Bien joué Aléa, ça a l'air de marcher.

# aléa (non optimisé) ----
n <- 250000
R <- 3; r <- 2
M <- function(theta){
  rbind(c(cos(theta), 0, -sin(theta)),
        c(0,1,0),
        c(sin(theta),0,cos(theta)))
}
out <- matrix(NA_real_, nrow=n, ncol=3)
j <- 0
while(j < n){
  X <- runif(1, R-r, R+r)
  Y <- runif(1 , -r, r)
  test1 <- (X-R)^2+Y^2 > r^2
  if(test1) next
  test2 <- runif(1, 0, R+r) > X
  if(test2) next
  j <- j + 1
  out[j,] <- c(M(runif(1,0,2*pi)) %*% c(X,Y,0))
}

Test: évaluation d'une intégrale:

> f <- function(x) exp(-x[1]) + exp(-x[2]) + exp(-x[3])
> mean(apply(out, 1, f)) * 2*pi^2*R*r^2
[1] 4476.286

Comparaison avec intégration numérique (kif-kif):

> library(cubature)
> g <- function(uvw){
+   u <- uvw[1]; v <- uvw[2]; w <- uvw[3]
+   f(c((R+w*cos(v))*cos(u), (R+w*cos(v))*sin(u), w*sin(v)))*(R+w*cos(v))*w
+ }
> pcubature(g, c(0,0,0), c(2*pi,2*pi,r))$integral
[1] 4472.645

Comparaison avec les simulations du package alphashape3d (pas trop kif-kif):

> sims <- alphashape3d::rtorus(n, r, R)
> mean(apply(sims, 1, f)) * 2*pi^2*R*r^2
[1] 4494.227

aléa · July 2018

On peut améliorer un peu la vitesse, pour rejeter moins souvent, en prenant
$U$ suivant la loi uniforme sur $[0,1]$, $\alpha$ la loi uniforme sur $[0,2\pi]$ et toujours $T$ suivant la loi uniforme sur $[0,h+r]$.
On pose alors $X=h+r\sqrt{U}\cos \alpha$, $Y=r\sqrt{U}\sin\alpha$ et on rejette si $T>X$.
La suite est inchangée.

Saturne · July 2018

Je vais y regarder.
Autre chose : on peut supprimer la 3ème colonne de ta matrice vu qu'on multiplie par $(X,Y,0)'$ (et supprimer ce $0$)..

Saturne · July 2018

Merci aléa, ton deuxième algo fonctionne aussi.

Pour compléter ce fil, voici un algo pour générer uniformément sur la surface d'un tore ; math.stackexchange.

mohanad1 · July 2018

Bonjour,
Voila une question qui a été déjà posée sous différentes formes
Il y a un terme qui a attiré mon attention : "répartition uniforme". J'ai déjà été confronté à ce problème. Début des années 80, je voulais dessiner des petits symboles, aléatoirement, pour figurer des espaces verts sur des plans. A cette époque, la fonction rand() était inconnue, en tout cas non disponible. Donc, j'ai bricolé un truc qui calculait un X et un Y aléatoire pour dessiner mes symboles. C'était donc une fonction uniforme.
A ma grande surprise, la répartition sur le dessin n'était pas uniforme. Alors mes notions de probabilités me sont revenues en mémoire et, évidemment, la répartition était gaussienne, comme le dit le TCL (dont je ne connaissais pas le nom).
Par ailleurs, on ne voit plus Egoroff depuis longtemps, et qu''est-il arrivé à Remarque, il n'est même plus inscrit ?
Bonne journée.

leon1789 · July 2018

Bonjour

mohanad1 écrivait:

> A cette époque, la fonction rand() était inconnue, en tout cas non disponible.

Bien sûr qu'il y avait des fonctions de tirage aléatoire à l'époque (sous Pascal par exemple).
Elles étaient peut-être moins sophistiquées qu'aujourd'hui...

> C'était donc une fonction uniforme.

Une fonction uniforme : quesako ??

la répartition était gaussienne,

Donc les petits symboles étaient surtout concentrés à un certain endroit dans l'espace vert !
Tu n'as pas réussi à obtenir une répartition uniforme. C'est pourtant pas si compliqué en général...
Il fallait demander à un paysagiste de faire le travail B-)