【随机算法】拒绝采样

摘要: 拒绝采样的算法原理和两道例题

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

---

在数值分析和计算统计中，拒绝采样是用于从分布生成观测值的基本技术，是一种精确的模拟方法。该方法适用于具有密度的任何分布。

拒绝采样可以这样直观理解：在一维中对随机变量进行采样，可以对二维笛卡尔图执行均匀随机采样，并将样本保持在其密度函数图下的区域中。

当然拒绝采样也有一定的问题，主要是下面两个：

1. 为了提高采样的接受率，需要使得容易实现的辅助分布尽可能贴近原分布，但这种分布实战中并不好找
2. 拒绝采样的接受率随着维度的增加而减小。因此拒绝采样不适合高维采样。

Metropolis Sampling 能够很好地解决这两个问题，以后我们再学习。下面我们通过两个力扣题目来看一下拒绝采样是如何应用的。

---

题目1

• 478. 在圆内随机生成点

给定圆的半径和圆心的 x、y 坐标，写一个在圆中产生均匀随机点的函数 randPoint 。

说明:

输入值和输出值都将是浮点数。
圆的半径和圆心的 x、y 坐标将作为参数传递给类的构造函数。
圆周上的点也认为是在圆中。
randPoint 返回一个包含随机点的x坐标和y坐标的大小为2的数组。

算法1: 直接采样

圆内部分的方程如下：

$$x^{2} + y^{2} = r^{2}$$

要在这个圆内等概率生成点，相当于有无数个圆环：

$$x^{2} + y^{2} = random * r^{2}$$

随机取到每个环的概率不可以相同，因为里面的环周长小，均匀取半径会造成里面的点密集，外面的点稀疏。

随机生成 $r1 \in [0, 1]$，然后随机生成的半径为 $r \times \sqrt{r1}$。

随机生成 $\theta \in [0, 2\pi)$。

代码 (C++)

class Solution {
public:
    Solution(double radius, double x_center, double y_center) {
        x = x_center;
        y = y_center;
        r = radius;
        dre = std::default_random_engine();
        d_theta = std::uniform_real_distribution<double>(0, 2 * PI);
        d_r = std::uniform_real_distribution<double>(0.0, std::nextafter(1.0, std::numeric_limits<double>::max()));
    }

    vector<double> randPoint() {
        double r = sqrt(d_r(dre)) * (this -> r);
        double theta = d_theta(dre);
        double x = r * cos(theta);
        double y = r * sin(theta);
        return {this -> x + x, this -> y + y};
    }

private:
    const long double PI = 3.14159265358979323846;
    double r, x, y;
    std::default_random_engine dre;
    std::uniform_real_distribution<double> d_theta;
    std::uniform_real_distribution<double> d_r;
};

算法2：拒绝采样

用边长 2R 的外接正方形随机生成点，若该点落在圆内，则返回该点，否则拒绝该点，重新生成点。

$$P(不被拒绝) = \frac{\pi R^{2}}{4R^{2}}$$

代码 (C++)

class Solution {
public:
    Solution(double radius, double x_center, double y_center) {
        r = radius;
        x = x_center;
        y = y_center;
        dre = std::default_random_engine();
        x1 = x_center - r;
        x2 = x_center + r;
        y1 = y_center - r;
        y2 = y_center + r;
        dr_x = std::uniform_real_distribution<double>(x1, std::nextafter(x2, std::numeric_limits<double>::max()));
        dr_y = std::uniform_real_distribution<double>(y1, std::nextafter(y2, std::numeric_limits<double>::max()));
    }

    vector<double> randPoint() {
        double i, j;
        while(true)
        {
            i = dr_x(dre);
            j = dr_y(dre);
            if(sqrt(pow((x - i), 2) + pow((y - j), 2)) <= r)
                break;
        }
        return vector<double>{i, j};
    }

private:
    double r, x, y;
    // 左下 (x1, y1) 右上 (x2, y2)
    double x1, y1, x2, y2;
    std::default_random_engine dre;
    std::uniform_real_distribution<double> dr_x;
    std::uniform_real_distribution<double> dr_y;
};

---

题目2

• 470. 用 Rand7() 实现 Rand10()

已有方法 rand7 可生成 1 到 7 范围内的均匀随机整数，试写一个方法 rand10 生成 1 到 10 范围内的均匀随机整数。

不要使用系统的 Math.random() 方法。

算法: 拒绝采样

调用两次 rand7()，得到 [1, 49] 之间的随机整数，其中 [1, 40] 接受，并映射到 [1, 10]，[41, 49] 这 9 个数字拒绝。

int x1 = rand7();
int x2 = rand7();
int x = x1 + (x2 - 1) * 7;

rand7() 期望调用次数

$$\begin{align} E &= 2 + 2 \times \frac{9}{49} + 2 \times (\frac{9}{49})^{2} + ... \\ &= 2 \times \sum\limits_{n=0}\limits^{\infty}(\frac{9}{49})^{n} \\ &= 2 \times \frac{1}{1 - \frac{9}{49}} \\ &= 2.45 \end{align}$$

代码 (C++)

class Solution {
public:
    int rand10() {
        while(true)
        {
            int x1 = rand7();
            int x2 = rand7();
            int x = x1 + (x2 - 1) * 7;
            if(x <= 40)
                return (x % 10) + 1;
        }
    }
};

优化接受率：利用被拒绝的数字减少 rand7() 调用次数

Ref: 官解。

调用两次 rand7() 得到 x1, x2，[41, 49] 中的数字 y 被拒绝。

再调用一次 x3 = rand7()，y 是 [1, 9] 随机数，x3 是 [1, 7] 随机数，x3 + (y - 1) * 7 可以得到 [1, 63] 的随机数，[61, 63] 的数字 z 被拒绝。

再调用一次 x4 = rand7()，z 是 [1, 3] 随机数，x4 是 [1, 7] 随机数，x4 + (z - 1) * 7 得到 [1, 21] 的随机数，[21, 21] 的数字 w 拒绝掉，w 没有利用价值了。

代码 (C++)

class Solution {
public:
    int rand10() {
        while (true) {
            int x1 = rand7();
            int x2 = rand7();
            int x = x2 + (x1 - 1) * 7;
            if (x <= 40)
                return 1 + x % 10;
            int y = x - 40;
            int x3 = rand7();
            // get uniform dist from 1 - 63
            x = x3 + (y - 1) * 7;
            if (x <= 60)
                return 1 + x % 10;
            int z = x - 60;
            int x4 = rand7();
            // get uniform dist from 1 - 21
            x = x4 + (z - 1) * 7;
            if (x <= 20)
                return 1 + x % 10;
        }
    }
};

---