拒绝采样的概率论应用：从零开始理解

创作时间:

作者:

@小白创作中心

拒绝采样的概率论应用：从零开始理解

引用

网易

等

来源

https://fuxi.163.com/database/1112

https://blog.csdn.net/qq_51320133/article/details/137910261

https://cloud.baidu.com/article/3329831

https://blog.csdn.net/libing_zeng/article/details/81842245

https://blog.csdn.net/jteng/article/details/54344766

https://qianfanmarket.baidu.com/article/detail/1192937

https://blog.csdn.net/Anne033/article/details/109841951

https://blog.csdn.net/u010159842/article/details/78959515

在统计学和机器学习中，我们经常需要从一个特定的概率分布中生成随机样本。然而，对于一些复杂的分布，直接采样往往是不可能的。这时，一种名为“拒绝采样”（Rejection Sampling）的技术就派上了用场。本文将带你从零开始理解这一强大的统计工具。

从一个简单的例子说起

让我们从一个简单的例子开始：假设你想估算圆周率π的值。一个常见的方法是使用蒙特卡洛模拟：在一个边长为2的正方形内随机投点，正方形内有一个半径为1的圆。通过计算落在圆内的点的比例，我们可以估算出π的值。

然而，如果我们只允许在正方形的某些区域内投点，比如只能在正方形的左下角1/4区域内投点，这时该怎么办？这时，拒绝采样就派上用场了。我们可以在整个正方形内随机投点，但只保留（接受）那些落在允许区域内的点，其他点则被拒绝。通过这种方式，我们仍然可以准确估算出π的值。

这个例子展示了拒绝采样的基本思想：通过一个简单的分布（在这里是整个正方形内的均匀分布）来近似一个更复杂的分布（允许投点的区域），并通过接受或拒绝样本，最终得到符合目标分布的样本集。

拒绝采样的数学原理

在更一般的场景中，拒绝采样用于从一个复杂的目标分布p(x)中生成样本。其核心思想是使用一个简单的提议分布q(x)来生成候选样本，并根据一定的概率接受或拒绝这些样本。

具体来说，拒绝采样的步骤如下：

选择一个易于采样的提议分布q(x)，并找到一个常数k，使得对于所有的x，都有p(x) ≤ k*q(x)。
重复以下步骤直到获得足够的样本：
- 从提议分布q(x)中生成一个样本x'。
- 生成一个[0,1]区间内的随机数u。
- 计算接受概率α = p(x') / (k*q(x'))。
- 如果u ≤ α，则接受x'作为目标分布的样本；否则拒绝x'。

这个过程可以类比为在一张纸上画出目标分布的形状，然后用一个更大的、容易画的形状（提议分布）覆盖它。我们随机向大形状内扔纸团，但只保留那些落在目标形状内的纸团。通过这种方式，我们最终收集到的纸团位置就符合目标分布的形状了。

实际应用案例

为了更好地理解拒绝采样，让我们看一个具体的例子：假设我们需要从一个双峰分布中生成样本。这个分布由两个高斯分布叠加而成，直接采样非常困难。

我们可以选择一个标准差较大的正态分布作为提议分布，它的分布宽度足以覆盖目标分布。然后按照拒绝采样的步骤进行：

import numpy as np
from scipy.stats import norm

def target_distribution(x):
    """
    目标分布：一个双峰分布，由两个高斯分布叠加而成
    """
    return 0.75 * norm.pdf(x, loc=-2, scale=1) + 0.25 * norm.pdf(x, loc=2, scale=1)

def proposal_distribution(x):
    """
    提议分布：一个标准差为3的正态分布
    """
    return norm.pdf(x, loc=0, scale=3)

def rejection_sampling(n_samples, k=1.5):
    samples = []
    while len(samples) < n_samples:
        x_prime = np.random.normal(0, 3)  # 从提议分布中采样
        u = np.random.uniform(0, 1)       # 生成一个[0,1]内的随机数
        alpha = target_distribution(x_prime) / (k * proposal_distribution(x_prime))
        if u <= alpha:
            samples.append(x_prime)
    return samples

# 生成1000个样本
n_samples = 1000
rejection_samples = rejection_sampling(n_samples)

通过这段代码，我们可以从复杂的双峰分布中生成符合要求的样本，而无需直接对双峰分布进行采样。