I have a distribution of data that follows approximately a skewed gaussian distribution (count rate vs time). I fit the distribution with the following function in python:

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, time, counrate, p0=[5.5, T0, 4., 0.0001, 0.01], sigma=err_flux)
y_fit = skew(time, popt[0], popt[1], popt[2], popt[3], popt[4])

With this, I can get the parameters sigmag, mu, alpha, c, a and their errors by calculating:

perr = np.sqrt(np.diag(pcov))

However, I want to find the peak of the distribution, that is NOT described by the mu parameter in the function, which actually returns the mean while the peak in the skewed gaussian is given by the mode.

I can calculate the mode by obtaining the time corresponding to the max(y_fit). The question is: how do I get the error on the such calculated mode, given the above-explained procedure to get the maximum of the distribution?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.