error-handling - 使用scipy.odr拟合曲线

Question

我试图通过依赖于两个变量的拟合函数拟合一组数据点，我们将它们称为xdata和sdata。问题是我的曲线相当平坦，我希望它或多或少地“遵循点”。

我尝试使用scipy.odr拟合曲线，除了曲线太平坦以外，它工作得很好:

import numpy as np
from math import pi
from math import sqrt
from math import log
from scipy import optimize
import scipy.optimize
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.odr import *

mudr=np.array([ 57.43708609,  46.26119205,  55.60688742,  33.21615894,
        28.27072848,  22.54649007,  21.80662252,  11.21483444,   5.80211921]) 
#xdata points
dme=array([ 128662.54890776,  105265.32915726,  128652.56835434,
         77968.67019573,   66273.56542068,   58464.58559543,
         54570.66624991,   27286.90038703,   19480.92689266]) #xdata error
dmss22=np.array([  4.90050000e+17,   4.90050000e+17,   4.90050000e+17,
         4.90050000e+17,   4.90050000e+17,   4.90050000e+17,
         4.90050000e+17,   4.90050000e+17,   4.90050000e+17]) #sdata points
dmse=np.array([  1.09777592e+21,   1.11512117e+21,   1.13381702e+21,
         1.15033267e+21,   1.14883089e+21,   1.27076265e+21,
         1.22637165e+21,   1.19237598e+21,   1.64539205e+21]) # sdata error
F=np.array([ 115.01944248,  110.24354867,  112.77812389,  104.81830088,
        104.35746903,  101.32016814,  100.54513274,   96.94226549,
         93.00424779]) #ydata points
dF=np.array([  72710.75386699,   72590.6256987 ,  176539.40403673,
        130555.27503081,  124299.52080164,  176426.64340597,
        143013.52848306,  122117.93022746,  157547.78395513])#ydata error

def Ffitsso(p,X,B=2.58,Fc=92.2,mu=770,Za=0.9468): #fitfunction
    temp1 = (2*B*X[0])/(4*pi*Fc)**2
    temp2 = temp1*(afij[0]+afij[1]*np.log((2*B*X[0])/mu**2))
    temp3 = temp1**2*(afij[2]+afij[3]*np.log((2*B*X[0])/mu**2)+\
                   afij[4]*(np.log((2*B*X[0])/mu**2))**2)
    temp4 = temp1**3*(afij[5]+afij[6]*np.log((2*B*X[0])/mu**2)+\
                      afij[7]*(np.log((2*B*X[0])/mu**2))**2+\
                      afij[8]*(np.log((2*B*X[0])/mu**2))**3)
    return Fc/Za*(1+p[0]*X[1])*(1+temp2+temp3+temp4)+p[1]

#fitting using scipy.odr
xtot=np.row_stack( (mudr, dmss22) )
etot=np.row_stack( (Ze, dmss22e) )
fitting = Model(Ffitsso)
mydata = RealData(xtot, F, sx=etot2, sy=dF)
myodr = ODR(mydata, fitting, beta0=[0, 100])
myoutput = myodr.run()
myoutput.pprint()
bet=myoutput.beta

plt.plot(mudr,F,"b^")
plt.plot(mudr,Ffitsso(bet,[mudr,dmss22]))

与1相比，fit函数中的p [0] * X [0]应该较小，但通过拟合，p [0]的值按e-18的顺序排列，而dmss22的值按e-17的顺序排列。还不够小。

更糟糕的是，它是负数，意味着函数会减少，这是不应该发生的，它应该像绘制的数据点一样增加。

编辑:我已修复，不知道它对初始beta值如此敏感，将beta [0] = 1.5 * 10 (-15)设置为有效!**

Answer 1

这是一个具有曲线拟合和ODR拟合的图形拟合器，它使用scipy的差分进化(DE)遗传算法来为非线性求解器提供初始参数估计。 DE的scipy实现使用Latin Hypercube算法来确保彻底搜索参数空间，这需要在其中搜索参数范围-在本示例中，这些范围取自数据的最大值和最小值。请注意，为初始参数估计指定界限要比为单个特定值容易得多。

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.odr
from scipy.optimize import differential_evolution
import warnings

xData = numpy.array([1.1, 2.2, 3.3, 4.4, 5.0, 6.6, 7.7, 0.0])
yData = numpy.array([1.1, 20.2, 30.3, 40.4, 50.0, 60.6, 70.7, 0.1])


def func(x, a, b, c, d, offset): # curve fitting function for curve_fit()
    return a*numpy.exp(-(x-b)**2/(2*c**2)+d) + offset


def func_wrapper_for_ODR(parameters, x): # parameter order for ODR
    return func(x, *parameters)


# 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([minX, maxX]) # search bounds for b
    parameterBounds.append([minX, maxX]) # search bounds for c
    parameterBounds.append([minY, maxY]) # search bounds for d
    parameterBounds.append([0.0, maxY]) # search bounds for Offset

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

geneticParameters = generate_Initial_Parameters()


##########################
# curve_fit section
##########################

fittedParameters_curvefit, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters curve_fit:', fittedParameters_curvefit)
print()

modelPredictions_curvefit = func(xData, *fittedParameters_curvefit) 

absError_curvefit = modelPredictions_curvefit - yData

SE_curvefit = numpy.square(absError_curvefit) # squared errors
MSE_curvefit = numpy.mean(SE_curvefit) # mean squared errors
RMSE_curvefit = numpy.sqrt(MSE_curvefit) # Root Mean Squared Error, RMSE
Rsquared_curvefit = 1.0 - (numpy.var(absError_curvefit) / numpy.var(yData))

print()
print('RMSE curve_fit:', RMSE_curvefit)
print('R-squared curve_fit:', Rsquared_curvefit)

print()


##########################
# ODR section
##########################
data = scipy.odr.odrpack.Data(xData,yData)
model = scipy.odr.odrpack.Model(func_wrapper_for_ODR)
odr = scipy.odr.odrpack.ODR(data, model, beta0=geneticParameters)

# Run the regression.
odr_out = odr.run()

print('Fitted parameters ODR:', odr_out.beta)
print()

modelPredictions_odr = func(xData, *odr_out.beta) 

absError_odr = modelPredictions_odr - yData

SE_odr = numpy.square(absError_odr) # squared errors
MSE_odr = numpy.mean(SE_odr) # mean squared errors
RMSE_odr = numpy.sqrt(MSE_odr) # Root Mean Squared Error, RMSE
Rsquared_odr = 1.0 - (numpy.var(absError_odr) / numpy.var(yData))

print()
print('RMSE ODR:', RMSE_odr)
print('R-squared ODR:', Rsquared_odr)

print()


##########################################################
# graphics output section
def ModelsAndScatterPlot(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 plots
    xModel = numpy.linspace(min(xData), max(xData))
    yModel_curvefit = func(xModel, *fittedParameters_curvefit)
    yModel_odr = func(xModel, *odr_out.beta)

    # now the models as line plots
    axes.plot(xModel, yModel_curvefit)
    axes.plot(xModel, yModel_odr)

    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
ModelsAndScatterPlot(graphWidth, graphHeight)