我需要在半径为R的圆内生成一个均匀随机点。
我意识到,只要在区间[0…2π),且在区间(0…R)我最终会有更多的点指向中心,因为对于两个给定的半径,小半径内的点会比大半径内的点彼此更接近。
我在这里找到了一篇关于这个的博客文章,但我不明白他的推理。我认为这是正确的,但我真的很想知道他从哪里得到(2/R2)×r以及他如何得到最终解。
更新:在发布这个问题7年之后,关于平方根算法背后的数学问题,我仍然没有得到一个令人满意的答案。所以我花了一天时间自己写了一个答案。链接到我的答案。
我需要在半径为R的圆内生成一个均匀随机点。
我意识到,只要在区间[0…2π),且在区间(0…R)我最终会有更多的点指向中心,因为对于两个给定的半径,小半径内的点会比大半径内的点彼此更接近。
我在这里找到了一篇关于这个的博客文章,但我不明白他的推理。我认为这是正确的,但我真的很想知道他从哪里得到(2/R2)×r以及他如何得到最终解。
更新:在发布这个问题7年之后,关于平方根算法背后的数学问题,我仍然没有得到一个令人满意的答案。所以我花了一天时间自己写了一个答案。链接到我的答案。
当前回答
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;
其他回答
这里有一个快速而简单的解决方案。
在(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)的三角形,所有这些变换都会均匀地改变密度。你所做的就是在三角形中随机取一个点然后反过来得到圆中的一个点。
首先我们生成一个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°}
朴素解不起作用的原因是它给了靠近圆中心的点更高的概率密度。换句话说,半径为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];
} // ()
这可能会帮助那些对选择速度算法感兴趣的人;最快的方法是(可能?)拒绝抽样。
只需在单位正方形内生成一个点,并拒绝它,直到它在圆内。如(伪代码),
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这样复杂的公式。