土豆不吃鱼

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

问题描述

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

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

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

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

逆变换法

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

计算概率分布函数

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

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

直角坐标到极坐标的变换，雅可比行列式为$r$，得到极坐标下的概率密度函数：

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

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

$p(r) = \int _0 ^{2\pi} p(r,\theta) \mathrm d\theta=2r$

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

$p(\theta) = \frac{p(r, \theta)}{p(r)} = \frac{1}{2\pi}$

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

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

$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,y$直接从均匀分布采样是正方形，需要添加拒绝条件才行，如果使用极坐标变换，$r,\theta$分别从均匀分布采样，得到也是单位圆。这样为什么不行呢？

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

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

总结

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