How to use standard deviation errors from curve fit to plot error bar

Question

I'm wondering if I understand how standard deviation works.

The curve fit function returns the parameters of a function and the their covariances using to find the standard deviation.

Can I use the the standard deviations to plot the fitting curve with the error bars.

Here is an example. I try to keep everything as simple as possible since I'm learning.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

x = np.arange(0,10)
y = np.array([5,4.8,4.7,4.5,4.4,4.5,4.8,4.9,5,5.1])

def function(x,a,b,c):
    return a*x**2+ b*x+c

plt.scatter(x,y)

parameters, pcov = curve_fit(function, x, y)

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

print(perr)

plt.scatter(x, function(x, parameters[0], parameters[1], parameters[2]))

Answer 1

The answer from @Bill is one way to do this. I think there is a simpler way to do this using lmfit (disclosure: lead author).

First, it must be noted that your problem does not necessarily need an iterative curve-fitting approach, as it is a linear problem and can be solved by regression, for example with numpy.polyfit (which still uses least-squares). The approach here is to explore confidence bands and uncertainties with a curve-fitting method that can handle non-linear problems too.

Using lmfit, your script might look like:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

x = np.arange(0,10)
y = np.array([5,4.8,4.7,4.5,4.4,4.5,4.8,4.9,5,5.1])

def function(x, a, b, c):
    return a*x**2+ b*x+c

# first, use numpy polyfit:
a0, b0, c0 = np.polyfit(x, y, 2)

# create lmfit Model from your function:
model = Model(function)

# create Parameters with initial values from polyfit
params = model.make_params(a=a0, b=b0, c=c0)

# run the fit, save the result, print the report of values and statistics
result = model.fit(y, params, x=x)
print(result.fit_report())

# plot data, polyfit result, and least-squares fit
plt.scatter(x, y, label='data')
plt.scatter(x, function(x, a0, b0, c0), marker='+', label='polyfit')
plt.plot(x, result.best_fit, label='least-squares fit')

# use the ModelResult to evaluate the 1-sigma and 2-sigma confidence bands
dy1 = result.eval_uncertainty(sigma=1)
dy2 = result.eval_uncertainty(sigma=2)

# plot confidence bands:
plt.fill_between(x, result.best_fit-dy1, result.best_fit+dy1,
                 color='#AAAAAA60', label='$1\sigma$')
plt.fill_between(x, result.best_fit-dy2, result.best_fit+dy2,
                 color='#AAAAAA30', label='$2\sigma$')
plt.legend()
plt.show()

Note that

1. Start with numpy.polyfit: this will give the solution from regression.
2. Use those values as starting Parameter values for the curve fitting with lmfit.Model.
3. Use the ModelResult.eval_uncertainty() method to calculate uncertainty bands for the best-fit function.

The fit report is:

[[Model]]
    Model(function)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 5
    # data points      = 10
    # variables        = 3
    chi-square         = 0.08307576
    reduced chi-square = 0.01186797
    Akaike info crit   = -41.9058744
    Bayesian info crit = -40.9981191
    R-squared          = 0.84054557
[[Variables]]
    a:  0.02689394 +/- 0.00474101 (17.63%) (init = 0.02689394)
    b: -0.21598485 +/- 0.04432277 (20.52%) (init = -0.2159848)
    c:  4.97545455 +/- 0.08565372 (1.72%) (init = 4.975455)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.9627
    C(b, c) = -0.8099
    C(a, c) = +0.6642

which shows that the polyfit values are not changed, and also shows (1-sigma) uncertainties and correlations between parameters. The plot looks like