在圆内随机生成一个点(均匀地)

我需要在半径R的圆内生成一个均匀随机点。

我意识到,只要在区间[0…2π),且在区间(0…R)我最终会有更多的点指向中心,因为对于两个给定的半径,小半径内的点将比大半径内的点彼此更接近。

我发现了一个关于这个在这里的博客条目,但我不理解他的推理。我想这是正确的,但是我真的很想了解他从哪里得到(2/R2 R以及他如何得到最终的解决方案。


更新:在发布这个问题7年后,关于平方根算法背后的数学问题,我仍然没有收到一个令人满意的答案。所以我花了一天时间自己写了一个答案。链接到我的答案

216488 次浏览

半径和“靠近”该半径的点数之间存在线性关系,因此他需要使用半径分布,这也使得半径r附近的数据点数与r成比例。

这样想。如果你有一个矩形,其中一个轴是半径,一个是角,你取这个矩形内半径为0的点。它们都离原点很近(在圆上很近)然而,半径R附近的点,它们都落在圆的边缘附近(也就是说,彼此相距很远)。

这可能会让你知道为什么你会有这种行为。

在这个链接上导出的因子告诉你,矩形中有多少对应的区域需要调整,以便在映射到圆后不依赖于半径。

编辑:所以他在你分享的链接中写道,“通过计算累积分布的倒数,这很容易做到,我们得到r:”。

这里的基本前提是,通过将均匀分布映射为期望概率密度函数的累积分布函数的逆函数,可以从均匀分布创建一个具有期望分布的变量。为什么?现在把它当做理所当然,但这是事实。

这是我对数学的一些直观解释。密度函数f(r)关于r必须与r本身成比例。理解这个事实是任何微积分基础书的一部分。请参阅有关极区元素的部分。其他一些海报也提到了这一点。

我们记作f(r) = C*r;

这就是大部分的工作。现在,由于f(r)应该是一个概率密度,你可以很容易地看到,通过对f(r)在区间(0,r)上积分,你可以得到C = 2/ r ^2(这是给读者的练习)。

因此,f(r) = 2*r/ r ^2

好,这就是如何得到链接中的公式。

然后,最后一部分是从(0,1)中的均匀随机变量u你必须从这个期望密度f(r)映射到累积分布函数的逆函数。要理解为什么会这样,你可能需要找到像Papoulis这样的高级概率文本(或者自己推导)。

对f(r)积分得到f(r) = r^2/ r^2

为了求出它的反函数你设u = r^2/ r^2然后解出r,得到r = r *√(u)

直观上讲,u = 0映射到r = 0。同样,u = 1应该映射到r = r。同样,它通过平方根函数,这是有意义的,与链接匹配。

圆中的面积元是dA=rdr*dphi。这个额外的因子r破坏了你随机选择r和的想法。虽然phi分布平坦,但r不是,而是在1/r内平坦(也就是说,你更有可能击中边界而不是“靶心”)。

为了生成在圆上均匀分布的点从平面分布中选取r从1/r分布中选取。

或者使用Mehrdad提出的蒙特卡罗方法。

编辑

要在1/r中选择一个随机的r,你可以从区间[1/ r,无穷]中选择一个随机的x,并计算r=1/x。R以1/ R为单位平坦分布。

为了计算一个随机的,从区间[0,1]中选择一个随机的x,并计算=2*pi*x。

这取决于你对"均匀随机"的定义。这是一个微妙的点,你可以在这里的维基页面上阅读更多关于它的内容:http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29,在这里,同样的问题,对“均匀随机”给出不同的解释会给出不同的答案!

根据你如何选择这些点,分布可能会有所不同,即使它们在一些意义上是均匀随机的。

这篇博客似乎试图让它在以下意义上统一随机:如果你取圆的一个子圆,具有相同的中心,那么该点落在该区域的概率与该区域的面积成正比。我相信,这是试图遵循现在对2D区域的上面定义的区域的“均匀随机”的标准解释:一个点落在任何区域(具有明确定义的区域)的概率与该区域的面积成正比。

我认为在这种情况下,使用极坐标是一种使问题复杂化的方式,如果你选择随机点到一个边长为2R的正方形中,然后选择点(x,y),使x^2+y^2<=R^2

我用过这个方法: 这可能是完全未优化的(即它使用了一个点数组,所以它不能用于大圆圈),但它提供了足够的随机分布。如果你愿意,你可以跳过矩阵的创建,直接绘制。该方法是随机化矩形中落在圆内的所有点

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);
}

