将模型函数拟合到定义的数据范围

Question

假设我有一组(x=times,y=observation)有多个时间间隔的数据。无论数据趋势是什么，我们都假设它是线性的。在时间间隔内，有一个衰减，使数据偏离纯线性趋势，直到再次开始观察并恢复线性趋势。时间上有多个间隔，但在这个例子中，我只报告了最短的快照来说明问题。时间间隔是（正）线性趋势之间没有可用观察值的时间，因此连续之间的差异x=times（远）大于平均值。我想将衰减建模为函数 ( y_decay = C -D*x)的一部分

from scipy.optimize import curve_fit

import matplotlib.pyplot as plt

def f(x, A, B, C, D):

    line = A*x + B if ((x>=1) & (x<=3) | (x>=5) & (x<=9) | (x>=23) & (x<=25)) else C-D*x

    return line

x=[1,2,3, 12,13,14, 23,24,25]

y=[2,4,6, 5, 7, 9, 8, 10,12]

popt, pcov = curve_fit(f, x, y) 

figure = plt.figure(figsize=(5.15, 5.15))

figure.clf()

plot = plt.subplot(111)

ax1 = plt.gca()

plot.scatter(x,y)

plt.show()

如何将decay变量建模为函数的一部分并获得其最佳拟合值？

Answer 1

在完全周期性的情况下，我会做这样的事情：import matplotlib.pyplot as pltimport numpy as npfrom scipy.optimize import leastsqdef data_func( x, x0, a, bc, off, p1, p2):&nbsp; &nbsp; """&nbsp; &nbsp; fit function that uses modulus to get periodicity&nbsp; &nbsp; two linear functions are combines piecewise by boxing them&nbsp; &nbsp; using the heaviside function with the periodic input&nbsp; &nbsp; over all slope is added.&nbsp; &nbsp; Note that the "decay" part maybe positive with this solution.&nbsp; &nbsp; """&nbsp; &nbsp; P1 = abs(p1)&nbsp; &nbsp; P2 = abs(p2)&nbsp; &nbsp; X = x - x0&nbsp; &nbsp; P= P1 + P2&nbsp; &nbsp; mod = X % P&nbsp; &nbsp; y0 = a * P1&nbsp; &nbsp; beta = y0 * P / P2&nbsp; &nbsp; slope = y0 / P2&nbsp; &nbsp; box1 =&nbsp; np.heaviside( +np.abs( ( X - P1 / 2. ) % P - 0.5 * P ) - 0.5 * P2, .5 )&nbsp; &nbsp; box2 =&nbsp; np.heaviside( -np.abs( ( X - P1 / 2. ) % P - 0.5 * P ) + 0.5 * P2, .5 )&nbsp; &nbsp; out = a * mod * box1&nbsp;&nbsp; &nbsp; out += (beta - slope * mod&nbsp; )* box2&nbsp; &nbsp; out += off + bc * X&nbsp; &nbsp; return outdef residuals( params, xl ,yl ):&nbsp; &nbsp; x0, a, bc, off, p1, p2 = params&nbsp; &nbsp; diff = np.fromiter( ( y - data_func( x, x0, a, bc, off, p1, p2 ) for x, y in zip( xl, yl )&nbsp; ), np.float )&nbsp; &nbsp; return difftheOff = 0.7theP1= 1.8869theP2 = 5.21163theP = theP1 + theP2xdata = np.linspace(-1, 26, 51 )xdata = np.fromiter( ( x for x in xdata if (x-theOff)%theP <= theP1 ),np.float )ydata = np.fromiter( ( data_func( x, theOff, .6, .1, 17, theP1, theP2) for x in xdata ),np.float )tl = np.linspace(-1, 26, 150 )yl = np.fromiter( ( data_func( x, theOff, .6, .1, 17, theP1, theP2) for x in tl ),np.float )guess= [0, 0.55, .1, 16 , 2, 5 ]sol, err = leastsq( residuals, guess, args = ( xdata, ydata ) )print sol### getting the real slopes out of the datas1 = sol[1]+ sol[2]&nbsp;s2 =&nbsp; - sol[1] * sol[4] / sol[5] + sol[2]print "real slope1 = {}".format( s1 )print "real slope2 = {}".format( s2 )fit = np.fromiter( ( data_func( x, *sol ) for x in tl ),np.float )fig = plt.figure()ax = fig.add_subplot( 1, 1, 1 )### original dataax.plot( tl, yl, ls='--')ax.plot( xdata, ydata, ls='', marker='+')### fitax.plot( tl, fit )### check the slopesshort = np.linspace(0, 3, 3)ax.plot( short, [ 17 + s1 * s for s in short ] )short = np.linspace(3, 10, 3)ax.plot( short, [ 18 + s2 * s for s in short ] )ax.grid()plt.show()这使：>> [ 0.39352332&nbsp; 0.59149625&nbsp; 0.10850375 16.78546632&nbsp; 1.85009228&nbsp; 5.35049099]>> real slope1 = 0.7>> real slope2 = -0.0960237685357和自然地，间隙中缺乏信息会导致那里的斜率相当糟糕。结果，在周期性中存在相应的误差。如果知道这一点，那么准确度当然会提高。您需要对起始参数进行合理的猜测！

