蒙特卡洛积分中如何校正对数尺度/极坐标的非均匀采样？

一、极坐标蒙特卡洛积分的错误修正 #

问题根源 #

你的极坐标采样代码存在两个核心错误：

• 权重计算缺失：未结合极坐标的雅可比行列式（r）与采样概率密度的归一化项
• 缩放因子误用：错误用圆面积πR²作为最终缩放值，而非基于采样分布的概率密度计算权重

正确推导 #

对于圆域（半径R）的极坐标积分 I = ∫₀²π ∫₀ᴿ f(r cosφ, r sinφ) · r dr dφ：

• 采样r ~ Uniform(0, R)，概率密度 p(r) = 1/R
• 采样φ ~ Uniform(0, 2π)，概率密度 p(φ) = 1/(2π)
• 蒙特卡洛积分的单样本权重应为：f(...) · r / (p(r)·p(φ)) = f(...) · r · R · 2π
• 最终积分结果为：(1/N) × Σ[所有样本的权重项]

修正后的极坐标采样片段 #

if polar:
    for _ in range(size):
        r = np.random.random() * integration_radius
        phi = np.random.random() * 2. * np.pi
        x = r * np.cos(phi)
        y = r * np.sin(phi)
        # 正确权重：雅可比r × 采样概率密度的倒数(R*2π)
        weight = r * integration_radius * 2 * np.pi
        integral += function_to_integrate(x, y) * weight
    integral = integral / size

二、对数尺度采样的错误修正 #

问题根源 #

对数采样时，你误用原始区间上限1e7作为缩放因子，完全忽略了变量替换的雅可比行列式和采样概率密度的归一化要求。

正确推导 #

对于积分 I = ∫ₐᵇ x dx（a=1e-2, b=1e7），做对数变换 t = log₁₀(x)，则x=10ᵗ，dx=10ᵗ · ln(10) dt：

• 采样t ~ Uniform(-2, 7)，概率密度 p(t) = 1/(7 - (-2)) = 1/9
• 蒙特卡洛积分的单样本权重应为：x(t) · |dx/dt| / p(t) = 10ᵗ · 10ᵗ · ln(10) · 9
• 最终积分结果为：(1/N) × Σ[所有样本的权重项]

修正后的对数采样片段 #

if log:
    a_log = np.log10(1e-2)
    b_log = np.log10(1e7)
    dt = b_log - a_log
    for _ in range(size):
        t = np.random.uniform(a_log, b_log)
        x = 10 ** t
        # 正确权重：函数值×雅可比×采样区间长度（概率密度的倒数）
        jacobian = x * np.log(10)  # dx/dt = 10^t * ln10
        weight = x * jacobian * dt
        integral += weight
    integral = integral / size

三、完整修正后的代码 #

import numpy as np


def function_to_integrate(x, y):
    return np.exp(-x**2 - y**2)


def polar_MC(polar):
    size = 100000
    integral = 0.
    integration_radius = 4.
    if polar:
        for _ in range(size):
            r = np.random.random() * integration_radius
            phi = np.random.random() * 2. * np.pi
            x = r * np.cos(phi)
            y = r * np.sin(phi)
            # 极坐标下的正确权重：雅可比r × 采样概率密度的倒数(R*2π)
            weight = r * integration_radius * 2 * np.pi
            integral += function_to_integrate(x, y) * weight
        integral = integral / size
    else:
        length = 2. * integration_radius
        for _ in range(size):
            x = np.random.random() * length - length/2.
            y = np.random.random() * length - length/2.
            integral += function_to_integrate(x, y)
        integral = integral * length**2 / size
    print(f'POLAR: True integral should be pi ≈ {np.pi:.4f} ; MC: {integral:.4f}, polar={polar}')


def log_MC(log):
    size = 10000
    integral = 0.
    a = 1e-2
    b = 1e7
    true_value = 0.5 * (b**2 - a**2)
    if log:
        a_log = np.log10(a)
        b_log = np.log10(b)
        dt = b_log - a_log
        for _ in range(size):
            t = np.random.uniform(a_log, b_log)
            x = 10 ** t
            # 对数采样的正确权重：函数值×雅可比×采样区间长度
            jacobian = x * np.log(10)
            weight = x * jacobian * dt
            integral += weight
        integral = integral / size
    else:
        for _ in range(size):
            x = np.random.uniform(a, b)
            integral += x
        integral = integral * (b - a) / size
    print(f'LOG: True integral should be {true_value/1e13:.4f}*10^13 ; MC: {integral/1e13:.4f}*10^13, log={log}')


polar_MC(polar=True)
polar_MC(polar=False)

log_MC(log=True)
log_MC(log=False)

四、极坐标+对数半径采样的扩展思路 #

如果要在极坐标中对半径做对数采样，只需将极坐标中的r采样替换为对数采样：

1. 令t = log₁₀(r)，采样t在[log₁₀(r_min), log₁₀(r_max)]的均匀分布
2. 计算r=10ᵗ，雅可比行列式dr/dt = 10ᵗ · ln(10)
3. 极坐标的总雅可比为r · dr/dt = 10ᵗ · 10ᵗ · ln(10)
4. 结合φ的采样概率密度，最终权重为：f(r cosφ, r sinφ) · r · dr/dt · (dt_range) · 2π（dt_range是对数采样区间长度）

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