enter image description here

让我们像阿基米德那样处理这个问题。

我们如何在三角形ABC中均匀地生成一个点,其中|AB|=|BC|?让我们把它扩展到平行四边形ABCD。在ABCD中很容易均匀地生成点。我们均匀地选择AB上的X点和BC上的Y点并选择Z使XBYZ是一个平行四边形。为了在原始三角形中得到一个均匀选择的点,我们只需将ADC中出现的任何点沿AC折叠回ABC。

现在考虑一个圆。在极限情况下,我们可以把它想象成无穷多个等腰三角形ABC, B在原点,A和C在周长上,彼此逐渐接近。我们可以从这些三角形中选择一个角。所以我们现在需要通过在ABC条上选择一点来生成到中心的距离。同样,延伸到ABCD, D现在是圆中心半径的两倍。

使用上述方法可以很容易地在ABCD中选择一个随机点。在AB上随机选一个点,在BC上随机选一个点。Ie。在[0,R]上取一对随机数字x和y,给出离中心的距离。三角形是一条细条AB和BC本质上是平行的。所以Z点到原点的距离是x+y。如果x+y >r我们向下折叠。

这是R=1的完整算法。我希望你同意这很简单。它使用三角函数,但你可以保证它需要多长时间,以及它需要多少random()调用,这与拒绝抽样不同。

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

这里是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]

enter image description here

这里有一个快速而简单的解决方案。

在(0,1)范围内选择两个随机数,即ab。如果b < a,则交换它们。你的观点是(b*R*cos(2*pi*a/b), b*R*sin(2*pi*a/b))

您可以这样考虑这个解决方案。如果你把圆切开,然后把它拉直,你会得到一个直角三角形。将这个三角形缩小,你会得到一个从(0, 0)(1, 0)(1, 1)再回到(0, 0)的三角形。所有这些变换都会均匀地改变密度。你所做的就是在三角形中随机取一个点然后反过来得到圆中的一个点。

我不知道这个问题是否还有新的答案,但我自己碰巧也遇到过同样的问题。我试着跟自己“讲道理”寻找解决办法,我找到了一个。这可能和一些人在这里提出的建议是一样的,但不管怎样,它是这样的:

为了使圆表面的两个元素相等,假设dr相等,我们必须有dtheta1/dtheta2 = r2/r1。将该元素的概率表达式写成P(r, theta) = P{r1<r<R1 + dr, theta1<theta<Theta + dtheta1} = f(r, Theta)*dr*dtheta1,并设置两个概率(r1和r2)相等,我们得到(假设r和Theta是独立的)f(r1)/r1 = f(r2)/r2 =常数,从而得到f(r) = c*r。剩下的部分,根据f(r)是PDF的条件来确定常数c。

1)在-1和1之间随机选择一个X。

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

2)利用圆公式,在X和半径为1的情况下,计算Y的最大值和最小值:

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

3)在这两个极端之间随机选择一个Y:

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

4)将您的位置和半径值合并到最终值中:

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

下面是我的Python代码,用于从半径rad的圆中生成num随机点:

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()

朴素解不起作用的原因是它给了靠近圆中心的点更高的概率密度。换句话说,半径为r/2的圆被选中点的概率为r/2,但它的面积(点的数量)为*r^2/4。

因此,我们希望半径概率密度具有以下性质:

选择半径小于或等于给定r的概率必须与半径为r的圆的面积成正比(因为我们希望在点上有一个均匀的分布,面积越大意味着点越多)。

换句话说,我们希望在[0,r]之间选择半径的概率等于它在圆的总面积中所占的份额。圆的总面积是*R^2,半径为R的圆的面积是*R^2。因此,我们希望在[0,r]之间选择半径的概率为(pi*r^2)/(pi* r^2) = r^2/ r^2。

现在来算算:

选择半径在[0,r]之间的概率是p(r) dr从0到r的积分(这只是因为我们把所有较小半径的概率相加)。因此,我们要求积分(p(r)dr) = r^2/ r^2。我们可以清楚地看到R^2是常数,所以我们要做的就是算出p(R)积分后会得到R^2。答案显然是r *常数。积分(r *常数dr) = r^2/2 *常数。它必须等于r^2/ r^2,因此常数= 2/ r^2。因此,你有概率分布p(r) = r * 2/ r ^2

