numpy 多项式拟合并求解系数

ccrfmcuu  于 2023-08-05  发布在  其他
关注(0)|答案(2)|浏览(103)

我已经绘制了'energy 1'与'amp'的图。我写了下面的代码,a)为数据集绘制一个四阶多项式拟合曲线B)得到多项式的系数
但是,函数f1打印为“f1= 29.16 x^4 + 9.809e-07 x^3 - 73.08 x^2 - 1.84e-06 x + 0.2619”
并且系数打印为“coefficient= [ 2.91632913e+01 9.80928536e-07 -7.30843650e+01 -1.84038399e-06 2. 61942115 e-01]”
1.我怎样才能写出只有“偶次幂”项的函数,然后得到它们的系数?
1.如何重新检查/评估系数是否为特定的“振幅”值提供了正确的“能量1”值(这将再次与我的“能量1”数据相匹配)?

amp=(['amplitude'])
energy1=(name['solar'])
p1 = np.polyfit(amp, energy1, 4) #polynomial
f1 = np.poly1d(p1) #function
x1= np.linspace(-1.6, 1.6, 100)
plt.plot(amp, energy1, 'b^', label="solar", markersize=12)
plt.plot(x1, f1(x1), 'b', linewidth=2.0)

print ("f1=", f1)
c1 = np.polyfit(x1, f1(x1), 4)
print ("coefficient=", c1)

字符串

snvhrwxg

snvhrwxg1#

当涉及到适配自定义函数时,我建议使用scipy.optimize.curve_fit。定义函数时,第一个参数是x,其余参数是希望scipy优化的变量。由于雅可比矩阵在这种情况下很容易(你可以很容易地对函数的每个参数求导),所以也可以添加它来帮助优化。返回的第一部分,称为popt,包含优化的变量。
拟合将被优化,但并不精确(特别是在处理真实世界的数据时),所以我不确定您希望如何比较结果。在下面的代码中,我计算均方误差,并显示实际数据与拟合的关系图。根据文档,您还可以使用第二个返回参数,称为pcov“来计算参数的一个标准差误差”。

import numpy as np
from scipy.optimize import curve_fit

rng = np.random.default_rng(42)
N = 20
x = np.linspace(0, 2, N)
a, b, c = rng.uniform(-5, 5, 3)
y = a*x**4 + b*x**2 + c
y += rng.uniform(-2, 2, N)

def func(x, a, b, c):
    return a*x**4 + b*x**2 + c

def jac(x, a, b, c):
    return np.stack([x**4, x**2, 1*np.ones_like(x)], axis=1)

popt, pcov = curve_fit(func, x, y, jac=jac)

print(f"Actual: {a}, {b}, {c}")
print(f"Fit: {popt[0]}, {popt[1]}, {popt[2]}")

MSE = np.mean((y - func(x, *popt))**2)
print(f"MSE: {MSE}")

perr = np.sqrt(np.diag(pcov))
print(f"perr: {perr}")

字符串
输出量:

Actual: 2.7395604855596334, -0.6112156024794766, 3.585979199113824
Fit: 2.9453817011693193, -1.26046899494815, 4.061262973609011
MSE: 1.2228536711027866
perr: [0.19976097 0.74404971 0.48846652]


的数据

wmomyfyw

wmomyfyw2#

我提供了@jared的(正确的)非线性拟合解决方案的替代方案。
对于某些类型的数据,线性化可能更快或产生更稳定的数值拟合。这个想法是形成一个系数矩阵,其中自变量已经提升到您关心的幂(0,2,4),然后执行最小二乘线性矩阵求解。我鼓励您在真实的数据上比较此方法与非线性方法的运行时间和准确性。

import matplotlib.pyplot as plt
import numpy as np
from numpy.polynomial import Polynomial
from numpy.random._generator import default_rng

rand = default_rng(seed=0)
amp = np.sort(rand.uniform(low=-1.6, high=1.6, size=100))
artificial_coeffs = np.array([  # highest to lowest
     2.91632913e+01, 0,
    -7.30843650e+01, 0,
     2.61942115e-01,
])
artificial_poly = Polynomial(artificial_coeffs[::-1])  # lowest to highest
print('Target:', artificial_poly)
energy1 = artificial_poly(amp) + rand.normal(loc=0, scale=3, size=amp.size)

#  100x3           3x1          100x1
# [1 amp^2 amp^4] [f0 f2 f4] = [energy]
# [...          ]              [...   ]
fit = np.zeros_like(artificial_coeffs)
fit[::2], *_ = np.linalg.lstsq(
    a=amp[:, np.newaxis] ** ((0, 2, 4),),
    b=energy1, rcond=None,
)
polyfit = Polynomial(coef=fit)
print('  Fit:', polyfit)

fig, ax = plt.subplots()
ax.scatter(amp, energy1, label='experiment')
ax.plot(amp, polyfit(amp), label='fit', color='orange')
ax.legend()
plt.show()

个字符
x1c 0d1x的数据

相关问题