Как подогнать гауссиану с помощью Astropy

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

Если я подгоняю гауссиан, мне нравится давать исходной модели некоторые начальные параметры, основанные на вычислительном «зрении» их так (здесь я назвал поток ваших реальных данных и длину волны как orig_flux и orig_wavelength соответственно):

>>> an_amplitude = orig_flux.min()
>>> an_mean = orig_wavelength[orig_flux.argmin()]
>>> an_stddev = np.sqrt(np.sum((orig_wavelength - an_mean)**2) / (len(orig_wavelength) - 1))
>>> print(f'mean: {an_mean}, stddev: {an_stddev}, amplitude: {an_amplitude}')
mean: 5737.979797979798, stddev: 42.768052162734605, amplitude: 84.73925092448636

где для стандартного отклонения я использовал непредвзятую оценку стандартного отклонения.

Нанесение этого на мои поддельные данные показывает, что это разумные значения, которые я мог бы выбрать, если бы я также вручную просматривал данные:

>>> plt.plot(orig_wavelength, orig_flux, '.k', zorder=1)
>>> plt.scatter(an_mean, an_amplitude, color='red', s=100, zorder=2)
>>> plt.vlines([an_mean - an_stddev, an_mean + an_stddev], orig_flux.min(), orig_flux.max(),
...            linestyles='dashed', colors='gg', zorder=2)

Одна функция, которую я хотел добавить в astropy.modeling в прошлом, — это необязательные методы, которые можно прикрепить к некоторым моделям, чтобы дать разумные оценки их параметров на основе некоторых данных. Таким образом, для гауссов такой метод вернет почти то же самое, что я только что вычислил выше. Я не знаю, было ли это когда-либо реализовано.

Также стоит отметить, что ваша гауссова функция будет инвертирована (с отрицательной амплитудой) и что она смещена по оси потока примерно на 120 точек, поэтому я добавил Const1D к моей модели, чтобы учесть это, и вычесть смещение из амплитуды:

>>> an_disp = orig_flux.max()
>>> g_init = (
...     models.Const1D(an_disp) +
...     models.Gaussian1D(amplitude=(an_amplitude - an_disp), mean=an_mean, stddev=an_stddev)
... )
>>> fit_g = fitting.LevMarLSQFitter()
>>> g = fit_g(g_init, orig_wavelength, orig_flux)

В результате получается следующая подгонка, которая уже выглядит намного лучше:

>>> plt.plot(orig_wavelength, orig_flux, '.k')
>>> plt.plot(orig_wavelength, g(orig_wavelength), 'r-')

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