改进曲线拟合日志
Improve curve fitting log
我试着让我的曲线更合适。我的原始数据在一个 xlsx 文件中。我使用 pandas 提取它们。我想做两种不同的拟合，因为 Ra = 1e6 的行为发生了变化。我们知道 Ra 与 Nu**a 成正比。 a = 0.25 对于 Ra <1e6 如果不是 a = 0.33.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from scipy.optimize import curve_fit
import lmfit
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(data['Ra'].values, data['Nu_top'].values, label='Nu_top')
plt.scatter(data['Ra'].values, data['Nu_bottom'].values, label='Nu_bottom')
plt.errorbar(data['Ra'].values, data['Nu_top'].values , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(data['Ra'].values, data['Nu_bottom'].values , yerr=data['Ecart type bot'].values, linestyle="None")
def func(x,a):
return 10**(np.log10(x)/a)
"""maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]"""
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(data['Nu_top'].head(10).values, params, x=data['Ra'].head(10).values)
print(ret.fit_report())
popt, pcov = curve_fit(func, data['Ra'].head(10).values,
data['Nu_top'].head(10).values, sigma=data['Ecart type top'].head(10).values,
absolute_sigma=True, p0=[0.25])
plt.plot(data['Ra'].head(10).values, func(data['Ra'].head(10).values, *popt),
'r-', label='fit: a=%5.3f' % tuple(popt))
popt, pcov = curve_fit(func, data['Ra'].tail(4).values, data['Nu_top'].tail(4).values,
sigma=data['Ecart type top'].tail(4).values,
absolute_sigma=True, p0=[0.33])
plt.plot(data['Ra'].tail(4).values, func(data['Ra'].tail(4).values, *popt),
'b-', label='fit: a=%5.3f' % tuple(popt))
print(pcov)
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()
所以我使用日志：logRa = a * logNu。
Ra = x 轴
Nu = y 轴
这就是我以这种方式定义函数 func 的原因。
如您所见，我的两次合身并不完全正确。我的协方差等于 [0.00010971]。所以我不得不做错事，但我没有看到。我需要帮助。
这里的数据文件：
data.xlsx
我注意到 Ra 的数据值很大，在对它们进行缩放后，我执行了方程式搜索 - 这是我的代码结果。我使用标准 scipy 遗传算法模块 differential_evolution 来确定 curve_fit() 的初始参数值，并且该模块使用拉丁超立方体算法来确保彻底搜索参数 space 这需要搜索范围。给出初始参数估计的范围比找到特定值要容易得多。此等式适用于 nu_top 和 nu_bottom，请注意，图中未按对数比例缩放，因为在本例中没有必要。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import pandas
import warnings
filename = 'data.xlsx'
data=pandas.read_excel(filename,sheet_name='Sheet2',index=False,dtype={'Ra': float})
# notice the Ra scaling by 10000.0
xData = data['Ra'].values / 10000.0
yData = data['Nu_bottom']
def func(x, a, b, c): # "Combined Power And Exponential" from zunzun.com
return a * numpy.power(x, b) * numpy.exp(c * x)
# 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([0.0, 10.0]) # search bounds for a
parameterBounds.append([0.0, 10.0]) # search bounds for b
parameterBounds.append([0.0, 10.0]) # search bounds for c
# "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)
这里我把我的数据x和y放在log10()中。该图采用对数刻度。所以通常我应该有两个仿射函数，系数分别为 0.25 和 0.33。我改变了你的程序 James 中的函数 func 和 b 和 c 的界限，但我没有得到好的结果。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
def func(x,a):
return a*x
maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
popt, pcov = curve_fit(func, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), func(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
print(pcov)
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()
使用 polyfit 我得到了更好的结果。
使用我的代码，我打开文件并计算 log (Ra) 和 log (Nu)，然后以对数刻度绘制 (log (Ra)、log (Nu))。
对于 Ra <1e6，我应该有 a = 0.25，如果不是 a = 0.33
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from numpy import polyfit
import numpy.polynomial.polynomial as poly
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
x=np.log10(data['Ra'].values)
y1=np.log10(data['Nu_top'].values)
y2=np.log10(data['Nu_bottom'].values)
x2=np.log10(data['Ra'].head(11).values)
y4=np.log10(data['Nu_top'].head(11).values)
x3=np.log10(data['Ra'].tail(4).values)
y5=np.log10(data['Nu_top'].tail(4).values)
plt.xscale('log')
plt.yscale('log')
plt.scatter(x, y1, label='Nu_top')
plt.scatter(x, y2, label='Nu_bottom')
plt.errorbar(x, y1 , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(x, y2 , yerr=data['Ecart type bot'].values, linestyle="None")
"""a=np.ones(10, dtype=np.float)
weights = np.insert(a,0,1E10)"""
coefs = poly.polyfit(x2, y4, 1)
print(coefs)
ffit = poly.polyval(x2, coefs)
plt.plot(x2, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
absError = ffit - x2
SE = np.square(absError) # squared errors
MSE = np.mean(SE) # mean squared errors
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (np.var(absError) / np.var(x2))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
print('Predicted value at x=0:', ffit[0])
print()
coefs = poly.polyfit(x3, y5, 1)
ffit = poly.polyval(x3, coefs)
plt.plot(x3, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()
我的问题已解决，我设法用或多或少正确的结果拟合我的曲线
我试着让我的曲线更合适。我的原始数据在一个 xlsx 文件中。我使用 pandas 提取它们。我想做两种不同的拟合，因为 Ra = 1e6 的行为发生了变化。我们知道 Ra 与 Nu**a 成正比。 a = 0.25 对于 Ra <1e6 如果不是 a = 0.33.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from scipy.optimize import curve_fit
import lmfit
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(data['Ra'].values, data['Nu_top'].values, label='Nu_top')
plt.scatter(data['Ra'].values, data['Nu_bottom'].values, label='Nu_bottom')
plt.errorbar(data['Ra'].values, data['Nu_top'].values , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(data['Ra'].values, data['Nu_bottom'].values , yerr=data['Ecart type bot'].values, linestyle="None")
def func(x,a):
return 10**(np.log10(x)/a)
"""maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]"""
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(data['Nu_top'].head(10).values, params, x=data['Ra'].head(10).values)
print(ret.fit_report())
popt, pcov = curve_fit(func, data['Ra'].head(10).values,
data['Nu_top'].head(10).values, sigma=data['Ecart type top'].head(10).values,
absolute_sigma=True, p0=[0.25])
plt.plot(data['Ra'].head(10).values, func(data['Ra'].head(10).values, *popt),
'r-', label='fit: a=%5.3f' % tuple(popt))
popt, pcov = curve_fit(func, data['Ra'].tail(4).values, data['Nu_top'].tail(4).values,
sigma=data['Ecart type top'].tail(4).values,
absolute_sigma=True, p0=[0.33])
plt.plot(data['Ra'].tail(4).values, func(data['Ra'].tail(4).values, *popt),
'b-', label='fit: a=%5.3f' % tuple(popt))
print(pcov)
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()
所以我使用日志：logRa = a * logNu。
Ra = x 轴
Nu = y 轴
这就是我以这种方式定义函数 func 的原因。
如您所见，我的两次合身并不完全正确。我的协方差等于 [0.00010971]。所以我不得不做错事，但我没有看到。我需要帮助。
这里的数据文件：
data.xlsx
我注意到 Ra 的数据值很大，在对它们进行缩放后，我执行了方程式搜索 - 这是我的代码结果。我使用标准 scipy 遗传算法模块 differential_evolution 来确定 curve_fit() 的初始参数值，并且该模块使用拉丁超立方体算法来确保彻底搜索参数 space 这需要搜索范围。给出初始参数估计的范围比找到特定值要容易得多。此等式适用于 nu_top 和 nu_bottom，请注意，图中未按对数比例缩放，因为在本例中没有必要。
import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import pandas
import warnings
filename = 'data.xlsx'
data=pandas.read_excel(filename,sheet_name='Sheet2',index=False,dtype={'Ra': float})
# notice the Ra scaling by 10000.0
xData = data['Ra'].values / 10000.0
yData = data['Nu_bottom']
def func(x, a, b, c): # "Combined Power And Exponential" from zunzun.com
return a * numpy.power(x, b) * numpy.exp(c * x)
# 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([0.0, 10.0]) # search bounds for a
parameterBounds.append([0.0, 10.0]) # search bounds for b
parameterBounds.append([0.0, 10.0]) # search bounds for c
# "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)
这里我把我的数据x和y放在log10()中。该图采用对数刻度。所以通常我应该有两个仿射函数，系数分别为 0.25 和 0.33。我改变了你的程序 James 中的函数 func 和 b 和 c 的界限，但我没有得到好的结果。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
def func(x,a):
return a*x
maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
popt, pcov = curve_fit(func, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), func(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
print(pcov)
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()
使用 polyfit 我得到了更好的结果。 使用我的代码，我打开文件并计算 log (Ra) 和 log (Nu)，然后以对数刻度绘制 (log (Ra)、log (Nu))。 对于 Ra <1e6，我应该有 a = 0.25，如果不是 a = 0.33
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from numpy import polyfit
import numpy.polynomial.polynomial as poly
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
x=np.log10(data['Ra'].values)
y1=np.log10(data['Nu_top'].values)
y2=np.log10(data['Nu_bottom'].values)
x2=np.log10(data['Ra'].head(11).values)
y4=np.log10(data['Nu_top'].head(11).values)
x3=np.log10(data['Ra'].tail(4).values)
y5=np.log10(data['Nu_top'].tail(4).values)
plt.xscale('log')
plt.yscale('log')
plt.scatter(x, y1, label='Nu_top')
plt.scatter(x, y2, label='Nu_bottom')
plt.errorbar(x, y1 , yerr=data['Ecart type top'].values, linestyle="None")
plt.errorbar(x, y2 , yerr=data['Ecart type bot'].values, linestyle="None")
"""a=np.ones(10, dtype=np.float)
weights = np.insert(a,0,1E10)"""
coefs = poly.polyfit(x2, y4, 1)
print(coefs)
ffit = poly.polyval(x2, coefs)
plt.plot(x2, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
absError = ffit - x2
SE = np.square(absError) # squared errors
MSE = np.mean(SE) # mean squared errors
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (np.var(absError) / np.var(x2))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
print('Predicted value at x=0:', ffit[0])
print()
coefs = poly.polyfit(x3, y5, 1)
ffit = poly.polyval(x3, coefs)
plt.plot(x3, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()
我的问题已解决，我设法用或多或少正确的结果拟合我的曲线