天体测量的数据中，因变量通常是二维的，如位置偏差、自行、章动角、极移等，且这些二维观测量的分量存在一定的相关性。教科书上常用的最小二乘估计法一般是针对一维因变量型的数据，而缺乏基于二维因变量型的观测数据进行最小二乘估计的方案，尤其是需要考虑观测量分量之间相关性的情形。在这里我尝试推导过二维情形下的最小二乘估计的目标函数，并给出了相应的正规方程及其矩阵形式。

本文的PDF版本见这里。

引言

根据《天文数据处理》中的描述[1]，参数估计问题可以转换为下述的数学问题：若两个观测量 \(x\) 和 \(y\) 存在某种关于 \(m\) 个待定参数 \(c_k\)（\(k=1,2,\dots,m\)）的形式已知的函数关系，对于 \(N\) 对观测值（\(x_i\)，\(y_i\)）（\(i=1,2,\dots,N\)），确定参数 \(c_k\) 的某一估计值 \(\hat{c}_k\)，使某个目标函数 \begin{equation} d = d(x_1,x_2,\dots,x_N; y_1, y_2, \dots, y_N; c_1, c_2, \dots, c_m) \tag{1} \label{eq:sample} \end{equation} 因此，曲线拟合就是对目标函数进行最优化计算，寻求是目标函数 \(d\) 取极值的一组参数值。显然，这样的表述是针对观测量 \(y\) 为一维的情形给出的。 当观测量 \(\mathbfit{y}\)（为了区分我使用了黑体）为二维变量时，我们可以将上面的描述外推，得到当因变量为二维变量时的最小二乘参数估计形式。

最小二乘估计的目标函数

回顾一维情形

首先考虑双变量型数据，数据形式可以表述为 \begin{equation} (x_i,\sigma_{x,i}, y_i, \sigma_{y,i}), i=1,2,\dots,n \label{eq:data-1d} \tag{2} \end{equation} 其中，\(x_i\) 为自变量，\(y_i\) 为因变量，也就是观测量，\(\sigma_{x,i}\) 和 \(\sigma_{y,i}\) 分别为相应的测量误差。一般假设观测量彼此之间相互独立，即对于任意的 \(i \ne j\)，\(y_i\) 与 \(y_j\) 相互独立。一般假设自变量的测量误差为零，即 \begin{equation} \sigma_{x,i} = 0 \tag{3} \end{equation} 在天文观测资料处理过程中，对观测误差的估计通常与观测值本身之间相互独立，即 \(y_i\) 与 \(\sigma_{y,i}\) 相互独立。为了讨论方便，假设观测量的理论模型为线性模型，具体形式为 \begin{equation} F(x_i, \mathbfit{c}) = \sum_{j=1}^m c_j F_j (x_i) \tag{4} \end{equation} 其中模型中有 \(m\) 个未知参数需要估计，待估参数向量可记为 \begin{equation} \mathbfit{c} = (c_1, c_2, \dots, c_m)^{\rm T} \tag{5} \end{equation} 则第 \(i\) 次的观测值可写为 \begin{equation} y_i = F(x_i, \mathbfit{c}) + \epsilon_i \label{eq:model-1d} \tag{6} \end{equation} 其中 \(\epsilon_i\) 为第 \(i\) 次观测值与理论值之间的偏差，即残差。归一化残差（即加权后的残差）为 \begin{equation} X_i = \frac{\epsilon_i}{\sigma_{y,i}} \tag{7} \end{equation} 最小二乘意义下的参数估计就是选取各观测点的归一化残差平方和作为目标函数 \begin{equation} X^2 = \sum_{i=1} ^N X_i^2 = \sum_{i=1} ^N \frac{\epsilon_i^2}{\sigma_{y,i}^2} = \sum_{i=1} ^N \frac{\left(y_i - \sum_{j=1} ^m c_j F_j (x_i)\right)^2}{\sigma_{y,i}^2} \tag{8} \label{eq:x2-1d} \end{equation} \(X^2\) 也就是常说的 \(\chi^2\)。选取一组合适的参数值 \(\hat{c}_k\)，使得目标函数 \(X^2\) 取极小值，此时，\(X^2\) 相对于未知参数 \(c_k\) 的偏导为零，即 \begin{equation} 0 = \frac{\partial X^2}{\partial c_j} = \sum_{i=1}^N \frac{\partial X_i^2}{\partial c_j}, j=1,2,\dots,m \tag{9} \label{eq:dfunc-1d} \end{equation}

