Функция подгонки с curve_fit, но подгоночная кривая неверна

Question

def gaus(x,a,x0,sigma):
     return a*np.exp(-(x-x0)**2/(2*sigma**2))

times, amplitudes = openFile("../datafiles/number_of_counts.txt")

mean = sum(np.array(times)*np.array(amplitudes))/sum(amplitudes)            
sigma = np.sqrt(sum(np.array(amplitudes)*(np.array(times)-mean)**2)/sum(amplitudes))        

params,pcov = curve_fit(gaus,times, amplitudes,p0=[max(amplitudes),mean,sigma])


plt.plot(times, amplitudes)
plt.plot(times ,gaus(np.array(times),params[0],params[1],params[2]),'r',  label="fitted curve")

plt.ylabel("Coincidents")
plt.title("Coincident plot")
plt.legend()
plt.show()

Моя подгонка по Гауссу не работает должным образом, но выглядит как мягкая кривая, вместо того, чтобы подгоняться к острому пику, я предполагаю, что у меня есть какая-то очень глупая ошибка в моем сценарий, но не уверен, что. Кто-то, кто может это увидеть?

Answer 1

Ваши данные имеют постоянное смещение около 3750, но ваша функция модели gaus не может учесть это, поэтому вы подбираете нормальное распределение со смещением 0.

Требуется еще один параметр:

def gaus(x, a, x0, sigma, c):
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2)) + c

Тогда:

offset_guess = 3750  # maybe calculate it from the data as well
params, pcov = curve_fit(
    gaus, times, amplitudes,
    p0=[max(amplitudes), mean, sigma, offset_guess])

plt.plot(times, gaus(np.array(times), params[0], params[1], params[2], params[3]), ...)

Результат:

>>> print(params)
[1545.00193331  -20.45639132  -43.28484495 3792.41050636]

Answer 2

Я извлек данные из графика для анализа и обнаружил, что с извлеченными данными пик Вейбулла плюс смещение дал мне лучшее соответствие, чем гауссов пик с смещением. Вот графический сборщик Python с извлеченными данными и уравнением пика Вейбулла плюс смещение, вы сможете заменить действительные данные и запустить их напрямую. Обратите внимание на простое определение начальных оценок параметров.

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

xData = numpy.array([4.914e+03, 5.600e+03, 7.886e+03, 8.571e+03, 1.063e+04, 1.154e+04, 1.269e+04, 1.566e+04, 1.634e+04, 1.817e+04, 1.886e+04, 2.114e+04, 2.389e+04, 2.526e+04, 2.754e+04, 3.051e+04, 3.257e+04, 3.417e+04, 3.554e+04, 3.669e+04, 4.011e+04, 4.240e+04, 4.491e+04, 4.583e+04, 4.697e+04, 4.880e+04, 4.994e+04, 5.154e+04, 5.246e+04, 5.474e+04, 5.634e+04, 5.886e+04, 6.091e+04, 6.366e+04, 6.731e+04, 7.051e+04, 7.257e+04, 7.394e+04, 7.691e+04, 7.851e+04, 7.966e+04, 8.103e+04, 8.240e+04, 8.514e+04, 8.720e+04, 8.834e+04, 8.949e+04, 9.109e+04, 9.223e+04, 9.360e+04, 9.566e+04, 9.726e+04, 9.909e+04, 1.005e+05, 1.014e+05, 1.030e+05, 1.059e+05, 1.073e+05, 1.089e+05, 1.101e+05, 1.119e+05, 1.130e+05, 1.139e+05, 1.162e+05, 1.178e+05, 1.190e+05, 1.203e+05, 1.222e+05, 1.233e+05, 1.247e+05, 1.281e+05, 1.299e+05, 1.309e+05, 1.329e+05, 1.341e+05, 1.357e+05, 1.370e+05, 1.382e+05, 1.395e+05, 1.407e+05, 1.430e+05, 1.439e+05])
yData = numpy.array([3.300e+03, 8.100e+03, 6.100e+03, 1.010e+04, 9.700e+03, 7.300e+03, 7.500e+03, 6.900e+03, 8.100e+03, 3.900e+03, 5.700e+03, 4.900e+03, 4.500e+03, 8.300e+03, 4.100e+03, 8.100e+03, 5.300e+03, 8.100e+03, 6.700e+03, 1.130e+04, 9.300e+03, 6.300e+03, 9.500e+03, 8.900e+03, 1.490e+04, 6.300e+03, 1.190e+04, 6.300e+03, 7.700e+03, 1.310e+04, 9.500e+03, 1.590e+04, 1.050e+04, 1.930e+04, 4.890e+04, 7.350e+04, 5.230e+04, 5.130e+04, 2.350e+04, 1.950e+04, 1.010e+04, 1.510e+04, 9.500e+03, 9.500e+03, 6.900e+03, 6.900e+03, 1.210e+04, 6.300e+03, 7.700e+03, 5.700e+03, 1.410e+04, 8.700e+03, 1.390e+04, 4.900e+03, 7.500e+03, 4.900e+03, 9.500e+03, 5.300e+03, 9.300e+03, 6.300e+03, 1.250e+04, 4.300e+03, 7.700e+03, 6.900e+03, 9.700e+03, 8.500e+03, 1.130e+04, 5.300e+03, 5.100e+03, 1.700e+03, 8.700e+03, 7.300e+03, 6.300e+03, 2.100e+03, 3.100e+03, 7.100e+03, 4.900e+03, 6.100e+03, 3.700e+03, 9.300e+03, 5.500e+03, 5.700e+03])


def func(x, a, b, c, offset): # Weibull peak with offset from zunzun.com
    return a * numpy.exp(-0.5 * numpy.power(numpy.log(x/b) / c, 2.0)) + offset

a_est = max(yData)
b_est = max(yData)
c_est = 1.0
offset_est = min(yData)
initialParameters = numpy.array([a_est, b_est, c_est, offset_est])

# curve fit the test data
fittedParameters, pcov = curve_fit(func, xData, yData, initialParameters)

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print('Parameters:', fittedParameters)
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData), 500)
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)