算法概述
拒绝采样（Rejection Sampling）是一种从复杂概率分布中生成随机样本的基础蒙特卡洛方法。当目标分布的概率密度函数复杂、难以直接采样时，可以通过一个易于采样的提议分布（Proposal Distribution）配合接受-拒绝机制，间接获得服从目标分布的样本。
本文将用Java实现拒绝采样算法，涵盖均匀分布、正态分布、Beta分布等典型场景的采样实现，并分析其效率与适用范围。
核心原理
设目标分布的概率密度函数为 (f(x))，提议分布为 (g(x))，且存在常数 (M) 使得对所有 (x) 满足 (f(x) \leq M \cdot g(x))。拒绝采样的步骤如下：
- 从提议分布 (g(x)) 中采样得到一个候选值 (x)
- 从均匀分布 (U(0, 1)) 中采样得到 (u)
- 若 (u \leq \frac{f(x)}{M \cdot g(x)})，则接受 (x) 作为样本；否则拒绝并重复
接受概率为 (\frac{1}{M})，因此 (M) 越接近1，采样效率越高。
Java实现
项目结构
rejection-sampling/
├── RejectionSampler.java // 拒绝采样核心框架
├── UniformSampler.java // 均匀分布采样示例
├── NormalSampler.java // 正态分布采样示例
├── BetaSampler.java // Beta分布采样示例
└── Main.java // 测试入口
核心框架：RejectionSampler.java
import java.util.Random;
import java.util.function.DoubleFunction;
/**
* 拒绝采样通用框架
* 通过提议分布和接受-拒绝机制从目标分布中采样
*/
public class RejectionSampler {
private final Random random;
public RejectionSampler() {
this.random = new Random();
}
public RejectionSampler(long seed) {
this.random = new Random(seed);
}
/**
* 执行一次拒绝采样
* @param proposalSampler 提议分布采样器（返回候选值）
* @param targetDensity 目标分布概率密度函数
* @param proposalDensity 提议分布概率密度函数
* @param m 上界常数 M，需满足 target(x) <= M * proposal(x)
* @return 采样值
*/
public double sample(
java.util.function.Supplier<Double> proposalSampler,
DoubleFunction<Double> targetDensity,
DoubleFunction<Double> proposalDensity,
double m) {
if (m <= 0) {
throw new IllegalArgumentException("M必须为正数");
}
while (true) {
// 步骤1：从提议分布采样
double x = proposalSampler.get();
// 步骤2：生成均匀随机数
double u = random.nextDouble();
// 步骤3：计算接受概率
double acceptanceProb = targetDensity.apply(x) / (m * proposalDensity.apply(x));
// 接受或拒绝
if (u <= acceptanceProb) {
return x;
}
}
}
/**
* 批量采样
*/
public double[] sampleBatch(
java.util.function.Supplier<Double> proposalSampler,
DoubleFunction<Double> targetDensity,
DoubleFunction<Double> proposalDensity,
double m,
int n) {
double[] samples = new double[n];
for (int i = 0; i < n; i++) {
samples[i] = sample(proposalSampler, targetDensity, proposalDensity, m);
}
return samples;
}
/**
* 有界拒绝采样：在区间[a, b]内采样
* 使用均匀分布作为提议分布
*/
public double sampleBounded(
DoubleFunction<Double> targetDensity,
double a, double b,
double maxDensity) {
if (a >= b || maxDensity <= 0) {
throw new IllegalArgumentException("参数非法");
}
double interval = b - a;
// 均匀分布密度为 1/(b-a)，因此 M = maxDensity * (b-a)
double m = maxDensity * interval;
while (true) {
double x = a + random.nextDouble() * interval;
double u = random.nextDouble();
if (u <= targetDensity.apply(x) / maxDensity) {
return x;
}
}
}
}
正态分布采样：NormalSampler.java
/**
* 使用拒绝采样从标准正态分布 N(0,1) 中采样
* 提议分布使用柯西分布（Cauchy），因其尾部与正态分布相似
*/
public class NormalSampler {
private final RejectionSampler sampler;
private final java.util.Random random;
public NormalSampler() {
this.sampler = new RejectionSampler();
this.random = new java.util.Random();
}
/**
* 标准正态分布概率密度函数
*/
public static double normalPdf(double x) {
return Math.exp(-x * x / 2.0) / Math.sqrt(2.0 * Math.PI);
}
/**
* 柯西分布概率密度函数
*/
public static double cauchyPdf(double x) {
return 1.0 / (Math.PI * (1.0 + x * x));
}
/**
* 柯西分布采样：使用逆变换法
* CDF的逆函数为 tan(π * (u - 0.5))
*/
public double sampleCauchy() {
double u = random.nextDouble();
return Math.tan(Math.PI * (u - 0.5));
}
/**
* 通过拒绝采样获取标准正态分布样本
* 理论最优M = sqrt(2π/e) ≈ 1.52
*/
public double sampleStandardNormal() {
double m = Math.sqrt(2.0 * Math.PI / Math.E); // ≈ 1.52
return sampler.sample(
this::sampleCauchy,
NormalSampler::normalPdf,
NormalSampler::cauchyPdf,
m
);
}
/**
* 一般正态分布 N(mu, sigma^2)
*/
public double sampleNormal(double mu, double sigma) {
return mu + sigma * sampleStandardNormal();
}
public static void main(String[] args) {
NormalSampler sampler = new NormalSampler();
// 生成10000个样本并统计均值和方差
int n = 10000;
double sum = 0, sumSq = 0;
double[] samples = new double[n];
for (int i = 0; i < n; i++) {
samples[i] = sampler.sampleStandardNormal();
sum += samples[i];
sumSq += samples[i] * samples[i];
}
System.out.printf("样本均值: %.4f (理论: 0)%n", sum / n);
System.out.printf("样本方差: %.4f (理论: 1)%n", sumSq / n - (sum / n) * (sum / n));
}
}
Beta分布采样：BetaSampler.java
/**
* Beta分布拒绝采样实现
* Beta分布在[0,1]区间有界，适合使用均匀分布作为提议分布
*/
public class BetaSampler {
private final RejectionSampler sampler;
public BetaSampler() {
this.sampler = new RejectionSampler();
}
/**
* Beta分布概率密度函数
* f(x; alpha, beta) = x^(alpha-1) * (1-x)^(beta-1) / B(alpha, beta)
*/
public static double betaPdf(double x, double alpha, double beta) {
if (x <= 0 || x >= 1) return 0;
// 使用对数避免溢出
double logBeta = logGamma(alpha) + logGamma(beta) - logGamma(alpha + beta);
double logPdf = (alpha - 1) * Math.log(x) + (beta - 1) * Math.log(1 - x) - logBeta;
return Math.exp(logPdf);
}
/**
* 伽马函数的对数近似（Stirling近似）
* 对于正整数和半整数足够精确
*/
private static double logGamma(double z) {
if (z <= 0) throw new IllegalArgumentException("z必须为正");
// Lanczos近似简化版
if (z < 1) return logGamma(z + 1) - Math.log(z);
if (z == 1 || z == 2) return 0;
double x = z - 1;
double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
double a = p[0];
for (int i = 1; i < p.length; i++) {
a += p[i] / (x + i);
}
double t = x + p.length - 0.5;
return 0.5 * Math.log(2 * Math.PI) + Math.log(a) - t + (x + 0.5) * Math.log(t);
}
/**
* 采样Beta(alpha, beta)
* 使用均匀分布作为提议分布，M取PDF在众数处的值
*/
public double sampleBeta(double alpha, double beta) {
if (alpha <= 0 || beta <= 0) {
throw new IllegalArgumentException("alpha和beta必须为正");
}
// Beta分布的众数：(alpha-1)/(alpha+beta-2)，当alpha,beta>1时
double mode;
if (alpha > 1 && beta > 1) {
mode = (alpha - 1) / (alpha + beta - 2);
} else if (alpha <= 1 && beta <= 1) {
mode = 0.5; // 双凹，取中点
} else if (alpha <= 1) {
mode = 0.01; // 右偏，靠近0
} else {
mode = 0.99; // 左偏，靠近1
}
double maxDensity = betaPdf(mode, alpha, beta);
double m = maxDensity; // 均匀分布在[0,1]的密度为1
return sampler.sampleBounded(
x -> betaPdf(x, alpha, beta),
0, 1, maxDensity
);
}
public static void main(String[] args) {
BetaSampler sampler = new BetaSampler();
double alpha = 2.0, beta = 5.0;
int n = 5000;
double sum = 0;
for (int i = 0; i < n; i++) {
double s = sampler.sampleBeta(alpha, beta);
sum += s;
}
System.out.printf("Beta(%.1f, %.1f) 样本均值: %.4f (理论: %.4f)%n",
alpha, beta, sum / n, alpha / (alpha + beta));
}
}
测试入口：Main.java
public class Main {
public static void main(String[] args) {
System.out.println("=== 拒绝采样算法测试 ===\n");
// 测试1：标准正态分布
System.out.println("--- 标准正态分布 N(0,1) ---");
NormalSampler normalSampler = new NormalSampler();
int n = 10000;
double[] normalSamples = new double[n];
double sum = 0, sumSq = 0;
for (int i = 0; i < n; i++) {
normalSamples[i] = normalSampler.sampleStandardNormal();
sum += normalSamples[i];
sumSq += normalSamples[i] * normalSamples[i];
}
double mean = sum / n;
double var = sumSq / n - mean * mean;
System.out.printf("样本量: %d, 均值: %.4f, 方差: %.4f%n%n", n, mean, var);
// 测试2：Beta分布
System.out.println("--- Beta(2, 5) 分布 ---");
BetaSampler betaSampler = new BetaSampler();
sum = 0;
for (int i = 0; i < n; i++) {
sum += betaSampler.sampleBeta(2.0, 5.0);
}
System.out.printf("样本量: %d, 均值: %.4f (理论: 0.2857)%n%n", n, sum / n);
// 测试3：有界采样 - 半圆分布
System.out.println("--- 半圆分布（Wigner半圆律）---");
RejectionSampler sampler = new RejectionSampler();
sum = 0; sumSq = 0;
for (int i = 0; i < n; i++) {
double x = sampler.sampleBounded(
x -> Math.sqrt(Math.max(0, 1 - x * x)) * 2 / Math.PI,
-1, 1, 2.0 / Math.PI
);
sum += x;
sumSq += x * x;
}
System.out.printf("样本量: %d, 均值: %.4f, 二阶矩: %.4f%n", n, sum / n, sumSq / n);
}
}
复杂度分析
| 指标 | 复杂度 |
|---|---|
| 单次采样期望时间 | (O(M))，其中 (M) 为接受概率倒数 |
| 空间复杂度 | (O(1))，仅需存储常数和随机数生成器 |
| 批量采样 (n) 个 | (O(nM)) 期望时间 |
拒绝采样的核心瓶颈在于接受率。当目标分布与提议分布差异较大时（(M \gg 1)），大量候选样本被拒绝，效率急剧下降。在高维空间中，(M) 通常随维度指数增长，导致拒绝采样在高维场景下失效——这被称为”维度诅咒”。
适用场景与局限性
适用场景：
– 低维概率分布（1-3维）的随机采样
– 概率密度函数已知但累积分布函数难以求逆的情况
– 需要精确服从目标分布的蒙特卡洛模拟
局限性：
– 高维空间中效率极低
– 需要已知目标分布的密度上界 (M)
– 对于重尾分布，难以找到合适的提议分布
扩展方向
- 自适应拒绝采样：动态调整提议分布以提高接受率
- 切片采样：将多维采样转化为均匀分布采样
- Metropolis-Hastings算法：马尔可夫链蒙特卡洛方法的推广
总结
拒绝采样是理解蒙特卡洛方法的入门算法。通过Java实现，我们掌握了从复杂分布中生成样本的核心思想：用简单的提议分布包裹复杂的目标分布，以概率接受机制保证样本的正确性。尽管在高维场景下存在局限，但在低维概率建模、A/B测试、贝叶斯推断等领域，拒绝采样仍是一种简洁而有效的工具。
