蒙特卡洛积分效率翻倍：Halton与Hammersley低差异序列的Python实现详解

蒙特卡洛积分效率翻倍：Halton与Hammersley低差异序列的Python实现详解

在金融衍生品定价、物理仿真和机器学习等领域，高维积分计算一直是性能瓶颈。传统蒙特卡洛方法依赖伪随机采样，往往需要百万次迭代才能收敛。而低差异序列（Low-Discrepancy Sequences）通过数学构造的确定性采样，能在相同样本量下将误差降低40%-60%。本文将手把手教你用Python实现两种经典低差异序列——Halton与Hammersley，并揭示其在高维积分中的实战技巧。

1. 低差异序列的数学之美

1.1 从伪随机到确定性采样

伪随机数虽然满足统计均匀性，但在高维空间会出现聚类现象。下图对比了1000个二维采样点的分布：

import numpy as np import matplotlib.pyplot as plt # 伪随机采样 random_points = np.random.rand(1000, 2) # Halton序列 halton_points = np.array([halton_sequence(2, i) for i in range(1000)]) plt.figure(figsize=(12,5)) plt.subplot(121) plt.scatter(random_points[:,0], random_points[:,1], s=1) plt.title("伪随机采样") plt.subplot(122) plt.scatter(halton_points[:,0], halton_points[:,1], s=1) plt.title("Halton序列") plt.show()

关键差异：

• 伪随机点存在明显空白区域和聚集区
• Halton序列点间距均匀，覆盖率更高

1.2 基数逆序的魔法

Halton序列的核心是基数逆序函数。以基数b=3为例：

def radical_inverse(base, n): val, inv_base = 0.0, 1.0/base while n > 0: digit = n % base val += digit * inv_base n = n // base inv_base /= base return val

提示：基数应选择互质素数，常用前20个素数：[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71]

2. 高效实现技巧

2.1 维度灾难的破解之道

随着维度增加，Halton序列质量会下降。解决方案：

1. 素数筛选：跳过前k个素数（如跳过2,3,5）
2. 随机化：对序列施加随机平移

   def randomized_halton(dim, n, shift=None): if shift is None: shift = np.random.rand(dim) point = halton_sequence(dim, n) return (point + shift) % 1.0

3. 维度重排：交换高维坐标顺序

2.2 Hammersley的优化哲学

Hammersley序列在固定样本量时更高效：

def hammersley(dim, n, N): point = [n / float(N)] for i in range(1, dim): point.append(radical_inverse(PRIMES[i-1], n)) return point

性能对比（计算∫[0,1]^d sin(x1...xd) dx）：

3. 金融工程实战案例

3.1 亚式期权定价

使用Hammersley序列计算算术平均亚式看涨期权：

def asian_option(S0, K, r, sigma, T, d, N): dt = T/d discount = np.exp(-r*T) payoff = 0 for n in range(N): path = [] for i in range(d): z = hammersley(d, n, N)[i] phi = np.sqrt(2)*erfinv(2*z-1) # 逆正态变换 if i == 0: S = S0 * np.exp((r-0.5*sigma**2)*dt + sigma*np.sqrt(dt)*phi) else: S *= np.exp((r-0.5*sigma**2)*dt + sigma*np.sqrt(dt)*phi) path.append(S) avg = np.mean(path) payoff += max(avg - K, 0) return discount * payoff / N

参数调优指南：

• 当d>5时，建议使用随机化Halton
• 对于路径依赖期权，优先使用Hammersley
• 最优样本量N与维度d的关系：N ≈ 2^d

4. 超越均匀采样

4.1 球面均匀分布

将2D低差异序列映射到球面：

def spherical_hammersley(n, N): u, v = hammersley(2, n, N) theta = np.arccos(2*u - 1) # 极角 phi = 2 * np.pi * v # 方位角 x = np.sin(theta)*np.cos(phi) y = np.sin(theta)*np.sin(phi) z = np.cos(theta) return (x,y,z)

4.2 重要性采样结合

在概率密度高的区域增加采样：

def importance_halton(dim, n, pdf): point = halton_sequence(dim, n) jacobian = 1.0 for i in range(dim): point[i] = inverse_cdf(point[i], pdf) jacobian *= pdf(point[i]) return point, jacobian

在量子化学计算中，这种混合方法能将电子密度积分的收敛速度提升3倍以上。