Générer un point aléatoire dans un cercle (uniformément)

Abordons ceci comme Archimède l'aurait fait.

Comment Pouvons-nous générer un point uniformément dans un triangle ABC, où / AB / = / BC/? Rendons cela plus facile en étendant à un parallélogramme ABCD. Il est facile de générer des points uniformément dans ABCD. Nous choisissons uniformément un point aléatoire X sur AB et Y sur BC et choisissons Z de telle sorte que XBYZ soit un parallélogramme. Pour obtenir un point uniformément choisi dans le triangle d'origine, il suffit de plier tous les points qui apparaissent dans ADC vers le bas à ABC le long de AC.

Maintenant considérez un cercle. Dans la limite, nous pouvons le considérer comme infiniment nombreux triangles isocèles ABC avec B à l'origine et A et C sur la circonférence disparaissant les uns des autres. Nous pouvons choisir un de ces triangles simplement en choisissant un angle thêta. Nous devons donc maintenant générer une distance du centre en choisissant un point dans le ruban ABC. Encore une fois, étendre à ABCD, où D est maintenant deux fois le rayon du centre du cercle.

Choisir un point aléatoire dans ABCD est facile en utilisant la méthode ci-dessus. Choisir un point aléatoire sur AB. Choisissez uniformément un point aléatoire sur BC. IE. choisissez une paire de nombres aléatoires X et y uniformément sur [0, R] donnant des distances du centre. Notre triangle est un mince ruban donc AB et BC sont essentiellement parallèles. Donc, le point Z est simplement une distance x + y de l'origine. Si x + y > R nous replions vers le bas.

Voici l'algorithme complet pour R = 1. J'espère que vous êtes d'accord que c'est assez simple. Il utilise trig, mais vous pouvez donner une garantie sur combien de temps cela prendra, et combien d'appels random() il a besoin, contrairement à l'échantillonnage de rejet.

t = 2*pi*random()
u = random()+random()
r = if u>1 then 2-u else u
[r*cos(t), r*sin(t)]

Le voici dans Mathematica.

f[] := Block[{u, t, r},
  u = Random[] + Random[];
  t = Random[] 2 Pi;
  r = If[u > 1, 2 - u, u];
  {r Cos[t], r Sin[t]}
]

ListPlot[Table[f[], {10000}], AspectRatio -> Automatic]

Voici une solution simple et rapide.

Choisir deux nombres aléatoires dans l'intervalle (0, 1), à savoir a et b. Si b < a, échangez-les. Votre point est (b*R*cos(2*pi*a/b), b*R*sin(2*pi*a/b)).

, Vous pouvez penser à cette solution comme suit. Si vous prenez le cercle, le coupez, puis le redressez, vous obtiendrez un triangle rectangle. Mettez ce triangle à l'échelle, et vous auriez un triangle de (0, 0) à (1, 0) à (1, 1) et de nouveau à (0, 0). Toutes ces transformations changent la densité uniformément. Ce que vous avez fait est uniformément choisi un point au hasard dans le triangle inversé le processus pour obtenir un point dans le cercle.

Notez la densité de points proportionnelle au carré inverse du rayon, donc au lieu de choisir r à partir de [0, r_max], Choisissez à partir de [0, r_max^2], puis calculez vos coordonnées comme:

x = sqrt(r) * cos(angle)
y = sqrt(r) * sin(angle)

Cela vous donnera une distribution de points uniforme sur un disque.

Http://mathworld.wolfram.com/DiskPointPicking.html

Pensez-y de cette façon. Si vous avez un rectangle où un axe est Rayon et un est angle, et vous prenez les points à l'intérieur de ce rectangle qui sont proches du rayon 0. Ceux-ci tomberont tous très proches de l'origine (c'est-à-dire proches les uns des autres sur le cercle.) Cependant, les points proches du rayon R, ceux-ci tomberont tous près du bord du cercle (c'est-à-dire éloignés les uns des autres.)

