JointFitter#

在某些情况下,人们可能希望将多个数据集与共享参数的模型相匹配。这可以通过使用 astropy.modeling.fitting.JointFitter 。基本上,该修配器设置有一系列已定义的模型、不同模型之间通用的参数以及这些参数的初始值。然后,安装工被称为提供尽可能多的x和y数组,每个要适合的模型一个。拟合参数是将模型联合拟合到组合数据集的结果。

备注

JointFitter使用scipy.Optimize.leastsq。此外,它目前不支持固定参数、绑定参数或绑定参数。

示例:光谱线#

该示例用于具有不同光谱分辨率的两个光谱段,它们在两个段之间重叠的波长区域中具有相同的光谱线。

我们需要定义一个具有平均波长、面积和宽度参数的高斯函数。这是所需的 astropy.modeling.functional_models.Gaussian1D 函数具有平均波长、中心强度和宽度参数,但在不同的光谱分辨率下,一条线的中心强度会不同,但面积将是相同的。

首先,本例需要导入

>>> # imports
>>> import math
>>> import numpy as np
>>> from astropy.modeling import fitting, Fittable1DModel
>>> from astropy.modeling.parameters import Parameter
>>> from astropy.modeling.functional_models import FLOAT_EPSILON

现在,使用面积而不是强度作为参数定义AreaGaussian1D。这个新的模型是标准的Gaussian1D模型的修改和修剪版本。

>>> class AreaGaussian1D(Fittable1DModel):
...   """
...   One dimensional Gaussian model with area as a parameter.
...
...   Parameters
...   ----------
...   area : float or `~astropy.units.Quantity`.
...       Integrated area
...       Note: amplitude = area / (stddev * np.sqrt(2 * np.pi))
...   mean : float or `~astropy.units.Quantity`.
...       Mean of the Gaussian.
...   stddev : float or `~astropy.units.Quantity`.
...       Standard deviation of the Gaussian with FWHM = 2 * stddev * np.sqrt(2 * np.log(2)).
...   """
...   area = Parameter(default=1)
...   mean = Parameter(default=0)
...
...   # Ensure stddev makes sense if its bounds are not explicitly set.
...   # stddev must be non-zero and positive.
...   stddev = Parameter(default=1, bounds=(FLOAT_EPSILON, None))
...
...   @staticmethod
...   def evaluate(x, area, mean, stddev):
...       """
...       AreaGaussian1D model function.
...       """
...       return (area / (stddev * np.sqrt(2 * np.pi))) * np.exp(
...           -0.5 * (x - mean) ** 2 / stddev ** 2
...       )

模拟了拟合出的数据。第一个光谱段的光谱分辨率将比第二个段高2倍。第一段的波长从1到6,第二段的波长从4到10,从而产生从4到6的重叠波长区域。

>>> # Generate fake data
>>> mean = 5.1
>>> sigma1 = 0.2
>>> sigma2 = 0.4
>>> noise = 0.10
>>> # compute the central amplitudes so the lines in each segment have the
>>> # same area
>>> area = 1.5
>>> amp1 = area / np.sqrt(2.0 * math.pi * sigma1 ** 2)
>>> amp2 = area / np.sqrt(2.0 * math.pi * sigma2 ** 2)
>>> # segment 1
>>> rng = np.random.default_rng(147)
>>> x1 = np.linspace(1.0, 6.0, 200)
>>> y1 = amp1 * np.exp(-0.5 * (x1 - mean) ** 2 / sigma1 ** 2)
>>> y1 += rng.normal(0.0, noise, x1.shape)
>>> # segment 2
>>> x2 = np.linspace(4.0, 10.0, 200)
>>> y2 = amp2 * np.exp(-0.5 * (x2 - mean) ** 2 / sigma2 ** 2)
>>> y2 += rng.normal(0.0, noise, x2.shape)

现在定义适合的模型以及更适合使用的模型。然后对两个模拟数据集进行拟合。

>>> # define the two models to be fit
>>> gjf1 = AreaGaussian1D(area=1.0, mean=5.0, stddev=1.0)
>>> gjf2 = AreaGaussian1D(area=1.0, mean=5.0, stddev=1.0)
>>> # define the jointfitter specifying the parameters in common and their initial values
>>> fit_joint = fitting.JointFitter(
...    [gjf1, gjf2], {gjf1: ["area", "mean"], gjf2: ["area", "mean"]}, [1.0, 5.0]
... )
>>>
>>> # perform the fit
>>> g12 = fit_joint(x1, y1, x2, y2)

拟合参数表明,两个AreaGaussian1D模型的面积和平均波长完全相同,而宽度(Stddev)不同,反映了两段不同的光谱分辨率。

区域Gaussian1参数

>>> print(gjf1.param_names)
('area', 'mean', 'stddev')
>>> print(gjf1.parameters)
[1.49823951 5.10494811 0.19918164]

区域高斯2参数

>>> print(gjf1.param_names)
('area', 'mean', 'stddev')
>>> print(gjf2.parameters)
[1.49823951 5.10494811 0.39860539]

模拟数据和最佳拟合模型可以作图,表明两个AreaGaussian1D模型和两个光谱段之间有很好的一致性。

(png, svg, pdf)

../_images/jointfitter-1.png