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

1. 问题

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

2. 思路

简单说
(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)到这里就可以把问题解决了,下面是一些扩展知识。
(6)深入思考一下,为什么这样子做是有用的呢?下面详细阐述一下拒绝采样回到本题的一个思路历程。

拒绝采样(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。

时间复杂度:平均为O(1),最坏情况O(无穷)。
空间复杂度:O(1)。
rand7调用次数的期望值:2.45,计算如下

\[\begin{split}E &= 2 \cdot \frac{40}{49} + 4 \cdot \frac{9}{49} \cdot \frac{40}{49} + 6 \cdot (\frac{9}{49})^2 \cdot \frac{40}{49} + ... + 2k \cdot (\frac{9}{49})^{k-1} \cdot \frac{40}{49} \\ &= \frac{80}{49} ( 1 + 2 \cdot \frac{9}{49} + 3 \cdot (\frac{9}{49})^2 + ... + \ k \cdot (\frac{9}{49})^{k-1} ), \ k \to \infty \\ \frac{49}{80} E &= 1 + 2 \cdot \frac{9}{49} + 3 \cdot (\frac{9}{49})^2 + ... + \ k \cdot (\frac{9}{49})^{k-1} \\ ( \frac{49}{80} \cdot \frac{9}{49} ) \ E &= \frac{9}{49} + 2 \cdot (\frac{9}{49})^2 + 3 \cdot (\frac{9}{49})^3+ ... + \ (k-1) \cdot (\frac{9}{49})^{k-1} + k \cdot (\frac{9}{49})^{k} \\ ( \frac{49}{80} - \frac{49}{80} \cdot \frac{9}{49} ) \ E &= 1 + \frac{9}{49} + (\frac{9}{40})^2 + (\frac{9}{49})^3+ ... + (\frac{9}{49})^{k-1} - k \cdot (\frac{9}{49})^{k} \\ \frac{1}{2} \ E &= \frac{1 - (\frac{9}{49})^k}{1 - \frac{9}{49}} - k \cdot (\frac{9}{49})^{k}\end{split}\]

式子是一个等比数列乘以一个等差数列,使用错位相减,以及等比数列的求和公式得到最后的式子,k趋近于无穷时,\((\frac{9}{49})^k\)\(k \cdot (\frac{9}{49})^{k}\)都为0,最后得到\(E = 2 * \frac{49}{40} = 2.45\)

(方法二)利用被拒绝的数
在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时,只有一个数,这个时候就直接拒绝,继续按上面的一系列操作采样。

时间复杂度:平均为O(1),最坏情况O(无穷)。
空间复杂度:O(1)。
rand7调用次数的期望值:2.19

计算如下,方法和前面类似,使用错位相减和等比求和,以及极限求出,具体细节不多阐述。

\[\begin{split} E =\,& 2 \cdot \frac{40}{49} +3 \cdot \frac{9}{49} \cdot \frac{60}{63} + 4 \cdot \frac{9}{49} \cdot \frac{3}{63} \cdot \frac{20}{21} + \\ &\left(\frac{9}{49} \cdot \frac{3}{63} \cdot \frac{1}{21}\right) \times \left( 6 \cdot \frac{40}{49} + 7 \cdot \frac{9}{49} \cdot \frac{60}{63} + 8 \cdot \frac{9}{49} \cdot \frac{3}{63} \cdot \frac{20}{21} \right) + \\ &\left(\frac{9}{49} \cdot \frac{3}{63} \cdot \frac{1}{21}\right)^2 \times \left( 10 \cdot \frac{40}{49} + 11 \cdot \frac{9}{49} \cdot \frac{60}{63} + 12 \cdot \frac{9}{49} \cdot \frac{3}{63} \cdot \frac{20}{21} \right) + \\ &\ldots \\ &\left(\frac{9}{49} \cdot \frac{3}{63} \cdot \frac{1}{21}\right)^k \times \left( (2+4k) \cdot \frac{40}{49} + (3+4k) \cdot \frac{9}{49} \cdot \frac{60}{63} + (4+4k) \cdot \frac{9}{49} \cdot \frac{3}{63} \cdot \frac{20}{21} \right), k \to \infty \\ =\,& 2.19 \end{split}\]

3. 代码

# 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

4. 类似题目

478. Generate Random Point in a Circle
519. Random Flip Matrix

5. 参考资料

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

posted @ 2018-10-26 16:42 PilgrimHui 阅读(...) 评论(...) 编辑 收藏