延伸至二维情形

现在我们考虑当观测量 \(\mathbfit{y_i}\) 为二维变量时的情形，即

\begin{equation} \mathbfit{y_i} = \left( y^{(1)}_{i}, y^{(2)} _{i} \right)^{\mathrm T} \tag{10} \end{equation}

数据形式变为 \(\begin{equation} (x _i,\sigma_{x,i}, y^{(1)} _{i}, y^{(2)}_{i}, \sigma^{(1)} _{y,i}, \sigma^{(2)}_{y,i}, \rho _i), i=1,2,\dots,n \tag{11} \end{equation}\) 其中，\(\rho_i\) 表示第 \(i\) 次的观测值向量 \(\mathbfit{y_i}\) 两个分量 \(y^{(1)}_{i}\) 和 \(y^{(2)}_{i}\) 之间的（Pearson）相关性系数。 这里，我们再次假设任意两次观测的观测量分量之间相互独立，即对于任意的 \(i \ne j\)，\(y^{(1)}_{i}\) 和 \(y^{(2)}_{i}\)、\(y^{(1)}_{i}\) 和 \(y^{(2)}_{j}\)、\(y^{(1)}_{j}\) 和 \(y^{(2)}_{i}\)、\(y^{(1)}_{j}\) 和 \(y^{(2)}_{j}\) 均相互独立。

第 \(i\) 次的观测量 \(\mathbfit{y_i}\) 的两个分量对应的理论模型为 \(\begin{equation} \begin{aligned} F^{(1)}(x _i, \mathbfit{c}) = \sum_{j=1}^m c _j F^{(1)}_{j} (x _i) \\ F^{(2)}(x_i, \mathbfit{c}) = \sum_{j=1}^m c_j F^{(2)} _{j} (x_i) \end{aligned} \tag{12} \end{equation}\) 则第 \(i\) 次观测值的两个分量可以表示为 \(\begin{equation} \begin{aligned} y^{(1)} _{i} = F^{(1)} (x_i, \mathbfit{c}) + \epsilon^{(1)}_{i} \\ y^{(2)}_{i} = F^{(2)} (x_i, \mathbfit{c}) + \epsilon^{(2)}_{i} \end{aligned} \tag{13} \end{equation}\) 归一化的残差分别为 \(\begin{equation} X^{(1)} _{i} = \frac{\epsilon^{(1)}_{i}}{\sigma^{(1)} _{y,i}},\quad X^{(2)}_{i} = \frac{\epsilon^{(2)} _{i}}{\sigma^{(2)}_{y,i}} \tag{14} \end{equation}\) 借用Mignard等人[2]提出的 \(X\) 统计量，归一化的二维残差 \(X_i\) 可由下式导出

\[\begin{aligned} X _i^{2} &= { \left[\begin{array}{cc} X^{(1)} _{i} & X^{(2)} _{i} \end{array}\right] } \left[\begin{array}{cc} 1 & \rho _i \\\\ \rho _i & 1 \end{array}\right]^{-1} \left[\begin{array}{c} X^{(1)} _{i} \\\\ X^{(2)} _{i} \end{array}\right] \\\\ & = \frac{1}{1-\rho _i^2} \left[\begin{array}{cc} X^{(1)} _{i} & X^{(2)} _{i} \end{array}\right] \left[\begin{array}{cc} 1 & -\rho _i \\\\ -\rho _i & 1 \end{array}\right]\left[\begin{array}{c} X^{(1)} _{i} \\\\ X^{(2)} _{i} \end{array}\right] \\\\ & = \frac{1}{1-\rho _i^2}\left[\left( X^{(1)} _{i}\right)^2 + \left(X^{(2)} _{i}\right)^2 - 2\rho _i X^{(1)} _{i} X^{(2)} _{i} \right] \end{aligned} \tag{15}\]

