两个磁盘交叉区域的统一采样
给定2D均匀变量,我们可以在单元磁盘中生成均匀分布,如此处所述。
我的问题是类似的,我希望统一采样两个交叉磁盘的交集区域,其中一个磁盘始终是单元磁盘,另一个可以像这里一样自由移动和调整大小
我试图将该区域分成两个区域(如上所述),并根据受尊敬的磁盘对每个区域进行抽样。 我的方法是基于上面引用的统一磁盘算法。 为了采样中心线右侧的第一个区域,我将限制theta在两个交点内。 下一个r需要根据theta进行投影,使得这些点被推入我们中线与磁盘半径之间的区域。 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))
然而,这种方法相当昂贵,并且未能保持一致性。 只是为了澄清我不想使用拒绝抽样。
我不确定这是否优于拒收采样,但是这是一个求解包含反函数数值计算的圆段(中心角<= pi)的均匀采样的解决方案。 (两个圆的交点的均匀采样可以由分段,扇区和三角形的采样组成 - 取决于交点如何分解成更简单的图形。)
首先,我们需要知道如何生成具有给定分布F
的随机值Z
,即我们想要的
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
这意味着:如果Z
有要求的分布F
,那么F(Z)
是均匀分布的。 反过来说:
Z = F^-1(Y),
其中Y
在[0,1]
均匀分布,具有所要求的分布。
如果F
是形式
/ 0, x < a
F(x) = | (F0(x)-F0(a)) / (F0(b)-F0(a)), a <= x <= b
1, b < x
那么我们可以在[F(a),F(b)]
均匀选择一个Y0
[F(a),F(b)]
并设置Z = F0^-1(Y0)
。
我们选择通过(theta,r)
参数化分段,其中中心角theta
从一个分段侧测量。 当分段的中心角为alpha
时,分段的面积与从分段开始的角度theta
开始的扇区相交(对于单位圆, [0,alpha/2]
theta
)
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).
对于[0, 1-r_min]
的x
,分布必须是
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()
链接地址: http://www.djcxy.com/p/40875.html