470. Implement Rand10() Using Rand7() （拒绝采样Reject Sampling）

Posted 2021-01-15 liaohuiqiang

问题

已提供一个Rand7()的API可以随机生成1到7的数字，使用Rand7实现Rand10，Rand10可以随机生成1到10的数字。

思路

简单说：
（1）通过(Rand N - 1) % 10 + 1的方法，可以求出Rand10，当N是10的倍数的时候。
（2）用( Rand7 - 1 ) * 7 + Rand7可以随机生成1-49，记作Rand49。
（3）如果可以通过Rand49计算出Rand40，即随机生成1-40，就可以通过Rand40 % 10来取得Rand10。
（4）如何通过Rand49计算出Rand40呢？可以通过拒绝采样的方法来计算：用Rand49生成一个数，如果它位于41-49，则丢弃，继续生成，当生成一个1-40的数时返回。这样子做可以近似看成随机生成1-40之间的数。
（5）到这里就可以把代码写出来了，但是为什么这样子做是有用的呢？下面阐述一下从拒绝采样到本题的一个思路历程。

拒绝采样（Reject Sampling）
首先了解几个知识点：

连续分布中，累积分布函数CDF为(F(x) = P {X < x})，表示随机变量X落在区间((-infty, x))的概率。
概率密度函数PDF的意义在积分的时候得到体现，$ int_{a}^{b} f(x) dx = F(b) - F(a) = P { a le X lt b }$。

离散分布中，累积分布函数CDF一样有意义，(F(x) = P {X < x})表示小于x的那些(x_k)处的(p_k)之和。
概率质量函数PMF表示取得单个点的概率，(f(x) = P { X = x})。

拒绝采样基于这样一个前提： 对一个随机变量取样， 等价于从这个变量所服从的密度函数下方的区域均匀取样。

以下说的p(x)，q(x)指的是连续分布中的概率密度函数。我们用概率密度函数来表示一个概率分布。

假设我们想对p(x)分布采样，但是由于种种原因，我们很难对p(x)采样。但是另一个分布q(x)是可以采样的。那么我们可以根据q(x)的采样来得到一个在p(x)上的近似采样，具体做法如下。

（1）将q(x)和一个大于1的常数M相乘，使得M * q(x)可以把p(x)罩住，如下图所示（图源见参考文献）。

因为概率密度函数面积为1，所以q(x)肯定是有些地方比p(x)大，有些地方比q(x)小，如果全都比p(x)大，面积就超过1了。为了罩住p(x)，也就是处处比p(x)大，需要乘以一个大于1的常数M，这就是M * q(x)。

（2）从q(x)中获得一个采样点x(i)，对于x(i)计算一个接受概率 $ alpha = frac{p(x_i)}{Mq(x_i)}(，从Uniform(0,1)随机生成一个值)mu(，如果)alpha ge mu$，则接受x(i)作为一个来自p(x)的采样点。否则继续采样。

因为从q(x)分布取得这个点的概率是q(x)，并不是p(x)，（实际上对于连续分布，随机变量为单个数值的概率为0，但是这里我们讨论的是采样，从密度函数下方的区域均匀采样，这个场景下讨论时我们认为这个概率是有意义的值的而不是0），要把它变成p(x)需要多一个接受概率来进行选择，理想情况是选择p(x) / q(x)作为接收概率，因为q(x)乘以这个接受概率刚好就等于p(x)。

但是p(x) / q(x)是可能大于1的，正如前面所说p(x)有些地方比q(x)高，所以对q(x)乘上一个M，得到一个处于(0,1)的p(x) / (M * q(x))，可以作为接受概率，这个接受概率可以从一个(0,1)的均匀分布上取得。从图中可以看出，p(x)越接近M * q(x)的时候，接受概率越大，这很符合我们的采样直觉，因为M * q(x)正比于q(x)而M又是已经给出，所以p(x)越接近M * q(x)也表示p(x)越接近q(x)，这个时候它们的概率比较接近，当然要接受概率大一些。反过来，当p(x)和M * q(x)差很多时，p(x)和q(x)的差距也很大，这个时候根据q(x)采样到的样本很大概率要被拒绝。

最终采样的概率为(frac{p(x)}{M q(x)} q(x) = frac {1}{M} p(x))，其中(frac{p(x)}{M q(x)})可以从一个均匀分布采样，q(x)是已知的可以采样的分布，通过这两个分布采样可以得到一个近似p(x)的采样。显然，当我们取的M越接近1的时候，这个近似采样越接近p(x)上的采样。如果p(x)和q(x)完全一样，M就是1，在q上采样就等于在p上采样。

回到本题
Rand7以随机生成1-7，那么用( Rand7 - 1 ) * 7 + Rand7就可以随机生成1-49，记作Rand49()。
如果能够计算出Rand40，即随机生成1-40，就可以通过(Rand40 - 1) % 10 + 1的方法来取得Rand10。

那么知道Rand49的情况下如何计算出Rand40呢？这个场景就类似上面拒绝采样的场景：我们可以从q(x)分布上采样，如何模拟出一个p(x)上的近似采样。对离散分布当然也可以用拒绝采样，把上述的概率密度函数写成概率质量函数即可。

把Rand40对应的分布当作上述的p(x)，然后把Rand49对应的分布当作上述q(x)，然后M取49/40。
我们在q(x)上采样相当于用Rand49生成数字，每个数生成概率为 1/49，乘上M就是 1/40，对于每个随机变量，M * q(x)为1/40。p(x)是Rand40对应的分布，随机变量为1-40时，p(x) = 1/ 40，那么接受概率为p(x) / (M * q(x))就是1。随机变量为41-49时，p(x)为0，所以接受概率为0。

利用被拒绝的数
在Rand49取到41-49时，会被拒绝，那么可不可以在取到41-49时不拒绝，而是把这些数利用起来呢？

现在考虑取到的数是41-49之间的情况，减掉40规整为1-9，假设这个值为val，val处于1-9之间，考虑再用rand7()，构建一个[1, 63]的均匀分布，也就是 (val - 1) * rand7() + 7。再根据拒绝采样在[1, 63]的均匀分布中采1-60之间的值，然后用%10的方法得到Rand10。

接着上面的操作，如果采的数是61-63时，先减掉60规整为1-3，假设这个值为val，val处于1-3之间，考虑再用rand7()，构建出[1, 21]的均匀分布，也就是 (val - 1) * rand7() + 7。再根据拒绝采样采1-20之间的值，用%10的方法得到Rand10。

接着上面的操作，如果采的数是21时，只有一个数，这个时候就直接拒绝，继续按上面的一系列操作采样。

代码

# The rand7() API is already defined for you.
# def rand7():
# @return a random integer in the range 1 to 7

class Solution(object):
    def rand10(self):
        """
        :rtype: int
        """
        while True:
            val = (rand7() - 1) * 7 + rand7()
            if(val <= 40):
                break
        return (val-1) % 10 + 1

class Solution(object):
    def rand10(self):
        """
        :rtype: int
        """
        while True:
            val = (rand7() - 1) * 7 + rand7()
            if val <= 40:
                return (val-1) % 10 + 1
            val -= 40
            val = (val - 1) * 7 + rand7()
            if val <= 60:
                return (val-1) % 10 + 1
            val -= 60
            val = (val - 1) * 7 + rand7()
            if val <= 20:
                return (val-1) % 10 + 1

参考资料

蒙特卡洛采样之拒绝采样-csdn