Cela pourrait vous donner une idée de la raison pour laquelle vous obtenez ce comportement.

Le facteur qui est dérivé sur ce lien vous indique combien de zone correspondante dans le rectangle doit être ajustée pour ne pas dépendre du rayon Une fois qu'il est mappé au cercle.

Edit: donc ce qu'il écrit dans le lien que vous partagez est ,"c'est assez facile à faire en calculant l'inverse de la distribution cumulative, et nous obtenons pour r:".

La prémisse de base est ici que vous pouvez créer une variable avec une distribution souhaitée à partir d'un uniforme en mappant l'uniforme par la fonction inverse du cumulatif la fonction de distribution, de la fonction de densité de probabilité. Pourquoi? Il suffit de prendre pour acquis pour l'instant, mais c'est un fait.

Voici mon explication intuitive du calcul. La fonction de densité f(r) par rapport à r est proportionnelle à r lui-même. Comprendre ce fait fait partie de tous les livres de calcul de base. Voir les sections sur les éléments de la zone polaire. Certains autres affiches ont mentionné cette.

Donc nous l'appellerons f(r) = C * r;

Cela s'avère être la plupart de travail. Maintenant, puisque f(r) devrait être une densité de probabilité, vous pouvez facilement voir qu'en intégrant f (r) sur l'intervalle (0,R) Vous obtenez que C = 2/R^2 (c'est un exercice pour le lecteur.)

Ainsi, f (r) = 2 * r / r ^ 2

OK, c'est ainsi que vous obtenez la formule dans le lien.

Ensuite, la partie finale va de la variable aléatoire uniforme u dans (0,1) vous devez mapper par la fonction inverse de la fonction de distribution cumulative à partir de cette densité désirée f (r). Comprendre pourquoi c'est le cas, vous devez trouver un texte de probabilité avancé comme Papoulis probablement (ou dérivez-le vous-même.)

En intégrant f (r) vous obtenez F (r) = r^2 / R^2

Pour trouver la fonction inverse de ceci, vous définissez u = r^2 / R^2 puis résolvez pour r, ce qui vous donne r = r * sqrt (u)

Cela a tout à fait du sens intuitivement aussi, u = 0 devrait correspondre à r = 0. En outre, u = 1 shoudl map à r = R. En outre, il va par la fonction racine carrée, ce qui est logique et correspond au lien.

Cela dépend vraiment de ce que vous entendez par "uniformément aléatoire". C'est un point subtil et vous pouvez en savoir plus sur la page wiki ici: http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29 , où le même problème, donner des interprétations différentes à "uniformément aléatoire" donne des réponses différentes!

Selon la façon dont vous choisissez les points, la distribution peut varier, même si elles sont uniformément aléatoire, dans certains sens.

Il semble que le blog l'entrée essaie de la rendre uniformément aléatoire dans le sens suivant: si vous prenez un sous-cercle du cercle, avec le même centre, alors la probabilité que le point tombe dans cette région est proportionnelle à la zone de la région. Cela, je crois, tente de suivre l'interprétation maintenant standard de "uniformément aléatoire" pour les régions 2D avec des zones définies sur elles: la probabilité qu'un point tombe dans n'importe quelle région (avec une zone bien définie) est proportionnelle à la zone de cette région.

La raison pour laquelle la solution naïve ne fonctionne pas est qu'elle donne une densité de probabilité plus élevée aux points plus proches du centre du cercle. En d'autres termes, le cercle qui a le rayon r/2 a la probabilité r/2 d'obtenir un point sélectionné, mais il a la zone (nombre de points) pi*r^2/4.

Par conséquent, nous voulons qu'une densité de probabilité de rayon ait la propriété suivante:

La probabilité de choisir un rayon inférieur ou égal à un r donné doit être proportionnelle à l'aire de le cercle avec le rayon R. (parce que nous voulons avoir une distribution uniforme sur les points et les grandes surfaces signifient plus de points)

