I got a study of 210 samples and I tried fitting gamma distribution to them. I used method of moments and maximum likelihood estimation to calculate the parameters, but parameters came out quite different (that much that the shape of curve is different).
This is the histogram with green graph of gamma PDF, calculated from MLE and red graph calculated from MoM. The data is this.
The derivation of MoM and MLE are certainly correct (I tried them both on samples of 1000 generated by gamma distribution and they are very precise), but if you want to check it out:
Method of momements
For gamma distribution we have $\text{E}(\bar{X}) = \frac{a}{\lambda}$ and $\text{E}(\bar{X}^2) = \frac{a(a+1)}{\lambda^2}$. On the other hand we have formulas $\text{E}(\bar{X}) = \frac{1}{n} \sum_{i=1}^n x_i = \bar{x}$ and $\text{E}(\bar{X}^2) = \frac{1}{n} \sum_{i=1}^n x_i^2$. Combining these we get \begin{equation} \hat{a} = \hat{\lambda} \bar{x}~~~~\text{and}~~~~\hat{\lambda} = \frac{\bar{x}}{\frac{1}{n}\sum_{i=1}^n x_i^2 - \bar{x}^2}. \end{equation} From here I got $$a \approx 0.799174075139 ~~~~\text{and}~~~~ \lambda \approx 1.31877117345. $$
Maximum likelihood method
For gamma distribution we have
$$f_{X_i}(x_i) = \frac{\lambda^a}{\Gamma(a)} x_i^{a-1} e^{-\lambda x_i}.$$
Taking the product from $1$ to $n$ we get the likelihood function
$$L = \prod_{i=1}^n \frac{\lambda^a}{\Gamma(a)} x_i^{a-1} e^{-\lambda x_i}$$
and log-likelihood function
$$\ell = \log L = \sum_{i=1}^n \left(a \log\lambda - \log\Gamma(a) + (a-1)\log x_i -\lambda x_i\right) =\\
=na \log\lambda - n \log\Gamma(a) + (a-1)\sum_{i=1}^n \log(x_i) - \lambda \sum_{i=1}^n x_i.$$
Calculating partial derivatives with respect to $a$ and $\lambda$ and equating with $0$ we obtain
\begin{align}
\frac{\partial \ell}{\partial a} &= n \log \lambda - n \psi(a) + \sum_{i=1}^n \log(x_i) = 0,\\
\frac{\partial \ell}{\partial \lambda} &= \frac{na}{\lambda} - n \bar{x} = 0,
\end{align}
where $\psi(a)$ is a digamma function ($\psi(a) = \frac{\Gamma'(a)}{\Gamma(a)}$). Now I solved this system by plugging the second equation into first, threw everything on one side and using scipy.optimize.newton I found the root $a$. Plugging $a$ into second equation I got
$$a \approx 1.59540858049 ~~~~\text{in}~~~~ \lambda \approx 2.63269156405.$$
So what I would like to know is what are the reasons for this (e.g. small sample size, bad sample, etc.).
