Search code examples
rmontecarlogamma-distribution

Missing plot in gamma distribution with monte carlo approximation


i do not see why my plot is empty again, it worked and then i went over it to annotate it and it do not work anymore

thanks in advance


set.seed(1)
M_values=10^4
alpha <- 2
beta <- 3
a <- 1
b <- 3

# Generate a sequence of values for x
x <- seq(0, 6, length.out = 1000)

# Calculate the density of the Gamma distribution for each value of x
density <- dgamma(x, shape = alpha, rate = beta)

area_approximation <- numeric(length(M_values))#we have to create a vector to store our estimated areas for all the random points

# This is the monte carlo simulation, i used a loop to go over all the different random points
for (i in seq_along(M_values)) { #to iterates iover every m values
    M <- M_values[i]
    U <- runif(M, min = a, max = b) #random U generated within the interval [a,b]

    gamma_density <- dgamma(U, shape = alpha, rate = beta)#we have to calculate the density of the Gamma distribution at points U

    V <- runif(M, min = 0, max = max(density))#generates the random points on the y axis that does not go higher then the max density

    points_under_curve <- sum(V <= gamma_density)#number of points under the gamma curve
    

    area_approximation[i] <- (b - a) * max(density) * points_under_curve / M #approximation of the area under the Gamma curve
}

#plots estimated area vs. number of random points
plot(M_values, area_approximation, type = "l", col = "blue", lwd = 2, xlab = "Number of Random Points", ylab = "Estimated Area",main = "Convergence of Monte Carlo Estimation")

#true theoretical value
abline(h = theoretical_probability, col = "red", lty = 2)

legend("topright", legend = c("Estimated Area", "True Theoretical Value"),col = c("blue", "red"), lty = c(1, 2), lwd = c(2, 1), cex = 0.8)

empty plot

It worked then stop showing the plot and i do not know why I went over my entire code


Solution

  • If you want to iterate through all values from 1 to 10^4, you should assign

    M_values <- seq.int(10^4)
    

    rather than M_values <- 10^4

    and then you will obtain a plot like below

    enter image description here