每日算法 — 使用java实现蒙特卡洛积分：概率数值积分与重要性采样策略

引言

蒙特卡洛方法（Monte Carlo Method）是一类基于随机采样和概率统计的数值计算技术，以其简洁的实现和强大的高维扩展能力，在科学计算、金融工程、物理学、机器学习等领域占据重要地位。蒙特卡洛积分是其中最基本的应用，它通过在被积函数的定义域内随机采样，用样本均值估计积分值。本文将深入讲解蒙特卡洛积分的数学原理，实现基本采样、分层采样和重要性采样三种策略，并提供完整的Java代码，帮助读者理解概率数值方法的精髓。

蒙特卡洛积分的数学原理

基本思想

设需要计算定积分：

$$I = \int_a^b f(x) \, dx$$

根据大数定律，若 $X_1, X_2, …, X_n$ 是在 $[a, b]$ 上独立均匀分布的随机变量，则：

$$I \approx (b – a) \cdot \frac{1}{n} \sum_{i=1}^{n} f(X_i)$$

当 $n \to \infty$ 时，估计值以概率1收敛到真实积分值。

误差分析

蒙特卡洛积分的标准误差为：

$$\sigma_{MC} = \frac{b – a}{\sqrt{n}} \cdot \sigma_f$$

其中 $\sigma_f$ 是被积函数的标准差。关键特性：
– 误差以 $O(1/\sqrt{n})$ 衰减，与维度无关
– 传统网格积分在高维时误差以 $O(n^{-1/d})$ 衰减，维度灾难严重
– 当维度 $d > 4$ 时，蒙特卡洛方法通常优于传统数值积分

高维推广

对于 d 维积分：

$$I = \int_{[a,b]^d} f(\mathbf{x}) \, d\mathbf{x}$$

蒙特卡洛估计为：

$$I \approx V \cdot \frac{1}{n} \sum_{i=1}^{n} f(\mathbf{x}_i)$$

其中 $V = (b-a)^d$ 是积分区域体积，$\mathbf{x}_i$ 是 d 维均匀随机向量。

采样策略

1. 基本随机采样

直接使用均匀分布的伪随机数进行采样，实现最简单。

2. 分层采样（Stratified Sampling）

将积分区域划分为若干子区域（层），在每层内独立采样。可降低方差：

$$\text{Var}{\text{stratified}} = \sum{k} \frac{V_k^2}{n_k} \sigma_k^2$$

当层划分合理时，分层采样的方差不大于基本随机采样。

3. 重要性采样（Importance Sampling）

当函数在某些区域变化剧烈时，均匀采样效率低下。重要性采样引入一个提议分布 $p(x)$，从 $p(x)$ 中采样：

$$I = \int_a^b f(x) \, dx = \int_a^b \frac{f(x)}{p(x)} \cdot p(x) \, dx \approx \frac{1}{n} \sum_{i=1}^{n} \frac{f(x_i)}{p(x_i)}$$

理想情况下取 $p(x) \propto |f(x)|$，可大幅减小方差。

Java实现

import java.util.Random;
import java.util.function.DoubleFunction;

/**
 * 蒙特卡洛积分实现
 * 包含基本采样、分层采样和重要性采样三种策略
 */
public class MonteCarloIntegration {

    private Random random;
    private static final long DEFAULT_SEED = 42L;

    public MonteCarloIntegration() {
        this(DEFAULT_SEED);
    }

    public MonteCarloIntegration(long seed) {
        this.random = new Random(seed);
    }

    /**
     * 基本蒙特卡洛积分 - 均匀随机采样
     *
     * @param f       被积函数
     * @param a       积分下限
     * @param b       积分上限
     * @param n       采样点数
     * @return 积分估计值
     */
    public double integrateBasic(DoubleFunction<Double> f, double a, double b, int n) {
        double sum = 0.0;
        double range = b - a;

        for (int i = 0; i < n; i++) {
            double x = a + random.nextDouble() * range;
            sum += f.apply(x);
        }

        return range * sum / n;
    }