目标函数即归一化的二维残差平方和为 \(\begin{equation} X^2= \sum _{i=1} ^N X_i^2 = \sum_{i=1} ^N \frac{1}{1-\rho_i^2}\left[ \left( X^{(1)} _{i} \right)^2 + \left( X^{(2)} _{i} \right)^2 - 2\rho _i X^{(1)} _{i} X^{(2)} _{i} \right] \tag{16} \label{eq:x2-2d} \end{equation}\) 当观测量\(\mathbfit{y_i}\)的两个分量之间不相关时，即 \(\begin{equation} \rho_i = 0, \quad i = 1, 2, \dots, N \tag{17} \end{equation}\) 此时，目标函数变为 \(\begin{equation} X^2= \sum _{i=1} ^N X_i^2 = \sum_{i=1} ^N \left( X^{(1)}_{i} \right)^2 + \sum_{i=1} ^N \left( X^{(2)}_{i} \right)^2 = \left( X^{(1)} \right)^2 + \left( X^{(2)} \right)^2 \tag{18} \end{equation}\) 即将两个分量视为相互独立的一维观测量时得到的目标函数之和。

由极大似然估计间接导出

在观测量\(\mathbfit{y}\)为一维情形下，若观测残差\(\epsilon_i\) 只包含随机噪声，即观测残差服从高斯分布 \(\begin{equation} \epsilon _i \sim N(0, \sigma_{y,i}) \tag{19} \label{eq:nor-dist-1d} \end{equation}\) 此时，极大似然估计等价于最小二乘估计。这一结论也可以推广到二维甚至更高维情形。现在我先简单给出一维情形下的不严格证明，之后再从二维情形下的极大似然估计出发给出二维情形下的最小二乘估计的目标函数。

假设在方程\eqref{eq:model-1d}对应的模型下，在待估参数 \(\mathbfit{c}\) 处得到一对观测值（\(x_i\)，\(y_i\)）的概率密度为 \(f(x _i, y_i; \mathbfit{c})\)。现在进行了 \(N\) 次观测之后，得到一组如方程\eqref{eq:data-1d}所示的观测值这一事件 \(p\) 发生的概率密度为[3] \(\begin{equation} L(\mathbfit{c}) = \prod_{i=1}^{N} f(x_i, y _i; \mathbfit{c}) \tag{20} \end{equation}\) 上式中的\(L(\mathbfit{c})\) 即似然函数。由于事件\(p\)是真实发生的，因此在待估参数\(\mathbfit{c}\)取得最佳估计 \(\hat{\mathbfit{c}}\) 时，\(L(\mathbfit{c})\) 取得极大值。 但因为 \(L(\mathbfit{c})\) 是连乘形式，不利于计算，常选取其对数形式，即对数似然函数 \(l(\mathbfit{c})\) 来进行参数估计，其表达式为 \(\begin{equation} l(\mathbfit{c}) = \ln L(\mathbfit{c}) = \sum_{i=1}^{N} \ln f(x _i, y_i; \mathbfit{c}) \tag{21} \end{equation}\) 这就是极大似然估计法的目标函数。不同于最小二乘估计，这里取目标函数的极大值。 在方程\eqref{eq:nor-dist-1d}的假设下，有 \(\begin{equation} f(x_i, y_i; \mathbfit{c}) = f(\epsilon _i; \mathbfit{c}) = \frac{1}{\sqrt{2\pi} \sigma_i} \exp{\left(-\frac{\epsilon_i^2}{2\sigma^2_{y,i}}\right)} \tag{22} \end{equation}\) 由此可以计算得 \(\begin{equation} l(\mathbfit{c}) = \sum_{i=1}^{N} \left[ -\ln{\left( \sqrt{2\pi} \sigma _i \right)} -\frac{\epsilon _i^2}{2\sigma^2 _{y,i}}\right] = - \sum_{i=1}^{N} \frac{\epsilon_i^2}{2\sigma^2_{y,i}} + {\rm const} \tag{23} \end{equation}\) 上式等号右边最后一项为与任一待估参数\(c_k\)均无关的常数项。 可见，对数似然函数 \(l(\mathbfit{c})\) 取极大值等价于方程\eqref{eq:x2-1d}中的 \(X^2\) 取极小值，因此，此时的极大似然估计等价于最小二乘估计。 我不加证明地将这一结论外推到二维乃至高维情形，并利用这一结论反推出二维情形下最小二乘估计法的目标函数。

在二维情形下，若观测残差 \(\mathbfit{\epsilon_i}\) 服从二维高斯分布，则有

