如何使用 Astropy 拟合高斯分布
How to fit a Gaussian using Astropy
我正在尝试使用 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)")
我做错了什么或遗漏了什么？
我制作了一些类似于您的假数据，并在其上尝试了 运行 您的代码并获得了类似的结果。我认为问题在于，如果你不将模型的初始参数调整到至少与原始模型相似，否则 fitter 无论执行多少轮拟合都无法收敛。
如果我要拟合高斯分布，我喜欢根据计算 "eyeballing" 给初始模型一些初始参数，就像这样（在这里我将您的真实数据的通量和波长命名为 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 的一个功能是可以附加到某些模型的可选方法，以根据某些数据对其参数进行合理估计。所以对于 Gaussians 这样的方法会 return 很像我刚刚在上面计算的。不过，我不知道这是否曾经实施过。
同样值得注意的是，您的高斯分布将被反转（具有负振幅）并且它在通量轴上位移了大约 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-')
我不是建模或统计方面的专家，所以知识更渊博的人可能会对此有所改进。我添加了一个笔记本，其中包含我对问题的完整分析，包括我如何生成示例数据 here。
我正在尝试使用 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)")
我做错了什么或遗漏了什么？
我制作了一些类似于您的假数据，并在其上尝试了 运行 您的代码并获得了类似的结果。我认为问题在于，如果你不将模型的初始参数调整到至少与原始模型相似，否则 fitter 无论执行多少轮拟合都无法收敛。
如果我要拟合高斯分布，我喜欢根据计算 "eyeballing" 给初始模型一些初始参数，就像这样（在这里我将您的真实数据的通量和波长命名为 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 的一个功能是可以附加到某些模型的可选方法，以根据某些数据对其参数进行合理估计。所以对于 Gaussians 这样的方法会 return 很像我刚刚在上面计算的。不过，我不知道这是否曾经实施过。
同样值得注意的是，您的高斯分布将被反转（具有负振幅）并且它在通量轴上位移了大约 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-')
我不是建模或统计方面的专家，所以知识更渊博的人可能会对此有所改进。我添加了一个笔记本，其中包含我对问题的完整分析，包括我如何生成示例数据 here。