Monte_Carlo_integration_volume_of_n-ball

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
import time

# ========== 计算函数 ==========
def monte_carlo_hypersphere_volume(dim, num_samples):
    """蒙特卡洛体积计算函数"""
    points = np.random.uniform(-1, 1, (num_samples, dim))
    r_squared = np.sum(points**2, axis=1)
    count = np.sum(r_squared <= 1.0)
    return 2**dim * count / num_samples  # 立方体体积 × 内点比例

def theoretical_volume(dim):
    """理论体积公式"""
    return (np.pi ** (dim/2)) / gamma(1 + dim/2)

# ========== 计算参数 ==========
dims = np.arange(1, 11)   # 1-10维测试（更高维度需要极大样本）
num_samples = 10**6       # 每个维度采样数
mc_volumes = []           # 蒙特卡洛结果容器
true_volumes = []         # 理论值容器
durations = []            # 耗时监测

# ========== 执行计算 ==========
print(f"开始执行蒙特卡洛采样（维度 1-10，每维 {num_samples} 样本）")
for d in dims:
    start_time = time.time()

    # 计算理论与蒙特卡洛体积
    true_vol = theoretical_volume(d)
    est_vol = monte_carlo_hypersphere_volume(d, num_samples)

    # 记录结果
    true_volumes.append(true_vol)
    mc_volumes.append(est_vol)
    durations.append(time.time() - start_time)

    print(f"维度 {d:2d} | 理论值: {true_vol:.4f} | 蒙特卡洛: {est_vol:.4f} | 耗时: {durations[-1]:.2f}s")

# ========== 绘制结果 ==========
plt.figure(figsize=(10, 6))

# 理论值曲线（蓝色实线）
plt.plot(dims, true_volumes, 'b-o', linewidth=2, markersize=8,
         markerfacecolor='white', label='理论值')

# 蒙特卡洛结果（红色虚线）
plt.plot(dims, mc_volumes, 'r--s', linewidth=1.5, markersize=7,
         markerfacecolor='white', label='蒙特卡洛估计')

# 可视化装饰
plt.title("高维超球体积随维度变化", fontsize=14, pad=20)
plt.xlabel("维度 $d$", fontsize=12)
plt.ylabel("单位超球体积", fontsize=12)
plt.xticks(dims)
plt.yscale('log')  # 对数坐标展示指数衰减
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()

# 添加统计数据表
stats_text = (f"蒙特卡洛参数:\n"
             f"- 采样点/维度: {num_samples}\n"
             f"- 总耗时: {sum(durations):.1f}s\n"
             f"- 最大误差: {np.max(np.abs(np.array(mc_volumes)-np.array(true_volumes))/np.array(true_volumes))*100:.1f}%")
plt.gcf().text(0.72, 0.65, stats_text, fontsize=9,
               bbox=dict(facecolor='white', alpha=0.8))

plt.show()