\[\begin{aligned} & \quad f(\epsilon^{(1)} _i, \epsilon^{(2)} _i; \mathbfit{c}) \\\\ &= \frac{1}{2\pi \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} \sqrt{1-\rho^2 _i}} \exp{\left\{ -\frac{1}{2 (1-\rho _i^2)} \left[\left( \frac{\epsilon^{(1)} _i}{\sigma^{(1)} _{y,i}}\right)^2 + \left( \frac{\epsilon^{(2)} _i}{\sigma^{(2)} _{y,i}}\right)^2 -2\rho _i \frac{\epsilon^{(1)} _i}{\sigma^{(1)} _{y,i}}\frac{\epsilon^{(2)} _i}{\sigma^{(2)} _{y,i}} \right] \right \} } \\\\ &= \frac{1}{2\pi \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} \sqrt{1-\rho^2 _i}} \exp{\left \{ -\frac{1}{2 (1-\rho _i^2)} \left[\left( X _i^{(1)} \right)^2 + \left( X _i^{(2)}\right)^2 -2\rho _i X _i^{(1)}X _i^{(2)} \right] \right\} } \end{aligned} \tag{24}\]

由此可知对数似然函数 \(l(\mathbfit{c})\) 的表达式为

\[\begin{aligned} l(\mathbfit{c}) &= \sum _{i=1}^{N} \left\{ -\frac{1}{2 (1-\rho_i^2)} \left[\left( X _i^{(1)} \right)^2 + \left( X _i^{(2)}\right)^2 -2\rho _i X _i^{(1)}X _i^{(2)} \right] -\ln{ \left( 2\pi \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} \sqrt{1-\rho^2 _i} \right) } \right\} \\\\ & = - \sum _{i=1}^{N} -\frac{1}{2 (1-\rho _i^2)} \left[\left( X _i^{(1)} \right)^2 + \left( X _i^{(2)}\right)^2 -2\rho _i X _i^{(1)}X _i^{(2)} \right] + {\rm const} \end{aligned} \tag{25}\]

可见，上式的对数似然函数 \(l(\mathbfit{c})\) 取极大值等价于方程\eqref{eq:x2-2d}中的 \(X^2\) 取极小值。 这也反过来说明了方程\eqref{eq:x2-2d}的正确性。

二维情形下最小二乘估计法的实现

前面分别给出了最小二乘估计法在一维和二维情形下的目标函数，这里我将导出相应的正规方程，从而给出最小二乘估计法的具体实现。 我依然先推导一维情形的正规方程。 利用方程\eqref{eq:x2-1d}将方程\eqref{eq:dfunc-1d}展开，对任意的 \(j=1,2,\dots,m\) 均有 \(\begin{equation} 0 = \sum _{i=1} ^N \frac{\partial X^2_i}{\partial c_j} = \sum_{i=1} ^N 2X_i \frac{\partial X_i}{\partial c _j} = \sum_{i=1} ^N -2\left( y _i - \sum_k^m c_k F_k (x _i) \right) \frac{F_j (x_i)}{\sigma^2_{y,i}} \tag{26} \label{eq:par-eq-1d} \end{equation}\) 进一步展开，有 \(\begin{equation} \sum_{i=1} ^N \sum_k^m \frac{c _k F_j (x_i) F_k (x _i)}{\sigma^2_{y,i}} = \sum _{i=1} ^N \frac{F_j (x_i) y_i}{\sigma^2 _{y,i}} , j=1,2,\dots,m \tag{27} \label{eq:nor-eq-1d} \end{equation}\) 考虑一个模型理论值向量 \(\begin{equation} \mathbfit{F} = \left( F(x_1,\mathbfit{c}), F(x_2,\mathbfit{c}), \dots, F(x_N,\mathbfit{c}) \right)^{\rm T} \tag{28} \end{equation}\) 向量 \(\mathbfit{F}\) 相对于 \(\mathbfit{c}\) 的偏导构成了雅可比矩阵 \(\begin{equation} \mathcal{J} = \frac{\partial \mathbfit{F}}{\partial \mathbfit{c}}, \quad \mathcal{J} _{ij}=\frac{\partial F(x_i,\mathbfit{c})}{\partial c_j}=F_j(x _i,\mathbfit{c}) \tag{29} \end{equation}\) 向量 \(\mathbfit{F}\) 相对于任一待估参数 \(c_j\) 的偏导构成了向量 \(\begin{equation} \frac{\partial \mathbfit{F}}{\partial c_j} = \left( F_j(x _1,\mathbfit{c}), F_j(x_2,\mathbfit{c}), \dots, F_j(x _N,\mathbfit{c}) \right)^{\rm T} \tag{30} \end{equation}\) 可见，\(\dfrac{\partial \mathbfit{F}}{\partial c_j}\)是雅可比矩阵\(\mathcal{J}\)的第\(j\)列，而 \(\left( \dfrac{\partial \mathbfit{F}}{\partial c_j}\right)^{\rm T}\) 是雅可比矩阵的转置矩阵 \(\mathcal{J}^{\rm T}\) 的第 \(j\) 行。 引入一个加权矩阵 \(\mathcal{M}\)，它可以看作是观测向量 \(\mathbfit{y}\) 的协方差矩阵 \(\mathcal{V}\) 的求逆，其表达式分别为

