我需要在半径为R的圆内生成一个均匀随机点。
我意识到,只要在区间[0…2π),且在区间(0…R)我最终会有更多的点指向中心,因为对于两个给定的半径,小半径内的点会比大半径内的点彼此更接近。
我在这里找到了一篇关于这个的博客文章,但我不明白他的推理。我认为这是正确的,但我真的很想知道他从哪里得到(2/R2)×r以及他如何得到最终解。
更新:在发布这个问题7年之后,关于平方根算法背后的数学问题,我仍然没有得到一个令人满意的答案。所以我花了一天时间自己写了一个答案。链接到我的答案。
我需要在半径为R的圆内生成一个均匀随机点。
我意识到,只要在区间[0…2π),且在区间(0…R)我最终会有更多的点指向中心,因为对于两个给定的半径,小半径内的点会比大半径内的点彼此更接近。
我在这里找到了一篇关于这个的博客文章,但我不明白他的推理。我认为这是正确的,但我真的很想知道他从哪里得到(2/R2)×r以及他如何得到最终解。
更新:在发布这个问题7年之后,关于平方根算法背后的数学问题,我仍然没有得到一个令人满意的答案。所以我花了一天时间自己写了一个答案。链接到我的答案。
这样想。如果你有一个矩形,其中一个轴是半径,一个是角,你取这个矩形内半径为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提出的蒙特卡罗方法。
EDIT
要在1/r中选择一个随机的r,你可以从区间[1/ r,无穷]中选择一个随机的x,并计算r=1/x。R以1/ R为单位平坦分布。
为了计算一个随机的,从区间[0,1]中选择一个随机的x,并计算=2*pi*x。
这取决于你对"均匀随机"的定义。这是一个微妙的点,你可以在这里的wiki页面上阅读更多关于它的内容:http://en.wikipedia.org/wiki/Bertrand_paradox_%28probability%29,在这里同样的问题,对“均匀随机”给出不同的解释会给出不同的答案!
根据你如何选择这些点,分布可能会有所不同,即使它们在某种意义上是均匀随机的。
It seems like the blog entry is trying to make it uniformly random in the following sense: If you take a sub-circle of the circle, with the same center, then the probability that the point falls in that region is proportional to the area of the region. That, I believe, is attempting to follow the now standard interpretation of 'uniformly random' for 2D regions with areas defined on them: probability of a point falling in any region (with area well defined) is proportional to the area of that region.
我认为在这种情况下,使用极坐标是一种使问题复杂化的方法,如果你在一个边长为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);
}
让我们像阿基米德那样处理这个问题。
我们如何在三角形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]
这里有一个快速而简单的解决方案。
在(0,1)范围内选择两个随机数,即a和b。如果b < a,则交换它们。你的观点是(b * R * cos(2 *π* a / b), b * R * sin(2 *π* a / b))。
您可以这样考虑这个解决方案。如果你把圆切开,然后把它拉直,你会得到一个直角三角形。把这个三角形缩小,你会得到一个从(0,0)到(1,0)到(1,1)再回到(0,0)的三角形,所有这些变换都会均匀地改变密度。你所做的就是在三角形中随机取一个点然后反过来得到圆中的一个点。
我不知道这个问题是否还有新的答案,但我自己碰巧也遇到过同样的问题。我试着跟自己“讲道理”寻找解决办法,我找到了一个。这可能和一些人在这里提出的建议是一样的,但不管怎样,它是这样的:
in order for two elements of the circle's surface to be equal, assuming equal dr's, we must have dtheta1/dtheta2 = r2/r1. Writing expression of the probability for that element as P(r, theta) = P{ r1< r< r1 + dr, theta1< theta< theta + dtheta1} = f(r,theta)*dr*dtheta1, and setting the two probabilities (for r1 and r2) equal, we arrive to (assuming r and theta are independent) f(r1)/r1 = f(r2)/r2 = constant, which gives f(r) = c*r. And the rest, determining the constant c follows from the condition on f(r) being a PDF.
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。
现在来算算:
The probability of choosing a radius between [0,r] is the integral of p(r) dr from 0 to r (that's just because we add all the probabilities of the smaller radii). Thus we want integral(p(r)dr) = r^2/R^2. We can clearly see that R^2 is a constant, so all we need to do is figure out which p(r), when integrated would give us something like r^2. The answer is clearly r * constant. integral(r * constant dr) = r^2/2 * constant. This has to be equal to r^2/R^2, therefore constant = 2/R^2. Thus you have the probability distribution p(r) = r * 2/R^2
Note: Another more intuitive way to think about the problem is to imagine that you are trying to give each circle of radius r a probability density equal to the proportion of the number of points it has on its circumference. Thus a circle which has radius r will have 2 * pi * r "points" on its circumference. The total number of points is pi * R^2. Thus you should give the circle r a probability equal to (2 * pi * r) / (pi * R^2) = 2 * r/R^2. This is much easier to understand and more intuitive, but it's not quite as mathematically sound.
程序员解决方案:
创建一个位图(布尔值的矩阵)。你想要多大就有多大。 在位图中画一个圆。 创建一个圆的点查找表。 在这个查找表中选择一个随机索引。
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]中选择r,然后计算你的坐标:
x = sqrt(r) * cos(angle)
y = sqrt(r) * sin(angle)
这就得到了圆盘上均匀的点分布。
http://mathworld.wolfram.com/DiskPointPicking.html
我仍然不确定确切的“(2/R2)×r”,但显而易见的是,在给定的单位“dr”中需要分配的点的数量,即r的增加将与R2成正比,而不是r。
check this way...number of points at some angle theta and between r (0.1r to 0.2r) i.e. fraction of the r and number of points between r (0.6r to 0.7r) would be equal if you use standard generation, since the difference is only 0.1r between two intervals. but since area covered between points (0.6r to 0.7r) will be much larger than area covered between 0.1r to 0.2r, the equal number of points will be sparsely spaced in larger area, this I assume you already know, So the function to generate the random points must not be linear but quadratic, (since number of points required to be distributed in given unit 'dr' i.e. increase in r will be proportional to r2 and not r), so in this case it will be inverse of quadratic, since the delta we have (0.1r) in both intervals must be square of some function so it can act as seed value for linear generation of points (since afterwords, this seed is used linearly in sin and cos function), so we know, dr must be quadratic value and to make this seed quadratic, we need to originate this values from square root of r not r itself, I hope this makes it little more clear.
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://stackoverflow.com/a/5838055/5224246从@sigfpe
设ρ(半径)和φ(方位角)是两个随机变量,对应于圆内任意一点的极坐标。如果这些点是均匀分布的,那么ρ和φ的分布函数是什么?
对于任意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随机变量根号(X),其中X在[0,1]上是一致的,有这个确切的CDF(因为P[根号(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)
你会得到
首先我们生成一个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])
如何在半径为R的圆内随机生成一个点:
r = R * sqrt(random())
theta = random() * 2 * PI
(假设random()均匀地给出0到1之间的值)
如果你想把它转换成笛卡尔坐标,你可以做到
x = centerX + r * cos(theta)
y = centerY + r * sin(theta)
为什么sqrt(随机())?
让我们看看sqrt(random())之前的数学运算。为简单起见,假设我们是在单位圆上工作,即R = 1。
点与点之间的平均距离应该是相同的,不管我们看的距离中心有多远。这意味着,例如,观察一个周长为2的圆的周长,我们应该找到的点的数量是周长为1的圆周长上点的数量的两倍。
由于圆的周长(2πr)随r线性增长,因此随机点的数量应该随r线性增长。换句话说,期望的概率密度函数(PDF)线性增长。由于PDF的面积应该等于1,最大半径是1,我们有
所以我们知道随机值的理想密度应该是什么样的。 现在:当我们只有一个0到1之间的均匀随机值时,我们如何生成这样一个随机值?
我们用了一个叫做反变换采样的技巧
从PDF中创建累积分布函数(CDF) 沿着y = x镜像 将得到的函数应用于0到1之间的统一值。
听起来复杂吗?让我插入一段带有小侧轨的引语来传达直觉:
Suppose we want to generate a random point with the following distribution: That is 1/5 of the points uniformly between 1 and 2, and 4/5 of the points uniformly between 2 and 3. The CDF is, as the name suggests, the cumulative version of the PDF. Intuitively: While PDF(x) describes the number of random values at x, CDF(x) describes the number of random values less than x. In this case the CDF would look like: To see how this is useful, imagine that we shoot bullets from left to right at uniformly distributed heights. As the bullets hit the line, they drop down to the ground: See how the density of the bullets on the ground correspond to our desired distribution! We're almost there! The problem is that for this function, the y axis is the output and the x axis is the input. We can only "shoot bullets from the ground straight up"! We need the inverse function! This is why we mirror the whole thing; x becomes y and y becomes x: We call this CDF-1. To get values according to the desired distribution, we use 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(random()) =√random()
这就是我们要推导的:-)
你也可以用你的直觉。
圆的面积是*r^2
为r = 1
得到的面积是。让我们假设我们有某种函数f,它将N=10个点均匀分布在一个圆内。这个比值是10 /
现在我们把面积和点数翻倍
对于r=2和N=20
面积是4,比值是20/4或10/2。半径越大,比值就越小,因为它的增长是二次的,N是线性的。
为了解决这个问题,我们可以说
x = r^2
sqrt(x) = r
如果在极坐标下生成一个向量
length = random_0_1();
angle = random_0_2pi();
更多的点会落在中心周围。
length = sqrt(random_0_1());
angle = random_0_2pi();
长度不再是均匀分布的,但是向量现在是均匀分布的。
这可能会帮助那些对选择速度算法感兴趣的人;最快的方法是(可能?)拒绝抽样。
只需在单位正方形内生成一个点,并拒绝它,直到它在圆内。如(伪代码),
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这样复杂的公式。