    /**
     * 分层采样蒙特卡洛积分
     * 将区间划分为 numStrata 个等宽子区间，每个子区间均匀采样
     *
     * @param f          被积函数
     * @param a          积分下限
     * @param b          积分上限
     * @param n          总采样点数（将被均分到各层）
     * @param numStrata  层数
     * @return 积分估计值
     */
    public double integrateStratified(DoubleFunction<Double> f, double a, double b,
                                       int n, int numStrata) {
        double total = 0.0;
        double strataWidth = (b - a) / numStrata;
        int samplesPerStrata = n / numStrata;

        for (int k = 0; k < numStrata; k++) {
            double strataStart = a + k * strataWidth;
            double strataSum = 0.0;

            for (int i = 0; i < samplesPerStrata; i++) {
                double x = strataStart + random.nextDouble() * strataWidth;
                strataSum += f.apply(x);
            }

            total += strataWidth * strataSum / samplesPerStrata;
        }

        return total;
    }

    /**
     * 重要性采样蒙特卡洛积分
     * 使用提议分布进行有偏采样，通过权重修正估计值
     *
     * @param f           被积函数
     * @param sampler     从提议分布采样的函数
     * @param pdf         提议分布的概率密度函数
     * @param n           采样点数
     * @return 积分估计值
     */
    public double integrateImportance(DoubleFunction<Double> f,
                                       DoubleFunction<Double> sampler,
                                       DoubleFunction<Double> pdf,
                                       int n) {
        double sum = 0.0;

        for (int i = 0; i < n; i++) {
            double x = sampler.apply(random.nextDouble());
            double p = pdf.apply(x);
            if (Math.abs(p) < 1e-12) continue; // 避免除零
            sum += f.apply(x) / p;
        }

        return sum / n;
    }

    /**
     * 多维蒙特卡洛积分 - 基本采样
     *
     * @param f    d维被积函数
     * @param dim  维度
     * @param a    各维下限
     * @param b    各维上限
     * @param n    采样点数
     * @return 积分估计值
     */
    public double integrateMultiDimensional(MultiDimFunction f, int dim,
                                             double[] a, double[] b, int n) {
        double volume = 1.0;
        for (int d = 0; d < dim; d++) {
            volume *= (b[d] - a[d]);
        }

        double sum = 0.0;
        double[] x = new double[dim];

        for (int i = 0; i < n; i++) {
            for (int d = 0; d < dim; d++) {
                x[d] = a[d] + random.nextDouble() * (b[d] - a[d]);
            }
            sum += f.apply(x);
        }

        return volume * sum / n;
    }

    /**
     * 使用蒙特卡洛方法估算π值
     * 在单位正方形内随机撒点，统计落入四分之一圆的比例
     */
    public double estimatePi(int n) {
        int insideCircle = 0;

        for (int i = 0; i < n; i++) {
            double x = random.nextDouble();
            double y = random.nextDouble();
            if (x * x + y * y <= 1.0) {
                insideCircle++;
            }
        }

        return 4.0 * insideCircle / n;
    }

    /**
     * 并行蒙特卡洛积分（多线程加速）
     */
    public double integrateParallel(DoubleFunction<Double> f, double a, double b,
                                     int n, int numThreads) {
        long chunkSize = n / numThreads;

        return java.util.stream.IntStream.range(0, numThreads)
            .parallel()
            .mapToDouble(threadId -> {
                Random localRand = new Random(DEFAULT_SEED + threadId);
                double sum = 0.0;
                double range = b - a;
                for (long i = 0; i < chunkSize; i++) {
                    double x = a + localRand.nextDouble() * range;
                    sum += f.apply(x);
                }
                return range * sum / chunkSize;
            })
            .average()
            .orElse(0.0);
    }

    // 多维函数接口
    @FunctionalInterface
    interface MultiDimFunction {
        double apply(double[] x);
    }

    // ================== 测试与演示 ==================

