使用numpy的多项式外推

3hvapo4f  于 2023-06-23  发布在  其他
关注(0)|答案(2)|浏览(149)

我有一个信号，我希望首先使用多项式拟合来拟合，然后使用该多项式表达式进行外推。这意味着，最初，我的x轴范围是192.1 - 196THz。我有这个范围对应的y轴值。但我想绘制191.325 - 196.125 THz的信号，这意味着我需要预测原始x轴范围之外的y轴值。这是我的代码

import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly



y=[9.996651009,10.13327408,10.27003353,10.43754263,10.60523071,10.71136134,10.81753489,10.92299416,11.02853271,11.10093835,11.17349619,11.16565259,11.15787539,11.16191795,11.16626826,11.16087842,11.15561384,11.16940196,11.18346164,11.22434472,11.26537218,11.35981126,11.45427174,11.6102862,11.76633923,11.89058563,12.01486743,12.21133227,12.40782023,12.55682327,12.7059072,12.88910065,13.07234341,13.20139753,13.33048759,13.46960371,13.60875466,13.71956075,13.83042354,13.91157425,13.99277241,14.08548465,14.17824161,14.26091536,14.34364144,14.43504833,14.52651725,14.59548578,14.66449439,14.73751377,14.81058218,14.88687742,14.96322511,15.038278,15.11342222,15.1644832,15.21556914,15.31798775,15.42044285,15.51000031,15.59959896,15.68089263,15.76229354,15.86213633,15.96214484,16.0269976,16.0923806,16.16483146,16.23762336,16.24233076,16.24719576,16.26116812,16.27543692,16.28311024,16.29128992,16.23458771,16.17830522,16.12636211,16.07478778]

y.reverse()

# x axis for polyfit
i = 192.1e12
x=[]
while i <= 196e12:
  x.append(i)
  i += 50e9

plt.plot(x, y)

# x axis for extrapolation
i1= 191.325e12
x2=[]
while i1 <= 196.125e12:

  x2.append(i1)
  i1 += 50e9

z = poly.polyfit(x, y, 20)
f = poly.polyval(x, z)

plt.plot(x2, f(x2), linestyle='--')

问题：
1.我得到的错误：“numpy. ndarray”对象不可调用。（第37章）我找不到工作。
2.如果我把超过21作为系数的数量，那么我运行到下面的问题：SVD在线性最小二乘法中不收敛。（第34行）我基本上对这个问题一无所知。
有人能给我一些想法吗？谢谢你。

numpy

2条答案

按热度按时间
kpbpu008

kpbpu0081#

外推法是不准确的，但是你可以使用scipy.interpolate.interp1d来完成，告诉它不要给予边界错误（bounds_error=False），并通过外推法填充值（fill_value="extrapolation"）。在我的代码中，我使用了三次样条拟合。此外，使用numpy函数生成x数组。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

plt.close("all")

y = [9.996651009, 10.13327408, 10.27003353, 10.43754263, 10.60523071, 
     10.71136134, 10.81753489, 10.92299416, 11.02853271, 11.10093835, 
     11.17349619, 11.16565259, 11.15787539, 11.16191795, 11.16626826, 
     11.16087842, 11.15561384, 11.16940196, 11.18346164, 11.22434472, 
     11.26537218, 11.35981126, 11.45427174, 11.6102862, 11.76633923, 
     11.89058563, 12.01486743, 12.21133227, 12.40782023, 12.55682327, 
     12.7059072, 12.88910065, 13.07234341, 13.20139753, 13.33048759, 
     13.46960371, 13.60875466, 13.71956075, 13.83042354, 13.91157425, 
     13.99277241, 14.08548465, 14.17824161, 14.26091536, 14.34364144, 
     14.43504833, 14.52651725, 14.59548578, 14.66449439, 14.73751377, 
     14.81058218, 14.88687742, 14.96322511, 15.038278, 15.11342222, 
     15.1644832, 15.21556914, 15.31798775, 15.42044285, 15.51000031, 
     15.59959896, 15.68089263, 15.76229354, 15.86213633, 15.96214484, 
     16.0269976, 16.0923806, 16.16483146, 16.23762336, 16.24233076, 
     16.24719576, 16.26116812, 16.27543692, 16.28311024, 16.29128992, 
     16.23458771, 16.17830522, 16.12636211, 16.07478778]
