如何适合功能的外壳 [英] How to Fit to The Outer Shell of a Function

问题描述

我试图使高斯函数适合杂乱的函数.我只想适合外部外壳(这些不只是每个x的最大值，因为某些最大值也太小，因为样本量很小).

I am trying to make a gaussian fit on a function that is messy. I want to only fit the exterior outer shell (these are not just the max values at each x, because some of the max values will be too low too, because the sample size is low).

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')

推荐答案

这是一种类似于此帖子的迭代方法. .在图形的形状不允许使用凸包的意义上有所不同.因此，这个想法是创建一个成本函数，该函数试图最小化图形的面积，同时如果点在图形上方，则付出高昂的代价.根据OP中图形的类型，需要调整成本函数.还必须检查最终结果中所有点是否真的在图形下方.在这里，您可以摆弄成本函数的细节.例如，例如，我可以在tanh中包含一个与tanh( slope * ( x - offset) )类似的偏移量，以将解决方案推离数据更远.

This is an iterative approach similar to this post. It is different in the sense that the shape of the graph does not permit the use of convex hull. So the idea is to create a cost function that tries to minimize the area of the graph while paying high cost if a point is above the graph. Depending on the type of the graph in OP the cost function needs to be adapted. One also has to check if in the final result all points are really below the graph. Here one can fiddle with details of the cost function. One my, e.g., include an offset in the tanh like tanh( slope * ( x - offset) ) to push the solution farther away from the data.

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]

和

这篇关于如何适合功能的外壳的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！