两个磁盘交叉区域的统一采样

Question

给定2D均匀变量如所讨论的,我们可以产生在一个单元盘的均匀分布在这里.

我的问题是在我希望均匀地采样两个相交的盘的交叉区域相似,其中一个磁盘总是单元盘和另一个可自由移动和调整大小等在这里

我试图将该区域分成两个区域(如上所述),并根据受尊重的磁盘对每个区域进行单独采样.我的方法基于上面引用的统一磁盘算法.为了对中心线右边的第一个区域进行采样,我将θ限制在两个交叉点内.接下来需要根据θ进行投影,以便将点推到我们的中线和磁盘半径之间的区域.python示例代码可以在这里找到.

u = unifrom2D()
A;B; // Intersection points
for p in allPoints
    theta = u.x * (getTheta(A) - getTheta(B)) + getTheta(B)
    r = sqrt(u.y + (1- u.y)*length2(lineIntersection(theta)))  
    p = (r * cos(theta), r * sin(theta))

然而,这种方法相当昂贵,并且进一步无法保持均匀性.只是为了澄清我不想使用拒绝抽样.

Answer 1

I am not sure if this is better than rejection sampling, but here is a solution for uniform sampling of a circle segment (with center angle <= pi) involving the numerical computation of an inverse function. (The uniform sampling of the intersection of two circles can then be composed of the sampling of segments, sectors and triangles - depending on how the intersection can be split into simpler figures.)

First we need to know how to generate a random value Z with given distribution F, i.e. we want

P(Z < x) = F(x)                   <=>  (x = F^-1(y))
P(Z < F^-1(y)) = F(F^-1(y)) = y   <=>  (F is monotonous)
P(F(Z) < y) = y

This means: if Z has the requested distribution F, then F(Z) is distributed uniformly. The other way round:

Z = F^-1(Y), 

where Y is distributed uniformly in [0,1], has the requested distribution.

If F is of the form

       / 0,                             x < a
F(x) = | (F0(x)-F0(a)) / (F0(b)-F0(a)), a <= x <= b
       \ 1,                             b < x

then we can choose a Y0 uniformly in [F(a),F(b)] and set Z = F0^-1(Y0).

We choose to parametrize the segment by (theta,r), where the center angle theta is measured from one segment side. When the segment's center angle is alpha, the area of the segment intersected with a sector of angle theta starting where the segment starts is (for the unit circle, theta in [0,alpha/2])

F0_theta(theta) = 0.5*(theta - d*(s - d*tan(alpha/2-theta)))

其中s = AB/2 = sin(alpha/2)和d = dist(M,AB) = cos(alpha/2)（圆心到线段的距离）。（这种情况alpha/2 <= theta <= alpha是对称的，这里不考虑。）我们需要一个随机theta数P(theta < x) = F_theta(x)。的逆F_theta不能以符号方式计算 - 它必须由某种优化算法（例如 Newton-Raphson）确定。

一旦theta确定，我们需要一个r范围内的随机半径

[r_min, 1], r_min = d/cos(alpha/2-theta).

因为x在[0, 1-r_min]分布中必须是

F0_r(x) = (x+r_min)^2 - r_min^2 = x^2 + 2*x*r_min.

这里可以用符号方式计算逆：

F0_r^-1(y) = -r_min + sqrt(r_min^2+y)

以下是用于概念验证的 Python 实现：

from math import sin,cos,tan,sqrt
from scipy.optimize import newton

# area of segment of unit circle
# alpha: center angle of segment (0 <= alpha <= pi)
def segmentArea(alpha):
  return 0.5*(alpha - sin(alpha))

# generate a function that gives the area of a segment of a unit circle
# intersected with a sector of given angle, where the sector starts at one end of the segment. 
# The returned function is valid for [0,alpha/2].
# For theta=alpha/2 the returned function gives half of the segment area.
# alpha: center angle of segment (0 <= alpha <= pi)
def segmentAreaByAngle_gen(alpha):
  alpha_2 = 0.5*alpha
  s,d = sin(alpha_2),cos(alpha_2)
  return lambda theta: 0.5*(theta - d*(s - d*tan(alpha_2-theta)))

# generate derivative function generated by segmentAreaByAngle_gen
def segmentAreaByAngleDeriv_gen(alpha):
  alpha_2 = 0.5*alpha
  d = cos(alpha_2)
  return lambda theta: (lambda dr = d/cos(alpha_2-theta): 0.5*(1 - dr*dr))()

# generate inverse of function generated by segmentAreaByAngle_gen
def segmentAreaByAngleInv_gen(alpha):
  x0 = sqrt(0.5*segmentArea(alpha)) # initial guess by approximating half of segment with right-angled triangle
  return lambda area: newton(lambda theta: segmentAreaByAngle_gen(alpha)(theta) - area, x0, segmentAreaByAngleDeriv_gen(alpha))

# for a segment of the unit circle in canonical position
# (i.e. symmetric to x-axis, on positive side of x-axis)
# generate uniformly distributed random point in upper half
def randomPointInSegmentHalf(alpha):
  FInv = segmentAreaByAngleInv_gen(alpha)
  areaRandom = random.uniform(0,0.5*segmentArea(alpha))
  thetaRandom = FInv(areaRandom)
  alpha_2 = 0.5*alpha
  d = cos(alpha_2)
  rMin = d/cos(alpha_2-thetaRandom)
  secAreaRandom = random.uniform(0, 1-rMin*rMin)
  rRandom = sqrt(rMin*rMin + secAreaRandom)
  return rRandom*cos(alpha_2-thetaRandom), rRandom*sin(alpha_2-thetaRandom)

可视化似乎验证了均匀分布（具有中心角 的线段的上半部分pi/2）：

import matplotlib.pyplot as plot
segmentPoints = [randomPointInSegmentHalf(pi/2) for _ in range(500)]
plot.scatter(*zip(*segmentPoints))
plot.show()