Use R to derive the parameters of a Gamma distribution

Question

I'm not 100% sure if this question should be posted here or in Mathematics Stack Exchange.

I'm a newbie in R.

The probability distribution function of Gamma distribution is given by g(x) = 1/(b^a*Gamma(a))*x^(a-1)e^(-x/b)

I created a Monte Carlo simulation with the code below:

 # Number of simulations
 number_sims <- 10000

 # Samples of fixed size n = 1000
 n = 1000

 for (i in 1:n){
 gamma1 <- rgamma(10, shape=1)
 gamma2 <- rgamma(100, shape=1)
 gamma3 <- rgamma(1000, shape=1)
 gamma4 <- rgamma(5000, shape=1)
 }

 mean_gamma1 <- mean(gamma1)
 mean_gamma1
 mean_gamma2 <- mean(gamma2)
 mean_gamma2
 mean_gamma3 <- mean(gamma3)
 mean_gamma3
 mean_gamma4 <- mean(gamma4)
 mean_gamma4






   

I now intend to use the Method of Moments to derive empirically the values of a and b. I've seen gmm package, but it is not what I intended. I intended to do a step-by-step approach as if I were doing by hand. Do you know any package, or any method to do that? Or I would need to create my own function?

Answer 1

I think the "fitdistplus" package is what you are looking for:

library(fitdistrplus)
fitdist(simulation1, "gamma", method = "mme")



Fitting of the distribution ' gamma ' by matching moments 
Parameters:
      estimate
shape 9.223529
rate  9.217187