Python 拒绝采样算法优化与多峰分布模拟实战

自适应拒绝采样优化

针对基础方法中 M 值固定的问题，自适应版本在采样过程中根据已接受的样本重新估计分布均值和方差，动态调整提议分布，从而降低 M 值，提高接受率。

class AdaptiveRejectionSampler:
    def __init__(self, target_pdf, initial_proposal, domain=(-10, 10)):
        self.target_pdf = target_pdf
        self.proposal_pdf = initial_proposal['pdf']
        self.proposal_rv = initial_proposal['rv']
        self.domain = domain
        self.M_history = []

    def find_optimal_M(self, n_test=1000):
        test_points = np.linspace(self.domain[0], self.domain[1], n_test)
        ratios = [self.target_pdf(x) / self.proposal_pdf(x) for x in test_points]
        return np.max(ratios) * 1.1 if ratios else 2.0

    def adaptive_sampling(self, n_samples=1000, adapt_every=100):
        samples = []
        total_attempts = 0
        M = self.find_optimal_M()
        self.M_history.append(M)

        for batch in range(0, n_samples, adapt_every):
            batch_size = min(adapt_every, n_samples - len(samples))
            batch_samples, batch_attempts = self._sample_batch(batch_size, M)
            samples.extend(batch_samples)
            total_attempts += batch_attempts

            # 每接受一定数量样本后更新提议分布
            if len(samples) > 10:
                self._update_proposal(samples)
                M = self.find_optimal_M()
                self.M_history.append(M)

        acceptance_rate = len(samples) / total_attempts if total_attempts > 0 else 0
        return np.array(samples), acceptance_rate

    def _sample_batch(self, batch_size, M):
        samples = []
        attempts = 0
        while len(samples) < batch_size and attempts < batch_size * 100:
            attempts += 1
            x = self.proposal_rv()
            u = np.random.uniform(0, 1)
            if u < self.target_pdf(x) / (M * self.proposal_pdf(x)):
                samples.append(x)
        return samples, attempts

    def _update_proposal(self, accepted_samples):
        mean = np.mean(accepted_samples)
        std = np.std(accepted_samples) + 0.1
        self.proposal_pdf = lambda x: stats.norm.pdf(x, mean, std)
        self.proposal_rv = lambda: np.random.normal(mean, std)

分层与混合提议分布

对于多峰分布，单一的高斯提议往往难以兼顾所有峰值。分层采样将定义域切分，而混合提议则直接组合多个分布组件。以下是混合提议分布的实现逻辑：

class MixtureProposalSampler:
    def __init__(self, target_pdf, proposal_components):
        self.target_pdf = target_pdf
        self.components = proposal_components
        self.n_components = len(proposal_components)
        weights = np.array([comp['weight'] for comp in proposal_components])
        self.weights = weights / np.sum(weights)

    def sample(self, n_samples=1000):
        samples = []
        attempts = 0
        while len(samples) < n_samples and attempts < n_samples * 100:
            attempts += 1
            # 1. 根据权重选择提议分布组件
            comp_idx = np.random.choice(self.n_components, p=self.weights)
            component = self.components[comp_idx]
            # 2. 生成样本
            x = component['rv']()
            # 3. 计算混合提议分布的概率密度
            mixture_pdf = sum(w * c['pdf'](x) for w, c in zip(self.weights, self.components))
            # 4. 接受/拒绝决策
            if mixture_pdf > 0:
                ratio = self.target_pdf(x) / mixture_pdf
                M = 2.0
                u = np.random.uniform(0, 1)
                if u < ratio / M:
                    samples.append(x)
        acceptance_rate = len(samples) / attempts if attempts > 0 else 0
        return np.array(samples), acceptance_rate

性能对比与分析

运行上述代码对四种方法进行对比测试，结果如下表所示。实验基于 2000 个样本的生成任务。

从数据可以看出，混合提议分布在保持较高接受率的同时，显著降低了计算时间。这是因为混合分布能够更紧密地贴合双峰目标分布的形状，减少了无效采样。自适应方法虽然比基础方法好，但频繁更新分布参数的开销使其速度不如混合方法。分层方法由于需要遍历多个区间，在此场景下效率最低。

可视化结果显示，混合提议分布生成的样本直方图与目标分布曲线拟合度最高，尤其是在双峰位置的重合度上表现优异。

总结

在实际应用中，选择哪种拒绝采样策略取决于目标分布的特性：

• 单峰简单分布：基础拒绝采样通常足够，实现成本最低。
• 复杂多峰分布：推荐使用混合提议分布，它能灵活匹配多个峰值。
• 分布未知或动态变化：自适应拒绝采样能根据数据流自动调整，鲁棒性更强。

通过合理设计提议分布和优化 M 值，可以显著提升蒙特卡洛采样的效率，这对大模型训练中的采样环节尤为重要。