SciPy Curve_fit() 不适合曲线
SciPy Curve_fit() doesn't fit curve
我做了一个随机图，并尝试使用 SciPy curve_fit 来拟合最佳曲线，但是失败了。
首先，我生成了一个随机指数衰减图，其中 A, w, T2 是使用 numpy 随机生成的：
def expDec(t, A, w, 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, A, w, 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, A, w, 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, A, w, 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, t, 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'], x['Output'])
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, w, T2 = np.random.uniform(3, 6, size=3)
dom = np.linspace(-9, 9, 1000)
t = np.arange(-9., 10.)
hr = exp_dec(t, A, w, T2)
fit, _ = curve_fit(exp_dec, t, hr)
fig, ax = plt.subplots()
ax.plot(dom, exp_dec(dom, A, w, T2), 'g', label='target')
ax.scatter(t, hr, c='r', label='samples')
ax.plot(dom, exp_dec(dom, *fit), 'b', label='fit')
ax.plot(dom, exp_dec(dom, 1, 1, 1), 'k:', label='start')
ax.legend()
要解释最后绘制的项目，请查看 curve_fit 的文档。请注意，有一个参数 p0，如果您不提供该参数，则默认为所有参数。这是您的拟合开始猜测值的初始猜测。
看这张图，大概就知道问题出在哪里了。开始猜测的频率比您的数据低得多。因为采样频率非常接近振荡频率，所以在能够充分增加频率以获得正确函数之前，拟合会达到局部最小值。您可以通过几种不同的方式解决此问题。
一种方法是给 curve_fit 一个更好的初始猜测。如果您知道振幅、频率和衰减率的界限，请使用它们。振幅通常是一个简单的线性拟合。最困难的通常是频率，正如您在这里看到的，最好高估它。但是如果你高估了太多，你最终可能会得到原始数据的谐波。
这里有几个样本拟合，它们在优化中显示了不同的局部最小值。第二个显示了高估振荡频率的谐波情况：
一组合适的起始参数是您随机范围的上限：
fit, _ = curve_fit(exp_dec, t, hr, p0=[6, 6, 6])
绿色曲线与蓝色曲线如此接近，你看不到它：
>>> A, w, T2
(4.123620356542087, 5.852142919229749, 5.195981825434215)
>>> tuple(fit)
(4.123620356542086, 5.852142919229749, 5.195981825434215)
解决该问题的另一种方法是更频繁地对数据进行采样。更多的数据通常意味着在优化中达到错误的局部最小值的可能性更低。然而，在处理正弦函数时，由于匹配的工作方式，这并不总是有帮助。这是一个样本数量为 10 倍的示例（仅适合 2 倍且默认猜测完全失败）：
...
t = np.arange(-9., 10., 0.1)
...
我做了一个随机图，并尝试使用 SciPy curve_fit 来拟合最佳曲线，但是失败了。
首先，我生成了一个随机指数衰减图，其中 A, w, T2 是使用 numpy 随机生成的：
def expDec(t, A, w, 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, A, w, 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, A, w, 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, A, w, 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, t, 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'], x['Output'])
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, w, T2 = np.random.uniform(3, 6, size=3)
dom = np.linspace(-9, 9, 1000)
t = np.arange(-9., 10.)
hr = exp_dec(t, A, w, T2)
fit, _ = curve_fit(exp_dec, t, hr)
fig, ax = plt.subplots()
ax.plot(dom, exp_dec(dom, A, w, T2), 'g', label='target')
ax.scatter(t, hr, c='r', label='samples')
ax.plot(dom, exp_dec(dom, *fit), 'b', label='fit')
ax.plot(dom, exp_dec(dom, 1, 1, 1), 'k:', label='start')
ax.legend()
要解释最后绘制的项目，请查看 curve_fit 的文档。请注意，有一个参数 p0，如果您不提供该参数，则默认为所有参数。这是您的拟合开始猜测值的初始猜测。
看这张图，大概就知道问题出在哪里了。开始猜测的频率比您的数据低得多。因为采样频率非常接近振荡频率，所以在能够充分增加频率以获得正确函数之前，拟合会达到局部最小值。您可以通过几种不同的方式解决此问题。
一种方法是给 curve_fit 一个更好的初始猜测。如果您知道振幅、频率和衰减率的界限，请使用它们。振幅通常是一个简单的线性拟合。最困难的通常是频率，正如您在这里看到的，最好高估它。但是如果你高估了太多，你最终可能会得到原始数据的谐波。
这里有几个样本拟合，它们在优化中显示了不同的局部最小值。第二个显示了高估振荡频率的谐波情况：
一组合适的起始参数是您随机范围的上限：
fit, _ = curve_fit(exp_dec, t, hr, p0=[6, 6, 6])
绿色曲线与蓝色曲线如此接近，你看不到它：
>>> A, w, T2
(4.123620356542087, 5.852142919229749, 5.195981825434215)
>>> tuple(fit)
(4.123620356542086, 5.852142919229749, 5.195981825434215)
解决该问题的另一种方法是更频繁地对数据进行采样。更多的数据通常意味着在优化中达到错误的局部最小值的可能性更低。然而，在处理正弦函数时，由于匹配的工作方式，这并不总是有帮助。这是一个样本数量为 10 倍的示例（仅适合 2 倍且默认猜测完全失败）：
...
t = np.arange(-9., 10., 0.1)
...