title: Why does Scipy's fit() function for gamma distribution producing a totally dissimilar distribution?
body: <p>I'm trying to use <code>scipy</code>'s fit function to fit a gamma distribution to observed data. The histogram for the observed data is below:</p>
  
  <p><a href="https://i.stack.imgur.com/NNyTD.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NNyTD.png" alt="enter image description here"></a></p>
  
  <p>The mean and variance are</p>
  
  <pre><code>print("mean",np.mean(observed_data)) # mean 0.427611176580073
  print("Var",np.var(observed_data)) # Var 0.6898193689790143
  </code></pre>
  
  <p>However, if I use <code>scipy.stats.gamma.fit()</code> to fit this observed data to a gamma distribution and then sample again from that distribution, the mean and variance are totally different:</p>
  
  <p><a href="https://i.stack.imgur.com/DWW9b.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DWW9b.png" alt="enter image description here"></a></p>
  
  <p>I know that <code>scipy</code> is fitting based on MLE, but I'm not understanding the intuition behind why these key statistics are so off - the mean and variance are totally different. In fact, I can get much better results simply running this through my own solver:</p>
  
  <pre><code>from scipy.optimize import fsolve
  from typing import List
  def fit_gamma_distribution(data: List[float]):
      mean = np.mean(data)
      variance = np.var(data)
  
      def equations(p):
          k, theta = p
          return (k * theta - mean, k * theta **2 - variance)
      solved_k, solved_theta = fsolve(equations, (1,1))
  
      if np.isclose(np.array([solved_k * solved_theta]), np.array([mean]), rtol=0.01):
          return fsolve(equations, (1,1))
  
  k, theta = fit_gamma_distribution(observed_data)
  
  new_dist = np.random.gamma(shape=k, scale=
  theta, size=len(observed_data))
  plt.hist(new_dist, alpha= 0.5, bins=40)
  plt.hist(observed_data, alpha=0.2, bins=40)
  plt.xlim(0,5)
  plt.title(f"New sampled distribution: μ = {round(np.mean(new_dist),2)} observed μ = {round(np.mean(observed_data), 2)}")
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/8RZHa.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8RZHa.png" alt="enter image description here"></a></p>
  
  <p>It's a much better fit than the <code>scipy</code> one in my opinion. Why is this the case?</p>