\[\mathcal{V} = \left(\begin{array}{cccc} \sigma _{y,1}^{2} & O & \dots & O \\\\ O & \sigma _{y,2}^{2} & \dots & O \\\\ \vdots & \vdots & \ddots & \vdots \\\\ O & O & \dots & \sigma _{y,N}^{2} \end{array}\right), \mathcal{M} = \left(\begin{array}{cccc} \sigma _{y,1}^{-2} & O & \dots & O \\\\ O & \sigma _{y,2}^{-2} & \dots & O \\\\ \vdots & \vdots & \ddots & \vdots \\\\ O & O & \dots & \sigma _{y,N}^{-2} \end{array}\right) = \mathcal{V}^{-1} \tag{31} \label{eq:cov-wgt}\]

一组观测值构成一个观测值向量，即 \begin{equation} \mathbfit{y} = (y _1, y_2, \dots, y_N)^{\rm T} \tag{32} \end{equation} 那么，方程\eqref{eq:nor-eq-1d}中对 \(i\) 求和可以写为

\[\begin{aligned} \sum _{i=1} ^N \dfrac{F _j (x _i) F _k (x _i)}{\sigma^2 _{y,i}} c _k &= \left( \dfrac{\partial \mathbfit{F}}{\partial c _j} \right)^{\rm T} \mathcal{M} \dfrac{\partial \mathbfit{F}}{\partial c _k} c _k \\\\ \sum _{i=1} ^N \dfrac{F _j (x _i) y _i}{\sigma^2 _{y,i}} &= \left( \dfrac{\partial \mathbfit{F}}{\partial c _j} \right)^{\rm T} \mathcal{M} \mathbfit{y} \end{aligned} \tag{33} \label{eq:nor-eq-2d-k}\]

类似地，方程\eqref{eq:nor-eq-2d-k}再对 \(k\) 求和也可以写成矩阵乘积的形式，为 \(\begin{equation} \left( \frac{\partial \mathbfit{F}}{\partial c _j} \right)^{\rm T} \sum_{k=1} ^m \frac{\partial \mathbfit{F}}{\partial c _k} c_k = \left( \frac{\partial \mathbfit{F}}{\partial c_j} \right)^{\rm T} \mathcal{J} \mathbfit{c} \tag{34} \end{equation}\) 方程\eqref{eq:nor-eq-1d}可以写成矩阵相乘的形式，有 \(\begin{equation} \left( \frac{\partial \mathbfit{F}}{\partial c_j} \right)^{\rm T} \mathcal{M} \mathcal{J} \mathbfit{c} = \left( \frac{\partial \mathbfit{F}}{\partial c _j} \right)^{\rm T} \mathcal{M} \mathbfit{y}, j=1,2,\dots,m \tag{35} \label{eq:JVMat-1d} \end{equation}\) 将所有的 \(j=1,2,\dots,m\) 对应的上述方程叠加在一起，则有 \begin{equation} \mathcal{J}^T \mathcal{M} \mathcal{J} \mathbfit{c} = \mathcal{J}^T \mathcal{M} \mathbfit{y} \tag{36} \label{eq:nor-mat-eq} \end{equation} 从而，可以求出未知参数向量的最小二乘估计 \begin{equation} \hat{\mathbfit{c}} = \left(\mathcal{J}^T \mathcal{M} \mathcal{J}\right)^{-1} \mathcal{J}^T \mathcal{M} \mathbfit{y} \tag{37} \label{eq:pmt} \end{equation}

