Define skewed gaussian function that returns two parameters after fitting

Question

I know such questions have been asked a few times already, but still I can't get it: I want to define a function that returns more than one parameter.

I have written the following code to fit data to a skewed gaussian function:

def skew(x, sigmag, mu, alpha, c, a):
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2))))
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2)))))
    return 2*a*normpdf*normcdf + c

popt, pcov = curve_fit(skew, xdata, ydata, p0=(5.5, 57636., 4.5, 0.0001, 0.01))
y_fit= skew(xdata, popt[0], popt[1], popt[2], popt[3], popt[4])

However, my idea is to get the peak of the data distribution, and not the mean returned by the skew function as one of the best fit values. Therefore, I would need the mode of the distribution, which can be calculated as the maximum of the normpdf.

How do I get normpdf from my defined function and get its maximum over fitted data?

Answer 1

Your code is not a Minimal, Complete, and Verifiable example we can run, and you supplied no sample output, but I think I see the question. You appear to be asking about the comma ',' tuple unpacking operator used in the popt, pcov = curve_fit(...) line. We will leave that line unchanged, and recover mode from your function. Use this:

def skew2(x, sigmag, mu, alpha, c, a):
    normpdf = (1 / (sigmag * np.sqrt(2 * math.pi))) * np.exp(-(np.power((x - mu), 2) / (2 * np.power(sigmag, 2))))
    normcdf = (0.5 * (1 + sp.erf((alpha * ((x - mu) / sigmag)) / (np.sqrt(2)))))
    return 2 * a * normpdf * normcdf + c, max(normpdf)

def skew(x, sigmag, mu, alpha, c, a):
    return skew2(x, sigmag, mu, alpha, c, a)[0]

popt, pcov = curve_fit(skew, xdata, ydata, p0=(5.5, 57636., 4.5, 0.0001, 0.01))
y_fit, mode = skew2(xdata, *popt[:5])