使用numpy拟合数据

59

我有以下数据:

>>> x
array([ 3.08,  3.1 ,  3.12,  3.14,  3.16,  3.18,  3.2 ,  3.22,  3.24,
    3.26,  3.28,  3.3 ,  3.32,  3.34,  3.36,  3.38,  3.4 ,  3.42,
    3.44,  3.46,  3.48,  3.5 ,  3.52,  3.54,  3.56,  3.58,  3.6 ,
    3.62,  3.64,  3.66,  3.68])

>>> y
array([ 0.000857,  0.001182,  0.001619,  0.002113,  0.002702,  0.003351,
    0.004062,  0.004754,  0.00546 ,  0.006183,  0.006816,  0.007362,
    0.007844,  0.008207,  0.008474,  0.008541,  0.008539,  0.008445,
    0.008251,  0.007974,  0.007608,  0.007193,  0.006752,  0.006269,
    0.005799,  0.005302,  0.004822,  0.004339,  0.00391 ,  0.003481,
    0.003095])

现在,我想用一个4次多项式来拟合这些数据。所以我做了这个:

>>> coefs = np.polynomial.polynomial.polyfit(x, y, 4)
>>> ffit = np.poly1d(coefs)

现在我创建了一个新的网格以评估拟合函数 ffit 的 x 值:

>>> x_new = np.linspace(x[0], x[-1], num=len(x)*10)

当我使用以下命令绘制所有的图表(数据集和拟合曲线)时:

>>> fig1 = plt.figure()                                                                                           
>>> ax1 = fig1.add_subplot(111)                                                                                   
>>> ax1.scatter(x, y, facecolors='None')                                                                     
>>> ax1.plot(x_new, ffit(x_new))                                                                     
>>> plt.show()

我得到了以下结果:

fitting_data.png fitting_data.png

我期望拟合函数能够正确拟合数据(至少在数据的最大值附近)。我做错了什么?

3个回答

115

不幸的是,np.polynomial.polynomial.polyfit返回的系数与np.polyfitnp.polyval(或者您使用的np.poly1d)相反。为了说明:

In [40]: np.polynomial.polynomial.polyfit(x, y, 4)
Out[40]: 
array([  84.29340848, -100.53595376,   44.83281408,   -8.85931101,
          0.65459882])

In [41]: np.polyfit(x, y, 4)
Out[41]: 
array([   0.65459882,   -8.859311  ,   44.83281407, -100.53595375,
         84.29340846])

通常情况下:np.polynomial.polynomial.polyfit返回系数[A,B,C],表示A + Bx + Cx^2 + ...,而np.polyfit则返回:... + Ax^2 + Bx + C

因此,如果要使用这组函数的组合,则必须颠倒系数的顺序,如下所示:

ffit = np.polyval(coefs[::-1], x_new)

然而,文档明确说明要避免使用np.polyfitnp.polyvalnp.poly1d,而应仅使用新的多项式包。

最安全的做法是仅使用多项式包:

import numpy.polynomial.polynomial as poly

coefs = poly.polyfit(x, y, 4)
ffit = poly.polyval(x_new, coefs)
plt.plot(x_new, ffit)

或者,要创建多项式函数:

ffit = poly.Polynomial(coefs)    # instead of np.poly1d
plt.plot(x_new, ffit(x_new))

适合度和数据图


37
请注意,您可以直接使用 Polynomial 类进行拟合并返回一个 Polynomial 实例。
from numpy.polynomial import Polynomial

p = Polynomial.fit(x, y, 4)
plt.plot(*p.linspace())

p 使用缩放和移位的 x 值以保证数值稳定性。如果需要使用常规形式的系数,您需要进行以下操作:

pnormal = p.convert(domain=(-1, 1))

6
对于系数的转换可以使用+1,如果需要在默认域上与其他多项式进行计算,则非常有用。请注意,这可以通过fit()方法直接完成,使用相同的“domain”参数。 - BenC
通过在 domain 参数中使用值 [](自版本1.5起),也可以避免使用 convert() - Mark H
要获取未缩放和未移位的系数,也可以执行p.convert().coef - undefined

2

使用numpy和matplotlib拟合Chebyshev级数和多项式级数最小二乘拟合曲线

快速摘要

下面是你需要注意的关键代码行,用于执行数据拟合,例如:

