Finding shape and scale parameters of gamma distribution

I have a dataset [click here][1] and trying to find the shape and scale parameters of gamma distribution.

> library(fitdistrplus)
> H <- "US.txt"
> rate <- mean(H$$All) > n = length(H$$Year)
> descdist(H$All, discrete = FALSE) summary statistics ------ min: 0 max: 7 median: 1 mean: 1.69375 estimated sd: 1.449341 estimated skewness: 1.126422 [![possible distribution][2]][2] #ratio of the variance to the mean > round(var(H$All)/rate, 2)
[1] 1.24


#Is this unusual for a Poisson distribution?

#check by performing a Monte Carlo (MC) simulation #repeat this ( m=1000 ) times. Let ( n ) be the number of years in your data record and ( \lambda ) be the rate. For each sample you compute the ratio of the variance to the mean.

> set.seed(3042)
> ratio = numeric()
> m = 1000
for (i in 1:m) {
h = rpois(n = n, lambda = rate)
ratio[i] = var(h)/mean(h)
}


The vector ratio contains 1000 values of the ratio. So determine the proportion of ratios greater than 1.24

> sum(ratio > var(H$All)/rate)/m 0.028  Modify your MC simulation using the gamma distribution for the rate and then examine the ratio of variance to the mean from a set of Poisson counts with the variable rate. The gamma distribution describes the variability in the rate using the shape and scale parameters. > ratio = numeric() > set.seed(3042) > m = 1000 for (i in 1:m) { h = rpois(n = n, lambda = rgamma(m, shape = ???, scale = ???)) ratio[i] = var(h)/mean(h) } sum(ratio > var(H$All)/rate)/m


So I want to specify the shape to be ??? and the scale to be ??? so the product matches closely the long-term average count of the dataset. To find shape and scale I used below codes but received error.

> fitdistr(H$All, "gamma") Error in stats::optim(x = c(1L, 3L, 0L, 2L, 1L, 2L, 1L, 1L, 1L, 3L, 2L, : initial value in 'vmmin' is not finite  Tried weibull > fitdistr(H$All, "weibull")
Error in fitdistr(H$All, "weibull") : Weibull values must be > 0  Why am I not getting shape and scale factor here? Additionally, if I do it manually m= 1.449341 sd= 1.69375 > shape = (m/sd)^2 > scale = (sd)^2/m > shape [1] 0.7322216 > scale [1] 1.979375  • You have integers$0,1,2,3,4,5,6,7$so you might consider trying to fit a distribution supported on integers such as binomial, Poisson, negative binomial Jul 10, 2022 at 9:11 • @Henry please have a look at the updated question. Jul 10, 2022 at 10:37 2 Answers I'm going to break ranks with others complaining that you can't fit a continuous distribution to a discrete valued sample. We actually do this all the time, and it's an interesting problem in statistical computing and asymptotics to consider what happens when the distributional assumptions aren't exactly met. So... you know, think critically about your problems. We can give you the benefit of a doubt and treat this as an interesting, high level problem. Otherwise, it should be fairly obvious, if you want data to be generated from a given distribution with given moments, you must simulate data from that distribution with the given moments. Poisson distributions are a 1 parameter, skewed distribution on the positive integers, Gamma distributions are 2 parameter skewed distributions on the positive reals. For a Poisson $$\lambda$$ RV, the mean is $$\lambda$$ and the variance is $$\lambda$$, i.e. the variance is equal to the mean. For a Gamma ($$\alpha, \beta$$) (shape, scale) distribution, the mean is $$\alpha/\beta$$ and the variance is $$\alpha/\beta^2$$. In other words, if you want to have variance=mean, you need to constrain $$\beta=1$$. However, you can quickly see there is no 1-1 correspondence with the higher order moments. A Gamma distribution will have skewness $$2/\sqrt{\alpha}$$ whereas a Poisson distribution will have $$1/\sqrt{\lambda}$$. In other words, a Gamma having the same mean, variance as Poisson is always twice as skewed. The properties of MLEs are certainly "interesting" if not favorable in this case. A misspecified MLE will give a distributional estimate that doesn't even have an unbiased estimate of the mean. I can't go through the same exercise with a Weibull, which (although being a 2 parameter, skewed, positive valued distribution) is completely different from Gamma. The constraint between mean and variance is far more complicated, and skewness even more so. It shouldn't be surprising to see that the same issues can give a problem. First, I don't understand why you are aiming to model an integer (are those counts of an event?) with continuous distributions such as Gamma and Weibull. I think Poisson and Negative Binomial would be more natural choices. I tested your code and I was able to reproduce the problem. If I interpret correctly, you can't have a zero as a realization of a Gamma distribution (see here, "$$E[X] = kθ = α/β$$ is fixed and greater than zero", zero is not in the support). Concerning the Weibull attempt, I believe this answer could be helpful, the culprit is similar and again due to exact zeros. I hope this is helpful, I couldn't comment to ask for clarification as my reputation is still too low. • "are those counts of an event"-yes it is 'frequency of an event'- I did Poisson but the result is (variance is not equal to mean) 1.24, indicates that the variance of event counts is 24% larger than the mean. Isn't this unusual for a Poisson distribution? Jul 10, 2022 at 13:40 • Ok, now I see. Your MC simulation was an attempt to measure the frequency of this phenomenon, which seems to be fairly rare. When the count data has more variance than expected then we have "over-dispersion" wiki. In that case, replacing Poisson with Negative-binomial seems to be a viable strategy. Jul 10, 2022 at 13:56 • (0.28) 2.8% of the ratios are larger, so the answer from MC experiment is 'yes,' the variability in event counts is higher than expected (unusual) from a Poisson distribution with a constant rate. Therefore, MC simulation was modified using the gamma distribution for the rate and then examine the ratio of variance to the mean from a set of Poisson counts with the variable rate. Jul 10, 2022 at 13:58 • we use here scale and shape to find the rate i.e lamda" suppose, h = rpois(n = n, lambda = rgamma(m, shape = 5.6, scale = 0.3)) In this case, shape of 5.6 and scale of 0.3 gives a close value to the mean 1.69375. In order to find the value of shape and scale, I used fitdistr(H$All, "gamma") or fitdistr(H\$All, "weibull"). I know that value of shape is 5.6 and scale is 0.3, but how do we get it ? Jul 10, 2022 at 14:01
• I still think that the problem is that you can't fit a Gamma or a Weibull to your data as they are (counts). But you can fit an over-dispersed model and trace it back to the likely parameters of your Gamma distribution. From the Gamma–Poisson mixture section here you can see that shape = r and scale θ = p/(1 − p) . If you can get r and p fitting your data, you will have the shape and the scale of your Gamma. Type ?dnbinom in R and you will get additional information that will be helpful. Jul 10, 2022 at 14:36