    public static void main(String[] args) {
        MonteCarloIntegration mc = new MonteCarloIntegration();

        System.out.println("========== 蒙特卡洛积分演示 ==========\n");

        // 测试1：估算π值
        System.out.println("【π值估算】");
        for (int n : new int[]{1000, 10000, 100000, 1000000}) {
            double pi = mc.estimatePi(n);
            System.out.printf("  采样 %,d 次: π ≈ %.6f, 误差 = %.6f%n",
                n, pi, Math.abs(pi - Math.PI));
        }

        // 测试2：一维积分 ∫[0,1] x^2 dx = 1/3
        System.out.println("\n【一维积分：∫[0,1] x² dx = 1/3 ≈ 0.333333】");
        DoubleFunction<Double> f1 = x -> x * x;
        double exact1 = 1.0 / 3.0;

        for (int n : new int[]{1000, 10000, 100000}) {
            double basic = mc.integrateBasic(f1, 0, 1, n);
            double stratified = mc.integrateStratified(f1, 0, 1, n, 10);
            System.out.printf("  N=%,d | 基本采样: %.6f (误差%.6f) | 分层采样: %.6f (误差%.6f)%n",
                n, basic, Math.abs(basic - exact1),
                stratified, Math.abs(stratified - exact1));
        }

        // 测试3：重要性采样 ∫[0,1] e^x dx = e - 1 ≈ 1.718282
        System.out.println("\n【重要性采样：∫[0,1] e^x dx = e - 1 ≈ 1.718282】");
        DoubleFunction<Double> f3 = Math::exp;
        double exact3 = Math.E - 1;

        // 提议分布：使用线性递增分布 p(x) = 2x，CDF: F(x) = x²，逆变换: x = sqrt(u)
        DoubleFunction<Double> sampler = u -> Math.sqrt(u);
        DoubleFunction<Double> pdf = x -> 2 * x;

        for (int n : new int[]{1000, 10000, 100000}) {
            double basic = mc.integrateBasic(f3, 0, 1, n);
            double importance = mc.integrateImportance(f3, sampler, pdf, n);
            System.out.printf("  N=%,d | 基本采样: %.6f (误差%.6f) | 重要性采样: %.6f (误差%.6f)%n",
                n, basic, Math.abs(basic - exact3),
                importance, Math.abs(importance - exact3));
        }

        // 测试4：高维积分 ∫[0,1]^d (x1 + x2 + ... + xd) dx = d/2
        System.out.println("\n【高维积分：∫[0,1]^d Σx_i dx = d/2】");
        for (int dim : new int[]{2, 5, 10, 20}) {
            double[] a = new double[dim];
            double[] b = new double[dim];
            java.util.Arrays.fill(b, 1.0);
            double exact = dim / 2.0;

            MultiDimFunction multiDimF = x -> {
                double sum = 0;
                for (double v : x) sum += v;
                return sum;
            };

            double result = mc.integrateMultiDimensional(multiDimF, dim, a, b, 100000);
            System.out.printf("  维度 d=%2d: 估计值=%.4f, 精确值=%.4f, 误差=%.4f%n",
                dim, result, exact, Math.abs(result - exact));
        }

        // 测试5：并行加速
        System.out.println("\n【并行积分加速：∫[0,1] sin(x) dx = 1 - cos(1) ≈ 0.459698】");
        DoubleFunction<Double> f5 = Math::sin;
        double exact5 = 1 - Math.cos(1);
        int n = 10_000_000;
        long start = System.currentTimeMillis();
        double serial = mc.integrateBasic(f5, 0, 1, n);
        long serialTime = System.currentTimeMillis() - start;

        start = System.currentTimeMillis();
        double parallel = mc.integrateParallel(f5, 0, 1, n,
            Runtime.getRuntime().availableProcessors());
        long parallelTime = System.currentTimeMillis() - start;

        System.out.printf("  串行: %.6f (耗时 %d ms)%n", serial, serialTime);
        System.out.printf("  并行: %.6f (耗时 %d ms)%n", parallel, parallelTime);
        System.out.printf("  精确值: %.6f | 加速比: %.2fx%n", exact5,
            (double) serialTime / parallelTime);
    }
}

误差对比与收敛性

以下表格总结了不同采样策略在估算 $I = \int_0^1 e^x dx$ 时的表现：

应用场景

1. 金融工程：期权定价（Black-Scholes模型的蒙特卡洛模拟）、风险价值（VaR）计算
2. 物理学：统计力学配分函数、量子场论路径积分、中子输运模拟
3. 机器学习：贝叶斯推断中的后验采样、强化学习的策略评估
4. 计算机图形学：全局光照的路径追踪、体积渲染
5. 可靠性工程：复杂系统的失效概率分析
6. 生物信息学：蛋白质折叠的自由能计算

总结

蒙特卡洛积分是概率数值方法的典范，其核心优势在于维度无关的 $O(1/\sqrt{n})$ 收敛速度，使其成为高维积分的首选方法。本文实现了基本随机采样、分层采样和重要性采样三种策略，并展示了并行化加速技巧。理解蒙特卡洛方法不仅对于数值计算至关重要，更是深入理解随机算法、统计模拟和机器学习采样方法（如MCMC、变分推断）的基础。在实际应用中，选择合适的采样策略和方差缩减技术，往往能将计算效率提升一个数量级。