y = np.array(y)[::-1]

step = 50e9
x = np.arange(192.1e12, 196e12+step, step)
x2 = np.arange(191.325e12, 196.125e12+step, step)

f = interp1d(x, y, kind="cubic", bounds_error=False, fill_value="extrapolate")

fig, ax = plt.subplots()
ax.plot(x, y, ".", label="data")
ax.plot(x2, f(x2), "--", label="fit")
ax.legend()
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.show()

编辑：如果你确实需要多项式拟合，那么使用np.poly1d。这返回一个可以被计算的函数，不像polyval只返回值（因此您得到的错误）。

import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

y = [9.996651009, 10.13327408, 10.27003353, 10.43754263, 10.60523071, 
     10.71136134, 10.81753489, 10.92299416, 11.02853271, 11.10093835, 
     11.17349619, 11.16565259, 11.15787539, 11.16191795, 11.16626826, 
     11.16087842, 11.15561384, 11.16940196, 11.18346164, 11.22434472, 
     11.26537218, 11.35981126, 11.45427174, 11.6102862, 11.76633923, 
     11.89058563, 12.01486743, 12.21133227, 12.40782023, 12.55682327, 
     12.7059072, 12.88910065, 13.07234341, 13.20139753, 13.33048759, 
     13.46960371, 13.60875466, 13.71956075, 13.83042354, 13.91157425, 
     13.99277241, 14.08548465, 14.17824161, 14.26091536, 14.34364144, 
     14.43504833, 14.52651725, 14.59548578, 14.66449439, 14.73751377, 
     14.81058218, 14.88687742, 14.96322511, 15.038278, 15.11342222, 
     15.1644832, 15.21556914, 15.31798775, 15.42044285, 15.51000031, 
     15.59959896, 15.68089263, 15.76229354, 15.86213633, 15.96214484, 
     16.0269976, 16.0923806, 16.16483146, 16.23762336, 16.24233076, 
     16.24719576, 16.26116812, 16.27543692, 16.28311024, 16.29128992, 
     16.23458771, 16.17830522, 16.12636211, 16.07478778]
y = np.array(y)[::-1]

step = 50e9
x = np.arange(192.1e12, 196e12+step, step)
x2 = np.arange(191.325e12, 196.125e12+step, step)

z = np.polyfit(x, y, 20)
p = np.poly1d(z)

fig, ax = plt.subplots()
ax.plot(x, y, ".", label="data")
ax.plot(x2, p(x2), "--", label="fit")
ax.legend()
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.show()

编辑：使用更新的np.polynomial.Polynomial方法，我们使用5阶多项式拟合得到了类似的结果。

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

plt.close("all")

y = [9.996651009, 10.13327408, 10.27003353, 10.43754263, 10.60523071, 
     10.71136134, 10.81753489, 10.92299416, 11.02853271, 11.10093835, 
     11.17349619, 11.16565259, 11.15787539, 11.16191795, 11.16626826, 
     11.16087842, 11.15561384, 11.16940196, 11.18346164, 11.22434472, 
     11.26537218, 11.35981126, 11.45427174, 11.6102862, 11.76633923, 
     11.89058563, 12.01486743, 12.21133227, 12.40782023, 12.55682327, 
     12.7059072, 12.88910065, 13.07234341, 13.20139753, 13.33048759, 
     13.46960371, 13.60875466, 13.71956075, 13.83042354, 13.91157425, 
     13.99277241, 14.08548465, 14.17824161, 14.26091536, 14.34364144, 
     14.43504833, 14.52651725, 14.59548578, 14.66449439, 14.73751377, 
     14.81058218, 14.88687742, 14.96322511, 15.038278, 15.11342222, 
     15.1644832, 15.21556914, 15.31798775, 15.42044285, 15.51000031, 
     15.59959896, 15.68089263, 15.76229354, 15.86213633, 15.96214484, 
     16.0269976, 16.0923806, 16.16483146, 16.23762336, 16.24233076, 
     16.24719576, 16.26116812, 16.27543692, 16.28311024, 16.29128992, 
     16.23458771, 16.17830522, 16.12636211, 16.07478778]
y = np.array(y)[::-1]

step = 50e9
x = np.arange(192.1e12, 196e12+step, step)
x2 = np.arange(191.325e12, 196.125e12+step, step)

p = Polynomial.fit(x, y, 5)

