如何使用熵拟合高斯

Question

我正在尝试使用astropy.modeling软件包使高斯适合一组数据点，但我得到的只是一条平线。见下文：

这是我的代码：

%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)")

我在做什么错还是想念？

Answer 1

我制作了一些类似于您的伪造数据，并尝试在其上运行代码并获得了相似的结果。我认为问题在于，如果您不将模型的初始参数调整为至少类似于原始模型，否则无论拟合执行了多轮拟合，拟合器都将无法收敛。

[如果我要拟合高斯，我想给初始模型一些基于计算的初始参数，像这样（在这里，我将实际数据的通量和波长分别命名为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
标准偏差我用的是unbiased standard deviation estimate。

在我的虚假数据上显示此值表明，如果我也手动查看数据，这些值可能是我可能选择的合理值：

>>> 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来解决这个问题，并减去了位移从振幅：

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-')

我不是建模或统计方面的专家，因此具有较深知识的人可能会对此进行改进。我已经添加了一个笔记本，可以对问题进行全面分析，包括如何生成示例数据。