import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial

poly_series = Polynomial.fit(x, y, deg=5)  # 5th order polynomial curve fit
x_poly, y_poly = poly_series.linspace()
plt.plot(x_poly, y_poly)

以下是我完整程序的输出:

enter image description here

详情

即使阅读了这里所有的答案,由于Matplotlib的变化和我遇到的其他问题,我一直在苦苦挣扎,直到最终找到了这个演示文稿!:https://numpy.org/doc/stable/reference/routines.polynomials.classes.html

它包含一个非常好的示例,在“拟合”部分的底部!:

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Chebyshev as T

np.random.seed(11)
x = np.linspace(0, 2*np.pi, 20)
y = np.sin(x) + np.random.normal(scale=.1, size=x.shape)
p = T.fit(x, y, 5)
plt.plot(x, y, 'o')
xx, yy = p.linspace()
plt.plot(xx, yy, lw=2)
p.domain
p.window
plt.show()

这是一个增强版,包含我的修改和一些有用的注释。我从我的plot_best_fit_polynomial.py演示中修改了下面的代码,该演示位于我的eRCaGuy_hello_world存储库中:

import matplotlib.pyplot as plt
from numpy.polynomial import Chebyshev
from numpy.polynomial import Polynomial

# data to fit
x =   [0.        , 0.33069396, 0.66138793, 0.99208189, 1.32277585,
       1.65346982, 1.98416378, 2.31485774, 2.64555171, 2.97624567,
       3.30693964, 3.6376336 , 3.96832756, 4.29902153, 4.62971549,
       4.96040945, 5.29110342, 5.62179738, 5.95249134, 6.28318531]
y =   [ 0.17494547,  0.29609217,  0.5657562 ,  0.57183462,  0.9685718 ,
        0.96462136,  0.86211039,  0.76726418,  0.51805246,  0.05803429,
       -0.25321856, -0.52352074, -0.66675568, -0.85965411, -1.12713934,
       -1.08134779, -0.76348274, -0.45674931, -0.32780698, -0.06834466]

# Obtain a 5th degree (order) least-squares fit curve to the x, y data using a
# Chebyshev Series
cheby_series = Chebyshev.fit(x, y, deg=5)
# Lines-space: get evenly-spaced points to plot a line; see:
# https://numpy.org/doc/stable/reference/generated/numpy.polynomial.chebyshev.Chebyshev.linspace.html
x_cheby, y_cheby = cheby_series.linspace()

# Now do the same things with a 5th order Polynomial Series fit as well!
# see: https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.html
poly_series = Polynomial.fit(x, y, deg=5)
x_poly, y_poly = poly_series.linspace()

# plot all the data

f1 = plt.figure()

plt.plot(x, y, 'o')
plt.plot(x_cheby, y_cheby, linewidth=6, alpha=0.5,
   label="Chebyshev Series 5th degree\nleast squares best fit curve")
plt.plot(x_poly, y_poly, 'k', linewidth=1,
   label="Polynomial Series 5th degree\nleast squares best fit curve")

plt.legend()
plt.show()

相关:

  1. 我刚刚也学会了如何给图表加标签:在Matplotlib中如何添加图标题、图副标题、图脚注、绘图标题、坐标轴标签、图例标签和(x,y)点标签

参考资料:

  1. https://numpy.org/doc/stable/reference/routines.polynomials.classes.html - 在底部查看演示
  2. https://numpy.org/doc/stable/reference/routines.polynomials.html#documentation-for-the-polynomial-package
  3. https://numpy.org/doc/stable/reference/routines.polynomials.chebyshev.html
  4. https://numpy.org/doc/stable/reference/generated/numpy.polynomial.chebyshev.Chebyshev.linspace.html
  5. https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.html

你的第一个示例图并不代表代码:代码执行的是一个五次多项式拟合,而图表显示它是一个切比雪夫拟合,这是不正确的。 - undefined
1
@user115625,我能理解你的困惑。为了更清楚地表达,我已经更新了我的回答。现在我说:“这是我完整程序的输出结果下方(_不是_仅仅上面的部分程序片段)”。快速摘要部分的目标之一是为了让你知道在详细信息部分的完整代码中应该查找的位置。 - undefined
而在“详细信息”部分的完整代码是生成该图的原因。 - undefined

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接