我需要在半径为 R 的圆内生成一个均匀随机的点。
我意识到,只需在区间 [0 ... 2π) 中选择一个均匀随机的角度,在区间 (0 ... R) 中选择一个均匀随机的半径,我最终会得到更多指向中心的点,因为对于两个给定半径,较小半径中的点将比较大半径中的点更靠近彼此。
我在此 over here 上找到了一篇博客文章,但我不明白他的推理。我想这是正确的,但我真的很想了解他从哪里得到 (2/R2)×r 和他是如何得出最终解决方案的。
更新:在发布此问题 7 年后,我仍然没有收到关于平方根算法背后数学的实际问题的满意答案。所以我花了一天时间自己写一个答案。 Link to my answer。
如何在半径为 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 的圆周长上的点数的两倍。
https://i.stack.imgur.com/lAyvB.png
由于圆的周长 (2πr) 随 r 线性增长,因此随机点的数量应随 r 线性增长。换句话说,所需的 probability density function (PDF) 线性增长。由于 PDF 的面积应该等于 1,最大半径是 1,所以我们有
https://i.stack.imgur.com/WI9F7.png
所以我们知道我们的随机值的期望密度应该是什么样子。现在:当我们只有一个介于 0 和 1 之间的均匀随机值时,我们如何生成这样一个随机值?
我们使用了一个叫做 inverse transform sampling 的技巧
从 PDF 创建累积分布函数 (CDF) 沿 y = x 镜像此函数 将结果函数应用于 0 和 1 之间的统一值。
听起来很复杂?让我插入一个带有小旁道的块引用来传达直觉:
假设我们要生成一个具有以下分布的随机点:即均匀分布在 1 和 2 之间的点的 1/5,均匀分布在 2 和 3 之间的点的 4/5。顾名思义,CDF 是PDF 的累积版本。直观地说:虽然 PDF(x) 描述了 x 处的随机值的数量,但 CDF(x) 描述了小于 x 的随机值的数量。在这种情况下,CDF 看起来像: 要了解它的用途,假设我们以均匀分布的高度从左到右发射子弹。当子弹击中线时,它们会掉到地上:看看地上的子弹密度如何与我们想要的分布相对应!我们快到了!问题是对于这个函数,y 轴是输出,x 轴是输入。我们只能“从地上直射子弹”!我们需要反函数!这就是我们反映整个事情的原因; x 变成 y,y 变成 x:我们称之为 CDF-1。为了根据所需分布获取值,我们使用 CDF-1(random())。
…所以,回到生成我们的 PDF 等于 2x 的随机半径值。
第 1 步:创建 CDF:由于我们使用的是实数,因此 CDF 表示为 PDF 的积分。
CDF(x) = ∫ 2x = x2
第 2 步:沿 y = x 镜像 CDF:
从数学上讲,这归结为交换 x 和 y 并求解 y:
CDF:y = x2 交换:x = y2 求解:y = √x CDF-1:y = √x
第 3 步:将结果函数应用于 0 到 1 之间的统一值
CDF-1(随机()) = √随机()
这就是我们打算得出的:-)
让我们像阿基米德那样处理这个问题。
我们如何在三角形 ABC 中均匀地生成一个点,其中 |AB|=|BC|?让我们通过扩展到平行四边形 ABCD 来简化这个过程。在 ABCD 中很容易均匀地生成点。我们统一选择 AB 上的随机点 X 和 BC 上的 Y 并选择 Z 使得 XBYZ 是平行四边形。为了在原始三角形中获得一个统一选择的点,我们只需将 ADC 中出现的任何点沿 AC 折回 ABC 即可。
现在考虑一个圆圈。在极限中,我们可以将其视为无限多的等腰三角形 ABC,其原点为 B,圆周上的 A 和 C 彼此接近消失。我们可以简单地通过选择角度 theta 来选择其中一个三角形。所以我们现在需要通过在长条 ABC 中选择一个点来生成到中心的距离。同样,延伸到 ABCD,其中 D 现在是圆心半径的两倍。
使用上述方法很容易在 ABCD 中选择一个随机点。在 AB 上随机选取一个点。在 BC 上统一选择一个随机点。 IE。在 [0,R] 上均匀选取一对随机数 x 和 y,给出与中心的距离。我们的三角形是一条细条,所以 AB 和 BC 基本上是平行的。所以点 Z 只是距离原点 x+y 的距离。如果 x+y>R 我们向下折叠。
这是 R=1 的完整算法。我希望你同意这很简单。它使用 trig,但您可以保证它需要多长时间,以及它需要多少 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]
https://i.stack.imgur.com/Nty8O.gif
random()+random()+random()
与一些更复杂的折叠一起使用(即,将无限小的平行六面体折叠成四面体)。虽然不相信这是一个好方法。
这是一个快速简单的解决方案。
在 (0, 1) 范围内选取两个随机数,即 a
和 b
。如果b < a
,交换它们。你的观点是(b*R*cos(2*pi*a/b), b*R*sin(2*pi*a/b))
。
您可以按如下方式考虑此解决方案。如果你把圆剪下来,然后把它拉直,你会得到一个直角三角形。缩小该三角形,您将得到一个从 (0, 0)
到 (1, 0)
到 (1, 1)
并再次返回到 (0, 0)
的三角形。所有这些变换都会均匀地改变密度。你所做的就是在三角形中均匀地选择一个随机点,然后颠倒这个过程,在圆圈中得到一个点。
b < a
时不交换,我们可以实现这一点!例如在 javascript jsfiddle.net/b0sb5ogL/1 中
注意点密度与半径的平方成反比,因此不是从 [0, r_max]
中选择 r
,而是从 [0, r_max^2]
中选择,然后将坐标计算为:
x = sqrt(r) * cos(angle)
y = sqrt(r) * sin(angle)
这将为您提供磁盘上的均匀点分布。
http://mathworld.wolfram.com/DiskPointPicking.html
这样想吧。如果您有一个矩形,其中一个轴是半径,一个轴是角度,并且您在该矩形内取半径 0 附近的点。这些点都将非常靠近原点(在圆上靠得很近)。但是,半径 R 附近的点,它们都将落在圆的边缘附近(即彼此相距很远)。
这可能会让您了解为什么会出现这种行为。
在该链接上派生的因子告诉您矩形中的对应区域需要调整多少,以便在映射到圆后不依赖于半径。
编辑:所以他在您分享的链接中写道,“通过计算累积分布的倒数很容易做到这一点,我们得到 r:”。
这里的基本前提是,您可以通过使用所需概率密度函数的累积分布函数的反函数映射制服来从制服中创建具有所需分布的变量。为什么?现在只是认为这是理所当然的,但这是事实。
这是我对数学的直观解释。关于 r 的密度函数 f(r) 必须与 r 本身成比例。理解这一事实是任何基本微积分书籍的一部分。请参阅有关极地元素的部分。其他一些海报也提到了这一点。
所以我们称它为 f(r) = C*r;
事实证明这是大部分工作。现在,由于 f(r) 应该是一个概率密度,您可以很容易地看到,通过在区间 (0,R) 上对 f(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 * sqrt(u)
这在直观上也完全有意义,u = 0 应该映射到 r = 0。此外,u = 1 应该映射到 r = R。此外,它通过平方根函数,这很有意义并匹配链接。
朴素解决方案不起作用的原因是它为靠近圆心的点提供了更高的概率密度。换句话说,半径为 r/2 的圆有 r/2 的概率在其中选择一个点,但它的面积(点数)为 pi*r^2/4。
因此,我们希望半径概率密度具有以下属性:
选择小于或等于给定 r 的半径的概率必须与半径为 r 的圆的面积成正比。 (因为我们希望点分布均匀,面积越大意味着点越多)
换句话说,我们希望在 [0,r] 之间选择半径的概率等于它在圆的总面积中所占的份额。圆的总面积为 pi*R^2,半径为 r 的圆的面积为 pi*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
注意:考虑这个问题的另一种更直观的方法是想象你试图让每个半径为 ra 的圆的概率密度等于它在圆周上的点数的比例。因此,半径为 r 的圆在其圆周上将有 2 * pi * r 个“点”。总点数为 pi * R^2。因此,您应该给圆 ra 概率等于 (2 * pi * r) / (pi * R^2) = 2 * r/R^2。这更容易理解和更直观,但它在数学上并不那么合理。
设 ρ(半径)和 φ(方位角)为两个随机变量,对应于圆内任意点的极坐标。如果点是均匀分布的,那么 ρ 和 φ 的分布函数是什么?
对于任何 r:0 < r < R,半径坐标 ρ 小于 r 的概率为
P[ρ < 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 随机变量 sqrt(X),其中 X 在 [0, 1) 上是一致的,具有这个精确的 CDF(因为 P[sqrt(X) < y] = P[x < y**2] = y* *2 表示 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)
你会得到
https://i.stack.imgur.com/m92pe.png
这真的取决于你所说的“均匀随机”是什么意思。这是一个微妙的点,您可以在此处的 wiki 页面上阅读更多相关信息:http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29,同样的问题,对“均匀随机”给出不同的解释会给出不同的答案!
根据您选择点的方式,分布可能会有所不同,即使它们在某种意义上是均匀随机的。
似乎博客条目试图在以下意义上使其均匀随机:如果你取圆的一个子圆,具有相同的中心,那么该点落在该区域的概率与面积成正比该区域。我相信,这试图遵循现在对具有 areas defined on them 的 2D 区域的“均匀随机”的标准解释:点落在任何区域(区域定义明确)的概率与该区域的面积成正比。
这是我的 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 = np.sqrt(np.random.uniform(0.0, rad**2, num))
?
我认为在这种情况下使用极坐标是使问题复杂化的一种方式,如果您将随机点选入边长为 2R 的正方形中,然后选择点 (x,y)
使得 x^2+y^2<=R^2
会容易得多。
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);
}
https://i.stack.imgur.com/gszyo.png
基于来自@sigfpe 的先前解决方案 https://stackoverflow.com/a/5838055/5224246
我曾经使用过这种方法:这可能是完全未优化的(即它使用点数组,因此它不能用于大圆圈)但提供了足够的随机分布。如果您愿意,可以跳过矩阵的创建并直接绘制。该方法是随机化矩形中落在圆内的所有点。
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);
}
https://i.stack.imgur.com/CgssC.png
首先我们生成一个 cdf[x] 它是
一个点到圆心的距离小于 x 的概率。假设圆的半径为 R。
显然,如果 x 为零,则 cdf[0] = 0
显然,如果 x 是 R,那么 cdf[R] = 1
显然如果 x = r 那么 cdf[r] = (Pi r^2)/(Pi R^2)
这是因为圆上的每个“小区域”都有相同的被选中的概率,所以概率与所讨论的区域成正比。距离圆心 x 的区域是 Pi r^2
所以 cdf[x] = x^2/R^2 因为 Pi 相互抵消
我们有 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 deg}
圆中的面积元素为 dA=rdr*dphi。那个额外的因素 r 破坏了你随机选择 ar 和 phi 的想法。虽然 phi 是平坦分布的,但 r 不是平坦分布,而是 1/r 平坦分布(即,与“靶心”相比,您更有可能触及边界)。
因此,要生成均匀分布在圆上的点,从平面分布中选取 phi,从 1/r 分布中选取 r。
或者使用 Mehrdad 提出的 Monte Carlo 方法。
编辑
要在 1/r 中选择一个随机 r 平坦,您可以从区间 [1/R, infinity] 中选择一个随机 x 并计算 r=1/x。然后 r 以 1/r 平分布。
要计算随机 phi,请从区间 [0, 1] 中选择一个随机 x 并计算 phi=2*pi*x。
你也可以使用你的直觉。
圆的面积是pi*r^2
对于r=1
这给了我们一个 pi
的区域。让我们假设我们有某种函数 f
可以将N=10
点均匀分布在一个圆内。这里的比率是 10 / pi
现在我们将面积和点数加倍
对于 r=2
和 N=20
这给出了 4pi
的面积,并且比率现在是 20/4pi
或 10/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
不再均匀分布,但向量现在将均匀分布。
半径与“靠近”该半径的点数之间存在线性关系,因此他需要使用半径分布,该半径分布也使得半径 r
附近的数据点数与 r
成正比。
我不知道这个问题是否仍然适用于已经给出所有答案的新解决方案,但我碰巧自己也遇到了完全相同的问题。我试图与自己“推理”一个解决方案,我找到了一个。这可能与某些人已经在这里提出的建议相同,但无论如何它是:
为了使圆表面的两个元素相等,假设 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。
我仍然不确定确切的“(2/R2)×r”,但显而易见的是需要在给定单位“dr”中分布的点数,即 r 的增加将与 r2 成正比,而不是与 r 成正比。
以这种方式检查...如果您使用标准生成,则在某个角度 theta 和 r(0.1r 到 0.2r)之间的点数,即 r 的分数和 r(0.6r 到 0.7r)之间的点数将相等,因为两个间隔之间的差异只有 0.1r。但是由于点之间覆盖的区域(0.6r 到 0.7r)将比 0.1r 到 0.2r 之间覆盖的区域大得多,相同数量的点将在更大的区域中稀疏分布,我假设你已经知道了,所以函数生成随机点不能是线性的,而是二次的,(因为需要分布在给定单位“dr”中的点数,即 r 的增加将与 r2 而不是 r 成正比),所以在这种情况下,它将与二次的,因为我们在两个区间中的 delta (0.1r) 必须是某个函数的平方,所以它可以作为点的线性生成的种子值(因为后记,这个种子在 sin 和 cos 函数中线性使用),所以我们知道,dr 必须是二次值,并且为了使这个种子二次,我们需要从 r 的平方根而不是 r 本身得出这个值,我希望这让它更清楚一点。
这么有趣的问题。上面多次解释了随着与轴原点的距离增加而选择点的概率降低的基本原理。我们通过取 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])
https://i.stack.imgur.com/nc2VW.png
这可能会帮助有兴趣选择速度算法的人;最快的方法是(可能?)拒绝抽样。
只需在单位正方形内生成一个点并拒绝它,直到它在一个圆圈内。例如(伪代码),
def sample(r=1):
while True:
x = random(-1, 1)
y = random(-1, 1)
if x*x + y*y <= 1:
return (x, y) * r
尽管它有时可能运行不止一次或两次(而且它不是恒定时间或不适合并行执行),但它要快得多,因为它不使用像 sin
或 cos
这样的复杂公式。
sqrt(...)
的调用。 (当且仅当 x
是 <= 1
时,sqrt(x)
才会是 <= 1
)
程序员解决方案:
创建位图(布尔值矩阵)。它可以任意大。
在该位图中画一个圆圈。
创建圆点的查找表。
在此查找表中选择一个随机索引。
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];
} // ()
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;
不定期副业成功案例分享
random(min_radius², max_radius²)
时,您的意思是否等同于random() * (max_radius² - min_radius²) + min_radius²
,其中random()
返回一个介于 0 和 1 之间的统一值?