CODEWITHSUNDEEP
matplotlibscipycurve-fittinggaussian
new-to-coding
new-to-coding

Reputation: 51

How to fit a Gaussian best fit for the data

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 -

enter image description here

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 -

enter image description here

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

Upvotes: 1

Views: 856

Answers (2)

JJacquelin
JJacquelin

Reputation: 1705

By inspection one obweve that the curve isn't symmetrical. This is even more visible in logarithmic y-scale zoomed in the range close to the peak.

enter image description here

This draw to think that the Gaussian model (which is symetrical) cannot be fitted correctly.

Also one observe that the part of the curve around the peak isn't far from to be linear. Thus a better model might be made of the combination of two exponential functions, for example :

enter image description here

enter image description here

I suppose that you can code this function in your nonlinear regression sofware. The above rough values of the parameters can be used as starting values for the iterative calculus.

Upvotes: 1

blunova
blunova

Reputation: 2532

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()

enter image description here

Upvotes: 0

Related Questions