Follow

Keep Up to Date with the Most Important News

By pressing the Subscribe button, you confirm that you have read and are agreeing to our Privacy Policy and Terms of Use
Contact

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

MEDevel.com: Open-source for Healthcare and Education

Collecting and validating open-source software for healthcare, education, enterprise, development, medical imaging, medical records, and digital pathology.

Visit Medevel

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

Add a comment

Leave a Reply

Keep Up to Date with the Most Important News

By pressing the Subscribe button, you confirm that you have read and are agreeing to our Privacy Policy and Terms of Use

Discover more from Dev solutions

Subscribe now to keep reading and get access to the full archive.

Continue reading