复合模型光谱拟合#
用例： 拟合活跃星系NGC 5548光谱中Lyman-alpha周围的复杂连续谱。
数据： 带有单位的3列ECSV文件。
工具： specutils，numpy。
跨仪器： 所有仪器。
文档： 本笔记本是STScI更大后处理数据分析工具生态系统的一部分。
引言#
在这个例子中，我们拟合活跃星系（NGC 5548）光谱中Lyman-alpha周围的复杂连续谱。这涉及到幂律（powerlaw）、消光（extinction）、各种宽度的发射线（emission lines）和吸收线（absorption lines）。只有光谱的某些区域（远离强吸收线）被拟合。模型有一些固定参数和一些自由参数，以及相互关联的参数。我们使用Astropy复合模型（compound-model）机制来同时拟合所有组件。
该示例仅部分使用了specutils。它将数据读取到Spectrum1D数据结构中。然而，当我们实际拟合模型时，我们只是抓取numpy数组（没有单位，因为那导致了一些错误）。
import numpy as np # 导入NumPy库，用于数值计算
import specutils # 导入specutils库，用于光谱数据处理
import time # 导入time库，用于时间相关操作
import astropy.modeling.fitting as fitting # 导入Astropy的拟合模块
from astropy.table import Table, QTable # 从Astropy导入表格处理模块
from astropy.nddata import StdDevUncertainty # 导入标准偏差不确定性类
import astropy.units as u # 导入Astropy单位模块
from astropy.visualization import quantity_support # 导入量支持模块，用于可视化
# 导入astropy库
import astropy
# 导入matplotlib库
import matplotlib
# 打印Astropy的版本
print("Astropy Version: ", astropy.__version__)
# 打印Numpy的版本
print("Numpy Version: ", np.__version__)
# 打印Specutils的版本
print("Specutils Version: ", specutils.__version__)
# 打印Matplotlib的版本
print("Matplotlib Version: ", matplotlib.__version__)
import matplotlib.pyplot as plt # 导入matplotlib库中的pyplot模块，用于绘图
# inline -- 非交互式单元格，notebook -- 交互式单元格
%matplotlib inline # 设置为inline模式，以便在Jupyter Notebook中直接显示图像
# %matplotlib notebook # 可选：设置为notebook模式，以便在Jupyter Notebook中实现交互式图像
%config InlineBackend.figure_format ='retina' # 设置图像格式为retina，以优化Macbook的显示效果
数据输入#
光谱（spectrum）：简单的三列ECSV文件，每列都有单位
拟合中包含的波长区域：简单的两列ASCII文件，包含下限和上限
首先设置路径名称。
开发者注释#
这些数据文件是小型ASCII文件，因此它们与笔记本一起存放在github仓库中。
datafile = "./n5548_mean_g130mb4.ecsv" # 定义数据文件路径
regionsfile = "./n5548_lyalpha_sample.dat" # 定义区域文件路径
使用Astropy的QTable读取表格，以便保留单位。
开发者笔记#
如果这个示例能提供一个关于在表格中编码不确定性的推荐实践的例子，那就很好了。问题：
不确定性应该有单位吗？
“不确定性”这个名称通常被理解为标准差（或等效物）。因此，如果用户希望将其定义为方差，例如，用户应该将列命名为
variance。要么添加一个表格列使其成为Astropy的不确定性对象，要么将其作为单独的变量存储，这样做有点不美观。
# 从指定的文件中读取数据，格式为ASCII ECSV
data = QTable.read(datafile, format='ascii.ecsv')
# 显示数据的前3行
data[:3]
将光谱放入Spectrum1D对象中#
我们需要先将不确定性转换为astropy不确定性对象。
开发者笔记#
我认为将不确定性类型作为Spectrum1D调用的选项隐藏起来可能会更简单，默认设置为标准差（standard-deviation）。
例如……我最初尝试了我认为合理的方式：
data['stdev'] = StdDevUncertainty(data['uncertainty'])
spectrum = specutils.Spectrum1D(spectral_axis=data['wavelength'],
flux=data['flux'],
uncertainty=data['stdev'])
这会引发一个。
# 创建一个标准差不确定性对象，使用数据中的不确定性值
uncertainty = StdDevUncertainty(data['uncertainty'])
# 创建一个一维光谱对象，包含波长、光谱和不确定性
spectrum = specutils.Spectrum1D(spectral_axis=data['wavelength'],
flux=data['flux'], uncertainty=uncertainty)
# 打印光谱对象的信息
print(spectrum)
读取光谱区域#
将这些转换为 specutils 光谱区域。
# 从指定的文件中读取区域数据，格式为ASCII
regionstab = QTable.read(regionsfile, format='ascii')
# 初始化一个空列表，用于存储子区域
subregions = []
# 遍历区域表中的每一对起始和结束值
for x0, x1 in zip(regionstab['col1'], regionstab['col2']):
# 将起始和结束值转换为天文单位（Ångström），并添加到子区域列表中
subregions += [(x0 * u.AA, x1 * u.AA)]
# 创建光谱区域对象，包含所有子区域
regions = specutils.SpectralRegion(subregions)
# 输出光谱区域对象
regions
从区域创建掩膜#
开发者注释#
我们可能可以通过提取区域来完成整个工作流程，但我认为这样做最终会比使用掩膜更复杂。
def mask_from_regions(spectrum, regions):
# 创建一个与光谱轴形状相同的布尔数组，初始值为False
mask = np.zeros(spectrum.spectral_axis.shape, dtype=np.bool_)
# 遍历每个区域
for r in regions:
# 创建一个子掩码，标记在当前区域内的光谱轴值
submask = (spectrum.spectral_axis > r.lower) & (spectrum.spectral_axis <= r.upper)
# 更新主掩码，将当前区域的子掩码合并到主掩码中
mask = mask | submask
# 将生成的掩码赋值给光谱对象的mask属性
spectrum.mask = mask
# 从指定的区域生成掩膜
mask_from_regions(spectrum, regions)
# 打印掩膜的最小值和最大值
print(spectrum.mask.min(), spectrum.mask.max())
绘制光谱并突出显示掩模的便利例程#
开发者笔记#
我认为这种展示掩模的方式相当优雅。我有点更喜欢这种方式，而不是让阴影区域一直延伸到顶部。（在展示残差时，我会这样做以便进行比较）。
另一种方法可能是将掩模区域的线条颜色改为灰色，或者改变不透明度，但这会更难实现，并且可能会更慢。
这个例程展示了使用单位和 Spectrum1D 的一个好处——你可以自动在轴上添加单位。
这也暴露了一个限制：目前，虽然 quantity_support() 可以自动在轴上添加单位，但它并没有添加标签。标记轴时的标准惯例是给出
label (units)。这里使用get_label()和set_label()添加标签的技巧可能对大多数用户来说过于晦涩。我们是否会将所有可视化功能保留在 specviz 中，还是应该在 specutils 中提供一些可视化工具？
def plot_spectrum(spectrum,color='b',alpha=0.5,figsize=(15,6),
label=None,ax=None,plot_mask=True,
mask_color='g',mask_alpha=0.1,title=None):
with quantity_support(): # 使用量支持上下文管理器
if ax is None: # 如果没有提供轴对象
fig, ax = plt.subplots(figsize=figsize) # 创建一个新的图形和轴
ax.plot(spectrum.spectral_axis,spectrum.flux,color=color,alpha=alpha,label=label) # 绘制光谱数据
if plot_mask: # 如果需要绘制掩膜
if spectrum.mask is not None: # 检查光谱是否有掩膜
ax.fill_between(spectrum.spectral_axis,0,spectrum.flux*spectrum.mask,
alpha=mask_alpha,color=mask_color) # 填充掩膜区域
ax.set_xlabel(r"Wavelength (" + ax.get_xlabel() + ")") # 设置x轴标签为波长
ax.set_ylabel(r"Flux (" + ax.get_ylabel() + ")") # 设置y轴标签为通量
return ax # 返回轴对象
ax = plot_spectrum(spectrum) # 绘制光谱并返回坐标轴对象
ax.set_title('NGC 5548 spectrum and mask to be used for fitting the model'); # 设置图表标题
初始猜测模型
该模型被直接导入：
import n5548_models as models # 导入n5548_models模块，通常用于处理NGC 5548的模型数据
上述定义的 .py 模块构建了一个或多个特殊类型的函数实例，这些函数定义在 astropy.modeling.models 包中，称为“复合模型”（compound model）。
复合模型只是 astropy.modeling.models 函数的组合，使用加法、乘法等组合运算符。
示例：
compound_model = models.PowerLaw1D(1.,1.) + models.Gaussian1D(1.,1.,1.)
将创建一个具有两个组件的复合模型实例。
一个实际的、可导入的模型定义将如下所示：
from custom_models import gaussian, powerlaw, ccmext
model1 = \
powerlaw(name = 'powerlaw1', # 功率律模型
amp = 6.586200E-14, # 振幅
x_0 = 1000.0, # 参考点
alpha = 0.4819233, # 指数
bounds = {'amp': (0., 1.00E-11), # 振幅的边界
'x_0': (0., 1.00E-11), # 参考点的边界
'alpha': (-5., 5.)}, # 指数的边界
fixed = {'x_0': True} # 固定参数
) \
\
gaussian(name = ‘C III 1176’, # 高斯模型
norm = 2.000000E-14, # 归一化因子 mean = 1195.006, # 均值 fwhm = 861.4926, # 全宽半高 bounds = {'norm': (0., 1.00E-10), # 归一化因子的边界 'mean': (1000., 2000.), # 均值的边界 'fwhm': (1000., 2000.), # 全宽半高的边界 'skew': (1., 1.)}, # 偏斜的边界 fixed = {'norm': True, # 固定归一化因子 'mean': True, # 固定均值 'fwhm': True, # 固定全宽半高 'skew': True}, # 固定偏斜 ) \
对于这个练习，我们选择名为 ‘model1’ 的模型：
# 将模型1赋值给复合模型
compound_model = models.model1
该模块使用了一些特殊的函数类型，这些类型通过重写模块 custom_models 中的标准函数来定义，位于 asytropy.modeling.models 中。
这种重写是必要的，因为 specfit 中的光谱成分不符合 astropy.modeling.models 中定义的标准。例如，specfit 中的高斯（Gaussian）由幅度（amplitude）、中心波长（central wavelength）、以 km/s 为单位的全宽半高（FWHM）和偏度（skewness）参数定义。而在 astropy.modeling.models 中，高斯由幅度、中心波长和与中心波长单位一致的宽度（width）定义，并且没有偏度参数。这些不兼容性通过在 fit_functions 模块中定义的子类得以解决。
开发者笔记#
似乎将带有偏度的高斯作为 Astropy 模型是有用的，因此这应该考虑上游移植。
对于全宽半高（FWHM），当存在偏度时，将其作为宽度参数可能更清晰，因为这比标准差更直观。
以下复合模型的打印语句还不错，但不是很好。可能值得为笔记本提供一个美观的打印格式，使其更易读。
此打印语句未指明哪些参数是固定的或浮动的。
它未指明哪些参数与其他参数相关联以及如何关联。
表格末尾有一个省略号。
# 打印复合模型的内容
print(compound_model)
ax = plot_spectrum(spectrum, label='data') # 绘制光谱数据并设置标签为'data'
ax.plot(spectrum.spectral_axis, compound_model(spectrum.spectral_axis.value), 'r', label='initial model') # 绘制初始模型，使用红色线条
ax.legend() # 显示图例
ax.set_title("Data and initial model") # 设置图表标题为"数据和初始模型"
拟合#
我们有数据和模型，现在需要将它们进行拟合。我们可以通过实例化一个Astropy拟合器来实现这一点，在这种情况下使用LevMarLSQFitter，它采用Levenberg-Marquardt算法进行最小二乘拟合。
fitter = fitting.LevMarLSQFitter() # 创建一个Levenberg-Marquardt最小二乘拟合器
在这个例子中，我们可以访问数据点的误差，因此我们可以使用它们的倒数作为拟合的权重。为了包含掩码（mask），其中像素值为1表示我们想要拟合的像素，而0表示我们想要排除的像素，我们将 \(\rm {mask} / \sigma\) 作为赋予拟合器的权重。
开发者笔记#
由于复合模型拟合在直接使用 Specutils 的
spectral_axis和flux量时没有成功，我决定去掉单位。在讨论中，似乎有一些机制可以在fit_line中去掉单位并重新添加（这比仅仅拟合一条直线更通用），但我没有想到去查看那里。尽管它与 scipy 拟合器的约定一致，我认为将 \(w = 1/\sigma\) 称为 权重 是相当令人困惑的。我检查了代码，对于这个拟合器，它确实是将其平方以用作权重，因此计算是正确的。但我认为在最小二乘拟合中，正确的权重更常见的定义是 \(w = 1/\sigma^2\)。由于这个不寻常的“权重”一词的使用，我现在已经深入研究，试图再次确认 Astropy 是否正确地进行了加权。
Nadia Dencheva 将检查非最小二乘拟合器在其说明中说使用 \(1/\sigma\) 时是否在计算优度函数之前进行了平方处理。
# 获取光谱的波长值
wavelength = spectrum.spectral_axis.value
# 获取光谱的 flux 值
flux = spectrum.flux.value
# 计算逆标准差，使用掩膜和不确定性数组
inverse_sigma = spectrum.mask / spectrum.uncertainty.array
要进行拟合，调用拟合器实例并传入数据、权重以及一些控制参数（如果需要的话）。我们也来做一些计时：
start_time = time.time() # 记录开始时间
fit_result = fitter(compound_model, wavelength, flux, weights=inverse_sigma, acc=1.E-30, maxiter=6000) # 使用fitter函数拟合数据
end_time = time.time() # 记录结束时间
print("Elapsed time: ", end_time - start_time) # 打印耗时
开发者备注#
我不确定如何解读下面的 fit_info 消息
实际和预测的平方和相对减少量均至多为 0.0000001
# 打印拟合器的状态信息
print(fitter.fit_info['message']) # 输出拟合过程中的消息
结果是另一个复合模型的实例，拟合值被设置为参数值：
print(fit_result) # 打印拟合结果
让我们打印一些派生结果：
查看拟合结果#
创建一个小例程来估计卡方（chisq）。它需要计算自由度（degrees of freedom），考虑到掩膜（mask）和自由参数（free parameters）的数量。
def chisq(x, y, err, mask, model, nfree):
# 计算卡方值
chisq = (y - model(x))**2 / err**2 # 计算每个点的卡方贡献
chisq = np.sum(chisq * mask) # 仅对有效数据点求和，mask用于选择有效点
npoints = np.sum(mask) # 计算有效数据点的数量
return np.sqrt(chisq / (npoints - nfree - 1)) # 返回归一化后的卡方值
确定有多少个固定参数和自由参数。
开发者注释#
这个过程相当复杂且繁琐。拟合器知道有多少个数据点和自由参数，但我认为它并没有将这些信息作为“元数据”传递给拟合结果。
fit_info是从scipy.optimize.leastsq返回的一个字典，但它并不包含自由参数的数量。
if 'fixed' in fit_result.parameter_constraints: # 检查是否存在固定参数约束
fix = np.asarray(fit_result.fixed.values()) # 将固定参数的值转换为NumPy数组
n_fixed_parameters = np.sum(np.where(fix, 1, 0)) # 计算固定参数的数量
else:
n_fixed_parameters = 0 # 如果没有固定参数，则数量为0
if 'tied' in fit_result.parameter_constraints: # 检查是否存在绑定参数约束
tie = np.asarray(fit_result.tied.values()) # 将绑定参数的值转换为NumPy数组
n_tied_parameters = np.sum(np.where(tie, 1, 0)) # 计算绑定参数的数量
else:
n_tied_parameters = 0 # 如果没有绑定参数，则数量为0
n_free_par = len(fit_result.parameters) - n_fixed_parameters - n_tied_parameters # 计算自由参数的数量
打印出 \(\chi^2\) 及其他相关信息。
开发者备注#
这种输出应该是标准输出。
chisq_in = chisq(wavelength, flux, spectrum.uncertainty.array, spectrum.mask, compound_model, n_free_par) # 计算输入模型的卡方值
chisq_out = chisq(wavelength, flux, spectrum.uncertainty.array, spectrum.mask, fit_result, n_free_par) # 计算输出模型的卡方值
print("chisq from input model: %f" % chisq_in) # 打印输入模型的卡方值
print("chisq from output model: %f" % chisq_out) # 打印输出模型的卡方值
print("Total data points: %d" % len(wavelength)) # 打印总数据点数量
print("Data points in wavelength ranges: %d" % np.sum(spectrum.mask)) # 打印在波长范围内的数据点数量
print("Number of free parameters: %d" % n_free_par) # 打印自由参数的数量
print("Number of iterations: %d" % fitter.fit_info['nfev']) # 打印迭代次数
print ("Fit engine took %d elapsed seconds." % (end_time - start_time)) # 打印拟合引擎消耗的时间
与协方差矩阵相关的每个自由参数的误差。
开发者备注#
这里是空的。但实际上它在这个拟合器的
fit_info中，所以我们可以提取它。在这个笔记本的更高级版本中，展示如何在三角（角落）图中绘制误差椭圆将是很好的，这是一种可视化参数对之间相关性的好方法。
cov = fitter.fit_info['param_cov'] # 获取拟合信息中的参数协方差矩阵
param_errors = {} # 初始化一个字典来存储参数误差
i = 0 # 初始化索引
if cov is not None: # 检查协方差矩阵是否存在
# 从协方差矩阵中提取方差
fit_errors = {} # 初始化一个字典来存储拟合误差
for param_name in fit_result.param_names: # 遍历所有参数名称
fixed = fit_result.fixed[param_name] # 检查参数是否被固定
tied = fit_result.tied[param_name] # 检查参数是否被绑定
if not fixed and not tied: # 如果参数既没有被固定也没有被绑定
fit_errors[param_name] = math.sqrt(cov[i,i]) # 计算标准误差并存储
i += 1 # 更新索引
# 将误差映射到输入模型的组件和参数
for param_name in fit_errors.keys(): # 遍历拟合误差字典的键
index, target_param_name = fit_result._param_map[param_name] # 获取参数的索引和目标参数名称
component_name = fit_result._submodels_names[index] # 获取组件名称
param_errors[(component_name, target_param_name)] = fit_errors[param_name] # 将误差存储到参数误差字典中
print(param_errors) # 打印参数误差字典
图表
输入和输出复合模型#
开发者笔记#
这个图的一个华丽但非常实用的交互版本可能具有以下功能：
缩放和平移等功能（可以通过使用
%matplotlib notebook来实现，但使用时滚动会很麻烦）切换数据的曲线和直方图的能力
切换遮罩的开关，改变透明度等功能
切换数据上的误差条的开关，改变透明度等功能
悬停信息：数据、残差、该点对卡方的百分比贡献
切换复合模型的各个组成部分的开关
在这个例子中，这并不是特别有用，因为最好将它们显示在 powerlaw*extinction 之上，但这确实是一个相当特殊的用例。
ax = plot_spectrum(spectrum, label='data') # 绘制光谱数据并设置标签为'data'
ax.plot(wavelength, compound_model(wavelength), 'r', label='initial model') # 绘制初始模型，颜色为红色
ax.plot(wavelength, fit_result(wavelength), 'g', label='fitted model') # 绘制拟合模型，颜色为绿色
ax.set_xlim((1200, 1275)) # 设置x轴范围为1200到1275
ax.set_ylim((0, 7e-13)) # 设置y轴范围为0到7e-13
ax.legend() # 显示图例
ax.set_title("Data and models") # 设置图表标题为"Data and models"
绘制残差#
专注于最有趣的区域。
ylim = (-2.e-13, 2.e-13) # 设置y轴的范围
fig, ax = plt.subplots(figsize=(15, 6)) # 创建一个15x6英寸的图形和坐标轴
ax.plot(wavelength, flux - compound_model(wavelength), label='original') # 绘制原始数据与复合模型的差异
ax.plot(wavelength, flux - fit_result(wavelength), label='original') # 绘制原始数据与拟合结果的差异
ax.fill_between(wavelength, ylim[0] * spectrum.mask, ylim[1] * spectrum.mask, # 填充y轴范围内的区域
alpha=0.1, color='g') # 设置填充的透明度和颜色
ax.set_xlim((1185., 1270.)) # 设置x轴的范围
ax.set_ylim(ylim) # 设置y轴的范围
ax.legend() # 显示图例
ax.set_xlabel(r'Wavelength ($\rm \AA$)', fontsize='large') # 设置x轴标签
ax.set_ylabel(r'Flux ($\rm erg\, cm^{-2}\, s^{-1}\, \AA^{-1}$)', fontsize='large') # 设置y轴标签
ax.set_title('Residuals', fontsize='large') # 设置图形标题
绘制单个组件#
在这里，我们将它们绘制在基本的 powerlaw*extinction 之上。
开发者笔记#
我试图标记或以其他方式指示哪些组件是相互关联的，但似乎这些信息会丢失，因为关联参数只是通过传递一个函数来完成的。可能更好的做法是让 Astropy 有一个基础类，供人们使用，该类将被引用参数的名称作为属性，并具有一个
__call__方法。这样至少可以建议人们在关联参数时使用这个方法，以便于后续的检查。
ax = plot_spectrum(spectrum, label='data') # 绘制光谱数据并设置标签为'data'
ax.plot(wavelength, fit_result(wavelength), 'g', label='fitted model') # 绘制拟合模型，颜色为绿色
plext = fit_result['powerlaw1'] * fit_result['extinction'] # 计算功率律与消光的乘积
ax.plot(wavelength, plext(wavelength), '--', alpha=0.5, linewidth=5, label="powerlaw + extinction") # 绘制功率律与消光的组合，虚线
for component in fit_result: # 遍历拟合结果中的每个组件
if component.name != 'powerlaw1' and component.name != 'extinction': # 排除功率律和消光
ax.plot(wavelength, (plext + component)(wavelength), label=component.name, alpha=0.5) # 绘制其他组件的组合
ax.set_xlim((1200, 1275)) # 设置x轴范围为1200到1275
ax.set_ylim((0, 9e-13)) # 设置y轴范围为0到9e-13
ax.legend(loc='upper right', fontsize='small') # 添加图例，位置在右上角，字体小
ax.set_title("Data and models"); # 设置图表标题为"Data and models"
看起来您没有提供任何代码。如果您能提供需要注释的Python代码，我将很乐意为您添加中文注释。请将代码粘贴到这里。