如何适应函数的外壳

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]

和

Answer 2

这是一种使用 scipy 的差分进化模块结合“砖墙”的不同方法，如果拟合期间的任何预测值大于相应的 Y 值，则拟合误差会非常大。我从@mikuszefski 的答案中无耻地挖了代码来生成本示例中使用的数据。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings

from scipy.optimize import differential_evolution

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

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 )


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


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = Gauss(xData, *parameterTuple)
    multiplier = 1.0
    for i in range(len(val)):
        if val[i] < yData[i]: # ****** brick wall ******
            multiplier = 1.0E10
    return np.sum((multiplier * (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)

    minData = min(minX, minY)
    maxData = max(maxX, maxY)

    parameterBounds = []
    parameterBounds.append([minData, maxData]) # parameter bounds for a
    parameterBounds.append([minData, maxData]) # parameter bounds for x0
    parameterBounds.append([minData, maxData]) # parameter bounds for sigma
    parameterBounds.append([minData, maxData]) # parameter bounds for offset

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

# generate initial parameter values
geneticParameters = generate_Initial_Parameters()

# create values for display of fitted function
y_fit = Gauss(xData, *geneticParameters)

plt.scatter(xData, yData, s=1 ) # plot the raw data
plt.plot(xData, y_fit) # plot the equation using the fitted parameters
plt.show()

print('parameters:', geneticParameters)