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

2024-06-07 07:17:11 发布

您现在位置:Python中文网/ 问答频道 /正文

这是一个与我之前的帖子类似的问题,但我真的不知道如何做到这一点(所以我没有包括我的代码)。假设我有一个二次函数

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

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

我注意到有一个名为^{}的工具,这似乎是一个合理的想法,但从它的文档中,我需要知道带的最大和最小边界,这在我的例子中并不清楚。有什么方法和工具可以让我画出乐队吗?谢谢你的帮助


Tags: 工具方法函数代码文档绘制符号数字
3条回答

对于给定的x,预期值将介于

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

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

另一种看待它的方式是,因为σ总是非负的,所以ax2 + bx + c的最大绝对偏移量总是

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

不管abc的符号是什么

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

不确定是否有工具可以为任何函数和系数组合构建频带。你可以通过分析计算出最大值和最小值(如Mad物理学家的答案所示),或者如果你不确定,可以计算出所有可能的函数输出(如下所示)

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

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

下面我用熊猫,因为我对它很舒服 但这并不是必须的,您可以使用纯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'])

enter image description here

我会沿着这条路线做一些事情(不确定我是否错过了任何标志或SMT)

编辑:添加通用多项式

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

bars

相关问题 更多 >