当观测量\(\mathbfit{y}\)为二维时，情况要复杂得多。 类似于方程\eqref{eq:par-eq-1d}，有

\[\begin{aligned} \frac{\partial X^2 _i}{\partial c _j} &= \frac{1}{1-\rho _i^2} \left( 2 X^{(1)} _{i} \frac{\partial X^{(1)} _{i}}{\partial c _j} + 2 X^{(2)} _{i} \frac{\partial X^{(2)} _{i}}{\partial c _j} - 2\rho _i \frac{\partial X^{(1)} _{i}}{\partial c _j} X^{(2)} _{i} - 2\rho _i X^{(1)} _{i} \frac{\partial X^{(2)} _{i}}{\partial c _j} \right) \\\\ &= \frac{2}{1-\rho _i^2} \left[ -\left( y^{(1)} _{i} - \sum _k^m c _k F^{(1)} _{k} (x _i) \right) \frac{F^{(1)} _{j} (x _i)}{(\sigma^{(1)} _{y,i})^2} - \left( y^{(2)} _{i} - \sum _k^m c _k F^{(2)} _{k} (x _i) \right) \frac{F^{(2)} _{j} (x _i)}{(\sigma^{(2)} _{y,i})^2} \right. \\\\ &\quad \left. + \rho _i \left( y^{(2)} _{i} - \sum _k^m c _k F^{(2)} _{k} (x _i) \right) \frac{F^{(1)} _{j} (x _i)}{\sigma^{(1)} _{y,i}\sigma^{(2)} _{y,i}} + \rho _i \left( y^{(1)} _{i} - \sum _k^m c _k F^{(1)} _{k} (x _i) \right) \frac{F^{(2)} _{j} (x _i)}{\sigma^{(1)} _{y,i}\sigma^{(2)} _{y,i}} \right], \\\\ & j=1,2,\dots,m \end{aligned} \tag{38} \label{eq:nor-eq-2d}\]

将上式对 \(i\) 求和，并令其等于0，有

\[\tag{39} \begin{aligned} &\sum _{i=1} ^N \sum _k^m \frac{ c _k}{1-\rho _i^2} \left[ \frac{F^{(1)} _{j} (x _i) F^{(1)} _{k} (x _i)}{(\sigma^{(1)} _{y,i})^2} + \frac{F^{(2)} _{j} (x _i) F^{(2)} _{k} (x _i)}{(\sigma^{(2)} _{y,i})^2} - \frac{\rho _i \left( F^{(1)} _{j} (x _i) F^{(2)} _{k} (x _i) + F^{(2)} _{j} (x _i) F^{(1)} _{k} (x _i) \right)}{\sigma^{(1)} _{y,i}\sigma^{(2)} _{y,i}} \right] \\\\ &= \sum _{i=1} ^N \sum _k^m \frac{1}{1-\rho _i^2} \left[ \frac{F^{(1)} _{j} (x _i) y^{(1)} _{i}}{(\sigma^{(1)} _{y,i})^2} + \frac{F^{(2)} _{j} (x _i) y^{(2)} _{i}}{(\sigma^{(2)} _{y,i})^2} - \frac{\rho _i \left( F^{(1)} _{j} (x _i) y^{(2)} _{i} + F^{(2)} _{j} (x _i) y^{(1)} _{i} \right)}{\sigma^{(1)} _{y,i}\sigma^{(2)} _{y,i}} \right], \\\\ & j=1,2,\dots,m \end{aligned} \label{eq:nor-eq-2d-new}\]

方程\eqref{eq:nor-eq-2d-new}等式左边求和子项可以写为

