从单位圆均匀采样

采样是从某种概率分布中获取一定的样本。在计算机如何产生随机数中讲了如何得到均匀分布的样本,也是一种均匀采样的方法;如何得到服从正态分布的样本中,讲了从正态分布采样的方法。前面也总结了常用的概率分布,其实利用逆变换法可以得到服从已知概率分布的样本,在前面的文章已经介绍过了。本文继续介绍逆变化法的应用,从单位圆均匀采样。

问题描述

之前我们从0-1区间内采样得到均匀的样本,这是一维问题,现在需要从单位圆盘上均匀采样,意味着单位圆内,每个点被采样的概率相等。如果两个坐标点x,yx,y分别从U(1,1)U(-1,1)均匀采样,其实得到是从正方形采样的样本。这肯定不满足要求,可能有的读者朋友觉得,加上一个判断条件,如果满足在圆内就接受,这就是我们需要的样本点,如果不在圆内,那就拒绝,继续采样。

如果你是这样的想的,恭喜你,你已经找到了一种非常简单的采样方法,称为“拒接采样”或者“接受-拒绝采样”。我们再试着计算采样效率,也就是落在单位圆内概率

r==π4=78.5%r = \frac{圆的面积}{正方形面积} = \frac{\pi}{4}=78.5\%

看起来采样的效率并不是很高,这也是拒绝采样的弊端,读者应该能实现拒绝采样的代码,如果有疑问可以在公众号后台留言。

逆变换法

这里重复一下逆变换定理:设XX为连续型随机变量,取值于区间(a,b)(a,b)(可包括\infty和端点), XX的概率密度在(a,b)(a,b)上取正值, F(X)F(X)的分布函数为F(X),UU(0,1)F(X),U \sim \text{U}(0,1)Y=F1(U)F()Y = F^{-1}(U) \sim F(\cdot),逆变换定理用语言来叙述就是,找到随机变量的概率分布函数,其反函数就是我们要找的变换,这个变换能把均匀分布映射到已知的概率分布上。

计算概率分布函数

假设概率密度函数为p(x,y)p(x,y),由于是从单位圆均匀采样,则概率密度函数应该是常数,每个点被采样的概率都是相等。根据概率归一化p(x,y)dxdy=1\int p(x, y) \mathrm dx \mathrm dy = 1,很容易得到p(x,y)=1πp(x,y)=\frac{1}{\pi},单位圆内每个点被采样的概率为1π\frac{1}{\pi}.

接下来做极坐标变换,由于极坐标中r,θr,\theta是相互独立的,比较容易分别得到对应的概率分布函数,根据逆变换定理,就能得到r,θr,\theta的变换函数,这个变换就是映射均匀采样到r,θr,\theta的采样。这里的思路需要理解清楚。

直角坐标到极坐标的变换,雅可比行列式为rr,得到极坐标下的概率密度函数:

p(r,θ)=rp(x,y)=rπp(r,\theta) = rp(x,y) = \frac{r}{\pi}

下面就是根据概率公式计算了,分别计算边缘概率:

p(r)=02πp(r,θ)dθ=2rp(r) = \int _0 ^{2\pi} p(r,\theta) \mathrm d\theta=2r

概率分布函数很容易得到Fr(r)=r2F_r(r) = r^2,根据逆变换定理极径rr的采样为r=Fr1(a)=ar = F_r^{-1}(a) = \sqrt a,其中aU(0,1)a\sim U(0,1),同理计算θ\theta的采样函数。

p(θ)=p(r,θ)p(r)=12πp(\theta) = \frac{p(r, \theta)}{p(r)} = \frac{1}{2\pi}

容易得到θ\theta的概率分布函数Fθ=θ2πF_\theta = \frac{\theta}{2\pi}, θ=Fθ1(b)=2πb\theta = F_\theta ^{-1} (b) = 2\pi b,其中bU(0,1)b\sim U(0, 1).

上面已经分别得到了r,θr,\theta的变换函数,根据极坐标公式于是得到单位圆的均匀采样公式:

x=rcosθ=acos(2πb)y=rsinθ=asin(2πb) x=r\cos\theta = \sqrt a \cos (2\pi b)\\ y=r\sin \theta = \sqrt a \sin(2\pi b)

def sampling_circle(sample_size):
    a = np.random.uniform(size=sample_size)
    b = np.random.uniform(size=sample_size)
    x = np.sqrt(a) * np.cos(2 * np.pi * b)
    y = np.sqrt(a) * np.sin(2 * np.pi * b)
    return x, y

下图展示几种采样图,直观来看分布是比较均匀的。

为什么不直接采样

x,yx,y直接从均匀分布采样是正方形,需要添加拒绝条件才行,如果使用极坐标变换,r,θr,\theta分别从均匀分布采样,得到也是单位圆。这样为什么不行呢?

先来直观看看,如果rU(0,1),θU(0,2π)r\sim U(0,1),\theta\sim U(0, 2\pi)从均匀分布采样,结果如何。

从上面的图可以看出,中心点更加密集,说明越靠近中心采样的概率更大一些。这是因为我们是均匀的采样rr,对于采样得到的r0r_0来说,其周长为2πr02\pi r_0,意味着不同周长采样概率是相等的。很明显周长越长,应该采样更多的点,这样才能保证不同周长的圆,被采样的概率相同。采样概率与极径有关。所以r,θr,\theta不能直接从均匀分布采样。

总结

本文介绍了如何从单位圆均匀采样的方法,提到了拒绝采样法,这是一种自然直观的方法,由于方法相对简单,就不做过多介绍了。重点介绍逆变换法,应用逆变换定理的只需要找到概率分布函数,其反函数就是我们需要逆变换。掌握了此方法,我们也能从单位球体均匀采样,球面均匀采样,三角形中均匀采样等等。