R 远离理论值的伽马密度曲线下面积的近似

您的代码有两个问题:

1. 曲线下的点数不能计算为sum(U >= a & U <= b)(因为按照您生成U的方式，该条件等于真).相反，它必须计算为sum(V <= gamma_density)(即y坐标位于曲线下方的点数)；
2. 面积近似必须考虑用于生成随机点的矩形的高度(因为这并不总是1，您应该应用类似于area-of-the-rectangle乘以proportion of points under the curve的公式).因此，您应该计算(b - a) * max(density) * points_under_curve / M.

Reprex:

set.seed(1)

M=10^5
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 = ".", cex = 0.5)

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

points_under_curve <- sum(V <= gamma_density)

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

print(area_approximation)
#> [1] 0.1975212

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

print(theoretical_probability)
#> [1] 0.1979142

^{创建于2024-02-18年第reprex v2.0.2页}