Answer 2

当假设所有数据都具有相同的斜率m并且所有“衰减”斜率时，那将是我的 AnsatzDimport matplotlib.pyplot as plt ### only for plotting; not important for the OPimport numpy as np ### for easy data manipulationfrom scipy.optimize import leastsq ### one of many possibilities for optimizationdef generic_data( m, D, n ): ### just generating data; not important for OP&nbsp; &nbsp; alpha0 = 0&nbsp; &nbsp; timel = [ 0 ] ### to avoid errer, remove at the end&nbsp; &nbsp; dl = list()&nbsp; &nbsp; for gaps in range( n + 1 ):&nbsp;&nbsp; &nbsp; &nbsp; &nbsp; for pnts in range( 3 + np.random.randint( 2 ) ):&nbsp;&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; timel.append ( timel[-1] +&nbsp; 0.5 + np.random.rand() )&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; dl.append( m * timel[ -1 ] + alpha0 )&nbsp; &nbsp; &nbsp; &nbsp; ###now the gap&nbsp; &nbsp; &nbsp; &nbsp; dt =&nbsp; 2. + 2 * np.random.rand()&nbsp; &nbsp; &nbsp; &nbsp; timel.append ( timel[-1] + dt )&nbsp; &nbsp; &nbsp; &nbsp; dl.append( dl[-1] - D * dt )&nbsp; &nbsp; &nbsp; &nbsp; alpha0 = dl[-1] - m * timel[-1]&nbsp; &nbsp; del timel[0]&nbsp; &nbsp; ### remove jump of last gap&nbsp; &nbsp; del timel[-1]&nbsp; &nbsp; del dl[-1]&nbsp; &nbsp; dl = np.fromiter( ( y + np.random.normal( scale=0.1 ) for y in dl ), np.float )&nbsp; &nbsp; return np.array( timel ), dldef split_into_blocks( tl, dl ):&nbsp; &nbsp; """&nbsp; &nbsp; identifying the data blocks by detecting positions of negative slope.&nbsp; &nbsp; first subtract a shifted version of the data from itself&nbsp; &nbsp; find the negatives and get their positions&nbsp; &nbsp; get sub-lists based on this positins&nbsp;&nbsp; &nbsp; """&nbsp; &nbsp; mask = np.where( dl[1::] - dl[:-1:] < 0, 1, 0 )&nbsp; &nbsp; where = np.argwhere( mask )&nbsp; &nbsp; where = where.reshape( 1, len( where ) )[0]&nbsp; &nbsp; where = np.append( where, len( dl ) - 1 )&nbsp; &nbsp; where = np.insert( where, 0, -1 )&nbsp; &nbsp; tll = list()&nbsp; &nbsp; dll = list()&nbsp; &nbsp; for s, e in zip( where[ :-1:], where[ 1:: ] ):&nbsp; &nbsp; &nbsp; &nbsp; dll.append( dl[ s + 1 : e + 1 ] )&nbsp; &nbsp; &nbsp; &nbsp; tll.append( tl[ s + 1 : e + 1 ] )&nbsp; &nbsp; return np.array( tll ), np.array( dll )def residuals( params, tblocks, dblocks ):&nbsp; &nbsp; """&nbsp; &nbsp; typical residual function to be called by leastsq&nbsp; &nbsp; """&nbsp; &nbsp; ### split parameters&nbsp; &nbsp; nob = len( tblocks )&nbsp; &nbsp; m = params[0] ### all data with same slope&nbsp; &nbsp; D = params[1] ### all decay with same slope&nbsp; &nbsp; alphal = params[2:2+nob] ### off sets differ&nbsp; &nbsp; betal = params[-nob+1::]&nbsp; &nbsp; out = list()&nbsp; &nbsp; ### evaluate diefference between data and test function&nbsp; &nbsp; for i, (tl, dl) in enumerate( zip(tblocks, dblocks ) ):&nbsp; &nbsp; &nbsp; &nbsp; diff = [ d - ( m * t + alphal[i] ) for t, d in zip( tl, dl ) ]&nbsp; &nbsp; &nbsp; &nbsp; out= out + diff&nbsp; &nbsp; ### evaluate differences for the gapfunction; this could be done&nbsp; &nbsp; ### completely independent, but I do it here to have all in one shot.&nbsp; &nbsp; for j in range( nob -1 ):&nbsp; &nbsp; &nbsp; &nbsp; out.append( dblocks[ j][-1] - ( betal[j] + D * tblocks[j][-1] ) ) ###left point gap&nbsp; &nbsp; &nbsp; &nbsp; out.append( dblocks[ j+1][0] - ( betal[j] + D * tblocks[j+1][0] ) ) ###right point gap&nbsp; &nbsp; return out### create generic datatl, dl =&nbsp; generic_data( 1.3, .3, 3 )tll, dll = split_into_blocks( tl, dl )### and fitnob = len(dll)m0 = +1.0D0 = -0.1guess = [m0, D0 ]+ nob * [-3] + ( nob - 1 ) * [ +4 ]sol, err = leastsq( residuals, x0=guess, args=( tll, dll ) )mf = sol[0]Df = sol[1]print mf, Dfalphal = sol[2:2+nob]betal = sol[-nob+1::]### ... and plotfig = plt.figure()ax = fig.add_subplot( 1, 1, 1 )###original dataax.plot( tl, dl, ls='', marker='o')###identify blocksfor a,b in zip( tll, dll ):&nbsp; &nbsp; ax.plot( a, b, ls='', marker='x')###fit resultsfor j in range(nob):&nbsp; &nbsp; tloc = np.linspace( tll[j][0] - 0.3, tll[j][-1] + 0.3 , 3 )&nbsp; &nbsp; ax.plot( tloc, [ mf * t + alphal[j] for t in tloc ] )for j in range(nob - 1):&nbsp; &nbsp; tloc = np.linspace( tll[j][-1] - 0.3, tll[j+1][0] + 0.3 , 3 )&nbsp; &nbsp; ax.plot( tloc, [ Df * t +betal[j] for t in tloc ] )plt.show()这导致>> 1.2864142170851447 -0.2818180721270913和然而，该模型可能要求衰减线不与数据范围内的数据线相交。这使得额外的摆弄成为必要，因为我们需要检查某种类型的边界。另一方面，可以只拟合数据并寻找具有满足前面提到的边界的最小可能斜率的衰减曲线。因此，在这种情况下，我会D从残差中删除拟合部分，然后再计算它。