Scipy curve_fit 给出错误答案
Scipy curve_fit gives wrong answer
我有一个如下图所示的振荡数据,我想为其拟合一条正弦曲线。但是,我的结果不正确。
我想拟合这条曲线的函数是:
def radius (z,phi, a0, k0,):
Z = z.reshape(z.shape[0],1)
k = np.array([k0,])
a = np.array([a0,])
r0 = 110
rs = r0 + np.sum(a*np.sin(k*Z +phi), axis=1)
return rs
正确的解决方案可能如下所示:
r_fit = radius(z, phi=np.pi/.8, a0=10,k0=0.017)
plt.plot(z, r, label='data')
plt.plot(z, r_fit, label='fitted curve')
plt.legend()
但是我拟合曲线的结果看起来是:
from scipy.optimize import curve_fit
popt, pcov = curve_fit(radius, xdata=z, ydata=r)
r_fit = radius(z, *popt)
plt.plot(z, r, label='data')
plt.plot(z, r_fit, label='fitted curve')
plt.legend()
我的数据也如下:
r = np.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
113.12440802, 119.30644014, 119.86570527, 119.75184665,
117.12160143, 101.55081608, 100.07280857, 100.12880236,
100.39251753, 103.05404178, 117.15257288, 119.74048706,
119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
100.05521567])
z = np.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
-224.36240585, -188.55933945, -155.21242348, -122.02778866,
-87.84335638, -47.0274899 , 0. , 47.54559191,
94.97469981, 141.33801462, 181.59490575, 215.77219256,
248.95956379, 282.28027286, 318.16440024, 360.7246922 ,
407.940799 ])
因为我的函数只代表傅里叶级数,所以我也尝试了 scipy.fftpack.fft(r) 但我无法重现与我计算 fft 的信号接近的信号。
问题是如果不提供初始猜测,解就无法收敛。尝试添加一个合理的初始猜测:
p0 = [np.pi/.8, 10, 0.017]
popt, pcov = curve_fit(radius, xdata=z, ydata=r, p0=p0)
请注意,如果您要使用其他方法之一,例如 trf
或 dogbox
,那么在没有最初猜测的情况下,这更有可能 return 由于运行时错误参数无法收敛。
这是一个带有正弦方程的图形化 Python 拟合器,您的数据使用 scipy.optimize 差分进化遗传算法模块来确定 curve_fit 非线性的初始参数估计求解器。该 scipy 模块使用拉丁超立方体算法来确保对参数 space 的彻底搜索需要搜索范围。在此示例中,这些边界取自数据的最大值和最小值。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
r = numpy.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
113.12440802, 119.30644014, 119.86570527, 119.75184665,
117.12160143, 101.55081608, 100.07280857, 100.12880236,
100.39251753, 103.05404178, 117.15257288, 119.74048706,
119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
100.05521567])
z = numpy.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
-224.36240585, -188.55933945, -155.21242348, -122.02778866,
-87.84335638, -47.0274899 , 0. , 47.54559191,
94.97469981, 141.33801462, 181.59490575, 215.77219256,
248.95956379, 282.28027286, 318.16440024, 360.7246922 ,
407.940799 ])
# rename data to match previous example code
xData = z
yData = r
def func (x, amplitude, center, width, offset): # equation sine[radians] + offset from zunzun.com
return amplitude * numpy.sin(numpy.pi * (x - center) / width) + offset
# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
val = func(xData, *parameterTuple)
return numpy.sum((yData - val) ** 2.0)
def generate_Initial_Parameters():
# min and max used for bounds
maxX = max(xData)
minX = min(xData)
maxY = max(yData)
minY = min(yData)
diffY = maxY - minY
diffX = maxX - minX
parameterBounds = []
parameterBounds.append([0.0, diffY]) # search bounds for amplitude
parameterBounds.append([minX, maxX]) # search bounds for center
parameterBounds.append([0.0, diffX]) # search bounds for width
parameterBounds.append([minY, maxY]) # search bounds for offset
# "seed" the numpy random number generator for repeatable results
result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
return result.x
# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()
# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()
modelPredictions = func(xData, *fittedParameters)
absError = modelPredictions - yData
SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))
print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
axes = f.add_subplot(111)
# first the raw data as a scatter plot
axes.plot(xData, yData, 'D')
# create data for the fitted equation plot
xModel = numpy.linspace(min(xData), max(xData))
yModel = func(xModel, *fittedParameters)
# now the model as a line plot
axes.plot(xModel, yModel)
axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label
plt.show()
plt.close('all') # clean up after using pyplot
graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)
我有一个如下图所示的振荡数据,我想为其拟合一条正弦曲线。但是,我的结果不正确。
我想拟合这条曲线的函数是:
def radius (z,phi, a0, k0,):
Z = z.reshape(z.shape[0],1)
k = np.array([k0,])
a = np.array([a0,])
r0 = 110
rs = r0 + np.sum(a*np.sin(k*Z +phi), axis=1)
return rs
正确的解决方案可能如下所示:
r_fit = radius(z, phi=np.pi/.8, a0=10,k0=0.017)
plt.plot(z, r, label='data')
plt.plot(z, r_fit, label='fitted curve')
plt.legend()
但是我拟合曲线的结果看起来是:
from scipy.optimize import curve_fit
popt, pcov = curve_fit(radius, xdata=z, ydata=r)
r_fit = radius(z, *popt)
plt.plot(z, r, label='data')
plt.plot(z, r_fit, label='fitted curve')
plt.legend()
我的数据也如下:
r = np.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
113.12440802, 119.30644014, 119.86570527, 119.75184665,
117.12160143, 101.55081608, 100.07280857, 100.12880236,
100.39251753, 103.05404178, 117.15257288, 119.74048706,
119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
100.05521567])
z = np.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
-224.36240585, -188.55933945, -155.21242348, -122.02778866,
-87.84335638, -47.0274899 , 0. , 47.54559191,
94.97469981, 141.33801462, 181.59490575, 215.77219256,
248.95956379, 282.28027286, 318.16440024, 360.7246922 ,
407.940799 ])
因为我的函数只代表傅里叶级数,所以我也尝试了 scipy.fftpack.fft(r) 但我无法重现与我计算 fft 的信号接近的信号。
问题是如果不提供初始猜测,解就无法收敛。尝试添加一个合理的初始猜测:
p0 = [np.pi/.8, 10, 0.017]
popt, pcov = curve_fit(radius, xdata=z, ydata=r, p0=p0)
请注意,如果您要使用其他方法之一,例如 trf
或 dogbox
,那么在没有最初猜测的情况下,这更有可能 return 由于运行时错误参数无法收敛。
这是一个带有正弦方程的图形化 Python 拟合器,您的数据使用 scipy.optimize 差分进化遗传算法模块来确定 curve_fit 非线性的初始参数估计求解器。该 scipy 模块使用拉丁超立方体算法来确保对参数 space 的彻底搜索需要搜索范围。在此示例中,这些边界取自数据的最大值和最小值。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
r = numpy.array([100.09061214, 100.17932773, 100.45526772, 102.27891728,
113.12440802, 119.30644014, 119.86570527, 119.75184665,
117.12160143, 101.55081608, 100.07280857, 100.12880236,
100.39251753, 103.05404178, 117.15257288, 119.74048706,
119.86955437, 119.37452005, 112.83384329, 101.0507198 ,
100.05521567])
z = numpy.array([-407.90074345, -360.38004677, -312.99221012, -266.36934609,
-224.36240585, -188.55933945, -155.21242348, -122.02778866,
-87.84335638, -47.0274899 , 0. , 47.54559191,
94.97469981, 141.33801462, 181.59490575, 215.77219256,
248.95956379, 282.28027286, 318.16440024, 360.7246922 ,
407.940799 ])
# rename data to match previous example code
xData = z
yData = r
def func (x, amplitude, center, width, offset): # equation sine[radians] + offset from zunzun.com
return amplitude * numpy.sin(numpy.pi * (x - center) / width) + offset
# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
val = func(xData, *parameterTuple)
return numpy.sum((yData - val) ** 2.0)
def generate_Initial_Parameters():
# min and max used for bounds
maxX = max(xData)
minX = min(xData)
maxY = max(yData)
minY = min(yData)
diffY = maxY - minY
diffX = maxX - minX
parameterBounds = []
parameterBounds.append([0.0, diffY]) # search bounds for amplitude
parameterBounds.append([minX, maxX]) # search bounds for center
parameterBounds.append([0.0, diffX]) # search bounds for width
parameterBounds.append([minY, maxY]) # search bounds for offset
# "seed" the numpy random number generator for repeatable results
result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
return result.x
# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()
# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()
modelPredictions = func(xData, *fittedParameters)
absError = modelPredictions - yData
SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))
print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
axes = f.add_subplot(111)
# first the raw data as a scatter plot
axes.plot(xData, yData, 'D')
# create data for the fitted equation plot
xModel = numpy.linspace(min(xData), max(xData))
yModel = func(xModel, *fittedParameters)
# now the model as a line plot
axes.plot(xModel, yModel)
axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label
plt.show()
plt.close('all') # clean up after using pyplot
graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)