python - 如何适应函数的外壳

Question

我正在尝试对一个困惑的函数进行高斯拟合。我只想适合外部外壳(这些不仅仅是每个 x 处的最大值，因为一些最大值也会太低，因为样本量很小)。

from scipy.optimize import curve_fit
def Gauss(x, a, x0, sigma, offset):
        return a * np.exp(-np.power(x - x0,2) / (2 * np.power(sigma,2))) + offset

def fitNormal(x, y):
    popt, pcov = curve_fit(Gauss, x, y, p0=[np.max(y), np.median(x), np.std(x), np.min(y)])
    return popt

plt.plot(xPlot,yPlot, 'k.')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Y(x)')

x,y = xPlot,yPlot
popt = fitNormal(x, y)
minx, maxx = np.min(x), np.max(x)
xFit = np.arange(start=minx, stop=maxx, step=(maxx-minx)/1000)
yFitTest = Gauss(xPlot, popt[0], popt[1], popt[2], popt[3])

print('max fit test: ',np.max(yFitTest))
print('max y: ',np.max(yPlot))

maxIndex = np.where(yPlot==np.max(yPlot))[0][0]
factor = yPlot[maxIndex]/yFitTest[maxIndex]
yFit = Gauss(xPlot, popt[0], popt[1], popt[2], popt[3]) * factor

plt.plot(xFit,yFit,'r')

Answer 1

这是一种类似于 this post 的迭代方法。它的不同之处在于图的形状不允许使用凸包。因此，我们的想法是创建一个成本函数，尝试最小化图形的面积，同时如果某个点位于图形上方则付出高昂的成本。根据 OP 中图形的类型，需要调整成本函数。还必须检查最终结果中所有点是否确实位于图表下方。在这里，人们可以修改成本函数的细节。例如，我可以在 tanh 中包含一个偏移量，例如 tanh(slope * ( x - offset) ) 以使解远离数据。

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import leastsq

def g( x, a, s ): 
    return a * np.exp(-x**2 / s**2 )

def cost_function( params, xData, yData, slope, val ):
    a,s = params
    area = 0.5 * np.sqrt( np.pi ) * a * s
    diff = np.fromiter ( ( y - g( x, a, s) for x, y in zip( xData, yData ) ), np.float )
    cDiff = np.fromiter( ( val * ( 1 + np.tanh( slope * d ) ) for d in diff ), np.float )
    out = np.concatenate( [ [area] , cDiff ] )
    return out

xData = np.linspace( -5, 5, 500 )
yData = np.fromiter( (  g( x, .77, 2 ) * np.sin( 257.7 * x )**2 for x in xData ), np.float )


sol=[ [ 1, 2.2 ] ]
for i in range( 1, 6 ):
    solN, err = leastsq( cost_function, sol[-1] , args=( xData, yData, 10**i, 1 ) )
    sol += [ solN ]
    print sol

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1)
ax.scatter( xData, yData, s=1 ) 
for solN in sol:
    solY = np.fromiter( (  g( x, *solN ) for x in xData ), np.float )
    ax.plot( xData, solY ) 
plt.show()

给予

>> [0.8627445  3.55774814]
>> [0.77758636 2.52613376]
>> [0.76712184 2.1181137 ]
>> [0.76874125 2.01910211]
>> [0.7695663  2.00262339]

和