CODEWITHSUNDEEP
pythoncurve-fittinggaussian
freja
freja

Reputation: 59

Asymmetric Gaussian Fit in Python

I'm trying to fit an asymmetric Gaussian to this data: http://ge.tt/99iNaL53 (csv file).

I have tried to use a skewed Gaussian model from lmfit, and also a spline, but I'm not able to get the Gaussian model to fit well and the splines are not what I'm looking for (I don't want the spline to fit the data exactly as shown below, and altering the level of smoothing isn't helping).

Here is code using the above data that produces the plot below. The second figure is an example of what I'm trying to achieve with the goal of reading the rise and decay time from the fit.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from scipy.interpolate import UnivariateSpline
from lmfit.models import SkewedGaussianModel

data = np.loadtxt('data.csv', delimiter=',')

x = data[:,0]
y = data[:,1]

# Skewed Gaussian fit
model = SkewedGaussianModel()
params = model.make_params(amplitude=400, center=3, sigma=7, gamma=1)
result = model.fit(y, params, x=x)

# Cubic Spline
cs = CubicSpline(x, y)
x_range = np.arange(x[0], x[-1], 0.1)

# Univariate Spline
us = UnivariateSpline(x, y, k = 1)

# Univariate Spline (smoothed)
us2 = UnivariateSpline(x, y, k = 5)

plt.scatter(x, y, marker = '^', color = 'k', linewidth = 0.5, s = 10, label = 'data')
plt.plot(x_range, cs(x_range), label = 'Cubic Spline')
plt.plot(x_range, us(x_range), label = 'Univariate Spline, k = 1')
plt.plot(x_range, us2(x_range), label = 'Univariate Spline, k = 5')
plt.plot(x, result.best_fit, color = 'red', label = 'Skewed Gaussian Attempt')
plt.xlabel('x')
plt.ylabel('y')
plt.yscale('log')
plt.ylim(1,500)
plt.legend()
plt.show()

Data and fit attempts Example

Upvotes: 3

Views: 8728

Answers (2)

Michael Green
Michael Green

Reputation: 810

I wrote something for J. Chem. Ed. [1] that involved fitting asymmetric Gaussian functions to data, you can find the core repo here [2] but below is a snippet on how I went about fitting a data set where x = data[:,0] and y = data[:,1] to the type of function you're working with:

import numpy as np
from scipy.optimize import leastsq
from scipy.special import erf

initials = [6.5, 13, 1, 0] # initial guess

def asymGaussian(x, p):
    amp = (p[0] / (p[2] * np.sqrt(2 * np.pi)))
    spread = np.exp((-(x - p[1]) ** 2.0) / (2 * p[2] ** 2.0))
    skew = (1 + erf((p[3] * (x - p[1])) / (p[2] * np.sqrt(2))))
    return amp * spread * skew

def residuals(p,y,x):
    return y - asymGaussian(x, p)

# executes least-squares regression analysis to optimize initial parameters
cnsts = leastsq(
    residuals, 
    initials, 
    args=(
        data_set[:,1], # y value
        data_set[:,0]  # x value
        ))[0]

y = asymGaussian(data[:,0], cnsts)

finally just plot (y, data[:,0]). Hope this helps!

[1] https://pubs.acs.org/doi/10.1021/acs.jchemed.9b00818
[2] https://github.com/1mikegrn/pyGC

Upvotes: 5

M Newville
M Newville

Reputation: 7862

Is there a question here? I don't see one, actually.

That result from lmfit is the best fit to a skewed Gaussian model. You've chosen to plot the result on a log-scale. That completely changes the view of the quality of the fit or what is not fit well.

It seems like you're expecting a better fit, but not *too good. Well, it looks like your data is not perfectly represented by a single skewed Gaussian. It seems like you were not expecting it to be. You could try different forms for the model function, say a skewed Lorentzian or something. But your data has that low x shoulder that definitely does not look like your uncited figure.

Upvotes: 3

Related Questions