How can I plot a function with errors in coefficients?

Question

This is a question similar to my earlier post, but I really have no idea how to do that (so I didn't include my code). Suppose I have a quadratic function

y = 0.06(+/-0.16)x**2 - 0.65(+/-0.04)x + 1.2(+/-0.001)

The numbers after +/- sign represent the errors of each coefficient. I wonder how can I plot the exact function together with an error band?

I noticed there's a tool called plt.fill_between which seems like a plausible idea, but from its documentation, I need to know the maximum and minimum boundaries of the band, which is not clear in my case. Is there a way and tool I can plot the band? Thanks for the help!

Answer 1

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()

Answer 2

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.