问题描述
我制作了一个随机图,并尝试使用 SciPy curve_fit 将最佳曲线拟合到绘图中,但它失败了。
首先,我生成了一个随机指数衰减图,其中 A,w,T2
是使用 numpy 随机生成的:
def expDec(t,A,T2):
return A * np.cos(w * t) * (2.718**(-t / T2))
现在我让 SciPy 猜测最佳拟合曲线:
t = x['Input'].values
hr = x['Output'].values
c,cov = curve_fit(bpm,t,hr)
然后我绘制曲线
for i in range(n):
y[i] = bpm(x['Input'][i],c[0],c[1],c[2])
plt.plot(x['Input'],x['Output'])
plt.plot(x['Input'],y)
就是这样。这是合身看起来有多糟糕:
如果有人能帮忙,那就太好了。
MWE(也可交互使用 here)
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
inputs = []
outputs = []
# THIS GIVES THE DOMAIN
dom = np.linspace(-5,5,100)
# FUNCTION & ParaMETERS (RANDOMLY SELECTED)
A = np.random.uniform(3,6)
w = np.random.uniform(3,6)
T2 = np.random.uniform(3,6)
y = A * np.cos(w * dom) * (2.718**(-dom / T2))
# DEFInes EXPONENTIAL DECAY FUNCTION
def expDec(t,T2):
return A * np.cos(w * t) * (2.718**(-t / T2))
# SETS UP figURE FOR PLottING
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
# PLOTS THE FUNCTION
plt.plot(dom,y,'r')
# SHOW THE PLOT
plt.show()
for i in range(-9,10):
inputs.append(i)
outputs.append(expDec(i,T2))
# PUT IT DIRECTLY IN A PANDAS DATAFRAME
points = {'Input': inputs,'Output': outputs}
x = pd.DataFrame(points,columns = ['Input','Output'])
# FUNCTION WHOSE ParaMETERS PROGRAM SHOULD BE GUESSING
def bpm(t,T2):
return A * np.cos(w * t) * (2.718**(-t / T2))
# INPUT & OUTPUTS
t = x['Input'].values
hr = x['Output'].values
# USE SCIPY CURVE FIT TO USE NONLINEAR LEAST SQUARES TO FIND BEST ParaMETERS. TRY 1000 TIMES BEFORE STOPPING.
constants = curve_fit(bpm,hr,maxfev=1000)
# GET CONSTANTS FROM CURVE_FIT
A_fit = constants[0][0]
w_fit = constants[0][1]
T2_fit = constants[0][2]
# CREATE ARRAY TO HOLD FITTED OUTPUT
fit = []
# APPEND OUTPUT TO FIT=[] ARRAY
for i in range(-9,10):
fit.append(bpm(i,A_fit,w_fit,T2_fit))
# PLOTS BEST ParaMETERS
plt.plot(x['Input'],fit,"ro-")
解决方法
作为第一步,我想重写您的 MCVE 以使用矢量化操作和函数计算的单个实例。这会将所有内容减少到几行。我建议您在进行测试时也使用种子以提高可重复性:
def exp_dec(t,A,w,T2):
return A * np.cos(w * t) * np.exp(-t / T2)
np.random.seed(42)
A,T2 = np.random.uniform(3,6,size=3)
dom = np.linspace(-9,9,1000)
t = np.arange(-9.,10.)
hr = exp_dec(t,T2)
fit,_ = curve_fit(exp_dec,t,hr)
fig,ax = plt.subplots()
ax.plot(dom,exp_dec(dom,T2),'g',label='target')
ax.scatter(t,hr,c='r',label='samples')
ax.plot(dom,*fit),'b',label='fit')
ax.plot(dom,1,1),'k:',label='start')
ax.legend()
要解释最后绘制的项目,请查看 curve_fit
的文档。请注意,有一个参数 p0
,如果您不提供它,则默认为所有参数。这是您的拟合开始猜测值的初始猜测。
看看这张图片,你几乎可以看出问题所在。起始猜测的频率比您的数据低得多。因为采样频率非常接近振荡频率,所以在能够充分增加频率以获得正确的函数之前,拟合会达到局部最小值。您可以通过几种不同的方式解决此问题。
一种方法是给 curve_fit
一个更好的初始猜测。如果您知道幅度、频率和衰减率的界限,请使用它们。幅度通常是直接的线性拟合。最难的通常是频率,正如您在此处看到的那样,最好高估它。但如果你过高估计它,你最终可能会得到原始数据的谐波。
这里有几个样本拟合,它们显示了优化中的不同局部最小值。第二个显示了高估振荡频率的谐波情况:
一组合适的起始参数是您的随机范围的上限:
fit,p0=[6,6])
绿色曲线与蓝色非常接近,你看不到:
>>> A,T2
(4.123620356542087,5.852142919229749,5.195981825434215)
>>> tuple(fit)
(4.123620356542086,5.195981825434215)
解决问题的另一种方法是更频繁地对数据进行采样。更多的数据通常意味着在优化中达到错误的局部最小值的机会更低。但是,在处理正弦函数时,由于匹配的工作方式,这并不总是有帮助。这是一个示例数量为 10 倍的示例(仅使用 2 倍拟合并且默认猜测完全失败):
...
t = np.arange(-9.,10.,0.1)
...