Как соответствовать внешней оболочке функции

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

Это итеративный подход, подобный этому посту .Это отличается в том смысле, что форма графика не позволяет использовать выпуклую оболочку.Таким образом, идея состоит в том, чтобы создать функцию стоимости, которая пытается минимизировать площадь графика, оплачивая при этом высокую стоимость, если точка находится над графиком.В зависимости от типа графика в 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 Differental Evolution в сочетании с «кирпичной стеной», где, если любое прогнозируемое значение во время подгонки больше соответствующего значения 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)