Как подогнать гауссову с помощью Astropy - PullRequest
Новая
       67

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

3 голосов
/ 18 апреля 2020

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

https://i.stack.imgur.com/H37VC.png

https://i.stack.imgur.com/DOKSd.png

Вот мой код: 

%pylab inline
from astropy.modeling import models,fitting
from astropy import modeling

#Fitting a gaussian for the absorption lines
wavelength= linspace(galaxy1_wavelength_extracted_1.min(),galaxy1_wavelength_extracted_1.max(),200)
g_init = models.Gaussian1D(amplitude=1., mean=5000, stddev=1.)
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, galaxy1_wavelength_extracted_1, galaxy1_flux_extracted_1)

#Plotting 
plot(galaxy1_wavelength_extracted_1,galaxy1_flux_extracted_1,".k")
plot(wavelength, g(wavelength))
xlabel("Wavelength ($\\AA$)")
ylabel("Flux (counts)")

Что я делаю неправильно или отсутствует?

1 Ответ

2 голосов
/ 19 апреля 2020

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

Если я подгоняю гауссиан, я хотел бы дать исходной модели некоторые начальные параметры, основанные на вычислительной обработке их "глазного яблока", вот так (здесь я назвал поток и длину ваших реальных данных как 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)

enter image description here

Одна функция, которую я хотел добавить к 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-')

enter image description here

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

...