En d'autres termes,nous voulons que la probabilité de choisir un rayon entre [0, r] soit égale à sa part de l'aire globale du cercle. L'aire totale du cercle est pi * r^2, et l'aire du cercle de rayon r est pi * r^2. Ainsi, nous aimerions que la probabilité de choisir un rayon entre [0, r] soit (pi * r ^ 2) / (pi*r^2) = r ^ 2 / R^2.

Maintenant vient le calcul:

La probabilité de choisir un rayon entre [0, r] est l'intégrale de P (r) dr de 0 à r (c'est juste parce que nous ajoutons toutes les probabilités des rayons plus petits). Ainsi, nous voulons intégrale (p (r)dr) = r ^ 2 / R ^ 2. Nous pouvons clairement voir que R ^ 2 est une constante, donc tout ce que nous devons faire est de comprendre quel p (r), Une fois intégré, nous donnerait quelque chose comme r^2. La réponse est clairement R * constante. intégrale (R * constante dr) = R ^ 2/2 * constante. Ce doit être égal à r^2/R^2, donc constante = 2/R^2. Ainsi, vous avez la distribution de probabilité p (r) = r * 2 / R^2

Note: une autre façon plus intuitive de penser au problème est d'imaginer que vous essayez de donner à chaque cercle de rayon r une densité de probabilité égale à la proportion du nombre de points qu'il a sur sa circonférence. Donc un cercle de rayon r aurez 2 * pi * r "points" sur sa circonférence. Le nombre total de points est pi * r ^ 2. Ainsi, vous devez donner au cercle r une probabilité égale à (2 * pi * r) / (pi * R^2) = 2 * r/R^2. C'est beaucoup plus facile à comprendre et plus intuitif, mais ce n'est pas aussi mathématique.

Voici mon code Python pour générer num points aléatoires à partir d'un cercle de rayon rad:

import matplotlib.pyplot as plt
import numpy as np
rad = 10
num = 1000

t = np.random.uniform(0.0, 2.0*np.pi, num)
r = rad * np.sqrt(np.random.uniform(0.0, 1.0, num))
x = r * np.cos(t)
y = r * np.sin(t)

plt.plot(x, y, "ro", ms=1)
plt.axis([-15, 15, -15, 15])
plt.show()

Soit ρ (rayon) et φ (Azimut) deux variables aléatoires correspondant aux coordonnées polaires d'un point arbitraire à l'intérieur du cercle. Si les points sont répartis uniformément, Quelle est la fonction de répartition de ρ Et φ?

Pour tout r: 0

