使用Rejection Method从分布函数模拟数据的代码故障修复求助

首先，我得先假设你可能遇到的常见问题——比如目标分布未正确定义、接受条件逻辑错误、循环未累计足够样本，或者常数M的选取不合理。下面我会给出一个完整可运行的Python示例，同时拆解关键步骤帮你排查问题：

拒绝采样核心步骤回顾 #

拒绝采样的核心逻辑是：

• 从提议分布（这里是矩形内的均匀分布：x~U(0,2), y~U(0,2)）采样
• 计算接受概率：accept_prob = 目标分布f(x) / (M * 提议分布g(x))
  这里提议分布g(x)的密度是1/2（因为x在0-2均匀分布），如果我们取矩形上限y=2作为M*g(x)的上限，那么M*g(x)=2，接受概率可简化为f(x)/2（前提是f(x)在[0,2]内的最大值≤2）

修复后的完整代码示例 #

假设你的目标分布是归一化的概率密度函数，比如f(x) = 3x²/8（在[0,2]区间积分结果为1），代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 定义目标概率密度函数f(x)
def target_dist(x):
    # 这里替换成你实际的目标分布，确保在x∈[0,2]时有定义
    return (3 * x**2) / 8

# 拒绝采样函数
def rejection_sampling(n_samples):
    samples = []
    # 提议分布：x~U(0,2), y~U(0,2)
    while len(samples) < n_samples:
        # 从矩形区域采样
        x = np.random.uniform(0, 2)
        y = np.random.uniform(0, 2)
        # 接受条件：y ≤ 目标分布在x处的密度值
        if y <= target_dist(x):
            samples.append(x)
    return np.array(samples)

# 生成10000个样本
samples = rejection_sampling(10000)

# 可视化验证
plt.hist(samples, bins=50, density=True, alpha=0.7, label='采样结果')
x_range = np.linspace(0, 2, 100)
plt.plot(x_range, target_dist(x_range), 'r-', label='目标分布')
plt.legend()
plt.xlabel('x')
plt.ylabel('密度')
plt.title('拒绝采样结果验证')
plt.show()

常见故障排查点 #

• 目标分布未归一化/定义错误：如果你的目标分布不是概率密度函数，或者在[0,2]内的最大值超过了2（矩形的y上限），会导致接受概率大于1，或者几乎没有样本被接受。这时候你需要调整M的值，确保M*g(x)始终大于等于目标分布的最大值。
• 循环逻辑错误：比如没有用while循环直到收集到足够样本，而是固定次数采样，导致样本数量不足。
• 采样范围错误：比如x或y的采样范围不是[0,2]，导致无法覆盖目标分布的支持域。
• 接受条件写反：比如写成y >= target_dist(x)，这会导致错误的样本被接受。

你可以根据自己实际的目标分布，替换target_dist函数里的内容，确保逻辑符合上述步骤即可。

内容的提问来源于stack exchange，提问作者Apple_Tree