如何确定我的曲线应该是哪个模型并对其进行近似？
How to determine by which model my curve should be and approximate it?
我有一条曲线：enter image description here
及其数字化数据。 https://drive.google.com/open?id=1ZB39G3SmtamjVjmLzkC2JefloZ9iShpO
我应该如何为最小二乘法选择合适的函数以及如何在 Python 上实现这种近似？
我试过那样做
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym
x = np.loadtxt("x_data.txt", delimiter='\t', dtype=np.float)
y = np.loadtxt("y_data.txt", delimiter='\t', dtype=np.float)
plt.plot(x, y, 'ro',label="Original Data")
x = np.array(x, dtype=float)
y = np.array(y, dtype=float)
def func(x, a, b, c, d):
return a*x**3 + b*x**2 +c*x + d
popt, pcov = curve_fit(func, x, y)
xs = sym.Symbol('\lambda')
tex = sym.latex(func(xs,*popt)).replace('$', '')
plt.title(r'$f(\lambda)= %s$' %(tex),fontsize=16)
plt.plot(x, func(x, *popt), label="Fitted Curve")
plt.legend(loc='upper left')
plt.show()
谢谢
我从你的图中提取了数据进行分析，这是我对这个问题的第一次切入。由于您的绘图使用了十年对数缩放，因此我采用了 "left-side" Y 数据的反十年对数进行拟合。我对这个 "extracted" 数据的峰值方程搜索出现了一个对数正态方程，这里是一个图形化的 python 拟合器读取数据文件并拟合这个方程。此拟合器使用 scipy differential_evolution 遗传算法模块来确定非线性拟合器的初始参数估计值，这需要在其中进行搜索的参数范围。在此代码中，数据最大值和最小值与我的范围估计值一起使用。在初始参数值上估计 ranges 比具体值要容易得多。这个钳工应该能够直接读取您的数据文件。如果您可以 post 或 link 到实际数据，我可能会比这里显示的更合适。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
with open('./x_data.txt', 'rt') as f:
x_file = f.read()
with open('./y_data.txt', 'rt') as f:
y_file = f.read()
xlist = []
for line in x_file.split('\n'):
if line: # this allows blank lines in file
xlist.append(float(line.strip()))
ylist = []
for line in y_file.split('\n'):
if line: # this allows blank lines in file
ylist.append(float(line.strip()))
if len(xlist) != len(ylist):
print(len(xlist), len(ylist))
raise Exception('X and Y habe different length')
xData = numpy.array(xlist)
yData = numpy.array(ylist)
def func(t, a, b, c, d): # Log-Normal Peak A Shifted from zunzun.com
return a * numpy.exp(-0.5 * numpy.power((numpy.log(t-d)-b) / c, 2.0))
# 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)
parameterBounds = []
parameterBounds.append([minY, maxY]) # search bounds for a
parameterBounds.append([0.0, 2.0]) # search bounds for b
parameterBounds.append([-1.0, 0.0]) # search bounds for c
parameterBounds.append([-maxX, 0.0]) # search bounds for d
# "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), 1000)
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)
编辑：使用实际数据
根据现在可用的实际数据，我发现脉冲峰值方程可以很好地拟合数据，这是一个更新的拟合器。如果可能的话，我建议从时间 0 到 5 获取额外的数据，这将产生更好地表征峰值开始时区域的数据。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
with open('./x_data.txt', 'rt') as f:
x_file = f.read()
with open('./y_data.txt', 'rt') as f:
y_file = f.read()
xlist = []
for line in x_file.split('\n'):
if line: # this allows blank lines in file
xlist.append(float(line.strip()))
ylist = []
for line in y_file.split('\n'):
if line: # this allows blank lines in file
ylist.append(float(line.strip()))
if len(xlist) != len(ylist):
print(len(xlist), len(ylist))
raise Exception('X and Y have different length')
xData = numpy.array(xlist)
yData = numpy.array(ylist)
def func(t, a, b, c, Offset): # Pulse Peak With Offset from zunzun.com
return 4.0 * a * numpy.exp(-1.0 * (t-b) / c) * (1.0 - numpy.exp(-1.0 * (t-b) / c)) + 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)
parameterBounds = []
parameterBounds.append([minY, maxY]) # search bounds for a
parameterBounds.append([-5.0, 0.0]) # search bounds for b
parameterBounds.append([1.0, 10.0]) # search bounds for c
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), 1000)
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)
我有一条曲线：enter image description here
及其数字化数据。 https://drive.google.com/open?id=1ZB39G3SmtamjVjmLzkC2JefloZ9iShpO
我应该如何为最小二乘法选择合适的函数以及如何在 Python 上实现这种近似？
我试过那样做
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym
x = np.loadtxt("x_data.txt", delimiter='\t', dtype=np.float)
y = np.loadtxt("y_data.txt", delimiter='\t', dtype=np.float)
plt.plot(x, y, 'ro',label="Original Data")
x = np.array(x, dtype=float)
y = np.array(y, dtype=float)
def func(x, a, b, c, d):
return a*x**3 + b*x**2 +c*x + d
popt, pcov = curve_fit(func, x, y)
xs = sym.Symbol('\lambda')
tex = sym.latex(func(xs,*popt)).replace('$', '')
plt.title(r'$f(\lambda)= %s$' %(tex),fontsize=16)
plt.plot(x, func(x, *popt), label="Fitted Curve")
plt.legend(loc='upper left')
plt.show()
谢谢
我从你的图中提取了数据进行分析，这是我对这个问题的第一次切入。由于您的绘图使用了十年对数缩放，因此我采用了 "left-side" Y 数据的反十年对数进行拟合。我对这个 "extracted" 数据的峰值方程搜索出现了一个对数正态方程，这里是一个图形化的 python 拟合器读取数据文件并拟合这个方程。此拟合器使用 scipy differential_evolution 遗传算法模块来确定非线性拟合器的初始参数估计值，这需要在其中进行搜索的参数范围。在此代码中，数据最大值和最小值与我的范围估计值一起使用。在初始参数值上估计 ranges 比具体值要容易得多。这个钳工应该能够直接读取您的数据文件。如果您可以 post 或 link 到实际数据，我可能会比这里显示的更合适。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
with open('./x_data.txt', 'rt') as f:
x_file = f.read()
with open('./y_data.txt', 'rt') as f:
y_file = f.read()
xlist = []
for line in x_file.split('\n'):
if line: # this allows blank lines in file
xlist.append(float(line.strip()))
ylist = []
for line in y_file.split('\n'):
if line: # this allows blank lines in file
ylist.append(float(line.strip()))
if len(xlist) != len(ylist):
print(len(xlist), len(ylist))
raise Exception('X and Y habe different length')
xData = numpy.array(xlist)
yData = numpy.array(ylist)
def func(t, a, b, c, d): # Log-Normal Peak A Shifted from zunzun.com
return a * numpy.exp(-0.5 * numpy.power((numpy.log(t-d)-b) / c, 2.0))
# 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)
parameterBounds = []
parameterBounds.append([minY, maxY]) # search bounds for a
parameterBounds.append([0.0, 2.0]) # search bounds for b
parameterBounds.append([-1.0, 0.0]) # search bounds for c
parameterBounds.append([-maxX, 0.0]) # search bounds for d
# "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), 1000)
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)
编辑：使用实际数据
根据现在可用的实际数据，我发现脉冲峰值方程可以很好地拟合数据，这是一个更新的拟合器。如果可能的话，我建议从时间 0 到 5 获取额外的数据，这将产生更好地表征峰值开始时区域的数据。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings
with open('./x_data.txt', 'rt') as f:
x_file = f.read()
with open('./y_data.txt', 'rt') as f:
y_file = f.read()
xlist = []
for line in x_file.split('\n'):
if line: # this allows blank lines in file
xlist.append(float(line.strip()))
ylist = []
for line in y_file.split('\n'):
if line: # this allows blank lines in file
ylist.append(float(line.strip()))
if len(xlist) != len(ylist):
print(len(xlist), len(ylist))
raise Exception('X and Y have different length')
xData = numpy.array(xlist)
yData = numpy.array(ylist)
def func(t, a, b, c, Offset): # Pulse Peak With Offset from zunzun.com
return 4.0 * a * numpy.exp(-1.0 * (t-b) / c) * (1.0 - numpy.exp(-1.0 * (t-b) / c)) + 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)
parameterBounds = []
parameterBounds.append([minY, maxY]) # search bounds for a
parameterBounds.append([-5.0, 0.0]) # search bounds for b
parameterBounds.append([1.0, 10.0]) # search bounds for c
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), 1000)
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)