Gamma密度曲线下面积蒙特卡洛近似值与理论值偏差过大求助

我尝试用蒙特卡洛算法近似计算Gamma分布密度曲线在区间[a,b]下的面积：先绘制Gamma密度曲线，再用红点标记曲线下方的模拟点、蓝点标记曲线上方的点，最终估算曲线下面积。但近似结果和理论值偏差很大，试过增加模拟次数、更换理论值计算方式都没解决问题。

已知Gamma分布随机变量X满足$P(a≤X≤b)=P(X≤b)−P(X≤a)$，近似公式为$P(a≤X≤b)≈(b−a)×(红点数量/总点数)$，我的R代码如下：

set.seed(1)

M=10^3
alpha <- 2 # Shape parameter of the Gamma distribution
beta <- 3  # Rate parameter of the Gamma distribution
a <- 1     # Lower bound of the interval
b <- 3     # Upper bound of the interval

prob_a <- pgamma(a, shape = alpha, rate = beta)
prob_b <- pgamma(b, shape = alpha, rate = beta)

probability <- prob_b - prob_a

x <- seq(0, 6, length.out=1000)

density <- dgamma(x, shape = alpha, rate = beta)

plot(x, density, type = "l", col = "black", lwd = 2,
     xlab = "x", ylab = "Density",
     main = "Density of Gamma Distribution")

abline(v = a, col = "black", lwd = 1)
abline(v = b, col = "black", lwd = 1)

U <- runif(M, min = a, max = b)

gamma_density <- dgamma(U, shape = alpha, rate = beta)

V <- runif(M, min = 0, max = max(density))

points(U[V <= gamma_density], V[V <= gamma_density], col = "red", pch = 20, cex = 0.5)

points(U[V > gamma_density], V[V > gamma_density], col = "blue", pch = 20, cex = 0.5)

points_under_curve <- sum(U >= a & U <= b)

area_approximation <- (b - a) * points_under_curve / M

print(area_approximation)

theoretical_probability <- pgamma(b, shape = alpha, rate = beta) - pgamma(a, shape = alpha, rate = beta)

print(theoretical_probability)

问题根源 #

代码里的points_under_curve计算完全错误：U本身就是在[a,b]区间生成的均匀随机数，sum(U >= a & U <= b)的结果永远等于总点数M，导致近似面积变成(b-a)*1，和实际曲线下面积完全无关。正确的计算应该是统计落在曲线下方的点的数量，也就是满足V <= gamma_density的点数。

修正后的代码 #

set.seed(1)

M=10^4 # 可适当增加模拟次数提升精度
alpha <- 2
beta <- 3
a <- 1
b <- 3

# 计算理论概率
theoretical_probability <- pgamma(b, shape = alpha, rate = beta) - pgamma(a, shape = alpha, rate = beta)

# 绘制Gamma密度曲线
x <- seq(0, 6, length.out=1000)
density <- dgamma(x, shape = alpha, rate = beta)
plot(x, density, type = "l", col = "black", lwd = 2,
     xlab = "x", ylab = "Density",
     main = "Gamma Distribution & Monte Carlo Points")
abline(v = a, col = "gray", lty = 2)
abline(v = b, col = "gray", lty = 2)

# 蒙特卡洛模拟
U <- runif(M, min = a, max = b)
gamma_density <- dgamma(U, shape = alpha, rate = beta)
max_y <- max(density)
V <- runif(M, min = 0, max = max_y)

# 标记点
points(U[V <= gamma_density], V[V <= gamma_density], col = "red", pch = 20, cex = 0.5)
points(U[V > gamma_density], V[V > gamma_density], col = "blue", pch = 20, cex = 0.5)

# 计算近似面积/概率
points_under_curve <- sum(V <= gamma_density)
area_approximation <- (b - a) * (points_under_curve / M)

# 输出结果
cat("近似概率：", area_approximation, "\n")
cat("理论概率：", theoretical_probability, "\n")

额外调整建议 #

• 增加模拟次数（比如M=10^5）可以进一步缩小近似值和理论值的偏差
• 绘制曲线时用灰色虚线标记区间[a,b]，更清晰区分目标区域
• 用cat替代print输出结果，格式更易读

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