x = (x[1:] + x[:-1]) / 2

def gauss(x, mu, sigma, A):
    return A * exp(-(x - mu) ** 2 / 2 / sigma ** 2)

def bimodal(x, mu1, sigma1, A1, mu2, sigma2, A2):
    return gauss(x, mu1, sigma1, A1) + gauss(x, mu2, sigma2, A2)
expected = 1, 0.2, 250, 2, 0.2, 125
params, cov = curve_fit(bimodal, x, y, expected)
sigma = sqrt(diag(cov))
plot(x, bimodal(x, *params), color='red', lw=3, label='model')
legend()
print((params, '\n', sigma))
fit two gaussians to a histogram from one set of data, python