I have a data set that I want to fit to a gamma distribution. I'm trying out the method that I found here: http://scipy-lectures.org/intro/scipy/auto_examples/plot_curve_fit.html but it does not seem to work with my data. I've tried other methods too (Scipy's stats.gamma.fit, Fitter from the package with the same name).
My data can be found here: https://www.filehosting.org/file/details/832322/data.csv. It contains both the x- and y-axis.
This is my current code:
import csv
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import integrate as scint
from scipy import optimize
def gengamma_fit(x, a, b):
return 2 / (b**(2*a) * math.gamma(a)) * x**(2 * a - 1) * math.exp(-x / b)**2
gen_gamma_fit = np.vectorize(gengamma_fit)
filename = r'data.csv'
x1 = np.zeros(200)
y1 = np.zeros(200)
output1 = np.zeros((200,2))
with open(filename, newline='\n') as file:
data = csv.reader(file, delimiter=',')
i = 0
for row in data:
x1[i] = row[0]
y1[i] = row[1]
output1[i][0] = row[0]
output1[i][1] = float(row[1])
i += 1
# check to see if data is normalized
int1 = scint.simps(y1, x1)
print(int1)
plt.figure(figsize=(15,10))
plt.plot(x1,y1, color="red")
params, params_covariance = optimize.curve_fit(gen_gamma_fit, x1, y1, p0=[2.0, 0.04])
plt.plot(x1, gen_gamma_fit(x1, params[0], params[1]))
plt.show()
plt.close()
When I tried to do the fit using Mathematica, it gives me a very close fit. When I do it in Python, I get this:
It gets weirder (for me): when I use the parameters that Mathematica gives, Python does not plot it the same. I plotted the parameters using this line of code:
plt.plot(x1, gen_gamma_fit(x1, 2.57465, 0.10187))
As you can see, a = 2.57465 and b = 0.10187.
Am I making some obvious mistake here?