P [ρ 2

Où S1 et S0 sont les zones de cercle de rayon R et R respectivement. Si Le CDF peut être donné comme:

          0          if r<=0
  CDF =   (r/R)**2   if 0 < r <= R
          1          if r > R

Et PDF:

PDF = d/dr(CDF) = 2 * (r/R**2) (0 < r <= R).

Notez que pour R=1 variable aléatoire sqrt(X) où X est uniforme sur [0, 1) a exactement le CDF (car P[sqrt(X)

La distribution de φ est évidemment uniforme de 0 à 2*pi. Maintenant, vous pouvez créer des coordonnées polaires aléatoires et les convertir en cartésiens en utilisant des équations trigonométriques:

x = ρ * cos(φ)
y = ρ * sin(φ)

Ne peut pas résister à poster du code python pour R = 1.

from matplotlib import pyplot as plt
import numpy as np

rho = np.sqrt(np.random.uniform(0, 1, 5000))
phi = np.random.uniform(0, 2*np.pi, 5000)

x = rho * np.cos(phi)
y = rho * np.sin(phi)

plt.scatter(x, y, s = 4)

Vous allez obtenir

Solution en Java et l'exemple de distribution (2000 points)

public void getRandomPointInCircle() {
    double t = 2 * Math.PI * Math.random();
    double r = Math.sqrt(Math.random());
    double x = r * Math.cos(t);
    double y = r * Math.sin(t);
    System.out.println(x);
    System.out.println(y);
}

Basé sur la solution previus https://stackoverflow.com/a/5838055/5224246 depuis @ sigfpe

Je pense que dans ce cas, l'utilisation de coordonnées polaires est un moyen de compliquer le problème, il serait beaucoup plus facile si vous choisissez des points aléatoires dans un carré avec des côtés de longueur 2R, puis sélectionnez les points (x,y) tels que x^2+y^2<=R^2.

Nous générons D'abord un cdf [x] qui est

Probabilité qu'un point soit inférieur à la distance x du centre du cercle. Supposons que le cercle a un rayon de R.

Évidemment si x est zéro alors cdf[0] = 0

Évidemment si x est R alors le cdf [R] = 1

, Évidemment, si x = r alors le cdf[r] = (Pi r^2)/(Pi R^2)

C'est parce que chaque "petite zone" sur le cercle a la même probabilité d'être choisi, donc la probabilité est proportionnellement à la zone dans question. Et la zone donnée par une distance x du centre du cercle est Pi r ^ 2

Donc cdf [x] = x ^ 2 / R^2 parce que les Pi s'annulent

, Nous avons cdf[x]=x^2/R^2 où x va de 0 à R

Donc nous résolvons pour x

R^2 cdf[x] = x^2

x = R Sqrt[ cdf[x] ]

Nous pouvons maintenant remplacer cdf avec un nombre aléatoire de 0 à 1

x = R Sqrt[ RandomReal[{0,1}] ]

Enfin

r = R Sqrt[  RandomReal[{0,1}] ];
theta = 360 deg * RandomReal[{0,1}];
{r,theta}

, Nous obtenons les coordonnées polaires {0.601168 R, 311.915 deg}

Il existe une relation linéaire entre le rayon et le nombre de points "proches" de ce rayon, il doit donc utiliser une distribution de rayon qui rend également le nombre de points de données proches d'un rayon r proportionnel à r.

J'ai utilisé une fois cette méthode: Cela peut être totalement non optimisé (c'est-à-dire qu'il utilise un tableau de points donc inutilisable pour les grands cercles) mais donne suffisamment de distribution aléatoire. Vous pouvez sauter la création de la matrice et dessiner directement si vous le souhaitez. La méthode consiste à randomiser tous les points d'un rectangle qui tombent à l'intérieur du cercle.

bool[,] getMatrix(System.Drawing.Rectangle r) {
    bool[,] matrix = new bool[r.Width, r.Height];
    return matrix;
}

void fillMatrix(ref bool[,] matrix, Vector center) {
    double radius = center.X;
    Random r = new Random();
    for (int y = 0; y < matrix.GetLength(0); y++) {
        for (int x = 0; x < matrix.GetLength(1); x++)
        {
            double distance = (center - new Vector(x, y)).Length;
            if (distance < radius) {
                matrix[x, y] = r.NextDouble() > 0.5;
            }
        }
    }

}

private void drawMatrix(Vector centerPoint, double radius, bool[,] matrix) {
    var g = this.CreateGraphics();

    Bitmap pixel = new Bitmap(1,1);
    pixel.SetPixel(0, 0, Color.Black);

    for (int y = 0; y < matrix.GetLength(0); y++)
    {
        for (int x = 0; x < matrix.GetLength(1); x++)
        {
            if (matrix[x, y]) {
                g.DrawImage(pixel, new PointF((float)(centerPoint.X - radius + x), (float)(centerPoint.Y - radius + y)));
            }
        }
    }

    g.Dispose();
}

private void button1_Click(object sender, EventArgs e)
{
    System.Drawing.Rectangle r = new System.Drawing.Rectangle(100,100,200,200);
    double radius = r.Width / 2;
    Vector center = new Vector(r.Left + radius, r.Top + radius);
    Vector normalizedCenter = new Vector(radius, radius);
    bool[,] matrix = getMatrix(r);
    fillMatrix(ref matrix, normalizedCenter);
    drawMatrix(center, radius, matrix);
}

L'élément de zone dans un cercle est dA = rdr * dphi. Ce facteur supplémentaire r détruit votre idée de choisir au hasard un r et phi. Alors que phi est distribué à plat, r ne l'est pas, mais à plat dans 1 / r (c'est-à-dire que vous êtes plus susceptible de frapper la limite que "l'œil du taureau").

Donc, pour générer des points répartis uniformément sur le cercle, choisissez phi à partir d'une distribution plate et r à partir d'une distribution 1/R.

Utilisez alternativement la méthode de Monte Carlo proposée par Mehrdad.

Modifier

À choisissez un r aléatoire dans 1 / R Vous pouvez choisir un X aléatoire de l'intervalle [1/ R, infini] et calculer r=1 / X. r est ensuite distribué à plat dans 1 / r.

Pour calculer un phi aléatoire, choisissez un X aléatoire dans l'intervalle [0, 1] et calculez phi=2 * pi * x.

Je ne sais pas si cette question est toujours ouverte pour une nouvelle solution avec toute la réponse déjà donnée, mais il m'est arrivé d'avoir fait face exactement à la même question moi-même. J'ai essayé de" raisonner " avec moi-même pour une solution, et j'en ai trouvé une. C'est peut-être la même chose que certains l'ont déjà suggéré ici, mais de toute façon c'est ici:

Pour que deux éléments de la surface du cercle soient égaux, en supposant des dr égaux, nous devons avoir dtheta1 / dtheta2 = r2 / r1. Expression écrite de la probabilité pour cet élément comme P(r, theta) = P{ r1

Une solution de programmeur:

• Créer une carte de bits (une matrice de valeurs booléennes). Il peut être aussi grand que vous le souhaitez.
• dessinez un cercle dans cette carte de bits.
• crée une table de recherche des points du cercle.
• Choisissez un index aléatoire dans cette table de recherche.

const int RADIUS = 64;
const int MATRIX_SIZE = RADIUS * 2;

bool matrix[MATRIX_SIZE][MATRIX_SIZE] = {0};

struct Point { int x; int y; };

Point lookupTable[MATRIX_SIZE * MATRIX_SIZE];

void init()
{
  int numberOfOnBits = 0;

  for (int x = 0 ; x < MATRIX_SIZE ; ++x)
  {
    for (int y = 0 ; y < MATRIX_SIZE ; ++y)
    {
      if (x * x + y * y < RADIUS * RADIUS) 
      {
        matrix[x][y] = true;

        loopUpTable[numberOfOnBits].x = x;
        loopUpTable[numberOfOnBits].y = y;

        ++numberOfOnBits;

      } // if
    } // for
  } // for
} // ()

Point choose()
{
  int randomIndex = randomInt(numberOfBits);

  return loopUpTable[randomIndex];
} // ()

Le bitmap n'est nécessaire que pour l'explication de la logique. C'est le code sans le bitmap:

const int RADIUS = 64;
const int MATRIX_SIZE = RADIUS * 2;

struct Point { int x; int y; };

Point lookupTable[MATRIX_SIZE * MATRIX_SIZE];

void init()
{
  int numberOfOnBits = 0;

  for (int x = 0 ; x < MATRIX_SIZE ; ++x)
  {
    for (int y = 0 ; y < MATRIX_SIZE ; ++y)
    {
      if (x * x + y * y < RADIUS * RADIUS) 
      {
        loopUpTable[numberOfOnBits].x = x;
        loopUpTable[numberOfOnBits].y = y;

        ++numberOfOnBits;
      } // if
    } // for
  } // for
} // ()

Point choose()
{
  int randomIndex = randomInt(numberOfBits);

  return loopUpTable[randomIndex];
} // ()

Je ne suis toujours pas sûr du '(2/R2)×R ' exact, mais ce qui est apparent, c'est le nombre de points à distribuer dans l'unité donnée 'dr', c'est-à-dire que l'augmentation de r sera proportionnelle à r2 et non à r.

Vérifiez de cette façon...nombre de points à un certain angle thêta et entre r (0.1 r à 0.2 r) c'est-à-dire fraction du R et nombre de points entre r (0.6 r à 0.7 r) serait égal si vous utilisez la génération standard, puisque la différence est seulement 0.1 r entre deux intervalles. mais depuis la zone couverte entre les points (0.6 R à 0.7 r) sera beaucoup plus grande que la surface couverte entre 0.1 r à 0.2 r, le nombre égal de points sera peu espacé dans une plus grande surface, je suppose que vous le savez déjà, donc la fonction pour générer les points aléatoires ne doit pas être linéaire mais quadratique, (puisque le nombre de points requis pour être distribué dans une unité donnée ' dr ' c'est-à-dire l'augmentation de r sera proportionnelle à r2 et non r), donc dans ce cas ce sera inverse du quadratique, puisque le delta que nous avons (0.1 r) dans les deux intervalles doit être carré d'une fonction afin qu'il puisse agir comme valeur de graine pour la génération linéaire de points (puisque les mots après, cette graine est utilisée linéairement dans la fonction sin et cos), donc nous savons, dr doit être une valeur quadratique et pour rendre cette graine quadratique, nous devons provenir de ces valeurs à partir de la racine carrée de r PAS r lui-même, j'espère que cela le rend un peu plus clair.

Un tel problème amusant.
La logique de la probabilité qu'un point soit choisi s'abaissant à mesure que la distance par rapport à l'origine de l'axe augmente est expliquée plusieurs fois ci-dessus. Nous tenons compte de cela en prenant la racine de U [0,1]. Voici une solution générale pour un r positif en Python 3.

import numpy
import math
import matplotlib.pyplot as plt

def sq_point_in_circle(r):
    """
    Generate a random point in an r radius circle 
    centered around the start of the axis
    """

    t = 2*math.pi*numpy.random.uniform()
    R = (numpy.random.uniform(0,1) ** 0.5) * r

    return(R*math.cos(t), R*math.sin(t))

R = 200 # Radius
N = 1000 # Samples

points = numpy.array([sq_point_in_circle(R) for i in range(N)])
plt.scatter(points[:, 0], points[:,1])

Vous pouvez aussi utiliser votre intuition.

L'aire d'un cercle est pi*r^2

Pour r=1

Ceci nous donne une superficie de pi. Supposons que nous ayons une sorte de fonction f qui diviserait uniformément les points N=10 à l'intérieur d'un cercle. Le rapport ici est 10 / pi

Maintenant, nous doublons la surface et le nombre de points

Pour r=2 et N=20

Cela donne une aire de 4pi et le rapport est maintenant 20/4pi ou 10/2pi. Le rapport deviendra plus petit et plus le rayon est grand, plus sa croissance est quadratique et les N s'échelonnent linéairement.

Pour résoudre ce problème, nous pouvons simplement dire

x = r^2
sqrt(x) = r

Si vous générez un vecteur en coordonnées polaires comme ceci

length = random_0_1();
angle = random_0_2pi();

Plus de points atterriraient autour du centre.

length = sqrt(random_0_1());
angle = random_0_2pi();

length n'est plus uniformément distribué, mais le vecteur sera maintenant uniformément distribué.

1) Choisissez un X aléatoire entre -1 et 1.

var X:Number = Math.random() * 2 - 1;

2) en utilisant la formule du cercle, calculer les valeurs maximales et minimales de Y étant donné que X et un rayon de 1:

var YMin:Number = -Math.sqrt(1 - X * X);
var YMax:Number = Math.sqrt(1 - X * X);

3) Choisissez un y aléatoire entre ces extrêmes:

var Y:Number = Math.random() * (YMax - YMin) + YMin;

4) incorporez vos valeurs de localisation et de rayon dans la valeur finale:

var finalX:Number = X * radius + pos.x;
var finalY:Number = Y * radois + pos.y;