\[\begin{aligned} &\quad \left( \begin{array}{cc} F^{(1)} _{j} (x _i) &F^{(2)} _{j} (x _i) \end{array} \right) \left( \begin{array}{cc} (\sigma^{(1)} _{y,i})^{-1} &0 \\\\ 0 &(\sigma^{(2)} _{y,i})^{-1} \end{array} \right) \frac{1}{1-\rho _i^2} \left( \begin{array}{cc} 1 &-\rho _i \\\\ -\rho _i &1 \end{array} \right) \left( \begin{array}{cc} (\sigma^{(1)} _{y,i})^{-1} &0 \\\\ 0 &(\sigma^{(2)} _{y,i})^{-1} \end{array} \right) \left( \begin{array}{c} F^{(1)} _{k} (x _i) \\\\ F^{(2)} _{k} (x _i) \end{array} \right) \\\\ &= \left( \begin{array}{cc} F^{(1)} _{j} (x _i) &F^{(2)} _{j} (x _i) \end{array} \right) \left( \begin{array}{cc} (\sigma^{(1)} _{y,i})^2 &\rho _i \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} \\\\ \rho _i \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} &(\sigma^{(2)} _{y,i})^2 \end{array} \right)^{-1} \left( \begin{array}{c} F^{(1)} _{k} (x _i) \\\\ F^{(2)} _{k} (x _i) \end{array} \right) \end{aligned} \tag{40}\]

上式第二项即为第 \(i\) 次观测量 \(\mathbfit{y _i}\) 的协方差矩阵。 同样地，方程\eqref{eq:nor-eq-2d-new}等号右边求和子项可以写为

\[\left( \begin{array}{cc} F^{(1)} _{j} (x _i) &F^{(2)} _{j} (x _i) \end{array} \right) \left( \begin{array}{cc} (\sigma^{(1)} _{y,i})^2 &\rho _i \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} \\\\ \rho _i \sigma^{(1)} _{y,i} \sigma^{(2)} _{y,i} &(\sigma^{(2)} _{y,i})^2 \end{array} \right)^{-1} \left( \begin{array}{c} y^{(1)} _{i} \\\\ y^{(2)} _{i} \end{array} \right) \tag{41}\]

我们可以构建合适的矩阵 \(\mathcal{J}\)、\(\mathcal{M}\) 和 \(\mathbfit{y}\)，来将方程\eqref{eq:nor-eq-2d-new}改写成如方程\eqref{eq:JVMat-1d}的矩阵形式。 至少有两种方案可以实现上述目的。 第一种方案的矩阵 \(\mathcal{J}\)、\(\mathcal{V}\) 和 \(\mathbfit{y}\) 如下所示：

\[\begin{aligned} \mathcal{J} &= \left( \begin{array}{cc} \dfrac{\partial \mathbfit{F}^{(1)}}{\partial \mathbfit{c}} &\dfrac{\partial \mathbfit{F}^{(2)}}{\partial \mathbfit{c}} \end{array} \right), \\\\ \mathcal{J} _{ij} &= \frac{\partial F^{(1)}(x _i,\mathbfit{c})}{\partial c _j}=F^{(1)} _j(x _i,\mathbfit{c}) \\\\ \mathcal{J} _{(i+N),j} &= \frac{\partial F^{(2)}(x _i,\mathbfit{c})}{\partial c _j}=F^{(2)} _j(x _i,\mathbfit{c}) \end{aligned} \tag{42}\] \[\tag{43} \label{eq:covmat-2d-1} \mathcal{V} = \left(\begin{array}{cccccc} \left(\sigma^{(1)} _{y,1}\right)^2 &\dots &O &\rho _1 \sigma^{(1)} _{y,1} \sigma^{(2)} _{y,1} &\dots &O \\\\ \vdots &\ddots &\vdots &\vdots &\ddots &\vdots \\\\ O &\dots &\left(\sigma^{(1)} _{y,N}\right)^2 &O &\dots &\rho _N \sigma^{(1)} _{y,N} \sigma^{(2)} _{y,N} \\\\ \rho _1 \sigma^{(1)} _{y,1} \sigma^{(2)} _{y,1} &\dots &O & \left(\sigma^{(2)} _{y,1}\right)^2 &\dots &O \\\\ \vdots &\ddots &\vdots &\vdots & \ddots &\vdots \\\\ O &\dots &\rho _N \sigma^{(1)} _{y,N} \sigma^{(2)} _{y,N} &O &\dots &\left(\sigma^{(2)} _{y,N}\right)^2 \end{array}\right)\]

\begin{equation} \tag{44} \label{eq:obs-2d-1} \mathbfit{y} = \left( \begin{array}{c} \mathbfit{y}^{(1)} \
\mathbfit{y}^{(2)} \end{array} \right)^{\rm T} = (y^{(1)} _1, y^{(1)}_2, \dots, y^{(1)} _N, y^{(2)}_1, y^{(2)} _2, \dots, y^{(2)}_N)^{\rm T} \end{equation}