fig, ax = plt.subplots()
ax.plot(x, y, ".", label="data")
ax.plot(x2, p(x2), "--", label="fit")
ax.legend()
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.show()

赞(0）分享  回复(0）举报  2023-06-23
qncylg1j

qncylg1j2#

多项式会疯狂。然而，我想知道这个错误信息，我得到的一般，我也想做一个比较的精度为不同阶的多项式拟合
花点时间考虑一下多项式拟合代码在拟合多项式时会有哪些中间值是很有用的。你的X值大约为1.9 * 10^14。如果你再取它的22次方，你得到的是2.8 * 10^314。这是bigger than the maximum representable double-precision float。
有趣的是，即使你的X^20项符合浮点数的范围，它似乎仍然受到浮点精度问题的困扰。以下是它适合20次多项式的系数：

>>> print(z)
[ 7.01389718e+008 -7.23272471e-006 -6.18262251e-021  1.28113394e-034
  6.59034262e-049 -5.97742101e-066 -1.75197282e-077 -9.00617222e-092
  1.16972981e-106  3.58691274e-120 -9.24694178e-135  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000]

换句话说，它对超过X^10的所有值拟合的系数都为零。这是由于浮点精度问题造成的。
如何解决这个问题？你可以从机器学习从业者那里获得灵感，并将你的输入标准化。“Standardization”是指减去平均值并除以输入的标准差。由于这是线性变换，因此在数学上等同于没有标准化的拟合，同时避免了溢出问题。当我这样做时，它为X^22选择的系数是-0.11，这相当于没有标准化的-5.3 * 10^-316。如果没有标准化，这个数字就不能表示为浮点数。
下面是我使用的代码：

import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly
from sklearn.preprocessing import StandardScaler



y=[9.996651009,10.13327408,10.27003353,10.43754263,10.60523071,10.71136134,10.81753489,10.92299416,11.02853271,11.10093835,11.17349619,11.16565259,11.15787539,11.16191795,11.16626826,11.16087842,11.15561384,11.16940196,11.18346164,11.22434472,11.26537218,11.35981126,11.45427174,11.6102862,11.76633923,11.89058563,12.01486743,12.21133227,12.40782023,12.55682327,12.7059072,12.88910065,13.07234341,13.20139753,13.33048759,13.46960371,13.60875466,13.71956075,13.83042354,13.91157425,13.99277241,14.08548465,14.17824161,14.26091536,14.34364144,14.43504833,14.52651725,14.59548578,14.66449439,14.73751377,14.81058218,14.88687742,14.96322511,15.038278,15.11342222,15.1644832,15.21556914,15.31798775,15.42044285,15.51000031,15.59959896,15.68089263,15.76229354,15.86213633,15.96214484,16.0269976,16.0923806,16.16483146,16.23762336,16.24233076,16.24719576,16.26116812,16.27543692,16.28311024,16.29128992,16.23458771,16.17830522,16.12636211,16.07478778]

y.reverse()

# x axis for polyfit
i = 192.1e12
x=[]
while i <= 196e12:
    x.append(i)
    i += 50e9

plt.plot(x, y)

# x axis for extrapolation
i1= 191.325e12
x2=[]
while i1 <= 196.125e12:

    x2.append(i1)
    i1 += 50e9

x = np.array(x)
x2 = np.array(x2)

scaler = StandardScaler()
scaler.fit(x.reshape(-1, 1))
x_rescaled = scaler.transform(x.reshape(-1, 1))
x2_rescaled = scaler.transform(x2.reshape(-1, 1))

z = poly.polyfit(x_rescaled.reshape(-1), y, 20)
f = poly.polyval(x2_rescaled.reshape(-1), z)

plt.plot(x2, f, linestyle='--')
plt.ylim([7, 20])

我应该提到的是，在20次多项式上进行外推很可能在您拟合的区域之外是无用的。这种次数的多项式对微小的变化极为敏感。
例如，我添加了幅度为0.01的随机噪声，然后重新拟合多项式。这个过程我重复了九次。在视觉上，您几乎看不到噪声，但噪声对多项式拟合有巨大影响。结果多项式的一半向上，一半向下。

参见Why is the use of high order polynomials for regression discouraged?

赞(0）分享  回复(0）举报  2023-06-23
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备2022004421号-2 NULL123 https://www.null123.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com