Approximation of the area under the curve of the Gamma density far from the the theoretical value

I tried to set an approximation of the area under the Gamma density cuve, by first plot the density and add the points of the Monte Carlo algorithm (blue or red depending on their position to the Gamma curve). And finally evaluate the value of the area under the Gamma curve. But the results show that my approximation i far from the theoretical value and i do not know how to adjust it.

Knowing this : P(a≤X≤b)=P(X≤b)−P(X≤a) where X is a random variable with the Gamma distribution.

and this : P(a≤X≤b)≈(b−a)×(number of red points/total number of points)

Here is my code, i hope we can find the solution :)

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)

I tried increasing the number of simulations, look for other way to calculate the theoretical value....

There are two problems with your code:

1. The number of points under the curve cannot be computed as sum(U >= a & U <= b) (since that condition is equal true by the way you generate U). Instead, it must be computed as sum(V <= gamma_density) (i.e. the number of points for which the y coordinate lies under the curve);
2. The area approximation must take into account the height of the rectangle used to generate the random points (since that's not always 1 and you should apply a formula like area-of-the-rectangle times proportion of points under the curve). Therefore, you should compute (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

^{Created on 2024-02-18 with reprex v2.0.2}