另一种方案则是这样的：

\[\tag{45} \begin{aligned} \mathcal{J} _{(2i-1),j} &= \frac{\partial F^{(1)}(x _i,\mathbfit{c})}{\partial c _j}=F^{(1)} _j(x _i,\mathbfit{c}) \\\\ \mathcal{J} _{(2i),j} &= \frac{\partial F^{(2)}(x _i,\mathbfit{c})}{\partial c _j}=F^{(2)} _j(x _i,\mathbfit{c}) \end{aligned}\] \[\tag{46} \label{eq:covmat-2d-2} \mathcal{V} = \left(\begin{array}{ccccc} \left(\sigma^{(1)} _{y,1}\right)^2 &\rho _1 \sigma^{(1)} _{y,1} \sigma^{(2)} _{y,1} &\dots & O & O \\\\ \rho _1 \sigma^{(1)} _{y,1} \sigma^{(2)} _{y,1} &\left(\sigma^{(2)} _{y,1}\right)^2 &\dots &O &O \\\\ \vdots &\vdots &\ddots &\ddots &\vdots \\\\ \vdots &\vdots &\ddots &\left(\sigma^{(1)} _{y,N}\right)^2 &\rho _N \sigma^{(1)} _{y,N} \sigma^{(2)} _{y,N} \\\\ O &O &\dots & \rho _N \sigma^{(1)} _{y,N} \sigma^{(2)} _{y,N} &\left(\sigma^{(2)} _{y,N}\right)^2 \end{array}\right)\]

\begin{equation} \tag{47}\label{eq:obs-2d-2} \mathbfit{y} = (y^{(1)} _1, y^{(2)}_1, y^{(1)} _2, y^{(2)}_2, \dots, y^{(1)} _N, y^{(2)}_N)^{\rm T} \end{equation} \(\mathcal{M}\)可由方程\eqref{eq:covmat-2d-1}或\eqref{eq:covmat-2d-2}，结合方程\eqref{eq:cov-wgt}得出。利用方程\eqref{eq:nor-mat-eq}和\eqref{eq:pmt}，就可以得到未知参数的估计值。方程\eqref{eq:obs-2d-1}和\eqref{eq:obs-2d-2}实际上是将二维的观测值向量转换成了一维观测值向量。

当观测量 \(\mathbfit{y}\) 的两个分量不相关时，目标函数可以简化为 \(\begin{equation} \tag{48} X^2= \sum _{i=1} ^N X_i^2 = \sum_{i=1} ^N \left( X^{(1)}_{i} \right)^2 + \sum_{i=1} ^N \left( X^{(2)}_{i} \right)^2 = \left(X^{(1)}\right)^2 + \left(X^{(2)}\right)^2 \end{equation}\) 方程\eqref{eq:covmat-2d-1}和\eqref{eq:covmat-2d-2}给出的协方差矩阵 \(\mathcal{V}\) 和相应的加权矩阵 \(\mathcal{M}\) 均变为对角矩阵。

总结与讨论

借由上述推导，我们得到了在观测因变量为二维情况下的最小二乘估计的具体实现。对于更高维度的情形，可以进行类似的分析。对我个人而言，最容易理解的就是方程\eqref{eq:nor-mat-eq}到\eqref{eq:obs-2d-2}的矩阵形式。对于更加复杂的情况，如不同次的观测值 \(\mathbfit{y} _i\) 之间存在某种相关性且相关性系数（或者协方差）已知，只需将协方差矩阵 \(\mathcal{V}\) 的对应元素作出相应调整（从零调整为协方差）即可。观测量分量之间的相关性究竟如何影响到参数的最小二乘估计结果？对于这一问题，目前我还没有明确的结果和结论，后续我会利用一些典型的研究问题进行分析比较。

参考文献

• [1] 丁月蓉. 天文数据处理方法[M]. 南京大学出版社, 1998.

• [2] Mignard F, Klioner S, Lindegren L, et al. Gaia Data Release 1: Reference frame and optical properties of ICRF sources[J]. Astronomy & Astrophysics, 2016, 595: A5.

• [3] Feigelson E D, Babu G J. Modern statistical methods for astronomy: with R applications[M]. Cambridge University Press, 2012.