How can I plot a function with errors in coefficients?

I would do something along the lines (not sure if I missed any sign or smth).

EDIT: Add generic polynomial

• self.coefs: Coefficients for polynomial in the form of a0 + a1 * x + a2 * x ^ 2 + ...
• self.errs: Errors for each item in self.coefs in the form of e0, e1, e2, ...
• self.coefs_lower: Array of lower limit of coefficients
• self.coefs_upper: Array of upper limit of coefficients
• self.powers: Saved array of powers to raise x with np.power
• self._calc_terms: Function that takes as input x and the coefficients and returns each term to be added in an array: [a0, a1 * x + a2 * x ^ 2, ...]
• self.calc: Calculates the polynomial given x.
• self.calc_lower: Calculates the lowest value of polynomial given the coefficients and errors.
• self.calc_upper: Calculates the highest value of polynomial given the coefficients and errors.

calc_lower and calc_upper can be optimized with some coding & math for terms with odd power to find which one to use in each case (guess I am too lazy to calculate that programmatically)

import matplotlib.pyplot as plt
import numpy as np

class NthOrderPolynomial:
    def __init__(self, coefs, errs):
        assert len(coefs) == len(errs)
        self.coefs = np.array(coefs)[::-1]
        self.errs = np.array(errs)[::-1]
        self.coefs_lower = self.coefs - self.errs
        self.coefs_upper = self.coefs + self.errs
        self.powers = np.array(list(range(len(coefs))))

    @property
    def order(self):
        return len(self.coefs) - 1

    def _calc_terms(self, x, coefs):
        return coefs * np.power.outer(x, self.powers)

    def calc(self, x):
        terms = self._calc_terms(x, self.coefs)
        return np.sum(terms, axis=1)

    def calc_lower(self, x):
        terms_lower = self._calc_terms(x, self.coefs_lower)
        terms_upper = self._calc_terms(x, self.coefs_upper)
        min_terms = np.min([terms_lower, terms_upper], axis=0)
        return np.sum(min_terms, axis=1)

    def calc_upper(self, x):
        terms_lower = self._calc_terms(x, self.coefs_lower)
        terms_upper = self._calc_terms(x, self.coefs_upper)
        max_terms = np.max([terms_lower, terms_upper], axis=0)
        return np.sum(max_terms, axis=1)


coefs = [0.06, -0.65, 1.2]
errs = [0.16, 0.04, 0.01]
quadratic = NthOrderPolynomial(coefs, errs)

x = np.linspace(-10, 10, 21)
y = quadratic.calc(x)
y_lower = quadratic.calc_lower(x)
y_upper = quadratic.calc_upper(x)


fig, (ax1, ax2) = plt.subplots(2)
fig.tight_layout()

ax1.errorbar(x, y, yerr=(y - y_lower, y_upper - y), marker='o', markersize=4, linestyle='dotted')
ax2.fill_between(x, y_lower, y_upper)
plt.show()

For a given x, the expected value will be between

(a + σa)x2 + (b + σb)x + (c + σc),  x >= 0
(a + σa)x2 + (b - σb)x + (c + σc),  x < 0

and

(a - σa)x2 + (b - σb)x + (c - σc),  x >= 0
(a - σa)x2 + (b + σb)x + (c - σc),  x < 0

Another way to look at it is that since σ is always non-negative, the maximum absolute offset from ax2 + bx + c is always going to be

σax2 + σbx + σc,  x >= 0
σax2 - σbx + σc,  x < 0

regardless of the signs of a, b, and c.

This is a quick and dirty solution, but will likely be visually indistinguishable from the real result.

Not sure if there is a tool that can build the band for any function and combination of coefficients. You could figure out the max and min analytically (shown by Mad Physicist's answer) or calculate all possible function outputs if you are unsure (shown below).

Here is a way to determine the minimum and maximum function values using all possible combinations of the coefficients you specified.

import numpy as np
import pandas as pd


# write out the function so that if can take coefficient factors (+/- 1)

def y(x, coef1, coef2, coef3):
    return coef1 * x**2 - coef2 * x + coef3 

# compute each possible combination of coefficients
coef_signs = [1,-1]
# I use a list comprehension but there are several options
possible_coefs = [(0.06+c1*0.016, 0.65+c2*0.04,1.2+c3*0.001) for c1 in coef_signs for c2 in coef_signs for c3 in coef_signs]

Now we have a list of all possible function coefficient combinations.

possible_coefs

    [(0.076, 0.6900000000000001, 1.2009999999999998),
     (0.076, 0.6900000000000001, 1.199),
     (0.076, 0.61, 1.2009999999999998),
     (0.076, 0.61, 1.199),
     (0.044, 0.6900000000000001, 1.2009999999999998),
     (0.044, 0.6900000000000001, 1.199),
     (0.044, 0.61, 1.2009999999999998),
     (0.044, 0.61, 1.199)]

Below I use pandas because I am comfortable with it but it is not necessary, you could use pure numpy arrays just as easily to build the possible function values.

# specify some domain for the function
x = np.linspace(-10,10,100)

# build a Pandas dataframe from the combinations and resulting function values
df = pd.DataFrame(data={f'c0={c[0]};c1={c[1]}; c2={c[2]}':y(x, c[0],c[1],c[2]) for c in possible_coefs})

# calculate the min and max values from all which become the bands 
df['min'] = df.min(1)
df['max'] = df.max(1)

# plot
plt.fill_between(x, df['min'], df['max'])