Finding shape and scale parameters of gamma distribution

Question

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

jelsner tutorial [1]: https://drive.google.com/file/d/1JcVXCaxp74w3OJJfsdVcVvKBE0OeB4zo/view?usp=sharing [2]: https://i.stack.imgur.com/oapPG.png

Answer 1

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.

Answer 2

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.