另一种更直观的思考这个问题的方法是想象你试图给每个半径为r的圆一个概率密度,等于它的周长上的点的数量的比例。因此,半径为r的圆在其周长上有2 * π * r个“点”。点的总数是* R^2。因此,你应该给圆r一个概率等于(2 * * r) / (pi * r ^2) = 2 * r/ r ^2。这更容易理解,也更直观,但在数学上不太合理。

程序员解决方案:

  • 创建一个位图(布尔值的矩阵)。你想要多大就有多大。
  • 在位图中画一个圆。
  • 创建一个圆的点查找表。
  • 在这个查找表中选择一个随机索引。
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];
} // ()

位图仅用于解释逻辑。这是没有位图的代码:

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];
} // ()

注意点密度与半径的平方反比成正比,因此不是从[0, r_max]中选择r,而是从[0, r_max^2]中选择,然后计算你的坐标:

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

这就得到了圆盘上均匀的点分布。

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

我仍然不确定确切的“(2/R2)×r”,但显而易见的是,在给定的单位“dr”中需要分配的点的数量,即r的增加将与R2成正比,而不是r。

这边看……如果使用标准生成,在某个角度theta和r (0.1r到0.2r)之间的点数(即r的分数)和r (0.6r到0.7r)之间的点数将相等,因为两个间隔之间的差值仅为0.1r。但由于覆盖面积点之间(0.6到0.7 r)将远远大于0.1到0.2 r之间的覆盖面积,相同数量的点的距离将稀疏在较大的区域,我假设您已经知道,所以不能线性函数生成随机点,但二次(因为数量的点必须分布在单位“博士”即增加r2将成比例,而不是r),所以在这种情况下它将逆二次,因为我们在两个区间内的delta (0.1r)必须是某个函数的平方,所以它可以作为点的线性生成的种子值(因为后面,这个种子在sin和cos函数中线性使用),所以我们知道,dr必须是二次值,为了使这个种子成为二次的,我们需要从r的平方根中得到这个值,而不是r本身,我希望这让它更清楚一点。

Java解决方案和分发示例(2000分)

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);
}

Distribution 2000 points

基于之前的解决方案https://stackoverflow.com/a/5838055/5224246从@sigfpe

设ρ(半径)和φ(方位角)是两个随机变量,对应于圆内任意一点的极坐标。如果这些点是均匀分布的,那么ρ和φ的分布函数是什么?

对于任意r: 0 <r & lt;半径坐标ρ小于R的概率

P[ρ& lt;r] = P[点在半径r的圆内]= S1 / S0 =(r/ r)2

其中S1和S0分别为半径为r和r的圆的面积。 所以CDF可以给出:

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

和PDF格式:

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

请注意,对于R=1随机变量根号(X),其中X在[0,1]上是均匀的,有这个确切的CDF(因为P[根号(X) <y] = P[x <Y **2] = Y **2 for 0 <Y <= 1)。

φ在0 ~ 2*π范围内分布明显均匀。现在你可以创建随机极坐标,并使用三角方程将其转换为笛卡尔坐标:

x = ρ * cos(φ)
y = ρ * sin(φ)
无法抗拒发布R=1的python代码。

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)

你会得到

enter image description here

首先我们生成一个cdf[x]

一点到圆心的距离小于x的概率。假设圆的半径为R。

显然,如果x = 0,那么cdf[0] = 0

显然,如果x是R,则cdf[R] = 1

显然,如果x = r,则cdf[r] = (r^2)/(r^2)

这是因为圆上的每个“小区域”都有相同的被选中的概率,所以概率与问题区域成比例。距离圆心x的面积是r^2

所以cdf[x] = x^2/R^2因为两者相互抵消了

我们有cdf[x]=x^2/R^2其中x从0到R

我们解出x

R^2 cdf[x] = x^2


x = R Sqrt[ cdf[x] ]

现在我们可以用一个从0到1的随机数来替换cdf

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

最后

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

我们得到极坐标 {0.601168 R, 311.915°}

这样一个有趣的问题 一个点被选择的概率随着距离轴原点的增加而降低的基本原理在上面已经解释了多次。我们通过取U[0,1]的根来解释这一点。 下面是Python 3中正r的一般解决方案

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])

enter image description here

如何在半径为R的圆内生成一个随机点:

