# Find the error on the peak value (mode) of a skewed gaussian

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, popt, popt, popt, popt)


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?