如何绘制系数有误差的函数？

Question

这是一个与我之前的帖子类似的问题，但我真的不知道该怎么做（所以我没有包含我的代码）。假设我有一个二次函数

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

+/-符号后的数字代表每个系数的误差。我想知道如何将精确的函数与误差带一起绘制？

我注意到有一个名为plt.fill_between 的工具，这似乎是一个合理的想法，但从它的文档中，我需要知道波段的最大和最小边界，这在我的情况下并不清楚。有没有办法和工具我可以绘制乐队？感谢您的帮助！

Answer 1

对于给定的x，期望值将介于

(a + σ_a)x² + (b + σ_b)x + (c + σ_c ), x >= 0 (a + σ_a)x² + (b - σ_b)x + (c + σ_c ), x

和

(a - σ_a)x² + (b - σ_b)x + (c - σ_c ), x >= 0 (a - σ_a)x² + (b + σ_b)x + (c - σ_c ), x

另一种看待它的方式是，由于σ 始终为非负数，因此与ax<sup>2</sup> + bx + c 的最大绝对偏移量始终为

σ_ax² + σ_bx + σ_c, x >= 0 σ_ax² - σ_bx + σ_c, x

不管a、b和c的符号如何。

这是一个快速而肮脏的解决方案，但在视觉上可能与实际结果无法区分。

Answer 2

不确定是否有工具可以为任何函数和系数组合构建波段。您可以通过分析计算出最大值和最小值（如 Mad Physicist 的回答所示），或者如果您不确定（如下所示）计算所有可能的函数输出。

这是一种使用您指定的系数的所有可能组合来确定最小和最大函数值的方法。

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]

现在我们有了所有可能的函数系数组合的列表。

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

下面我使用 pandas 是因为我觉得它很舒服 但这不是必需的，您可以使用纯 numpy 数组来构建可能的函数值。

# 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'])

Answer 3

我会做一些类似的事情（不确定我是否错过了任何迹象或其他东西）。

编辑：添加通用多项式

• self.coefs：a0 + a1 * x + a2 * x ^ 2 + ... 形式的多项式系数
• self.errs：self.coefs 中每个项目的错误，格式为 e0、e1、e2、...
• self.coefs_lower：系数下限数组
• self.coefs_upper：系数上限数组
• self.powers：使用 np.power 提高 x 的保存能力数组
• self._calc_terms：将 x 和系数作为输入并返回要添加到数组中的每个项的函数：[a0, a1 * x + a2 * x ^ 2, ...]
• self.calc：计算给定 x 的多项式。
• self.calc_lower：根据系数和误差计算多项式的最小值。
• self.calc_upper：根据系数和误差计算多项式的最大值。

calc_lower 和 calc_upper 可以通过一些编码和数学来优化具有奇数能力的术语，以找到在每种情况下使用哪个术语（我想我懒得以编程方式计算）

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