蒙特卡洛法求无穷区域积分：分布选择与误差求解问询

我在《蒙特卡洛方法》课程作业中遇到了一个难题：需要用MC方法计算积分区域为$[0,+\infty) \times [0,+\infty) \times [0,+\infty)$的积分近似值，并求出置信度0.99对应的误差值。目前我用标准正态分布生成样本时误差极大，想请教应该选择何种分布生成样本，以及背后的核心逻辑是什么？

我的现有R代码如下：

n <- 100000 
alfa <- 0.01 # 对应99%置信度

# 目前用标准正态分布生成样本，原以为和被积函数形状类似
gen <- function(n){
  return(matrix(rnorm(3*n, 0, 1),ncol=3))
}

g <- function(x){
  # 计算抽样分布的密度（这里错误地用了rnorm而不是dnorm）
  tihedus <- dnorm(x[,1],0,1)*dnorm(x[,2],0,1)*dnorm(x[,3],0,1)
  return( (x[,1]+x[,2])*exp(-(x[,1]+x[,2]+2*x[,3]))/(x[,1]^2+x[,2]+x[,3]+1) / tihedus * ((x[,1]>=0) & (x[,2]>=0) & (x[,3]>=0)) )
}

MC(gen, g, n, alfa)

问题分析：为什么当前方法误差极大？ #

你的代码里有两个关键问题：

1. 负样本浪费：标准正态分布会生成大量负数值样本，但积分区域只包含非负半轴，这些负样本会被过滤，实际有效样本量远小于n，直接降低了估计效率。
2. 分布匹配度差：被积函数的核心衰减项是$\exp(-(x_1+x_2+2x_3))$，而标准正态分布的衰减是$\exp(-x^2/2)$，两者尾部衰减速度完全不匹配，导致g(x)的方差极大，蒙特卡洛估计的误差自然会很高。

推荐抽样分布：独立指数分布 #

针对这个积分，最合适的抽样分布是独立的指数分布，具体选择：

• $x_1 \sim \text{Exp}(1)$，概率密度为$f_1(x_1) = \exp(-x_1),\ x_1 \geq 0$
• $x_2 \sim \text{Exp}(1)$，概率密度为$f_2(x_2) = \exp(-x_2),\ x_2 \geq 0$
• $x_3 \sim \text{Exp}(2)$，概率密度为$f_3(x_3) = 2\exp(-2x_3),\ x_3 \geq 0$

核心逻辑：重要性抽样的“匹配性”原则 #

蒙特卡洛重要性抽样的核心是：让抽样分布的概率密度与被积函数的形状（尤其是尾部衰减趋势）尽可能匹配。这样被积函数除以抽样密度后的函数$g(x)$的方差会最小化，从而大幅降低估计误差。

对于你的积分，被积函数的衰减项刚好和指数分布的密度形式高度契合：

• $x_1,x_2$方向的衰减是$\exp(-x_1)$和$\exp(-x_2)$，对应选择速率参数为1的指数分布；
• $x_3$方向的衰减是$\exp(-2x_3)$，对应选择速率参数为2的指数分布，完美匹配其衰减速度。

修改后的R代码 #

n <- 100000 
alfa <- 0.01 

# 生成独立指数分布样本：x1~Exp(1), x2~Exp(1), x3~Exp(2)
gen <- function(n){
  x1 <- rexp(n, rate = 1)
  x2 <- rexp(n, rate = 1)
  x3 <- rexp(n, rate = 2)
  return(cbind(x1, x2, x3))
}

g <- function(x){
  # 计算抽样分布的联合密度
  f_x <- dexp(x[,1], rate=1) * dexp(x[,2], rate=1) * dexp(x[,3], rate=2)
  # 被积函数除以密度，无需再过滤负样本（指数分布只生成非负值）
  integrand <- (x[,1]+x[,2])*exp(-(x[,1]+x[,2]+2*x[,3]))/(x[,1]^2+x[,2]+x[,3]+1)
  return(integrand / f_x)
}

# 蒙特卡洛估计计算
mc_estimates <- g(gen(n))
mean_est <- mean(mc_estimates)
std_est <- sd(mc_estimates)
# 99%置信度对应的误差（双侧分位数z=2.576）
error_margin <- qnorm(1 - alfa/2) * (std_est / sqrt(n))

cat("积分近似值：", mean_est, "\n")
cat("99%置信度误差：", error_margin, "\n")

补充说明 #

• 指数分布本身只生成非负样本，所以不需要再做x>=0的判断，避免了样本浪费；
• 匹配的抽样分布会让g(x)的数值波动大幅减小，样本标准差std_est会远小于之前用正态分布的情况，最终误差也会显著降低；
• 置信度0.99的误差计算采用正态近似，对应双侧分位数为qnorm(0.995)=2.576，乘以标准误（样本标准差除以根号n）即可得到误差边际。

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