r = R * sqrt(random())
theta = random() * 2 * PI

(假设random()统一给出了0到1之间的值)

如果你想把它转换成笛卡尔坐标,你可以做到

x = centerX + r * cos(theta)
y = centerY + r * sin(theta)


为什么sqrt(random()) ?

让我们看看sqrt(random())之前的数学运算。为简单起见,假设我们处理的是单位圆,即R = 1。

点与点之间的平均距离应该是相同的,不管我们看的距离中心有多远。这意味着,例如,观察一个周长为2的圆的周长,我们应该找到的点的数量是周长为1的圆周长上点的数量的两倍。

< p >
,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/lAyvB.png " width = " 400 " / >

< / p >

由于圆的周长(2πr)随r线性增长,因此随机点的数量应该随r线性增长。换句话说,所需的概率密度函数 (PDF)线性增长。由于PDF的面积应该等于1,最大半径是1,我们有

< p >
,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/WI9F7.png " width = " 200 " / > < / p >

所以我们知道随机值的期望密度应该是什么样的。 现在:当我们只有一个0到1之间的均匀随机值时,我们如何生成这样一个随机值? < / p >

我们使用了一个叫做逆变换抽样的技巧

  1. 从PDF中创建累积分布函数 (CDF)
  2. 按照y = x进行镜像
  3. 将得到的函数应用于0到1之间的统一值。

听起来复杂吗?让我插入一段带有小侧轨的引语来传达直觉:

假设我们想要生成一个随机点,分布如下:

,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/58FAy.png " width = " 200 " / >

这是

  • 1和2之间的1/5点
  • 4/5的点在2和3之间。

顾名思义,CDF是PDF的累积版本。直观地:PDF(x)描述了随机值在x的数量,CDF(x)描述了随机值小于x的数量。

在这种情况下,CDF看起来像:

,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/Qe4u0.png " width = " 200 " / >

要了解这是如何有用的,想象我们从左向右在均匀分布的高度发射子弹。当子弹击中线时,它们掉到地上:

,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/eFWsP.png " width = " 200 " / >

看看地面上子弹的密度是如何与我们想要的分布相对应的!我们快到了!

问题是,对于这个函数,y轴是输出,而x轴是输入。我们只能“从地上直射子弹”!我们需要逆函数!

这就是为什么我们镜像整个事物;x变成yy变成x:

,,,,,,,,,,,,,,, < img src = " https://i.stack.imgur.com/UA2sT.png " width = " 200 " / >

我们称之为提供<一口> 1 > < /晚餐。要根据所需的分布获取值,我们使用CDF-1(random())。

所以,回到生成随机半径值,其中PDF等于2x

步骤1:创建CDF:

由于我们处理的是实数,CDF表示为PDF的积分。

__abc0 (__abc1) =∫2__abc1 = __abc1 __abc4

步骤2:按照y = x镜像CDF:

从数学上讲,这可以归结为交换xy并求解y:

< p > 提供:,,,,, y = x2
交换:,,,x = y2
解决:,,,y =·拉迪奇;x
CDF-1:  y = & radical

步骤3:将得到的函数应用于0到1之间的统一值

提供-1(random()) = & radical;random()

这就是我们要推导的:-)

你也可以用你的直觉。

圆的面积为pi*r^2

对于r=1

这给了我们一个pi的面积。让我们假设我们有某种函数__abc1,它将在一个圆内均匀分布N=10点。这里的比率是10 / pi

现在我们把面积和点数翻倍

对于r=2N=20

这给出了一个面积为4pi的比例现在是20/4pi10/2pi。半径越大,这个比值就会越来越小,因为它的增长是二次的,而N是线性缩放的。

为了解决这个问题,我们可以说

x = r^2
sqrt(x) = r

如果在极坐标下生成一个向量

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

更多的点会落在中心周围。

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

length不再是均匀分布的,但是向量现在将是均匀分布的。

这可能会帮助那些对选择速度算法感兴趣的人;最快的方法是(可能?)拒绝抽样。

只需在单位正方形内生成一个点,并拒绝它,直到它在圆内。如(伪代码),

def sample(r=1):
while True:
x = random(-1, 1)
y = random(-1, 1)
if x*x + y*y <= 1:
return (x, y) * r

虽然它有时可能会运行不止一次或两次(而且它不是常量时间或适合并行执行),但它要快得多,因为它不使用sincos这样复杂的公式。