How to fit a Gaussian best fit for the data

Question

I have a data set for which I am plotting a graph of time vs intensity at a particular frequency.

On the x axis is the time data set which is in a numpy array and on the y axis is the intensity array.

time = [ 0.3  1.3  2.3  3.3  4.3  5.3  6.3  7.3  8.3  9.3 10.3 11.3 12.3 13.3
 14.3 15.3 16.3 17.3 18.3 19.3 20.3 21.3 22.3 23.3 24.3 25.3 26.3 27.3
 28.3 29.3 30.3 31.3 32.3 33.3 34.3 35.3 36.3 37.3 38.3 39.3 40.3 41.3
 42.3 43.3 44.3 45.3 46.3 47.3 48.3 49.3 50.3 51.3 52.3 53.3 54.3 55.3
 56.3 57.3 58.3 59.3] 

intensity = [1.03587, 1.03187, 1.03561, 1.02893, 1.04659, 1.03633, 1.0481 ,
       1.04156, 1.02164, 1.02741, 1.02675, 1.03651, 1.03713, 1.0252 ,
       1.02853, 1.0378 , 1.04374, 1.01427, 1.0387 , 1.03389, 1.03148,
       1.04334, 1.042  , 1.04154, 1.0161 , 1.0469 , 1.03152, 1.22406,
       5.4362 , 7.92132, 6.50259, 4.7227 , 3.32571, 2.46484, 1.74615,
       1.51446, 1.2711 , 1.15098, 1.09623, 1.0697 , 1.06085, 1.05837,
       1.04151, 1.0358 , 1.03574, 1.05095, 1.03382, 1.04629, 1.03636,
       1.03219, 1.03555, 1.02886, 1.04652, 1.02617, 1.04363, 1.03591,
       1.04199, 1.03726, 1.03246, 1.0408 ]

When I plot this using matplotlib, using;

plt.figure(figsize=(15,6))
plt.title('Single frequency graph at 636 kHz', fontsize=18)
plt.plot(time,intensity)
plt.xticks(time[::3], fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Elapsed time (minutes:seconds)', fontsize=18)
plt.ylabel('Intensity at 1020 kHz', fontsize=18)
plt.savefig('WIND_Single_frequency_graph_1020_kHz')
plt.show()

The graph looks like -

I want to fit a gaussian curve for this data, and this is the code I used,

def Gauss(x, A, B):
    y = A*np.exp(-1*B*x**2)
    return y
parameters, covariance = curve_fit(Gauss, time, intensity_636)


fit_A = parameters[0]
fit_B = parameters[1]

fit_y = Gauss(time, fit_A, fit_B)

plt.figure(figsize=(15,6))
plt.plot(time, intensity, 'o', label = 'data')
plt.plot(time, fit_y, '-', label ='fit')
plt.legend()

And the best fit i obtain looks like this -

Where am I going wrong? How can I make the best fit curve fit the data better?

Answer 1

You need to define a more flexible model (more parameters) and to define reasonable initial values for them:

def f(x, a, b, mu, sigma):
    return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))

popt, pcov = curve_fit(f, time, intensity, p0=[1, 1, 30.3, 2])

plt.plot(time, intensity)
plt.plot(time, f(time, *popt))
plt.show()