x = arange(-10, 10, 0.01)

def f(p):
    mu, s = p
return exp(-(x - mu) ** 2 / (2 * s ** 2))

def fff(d):

    def ff(p):
        return d - f(p)
return ff

def noisy_data(p):
    return f(p) + normal(0, 0.1, len(x))

def fit(d):
    ff = fff(d)
p = leastsq(ff, [0, 1])[0]
return p

def bootstrap(d):
    p0 = fit(d)
residuals = f(p0) - d
s_residuals = std(residuals)
ps = []
for i in range(1000):
    new_d = d + normal(0, s_residuals, len(d))
ps.append(fit(new_d))
ps = array(ps)
mean_params = mean(ps, 0)
std_params = std(ps, 0)
return mean_params, std_params
data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)
print('95% confidence interval:')
print('mu: ', mean_params[0], ' +/- ', std_params[0] * 1.95996)
print('sigma: ', mean_params[1], ' +/- ', std_params[1] * 1.95996)
confidence interval with leastsq fit in scipy python